PseudoInverse for exact matrices
- To: mathgroup at smc.vnet.net
- Subject: [mg25154] PseudoInverse for exact matrices
- From: Gianluca Gorni <gorni at dimi.uniud.it>
- Date: Tue, 12 Sep 2000 02:58:43 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
Hello!
A mathematician that I know, who is not into programming,
and uses the computer as a calculator mainly, has been talking
of Mathematica with contempt and of a competing product with
praise, after trying to find the pseudoinverse of a 15 by 15
matrix of rational numbers. Mathematica run out of memory without
giving an answer, while the competition gave the result in
practically no time.
I was hurt in my pride, and looked into the matter.
Mathematica 4's documentation states that PseudoInverse[] works for both
numeric and symbolic matrices, and that it is based on singular
value decomposition. Now, I suppose that the singular value
decomposition is appropriate for floating point matrices, but it
looks crazy for exact matrices.
So I wrote the following algorithm that calculates the pseudoinverse
with only rational operations:
myPseudoInverse[m_List?MatrixQ /; Precision[m] === Infinity] :=
Module[{n = NullSpace[m]},
n = If[n === {}, {Table[0, {Length[First[m]]}]}, n];
Inverse[Transpose[m].m + Transpose[n].n].Transpose[m]];
and compared it with the built-in PseudoInverse on my (rather slow)
machine:
mat = Partition[Range[5*6], 5];
Timing[a = myPseudoInverse[mat];]
{0.0166667 Second, Null}
Timing[b = PseudoInverse[mat];]
{1.96667 Second, Null}
a == b
True
With the original 15 by 15 rational matrix, myPseudoInverse
gave an answer in less than a second, while PseudoInverse
crashed Mathematica and forced me to reboot.
In the meanwhile, I am afraid that Mathematica has made an outspoken and
influential detractor in the math department of a major university.
Gianluca Gorni