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