Re: Q: Conjugate
- To: mathgroup at smc.vnet.net
- Subject: [mg21410] Re: [mg21265] Q: Conjugate
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Thu, 6 Jan 2000 01:46:28 -0500 (EST)
- References: <199912200728.CAA01740@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
I think both symbolic PseudoInverse and QRDecomposition would
benefit from having a RealOnly->True/False option. In both cases the
internal efficiency and the result would both improve.
About a year ago I wrote simple Mathematica code (shown below) for
finding a symbolic pseudo-inverse that makes the assumption that all
symbols are real-valued. The gist is that we form a pair of weaker
generalized inverses using Gaussian elimination, and then put them
together to get the pseudo-inverse.
In a bit more detail, it goes like this. We use A+ to denote the
pseudo-inverse, also called the Moore-Penrose inverse, of the
matrix A. This is uniquely defined, although weaker notions of
a generalized inverse are not.
In the case where the matrix has full column rank, A+ can be
computed as
Inverse[Transpose[A].A]].Transpose[A]
To extend to the general case, we do the sort of computation one
would do to get this. Specifically, we augment the matrix
Transpose[A].A
with an identity matrix on the right, then perform row reduction,
and extract the augmented part. It turns out that this yields what
is called a {1,2,3}-generalized inverse, i1. But we can do similarly
with the transposed matrix, transpose the result, and get a
{1,2,4}-generalized inverse, i2. Then we can form i2.mat.i1 to get
the Moore-Penrose inverse.
The code below does not check whether the matrix has full column
rank. In that common case it will not be very efficient for the
simple reason that we can stop after finding the first generalized
inverse.
nearInverse[mat_?MatrixQ] := Module[
{len=Length[mat[[1]]], mattran=Transpose[mat], prod, aug, rowred},
prod = mattran.mat;
aug = Transpose[Join[Transpose[prod], IdentityMatrix[len]]];
rowred = RowReduce[aug, Method->OneStepRowReduction];
Together[Transpose[Take[Transpose[rowred], -len]] . mattran]
]
pseudoInverse[mat_?MatrixQ] := With[
{i1=nearInverse[mat], i2=Transpose[nearInverse[Transpose[mat]]]},
Together[i2.mat.i1]
]
Here is an example of rank 2.
mat = {{1,2*x,x^2,3},{1-4*x,-2+2*x,-6*x+x^2,-5}, {2*x,1,3*x,4}};
Now try the following:
InputForm[Timing[m1 = nearInverse[mat]]]
InputForm[Timing[m2 = Transpose[nearInverse[Transpose[mat]]]]]
The next few lines test that the two matrices satisfy certain
requirements of generalized inverses. In all cases we want matrices
of zeroes (of appropriate dimensions).
Together[m1.mat.m1 - m1]
Together[mat.m1.mat - mat]
Together[m2.mat.m2 - m2]
Together[mat.m2.mat - mat]
Together[Transpose[mat.m1] - mat.m1]
Together[Transpose[m2.mat] - m2.mat]
To get the pseudo-inverse, now do:
InputForm[Timing[m3 = pseudoInverse[mat]]]
You can check by comparing to
Together[ComplexExpand[PseudoInverse[mat]] - m3]
that m3 is correct. The appropriate pseudo-inverse tests are
indicated below. Again, we want matrices of zeroes in all cases.
Together[m3.mat.m3 - m3]
Together[mat.m3.mat - mat]
Together[Transpose[m3.mat] - m3.mat]
Together[Transpose[mat.m3] - mat.m3]
Here is a similar matrix, this time bivariate. Again, we indeed
get the Moore-Penrose inverse.
mat = {{1+y^2,2*x,x^2,3*y},
{1-4*x+y^2,-2+2*x-2*y,-6*x+x^2,3*y-8}, {2*x,1+y,3*x,4}};
In[30]:= InputForm[Timing[m3 = pseudoInverse[mat];]]
Out[30]//InputForm= {1.9000000000000004*Second, Null}
Much faster than using PseudoInverse this time.
In[31]:= InputForm[Timing[m4 =
Together[ComplexExpand[PseudoInverse[mat]]]]]
Out[31]//InputForm= {26.6*Second, Null}
I guess I should include the customary caveat to the effect that
one rarely needs a psuedo-inverse. For example there are better ways
to compute least-squares solutions to overdetermined or
underdetermined numeric linear systems using e.g. a QR or singular
values decomposition.
Daniel Lichtblau
Wolfram Research
Begin forwarded message:
From: Mark L Storch <mscc+ at andrew.cmu.edu>
To: mathgroup at smc.vnet.net
Subject: [mg21410] [mg21265] Q: Conjugate
Organization: Mse, Carnegie Mellon, Pittsburgh, PA
Hello,
I am using the PseudoInverse[] function to try to get some
analytic expressions for the solution to a relatively complex
system of equations. The mechanics of it works fine the problem
that I run into is in the output. I have a matrix equation of the
form A.B=C and I can solve for B by premultiplying both sides by the
inverse of A (in this case the pseudoInverse). My problem is that
all of the elements of A are simply numbers and the PseudoInverse
returns an expression that has a lot of Conjugate[] terms in it.
Is there any way to force Mathematica to recognize that Conjugate[x]
is simply x for all the terms other than setting a replacement rule
for each term individually?
Thank you,
Mark Storch