MPinv

57fc67e3-2320-4831-a059-30019495bfe6 56 4 5 false false 0 decimal

&[DATE] &[TIME] - &[FILENAME]

MPinv

Reference:

https://arxiv.org/ftp/arxiv/papers/0804/0804.4809.pdf

A T MPInvAux m A rows : n A cols : 12 mat A m n matrix ≡ 1 - 0 n m matrix 13 mat m n * 1 ≡ 2 - 1 A 1 el 1 -^11 mat 13 mat m n < flag 1 : M A A transpose * MEval : m n 12 mat n m 12 mat : 13 mat flag 0 : M A transpose A * MEval : 12 mat if ro 1 m range : co 1 n range : c 0 : L m n matrix : 14 mat k 1 n range r k n range : c1 1 c range : c c 1 + : 13 mat c 2 > L r c el M r k el L r c1 el L k c1 el transpose * (1 el - MEval : c 2 ≡ L r c el M r k el L r 1 el L k 1 el * - MEval : L r c el M r k el : if if L k c el T > L k c el L k c el sqrt MEval : k n < r k 1 + (n range : L r c el L r c el L k c el (1 -^* : 12 mat 0 if 21 line c c 1 - : if 31 line for flag c L co 1 c range el 13 mat 41 line if if 21 line :

Another tol's

t ∞ : k 1 n range M k k el 0 > (M k k el t < (& t M k k el : continue if for 21 line

v M abs vectorize 1 n n * (range el sort : t v v 0 ≡ vectorize augment eval 02 findrows 11 el T * eval : 21 line

A MRank A A length sqrt TOL * MPInvAux 2 el 11 line :

A MPInv flag c L 13 mat A A length sqrt TOL * MPInvAux : flag 0 < L M L transpose L * (1 -^MEval : flag 1 ≡ A transpose L * M * M * L transpose * MEval L M * M * L transpose * A transpose * MEval if 21 line if 21 line :

A MPInv.dn S U V' 13 mat A dn_LinAlgSVD transpose : Σ S : r 1 A rows A cols 12 mat min range eval : 12 mat V' transpose Σ r r el S r r el abs TOL > S r r el (1 -^0 if : (transpose * U transpose * 31 line :

A MRank.dn S A dn_LinAlgSVD 1 el : c S 1 S length range el : 12 mat c TOL > vectorize sum 21 line :

A MOrth.dn S U V' 13 mat A dn_LinAlgSVD transpose : c S 1 S length range el : c 1 c TOL > vectorize sum range eval : 12 mat U 1 A rows range c el 31 line :

m n MRand A m n matrix : r 1 m n * (range : A r el 109 -^109^random * eval : 31 line :

MPInv Examples

Pseudoinverse Summary

Dim

Rank

PseudoInverse

A n n el

A rank n ≡

A\002A\ A 1 -^≡

A m n el

A rank m ≡

A\002A\ A transpose A A transpose * (1 -^* ≡

A m n el

A rank n ≡

A\002A\ A transpose A * (1 -^A transpose * ≡

A m n el

A rank r ≡

A\002A\ G\002A\ F\002A\ * ≡

with

A F G * ≡

any full rank factorization of A

Examples

Setup

TOL 109 -^:

x MEval x :

A 211 mat :

A\002A\ A 1 -^: 0.5 11 mat

A MRank 1

A rank 1

G A MPInv : 0.5 11 mat

A\002A\ G - abs vectorize max 0

A 1212 mat :

A\002A\ A transpose A A transpose * (1 -^* : 0.2 0.4 21 mat

A MRank 1

A rank 1

G A MPInv : 0.2 0.4 21 mat

A\002A\ G - abs vectorize max 0

A 1221 mat :

A\002A\ A transpose A * (1 -^A transpose * : 0.2 0.4 12 mat

A MRank 1

A rank 1

G A MPInv : 0.2 0.4 12 mat

A\002A\ G - abs vectorize max 0

A 132422 mat :

A\002A\ A 1 -^: 2 - 1.5 1 0.5 - 22 mat

A MRank 2

A rank 2

G A MPInv : 2 - 1.5 1 0.5 - 22 mat

A\002A\ G - abs vectorize max 7.2 1014 -^*

A 020022 mat :

A\002A\ A MPInv.dn : 00 0.5 022 mat

A MRank 1

A rank 1

G A MPInv : 00 0.5 022 mat

A\002A\ G - abs vectorize max 0

A 12345678933 mat :

A\002A\ A MPInv.dn : 0.6389 - 0.1667 - 0.3056 0.0556 - 1.7283 1016 -^* 0.0556 0.5278 0.1667 0.1944 - 33 mat

A MRank 2

A rank 2

G A MPInv : 0.6389 - 0.1667 - 0.3056 0.0556 - 6.4056 1015 -^* - 0.0556 0.5278 0.1667 0.1944 - 33 mat

A\002A\ G - abs vectorize max 3.1 1013 -^*

A 1231456478376786444454 mat :

A\002A\ A MPInv.dn : 0.3818 - 0.1085 - 0.0833 - 0.0736 0.3643 0.6667 0.1667 0.3333 0.1667 - 0.6667 - 0.0039 0.062 0.1667 - 0.1008 0.0775 0.3818 - 0.1085 - 0.0833 - 0.0736 0.3643 45 mat

A MRank 3

A rank 3

G A MPInv : 0.3818 - 0.1085 - 0.0833 - 0.0736 0.3643 0.6667 0.1667 0.3333 0.1667 - 0.6667 - 0.0039 0.062 0.1667 - 0.1008 0.0775 0.3818 - 0.1085 - 0.0833 - 0.0736 0.3643 45 mat

A\002A\ G - abs vectorize max 1.7 1013 -^*

Rank

The code could be used for calculate the rank of a matrix.

n 28 :

a b 12 mat 10020012 mat :

c n 10 - :

A b n n MRand * a - :

P A dn_LinAlgSVD 2 el :

λ b c 1 MRand * a - :

λ n el 0 :

Matrix M is nxn, have rank c and eigenvalues λ

M P λ diag * P transpose * :

M 1.7294 2.1048 2.1048 14.9403 - 6.3433 4.0771 5.2002 - 4.6515 42 mat

SMath's rank have several bugs. For instance, can't calculate this

M rank #

M MRank 18

M MRank.dn 18

Timing

m n 12 mat 604012 mat :

A 200 m n MRand * 100 - :

A\002A\ A MPInv.dn :

A max A min 12 mat 99.9862 99.9942 - 12 mat

Id m n 12 mat min identity eval :

Slow but more accurate

x MEval x :

t.o 0 time :

G A MPInv :

0 time t.o - 26.3 s *

A\002A\ G - abs vectorize max 2.4 1016 -^*

Id m n > G A * A G * if :

Id Id - abs vectorize max 2.3 1013 -^*

Ido Id Id TOL sqrt > (* vectorize eval :

Ido 1000000000000001000000000000001000000000000001000000000000001000000000000001000000000000001000000000000001000000814 mat

Ido sum 40

equal to min(m,n)

Faster version

x MEval x eval :

t.o 0 time :

G A MPInv :

0 time t.o - 8.8 s *

A\002A\ G - abs vectorize max 1.1 1012 -^*

Id m n > G A * A G * if :

Id Id - abs vectorize max 7.5 1010 -^*

Reference:

https://arxiv.org/ftp/arxiv/papers/0804/0804.4809.pdf

Alvaro