svd nr

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

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

svd nr

SVD Decomposition, A = U S V', with A(m·n), U(m·n), Diag(S)(n), V(n·n) A K T N MSVD m A rows : (n A cols : (≥ K 1 ≡ μ m : μ n : if ro 1 m range : co 1 n range : 21 line U m μ matrix : U ro co el A ro co el : 21 line w n 1 matrix : t w : 12 mat σ 0 : g 0 : τ 0 : 13 mat 21 line 14 mat x sgn x sign x 0 ≡ (+ eval : a x y P a x y 12 mat * eval : 12 mat c s a b R c a * s b * + c b * s a * - 12 mat eval : k 1 n range κ k 2 + : t k el σ g * eval : g 0 : σ 0 : r k m range : 15 mat σ U r k el norm1 : (0 ≠ U r k el σ 1 -^U r k el * eval : s U r k el norme : 12 mat g U k k el sgn s * - eval : h U k k el g * s 2^- eval : 12 mat U k k el U k k el g - eval : κ 1 - n ≤ j κ 1 - (n range : s j el U r k el U r j el * eval : U r j el U r j el h 1 -^s j el * U r k el * + eval : 12 mat 21 line 0 if U r k el σ U r k el * eval : 51 line 0 if w k el σ g * eval : g 0 : σ 0 : c κ 1 - (κ 1 - n 12 mat max range eval : 14 mat k n < (σ U k c el norm1 : (0 ≠ (& U k c el σ 1 -^U k c el * eval : s U k c el norme : 12 mat g U k κ 1 - el sgn s * - eval : h U k κ 1 - el g * s 2^- eval : 12 mat U k κ 1 - el U k κ 1 - el g - eval : t c el h 1 -^U k c el * eval : 12 mat κ 1 - n ≤ j κ 1 - (m range : s j el U j c el transpose U k c el transpose * eval : U j c el U j c el s j el t c el * + eval : 12 mat 21 line 0 if U k c el σ U k c el * eval : 51 line 0 if τ τ w k el abs t k el abs + 12 mat norm1 : 51 line for τ τ T * eval : K 2 ≠ V n n matrix : S μ n matrix : 12 mat k n 1 range k n < c κ n range : g 0 ≠ V c k el g 1 -^U k k 1 + el 1 -^* U k c el * eval : κ n ≤ j κ n range : s j el U k c el transpose V c j el * eval : V c j el V c j el s j el V c k el * + eval : 12 mat 21 line 0 if 21 line 0 if V k c el 0 : V c k el 0 : 12 mat 31 line 0 if V k k el 1 : g t k el : κ k : 13 mat 21 line for k m 1 range r k m range : k n < U k k 1 + (n range el 0 : 0 if 12 mat k n ≤ w k el 0 ≠ k 1 + μ ≤ j k 1 + (μ range : s j el w k el 1 -^U k k el 1 -^* U r k el U r j el * (* eval : U r j el U r j el s j el U r k el * + eval : 31 line 0 if U r k el w k el 1 -^U r k el * eval : 21 line U r k el 0 : if U r k el 0 : if U k k el U k k el 1 + eval : 31 line for p n 1 range Iters 1 N range L p : L 1 ≥ L L 1 - : flag t L el abs τ < (L 1 ≡ (| eval : (break w L 1 - el abs τ < break continue if if 11 line for flag 1 ≠ c 0 : s 1 : 12 mat k L : k p ≤ k k 1 + : f t k el 12 mat t k el s c P : f abs τ < break g w k el : w k el f g 12 mat norme : (0 ≠ s c 12 mat w k el 1 -^f - g P : 0 if 12 mat j 1 m range U j L 1 - el U j k el 12 mat c s U j L 1 - el U j k el R : for 21 line if 21 line for 21 line 0 if z w p el : L p ≡ z 0 < w p el z - : V co p el V co p el - : 12 mat 0 if break 21 line Iters N < x w L el : y w p 1 - el : g t p 1 - el : h t p el : 14 mat f y z - (y z + (* g h - (g h + (* + 2 h * y * / eval : g f 112 mat norme : 12 mat f x z - (x z + (* x / h x / y f f sgn g * + / h - (* + eval : c 1 : s 1 : 12 mat k L : k p 1 - ≤ k k 1 + : h g 12 mat t k 1 + el s c P : t k el f h 12 mat norme : 12 mat t k el 0 ≠ s c 12 mat t k el 1 -^h f P : f g 12 mat c s x g R : 12 mat 0 if j 1 n range V j k el V j k 1 + el 12 mat c s V j k el V j k 1 + el R : for h y 12 mat w k 1 + el s c P : w k el f h 12 mat norme : 12 mat w k el 0 ≠ s c 12 mat w k el 1 -^h f P : f x 12 mat c s g y R : 12 mat 0 if j 1 m range U j k el U j k 1 + el 12 mat c s U j k el U j k 1 + el R : for 61 line for t L el 0 : t p el f : w p el x : 13 mat 61 line MaxIters N num2str rached. concat error if if 41 line for for ntx w abs vectorize co augment eval 1 csort : ntx ntx 2 col reverse eval : 21 line m n ≡ 0 r n 1 + (m range : ntx r el r : 21 line if 12 mat U ro ntx 1 μ range el el S co co el w ntx co el el : V co ntx co el el 13 mat 61 line p n 1 range Iters 1 N range L p : L 1 ≥ L L 1 - : flag t L el abs τ < (L 1 ≡ (| eval : (break w L 1 - el abs τ < break continue if if 11 line for flag 1 ≠ c 0 : s 1 : 12 mat k L : k p ≤ k k 1 + : f t k el 12 mat t k el s c P : f abs τ < break g w k el : w k el f g 12 mat norme : (0 ≠ s c 12 mat w k el 1 -^f - g P : 0 if 12 mat if 21 line for 21 line 0 if z w p el : L p ≡ z 0 < w p el z - : 0 if break 21 line Iters N < x w L el : y w p 1 - el : g t p 1 - el : h t p el : 14 mat f y z - (y z + (* g h - (g h + (* + 2 h * y * / eval : g f 112 mat norme : 12 mat f x z - (x z + (* x / h x / y f f sgn g * + / h - (* + eval : c 1 : s 1 : 12 mat k L : k p 1 - ≤ k k 1 + : h g 12 mat t k 1 + el s c P : t k el f h 12 mat norme : 12 mat t k el 0 ≠ s c 12 mat t k el 1 -^h f P : f g 12 mat c s x g R : 12 mat 0 if h y 12 mat w k 1 + el s c P : w k el f h 12 mat norme : 12 mat w k el 0 ≠ s c 12 mat w k el 1 -^h f P : f x 12 mat c s g y R : 12 mat 0 if 41 line for t L el 0 : t p el f : w p el x : 13 mat 61 line MaxIters N num2str rached. concat error if if 41 line for for w sort reverse eval 21 line if 61 line K 2 ≠ U S V 13 mat A transpose K T N MSVD : V S transpose U 13 mat 21 line A transpose K T N MSVD if if 11 line :

