# PseudoInverse of a variable matrix is very slow

I was trying to PseudoInverse a matrix of 4x4 size with variables as below and it took about 10 minutes and no results yet. Then I lost the patience and aborted it. I am wondering if it is a general case that PseudoInverse function of a variable matrix is very slow or is there any problem with the matrix?

Can someone give some insight into this?

Y = MatrixForm[{{gm1 + 1/RG + 1/ro1 + 1/Z1, -(1/Z1), -gm1 - 1/RG, -(1/ro1)},
{-(1/Z1), gm2 + 1/RG + 1/ro2 + 1/Z1, -(1/ro2), -gm2 - 1/RG}, {-(1/RG), -gm2 - 1/ro2, 1/RG + 1/ro2,
gm2}, {-gm1 - 1/ro1, -(1/RG), gm1, 1/RG + 1/ro1}}];
Z = PseudoInverse[Y];


UPDATE:

I ran it again and it finished after a long time with an output saying large output. Is there any way to speed up the calculation?

Also I found the discussion from the link and checked the function Gianluca Gorni wrote. It gave result almost immediately!

• try without MatrixForm as that is a wrapper. But why you are using PseudoInverse on symbolic matrix? Also, it is singular. Try Inverse and see. – Nasser Mar 11 '18 at 22:52
• I have just tried it and it doesn't seem to work. Still running and maybe 10 minutes or more... – anhnha Mar 11 '18 at 22:58
• Do not use the bugs tag until your observations have been confirmed to be a bug. – J. M. will be back soon Mar 11 '18 at 23:02
• @Nasser: I know that it is singular already. PseudoInverse  also works on singular matrix and I want to use that for some calculation. – anhnha Mar 11 '18 at 23:04
• @Nasser: I've just updated the post again and included the link to Gianluca Gorni's function. I checked his function and it gave the result almost immediately. Amazing! – anhnha Mar 11 '18 at 23:36

I'm writing this so that people have an idea of what they're up against.

First, let us temporarily replace the reciprocals with different symbols, to temporarily mitigate expression swell:

mn = {{gm1 + 1/RG + 1/ro1 + 1/Z1, -1/Z1, -gm1 - 1/RG, -1/ro1},
{-1/Z1, gm2 + 1/RG + 1/ro2 + 1/Z1, -1/ro2, -gm2 - 1/RG},
{-1/RG, -gm2 - 1/ro2, 1/RG + 1/ro2, gm2},
{-gm1 - 1/ro1, -1/RG, gm1, 1/RG + 1/ro1}} /.
{RG -> 1/rgi, Z1 -> 1/zi, ro1 -> 1/ri1, ro2 -> 1/ri2} // Simplify


When computing a complicated matrix expression, it helps to apply an appropriate decomposition to it. Let's use LUDecomposition[]:

{lu, piv, cond} = LUDecomposition[mn] // FullSimplify;


Note the LUDecomposition::sing warning; this is fine, since the proper factorization is computed anyway:

perm = IdentityMatrix[4][[piv]];
lf = LowerTriangularize[lu, -1] + IdentityMatrix[4];
uf = UpperTriangularize[lu];

perm.mn == lf.uf // Simplify
True


Note that things are arranged so that uf is the singular factor.

Conceptually, we can compute the pseudoinverse as

LeastSquares[uf, LinearSolve[lf, perm]]


and the LinearSolve[lf, perm] part is simple enough:

tmp = Simplify[LinearSolve[lf, perm]];


The problematic part is in computing the Moore-Penrose inverse of uf; structurally, it is similar to the following matrix:

cm = ReplacePart[UpperTriangularize[Array[C, {4, 4}]], {4, 4} -> 0]


$$\begin{pmatrix} \mathtt{C[1,1]} & \mathtt{C[1,2]} & \mathtt{C[1,3]} & \mathtt{C[1,4]} \\ 0 & \mathtt{C[2,2]} & \mathtt{C[2,3]} & \mathtt{C[2,4]} \\ 0 & 0 & \mathtt{C[3,3]} & \mathtt{C[3,4]} \\ 0 & 0 & 0 & 0 \\ \end{pmatrix}$$

which has a very complicated pseudoinverse (even after removing all the instances of Conjugate):

LeafCount[pinv = Simplify[PseudoInverse[cm] /. Conjugate -> Identity]]
1968


Then,

uinv = Simplify[pinv /.
DeleteCases[Thread[Flatten[cm] -> Flatten[uf]], HoldPattern[0 -> 0]]];


and finally

res = FullSimplify[uinv.tmp /. {rgi -> 1/RG, zi -> 1/Z1, ri1 -> 1/ro1, ri2 -> 1/ro2}]

• So does this mean Mathematica's PseudoInverse is not good? Did you see my update with the link to Gianluca Gorni 's function? I tried and it gave the result almost immediately. – anhnha Mar 12 '18 at 0:10
• I hadn't tried it yet, but looking at the code, his route seems to be more straightforward than mine. I would say that the problem with PseudoInverse[] is that it is trying to be too general (recall that Mathematica assumes everything is complex unless told otherwise). – J. M. will be back soon Mar 12 '18 at 0:18
• I would like to check how fast it is if all variables are real. How to add that assumption for the checking? – anhnha Mar 12 '18 at 0:25
• I don't really know; it doesn't seem like there's an appropriate option setting in PseudoInverse[]. – J. M. will be back soon Mar 12 '18 at 0:56

PseudoInverse has no option for assuming variables are real. This is a rare (I think) case where one can simply cheat using Assuming

mn = {{gm1 + 1/RG + 1/ro1 + 1/Z1, -1/Z1, -gm1 - 1/RG, -1/ro1}, {-1/Z1,
gm2 + 1/RG + 1/ro2 + 1/Z1, -1/ro2, -gm2 - 1/RG}, {-1/RG, -gm2 -
1/ro2, 1/RG + 1/ro2, gm2}, {-gm1 - 1/ro1, -1/RG, gm1,
1/RG + 1/ro1}};
mvars = Variables[mn];

AbsoluteTiming[
Assuming[Element[mvars, Reals], psi = PseudoInverse[mn];]]

(* Out[1711]= {9.957832, Null} *)

LeafCount[psi]

(* Out[1712]= 471786 *)


I would not claim this is good behavior, on the part of the code or the developer. It works because the internal code has some hooks to use Refine, and that catches the outer assumptions.

Some day we should make this a bona fide option.

• Wow, I learned to use Variables and Element to take all variables from the matrix and then assume they are reals. This method is so fast! – anhnha Mar 13 '18 at 13:29