MathGroup Archive 2000

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Q: Conjugate

  • To: mathgroup at
  • Subject: [mg21410] Re: [mg21265] Q: Conjugate
  • From: Daniel Lichtblau <danl at>
  • Date: Thu, 6 Jan 2000 01:46:28 -0500 (EST)
  • References: <>
  • Sender: owner-wri-mathgroup at

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


To extend to the general case, we do the sort of computation one
would do to get this. Specifically, we augment the matrix


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  

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]]]},

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 =
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>
To: mathgroup at
Subject: [mg21410] [mg21265] Q: Conjugate
Organization: Mse, Carnegie Mellon, Pittsburgh, PA

    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

  • Prev by Date: Re: Re: Problem with Rasters & Density Plots
  • Next by Date: double integrals
  • Previous by thread: Re: Unwanted definite-integral complex result from real integrand
  • Next by thread: double integrals