MSVD Decomposition, A = U S V', with A(m·n), U(m·n), Diag(S)(n), V(n·n)
MSVD(A,0) produces the "economy size" decomposition A K MSVD tol TOL A rows A cols 12 mat max * A norme * eval : K 2 ≠ U S V 13 mat A K 1015 -^30 MSVD : U tol MRZ S tol MRZ V tol MRZ 13 mat 21 line A K 1015 -^30 MSVD tol MRZ if 21 line :

MSVD Decomposition, A = U S V', with A(m·n), U(m·n), Diag(S)(n), V(n·n) A MSVD A 1 MSVD 11 line :

%Here is the code. To try it out, set your matrix to D in my code. %*************start code********************* D=A; cst=.000000001; %regularization parametr, similar to tol [r,c]=size(D); %get the number of row and columns of X Y=eye(r); Q(:,1)=D(:,1); I=eye(r); %generate Identity matrix of dim rxr sum=zeros([r,r]); %initialize sum with a rxc matrix with zeros for k=1:c-1, sum=sum + Q(:,k)*Q(:,k)'/(cst + D(:,k)'*Q(:,k)); Q(:,k+1)= (I-sum)*D(:,k+1); end for j=1:r, S=zeros([r,1]); M(c,j)=((Q(:,c)'*Y(:,j)))/(cst + D(:,c)'*Q(:,c)); for i=1:c-1, S=S + M(c-i+1,j)*D(:,c-i+1); dm= 1/(cst + D(:,c-i)'*Q(:,c-i)); M(c-i,j)=(dm*Q(:,c-i)')*(Y(:,j) - S); end end %*************end code*******

MPinv

A MRank w A 21015 -^30 MSVD : c w 1 w length range el : c TOL > vectorize sum 31 line :

Pseudonverse A T MPInv m A rows : n A cols : 12 mat A m n matrix ≡ m n matrix m n * 1 ≡ A 1 el 1 -^11 mat U S V 13 mat A 01015 -^30 MSVD : Σ S : r 1 m n 12 mat min range eval : 13 mat V Σ r r el S r r el abs T > S r r el (1 -^0 if : (transpose * U transpose * 21 line if if 21 line :

A MPInv A TOL A rows A cols 12 mat max * A norme * eval MPInv 11 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 1012 -^:

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 3.4 1015 -^*

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 3.8 1015 -^*

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 3.3 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 - 8.2757 1017 -^* 0.0556 0.5278 0.1667 0.1944 - 33 mat

A\002A\ G - abs vectorize max 61015 -^*

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.8 1014 -^*

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 8.9511 - 3.0757 - 3.0757 - 9.7087 2.0139 8.8169 - 7.8875 - 3.167 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 304012 mat :

A 200 m n MRand * 100 - :

A\002A\ A MPInv.dn :

A max A min 12 mat 99.7968 99.9745 - 12 mat

Id m n 12 mat min identity eval :

t.o 0 time :

G A MPInv :

0 time t.o - 50.5 s *

A\002A\ G - abs vectorize max 2.2 1015 -^*

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

Id Id - abs vectorize max 6.9 1013 -^*

Ido Id Id TOL sqrt > (* vectorize eval :

Ido 1000000000000001000000000000001000000000000001000000000000001000000000000001000000000000001000000000000001000000814 mat

Ido sum 30

equal to min(m,n)

Alvaro