MODULE LinEqSVD;
IMPORT Nbr := NbrRe, Vec := VecRe, Mtx := MtxRe, LinEq := LinEqRe, Errors := DataErrors, MathRe, Out := KernelLog ;
TYPE
Solver* = OBJECT
VAR u-, w-, vt-: Mtx.Matrix;
threshold*: Nbr.Real;
iterations*: LONGINT;
Reciprocal: Vec.Map;
zero, mag: Nbr.Real;
PROCEDURE & Initialize*( VAR A: Mtx.Matrix );
BEGIN
IF A # NIL THEN
u := A.Copy( ); LinEq.NormalizeMatrix( u, mag ); zero := 0; threshold := 1.0E-10; iterations := 30;
Reciprocal := reciprocal; decompose(); w := w*mag;
ELSE Errors.Error( "A NIL matrix was supplied." )
END
END Initialize;
PROCEDURE decompose;
VAR m, n: LONGINT;
rv1: POINTER TO ARRAY OF Nbr.Real;
anorm, scale, c, f, g, h, s, x, y, z: Nbr.Real; i, its, j, jj, k, l, nm: LONGINT; flag: BOOLEAN;
BEGIN
IF u = NIL THEN Errors.Error( "A NIL matrix was supplied." ); RETURN END;
m := u.rows; n := u.cols; NEW( rv1, n );
IF (w = NIL ) OR (w.rows # n) OR (w.cols # n) THEN NEW( w, 0, n, 0, n ) ELSE w.Multiply( zero ) END;
IF (vt = NIL ) OR (vt.rows # n) OR (vt.cols # n) THEN NEW( vt, 0, n, 0, n ) ELSE vt.Multiply( zero ) END;
g := 0; scale := 0; anorm := 0;
FOR i := 0 TO n - 1 DO
l := i + 1; rv1[i] := scale*g; g := 0; s := 0; scale := 0.0;
IF (i < m) THEN
FOR k := i TO m - 1 DO scale := scale + ABS( u.Get(k, i) ) END;
IF (scale # 0.0) THEN
FOR k := i TO m - 1 DO u.Set(k, i, u.Get(k, i)/scale); s := s + u.Get(k, i)*u.Get(k, i); END;
f := u.Get(i, i);
IF f >= 0 THEN g := -MathRe.Sqrt( s ) ELSE g := MathRe.Sqrt( s ) END;
h := f*g - s; u.Set(i, i, f - g);
FOR j := l TO n - 1 DO
s := 0.0;
FOR k := i TO m - 1 DO s := s + u.Get(k, i)*u.Get(k, j) END;
f := s/h;
FOR k := i TO m - 1 DO u.Set(k, j, u.Get(k, j) + f*u.Get(k, i)) END;
END;
FOR k := i TO m - 1 DO u.Set(k, i , u.Get(k, i)*scale) END;
END
END;
w.Set(i, i, scale*g); g := 0; s := 0; scale := 0.0;
IF (i < m) & (i # (n - 1)) THEN
FOR k := l TO n - 1 DO scale := scale + ABS( u.Get(i, k) ) END;
IF (scale # 0.0) THEN
FOR k := l TO n - 1 DO u.Set(i, k, u.Get(i, k)/scale); s := s + u.Get(i, k)*u.Get(i, k); END;
f := u.Get(i, l);
IF f >= 0 THEN g := -MathRe.Sqrt( s ) ELSE g := MathRe.Sqrt( s ) END;
h := f*g - s; u.Set(i, l, f - g);
FOR k := l TO n - 1 DO rv1[k] := u.Get(i, k)/h END;
FOR j := l TO m - 1 DO
s := 0.0;
FOR k := l TO n - 1 DO s := s + u.Get(j, k)*u.Get(i, k) END;
FOR k := l TO n - 1 DO u.Set(j, k, u.Get(j, k) + s*rv1[k] )END;
END;
FOR k := l TO n - 1 DO u.Set(i, k, u.Get(i, k)*scale) END;
END;
END;
anorm := Max( anorm, ABS( w.Get(i, i) ) + ABS( rv1[i] ) )
END;
FOR i := n - 1 TO 0 BY -1 DO
IF (i < (n - 1)) THEN
IF g # 0.0 THEN
FOR j := l TO n - 1 DO vt.Set(j, i, u.Get(i, j)/(u.Get(i, l)*g)) END;
FOR j := l TO n - 1 DO
s := 0;
FOR k := l TO n - 1 DO s := s + u.Get(i, k)*vt.Get(k, j) END;
FOR k := l TO n - 1 DO
IF (s # 0.0) THEN vt.Set(k, j, vt.Get(k, j) + s*vt.Get(k, i)) END;
END;
END;
END;
FOR j := l TO n - 1 DO vt.Set(i, j, 0.0); vt.Set(j, i, 0.0) END;
END;
vt.Set(i, i, 1.0); g := rv1[i]; l := i;
END;
FOR i := MinI( m, n ) - 1 TO 0 BY -1 DO
l := i + 1; g := w.Get(i, i);
FOR j := l TO n - 1 DO u.Set(i, j, 0.0) END;
IF (g # 0.0) THEN
g := 1.0/g;
FOR j := l TO n - 1 DO
s := 0.0;
FOR k := l TO m - 1 DO s := s + u.Get(k, i)*u.Get(k, j) END;
f := s*g/u.Get(i, i);
FOR k := i TO m - 1 DO
IF (f # 0.0) THEN u.Set(k, j, u.Get(k, j) + f*u.Get(k, i)) END;
END;
END;
FOR j := i TO m - 1 DO u.Set(j, i, u.Get(j, i)*g) END;
ELSE
FOR j := i TO m - 1 DO u.Set(j, i, 0.0) END;
END;
u.Set(i, i, u.Get(i, i) + 1.0) ;
END;
FOR k := n - 1 TO 0 BY -1 DO
its := 0;
LOOP
INC( its );
IF its > iterations THEN EXIT END;
flag := TRUE; l := k;
LOOP
nm := l - 1;
IF ((ABS( rv1[l] ) + anorm) = anorm) THEN flag := FALSE; EXIT
ELSIF ((ABS( w.Get(nm, nm) ) + anorm) = anorm) THEN EXIT
END;
DEC( l );
IF l < 0 THEN EXIT END;
END;
IF flag THEN
c := 0.0; s := 1.0; i := l;
LOOP
f := s*rv1[i]; rv1[i] := rv1[i]*c;
IF ((ABS( f ) + anorm) = anorm) THEN EXIT END;
g := w.Get(i, i); h := pythag( f, g ); w.Set(i, i, h); h := 1.0/h; c := g*h; s := -f*h;
FOR j := 0 TO m - 1 DO
y := u.Get(j, nm); z := u.Get(j, i); u.Set(j, nm, y*c + z*s); u.Set(j, i,z*c - y*s);
END;
IF i = k THEN EXIT END;
INC( i )
END;
END;
z := w.Get(k, k);
IF (l = k) THEN
IF (z < 0.0) THEN
w.Set(k, k,-z);
FOR j := 0 TO n - 1 DO vt.Set(j, k, -vt.Get(j, k)) END;
END;
EXIT;
END;
IF (its = iterations) THEN
Errors.Error( "Singular value decomposition iterations do not converge" );
u := NIL; w := NIL; vt := NIL; RETURN;
END;
x := w.Get(l, l); nm := k - 1; y := w.Get(nm, nm); g := rv1[nm]; h := rv1[k];
f := ((y - z)*(y + z) + (g - h)*(g + h))/(2.0*h*y); g := pythag( f, zero + 1 );
f := ((x - z)*(x + z) + h*((y/(f + sign( f )*ABS( g ))) - h))/x; c := 1.0; s := 1.0;
FOR j := l TO nm DO
i := j + 1; g := rv1[i]; y := w.Get(i, i); h := s*g; g := c*g; z := pythag( f, h ); rv1[j] := z; c := f/z;
s := h/z; f := x*c + g*s; g := g*c - x*s; h := y*s; y := y*c;
FOR jj := 0 TO n - 1 DO
x := vt.Get(jj, j); z := vt.Get(jj, i); vt.Set(jj, j, x*c + z*s); vt.Set(jj, i, z*c - x*s);
END;
z := pythag( f, h ); w.Set(j, j, z);
IF (z # 0.0) THEN z := 1.0/z; c := f*z; s := h*z; END;
f := c*g + s*y; x := c*y - s*g;
FOR jj := 0 TO m - 1 DO
y := u.Get(jj, j); z := u.Get(jj, i); u.Set(jj, j, y*c + z*s); u.Set(jj, i, z*c - y*s);
END;
END;
rv1[l] := 0.0; rv1[k] := f; w.Set(k, k, x);
END
END;
RETURN;
END decompose;
PROCEDURE PseudoInverse*( ): Mtx.Matrix;
VAR p, q, r, psinv: Mtx.Matrix;
BEGIN
q := w.Copy( ); q.MapAll( Reciprocal ); p := Mtx.Transpose( u ); r := vt*q; psinv := r*p;
RETURN psinv
END PseudoInverse;
PROCEDURE reciprocal( VAR x: Nbr.Real );
BEGIN
IF ABS( x ) < threshold THEN x := 0 ELSE x := 1.0/x END;
IF ABS( x ) < threshold THEN x := 0 END;
END reciprocal;
END Solver;
PROCEDURE pythag( a, b: Nbr.Real ): Nbr.Real;
VAR absa, absb, zero: Nbr.Real;
BEGIN
zero := 0; absa := ABS( a ); absb := ABS( b );
IF absa > absb THEN RETURN absa*MathRe.Sqrt( 1.0 + absb/absa*absb/absa );
ELSIF absb = 0 THEN RETURN zero
ELSE RETURN absb*MathRe.Sqrt( 1.0 + absa/absb*absa/absb )
END;
END pythag;
PROCEDURE Max( x, y: Nbr.Real ): Nbr.Real;
BEGIN
IF x > y THEN RETURN x ELSE RETURN y END;
END Max;
PROCEDURE MinI( i, j: LONGINT ): LONGINT;
BEGIN
IF i < j THEN RETURN i ELSE RETURN j; END;
END MinI;
PROCEDURE sign( x: Nbr.Real ): LONGINT;
BEGIN
IF x >= 0 THEN RETURN 1 ELSE RETURN -1 END;
END sign;
PROCEDURE Log( m: Mtx.Matrix );
VAR i, j: LONGINT;
BEGIN
FOR j := 0 TO m.rows - 1 DO
FOR i := 0 TO m.cols - 1 DO END;
Out.Ln;
END;
Out.Ln;
END Log;
PROCEDURE Test*;
VAR a, U, W, V, VT: Mtx.Matrix; zero: Nbr.Real; s: Solver;
BEGIN
Out.Ln; Out.String( "-------Singular Value Decomposition Test-------------" );
Out.String( " Singular Value Decomposition of Matrix A; " ); Out.Ln;
Out.String( " should yield eigenvalues w (in diagonal matrix) and " ); Out.Ln;
Out.String( " orthonormal matrices u & vT which, " ); Out.Ln;
Out.String( " multiplied by itself, should give unit matrices; " ); Out.Ln;
Out.String( " Finally, PseudoInverse of A is calculated, " ); Out.Ln;
Out.String( " and Pseudoinverse of Pseudoinverse should lead back to A (least square approximation of..) " ); Out.Ln;
NEW( a, 0, 4, 0, 5 );
zero := 0;
a.Set( 0, 0, zero + 1 ); a.Set( 1, 0, zero + 2 );
a.Set( 3, 0, zero + 3 ); a.Set( 0, 2, zero + 2 ); a.Set( 2, 1, zero + 3 ); a.Set( 3, 1, zero + 2 );
a.Set( 3, 4, zero + 4 );
NEW( s, a );
Out.String( "a:" ); Out.Ln; Log( a ); Out.String( "w:" ); Out.Ln; Log( s.w ); Out.String( "v,vt,v*vt:" ); Out.Ln;
Log( s.vt ); V := Mtx.Transpose( s.vt ); Log( V ); V := s.vt*V; Log( V ); Out.String( "u,ut,u*ut:" ); Out.Ln;
Log( s.u ); U := Mtx.Transpose( s.u ); Log( U ); U := s.u*U; Log( U ); Out.String( "u*w*v:" ); Out.Ln;
VT := Mtx.Transpose( s.vt ); W := s.w*VT; U := s.u*W; Log( U );
Out.String( "a, PseudoInverse(a), PseudoInverse(PseudoInverse(a)):" ); Out.Ln; Log( a );
a := s.PseudoInverse(); Log( a ); s.Initialize( a ); a := s.PseudoInverse(); Log( a );
END Test;
END LinEqSVD.
LinEqSVD.Test
System.Free LinEqSVD ~