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