Given matrix A[0..m)[0..n), this routine computes its singular value decomposition,
A = U W VT . The matrix U replaces A on output. The diagonal matrix of singular values W is output as a vector W[0..n). The matrix V (not the transpose VT) is output as V[0..n)[0..n). The number of rows, m, must be greater or equal to the number of columns, n; if it is smaller, then a should be filled up to square with zero rows.
~ {
Pythag:
PROC [a, b:
REAL]
RETURNS [
REAL]
Pythag computes SqRt[a*a + b*b] without destructive overflow or underflow.
~ {
at: REAL ~ ABS[a];
bt: REAL ~ ABS[b];
SELECT bt
FROM
< at => {ct: REAL ~ bt/at; RETURN[at*RealFns.SqRt[1.0 + ct*ct]]};
# 0.0 => {ct: REAL ~ at/bt; RETURN[bt*RealFns.SqRt[1.0 + ct*ct]]};
ENDCASE => RETURN[0.0];
};
Sign:
PROC [a, b:
REAL]
RETURNS [
REAL]
~ {IF b >= 0.0 THEN RETURN[ABS[a]] ELSE RETURN[-ABS[a]]};
nm, el, k, j, its, i, jj: INT;
z, y, x, scale, s, h, g, f, c, aNorm: REAL;
rv1: ColumnN ¬ NEW[LinearSystem.VecSeq[n]]; -- w, rv1 are diag, superdiag of reduced A
IF m < n THEN ERROR badShape["You must add enough zero rows to make A square."];
Householder reduction to bidiagonal form.
g ¬ scale ¬ aNorm ¬ 0.0;
FOR i
IN [0..n)
DO
el ¬ i+1; -- Invariant until start diagonalizing reduction
rv1[i] ¬ scale*g; -- Finish row reduction from bottom of this loop
Zero column i below diagonal
g ¬ s ¬ scale ¬ 0.0;
IF i < m
THEN {
FOR k IN [i..m) DO scale ¬ scale + ABS[a[k][i]]; ENDLOOP;
IF scale # 0.0
THEN {
Find norm of column i, avoiding underflow and overflow, for Householder reflection
FOR k IN [i..m) DO a[k][i] ¬ a[k][i]/scale; s ¬ s + a[k][i]*a[k][i]; ENDLOOP;
f ¬ a[i][i];
g ¬ -Sign[RealFns.SqRt[s], f];
h ¬ f*g - s;
Apply reflection I-2uuT to remaining columns
a[i][i] ¬ f - g;
FOR j
IN [el..n)
DO
s ¬ 0.0;
FOR k IN [i..m) DO s ¬ s + a[k][i]*a[k][j]; ENDLOOP;
f ¬ s/h;
FOR k IN [i..m) DO a[k][j] ¬ a[k][j] + f*a[k][i]; ENDLOOP;
ENDLOOP;
FOR k IN [i..m) DO a[k][i] ¬ scale*a[k][i]; ENDLOOP;
};
};
w[i] ¬ scale*g;
Zero row i beyond superdiagonal
g ¬ s ¬ scale ¬ 0.0;
IF i < m
AND i < n-1
THEN {
FOR k IN [el..n) DO scale ¬ scale + ABS[a[i][k]] ENDLOOP;
IF scale # 0.0
THEN {
Find norm of row i, avoiding underflow and overflow, for Householder reflection
FOR k IN [el..n) DO a[i][k] ¬ a[i][k]/scale; s ¬ s + a[i][k]*a[i][k]; ENDLOOP;
f ¬ a[i][el];
g ¬ -Sign[RealFns.SqRt[s], f];
h ¬ f*g - s;
Apply reflection I-2uuT to remaining rows
a[i][el] ¬ f - g;
FOR k IN [el..n) DO rv1[k] ¬ a[i][k]/h; ENDLOOP;
IF i < m-1
THEN {
FOR j
IN [el..m)
DO
s ¬ 0.0;
FOR k IN [el..n) DO s ¬ s + a[j][k]*a[i][k]; ENDLOOP;
FOR k IN [el..n) DO a[j][k] ¬ a[j][k] + s*rv1[k]; ENDLOOP;
ENDLOOP;
};
FOR k IN [el..n) DO a[i][k] ¬ scale*a[i][k]; ENDLOOP;
};
};
aNorm ¬ MAX[aNorm, ABS[w[i]] + ABS[rv1[i]]];
ENDLOOP;
Accumulation of right-hand transformations.
FOR i
DECREASING
IN [0..n)
DO
IF i < n-1
THEN {
IF g # 0.0
THEN {
-- Double division to avoid possible underflow:
FOR j IN [el..n) DO v[j][i] ¬ (a[i][j]/a[i][el])/g; ENDLOOP;
FOR j
IN [el..n)
DO
s ¬ 0.0;
FOR k IN [el..n) DO s ¬ s + a[i][k]*v[k][j]; ENDLOOP;
FOR k IN [el..n) DO v[k][j] ¬ v[k][j] + s*v[k][i]; ENDLOOP;
ENDLOOP;
};
FOR j IN [el..n) DO v[i][j] ¬ v[j][i] ¬ 0.0; ENDLOOP;
};
v[i][i] ¬ 1.0;
g ¬ rv1[i];
el ¬ i;
ENDLOOP;
Accumulation of left-hand transformations.
FOR i
DECREASING
IN [0..n)
DO
el ¬ i+1;
g ¬ w[i];
FOR j IN [el..n) DO a[i][j] ¬ 0.0; ENDLOOP;
IF g # 0.0
THEN {
g ¬ 1.0/g;
FOR j
IN [el..n)
DO
s ¬ 0.0;
FOR k IN [el..m) DO s ¬ s + a[k][i]*a[k][j]; ENDLOOP;
f ¬ (s/a[i][i])*g;
FOR k IN [i..m) DO a[k][j] ¬ a[k][j] + f*a[k][i]; ENDLOOP;
ENDLOOP;
FOR j IN [i..m) DO a[j][i] ¬ a[j][i]*g; ENDLOOP;
}
ELSE {FOR j IN [i..m) DO a[j][i] ¬ 0.0; ENDLOOP};
a[i][i] ¬ a[i][i] + 1.0;
ENDLOOP;
Diagonalization of the bidiagonal form.
Based on implicit QR reduction of symmetric tridiagonal BTB, for upper bidiagonal B
FOR k
DECREASING
IN [0..n)
DO
-- Loop over singular values.
FOR its
IN [1..30]
DO
-- Loop over allowed iterations.
FOR el
DECREASING
IN [0..k]
DO
-- Test for splitting:
nm ¬ el-1; -- Note that rv1[0] is always zero.
IF (ABS[rv1[el]] + aNorm) = aNorm THEN EXIT; -- Thus exit here, never fall through
IF (ABS[w[nm]] + aNorm) = aNorm THEN GO TO SVNegligible;
REPEAT
SVNegligible => {
-- Cancellation of rv1[el], if el > 0 :
Use Givens rotations to chase superdiag out of row
c ¬ 0.0;
s ¬ 1.0;
FOR i
IN [el..k]
DO
f ¬ s*rv1[i];
IF (
ABS[f] + aNorm) # aNorm
THEN {
g ¬ w[i];
h ¬ Pythag[f, g];
w[i] ¬ h;
h ¬ 1.0/h;
c ¬ g*h;
s ¬ -(f*h);
Accumulate rotation in left-hand transformations
FOR j
IN [0..m)
DO
y ¬ a[j][nm];
z ¬ a[j][i];
a[j][nm] ¬ y*c + z*s;
a[j][i] ¬ z*c - y*s;
ENDLOOP;
};
ENDLOOP;
};
ENDLOOP;
z ¬ w[k];
IF el = k
THEN {
-- Convergence.
IF z < 0.0
THEN {
-- Singular value is made nonnegative.
w[k] ¬ -z;
Accumulate reflection in right-hand transformations
FOR j IN [0..n) DO v[j][k] ¬ -v[j][k]; ENDLOOP;
};
EXIT; -- No more iterations needed for this singular value.
};
IF its = 30
THEN {SIGNAL noConvergence["No convergence in 30 SVDecomp iterations"]};
x ¬ w[el]; -- Shift from bottom 2-by-2 minor:
nm ¬ k-1;
y ¬ w[nm];
g ¬ rv1[nm];
h ¬ rv1[k];
f ¬ ((y-z)*(y+z) + (g-h)*(g+h))/(2.0*h*y);
g ¬ Pythag[f, 1.0];
f ¬ ((x-z)*(x+z) + h*((y/(f+Sign[g, f])) - h))/x;
Next QR transformation:
c ¬ 1.0;
s ¬ 1.0;
Use alternating right- and left-hand Givens rotations to do implicit shift QR
Rotate first two rows, then chase new off-bidiag non-zero back and forth 'til off matrix
FOR j
IN [el..nm]
DO
i ¬ j+1;
g ¬ rv1[i];
y ¬ w[i];
h ¬ s*g;
g ¬ c*g;
Rotate rows to zero below diagonal
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;
Accumulate Givens rotation in right-hand transformations
FOR jj
IN [0..n)
DO
x ¬ v[jj][j];
z ¬ v[jj][i];
v[jj][j] ¬ x*c + z*s;
v[jj][i] ¬ z*c - x*s;
ENDLOOP;
Rotate columns to zero above superdiagonal
z ¬ Pythag[f, h];
w[j] ¬ z; -- Rotation can be arbitrary if z = 0.
IF z # 0.0 THEN {z ¬ 1.0/z; c ¬ f*z; s ¬ h*z};
f ¬ c*g + s*y;
x ¬ c*y - s*g;
Accumulate Givens rotation in left-hand transformations
FOR jj
IN [0..m)
DO
y ¬ a[jj][j];
z ¬ a[jj][i];
a[jj][j] ¬ y*c + z*s;
a[jj][i] ¬ z*c - y*s;
ENDLOOP;
ENDLOOP;
rv1[el] ¬ 0.0;
rv1[k] ¬ f;
w[k] ¬ x;
ENDLOOP;
ENDLOOP;
};