SVD - Singular Value Decomposition 14-FEB-2017
This is a CAS program that runs on rev 8151 and computes the singular value decomposition of a matrix. It provides a few additional commands. (Although the SVD command exists on the HP Prime, it is quite limited and does not handle a number of cases).
Usage: svd2(A) where A is an \( m \times n \) matrix (real or complex).
Output: { U, S, V } where U and V are unitary \( m \times m \) and \( n \times n \) matrices, respectively, and S is a list (vector) of singular values ordered from largest to smallest, and U*diag2(S,m,n)*V = A.
Additional commands:
Usage: diag2(v,m,n) where \( v = [ v_1, v_2, \dotsm, v_r ] \)
Output: \( m \times n \) matrix whose diagonal entries are \( v_1, v_2, \dotsm, v_r \) with r=min(m,n) and 0 elsewhere.
Usage: pinv(A,eps) where A is an \( m \times n \) matrix (real or complex).
Output: \( A^\dagger \), the pseudo-inverse of A. The parameter eps is a cut-off value so that all singular values of A less than eps are considered 0.
Usage: qrp(A) where A is an \( m \times n \) matrix (real or complex).
Output: { Q, R, P } where P is a permutation matrix, Q is a unitary matrix, and R is upper triangular, and AP = QR. This is a QR factoring algorithm with pivoting so that the diagonal entries of R are non-increasing in absolute value.
Usage: lqp(A) where A is an \( m \times n \) matrix (real or complex).
Output: { L, Q, P } where P is a permutation matrix, Q is a unitary matrix, and L is lower triangular, and PA = LQ. This is an LQ factoring algorithm with pivoting so that the diagonal entries of L are non-increasing in absolute value.
Usage: svdmmult(A,B) where A and B are matrices.
Output: M where M = AB. The built-in matrix multiplication treats small values (approximately 1e-15?) as 0 and therefore destroys accuracy when working with tiny elements.
The code is given below. The SVD algorithm comes from a few papers by Golub, Kahan, and Demmel. Basically, it transforms A into A = U1 * B * V1 where U1 and V1 are unitary and B is upper bidiagonal through a sequence of Householder reflections. Then, the QR algorithm is applied to B to get it into diagonal form via Givens rotations. These rotations are also unitary so that if G is a Givens rotation, U1*B = U1*G^T*G*B. Then U and V are the accumulation of all the necessary rotations needed to transform B into diagonal form. Lastly, the singular values are sorted in decreasing order.
The program was started on the Prime, then moved over to xcas, and then finally put back onto the Prime. There are a few places where the original code (for xcas) had to be modified due to bugs in the Prime firmware with respect to CAS programming (e.g. downward loops are broken in rev 8151). Also, I am sure there are a few places where the code could be optimized for speed but I have not had time to test. In all the test cases I used, U*diag(S)*V - A is a matrix whose entries were in the order of 1E-11 or smaller, including rank deficient matrices.
Lastly, a big thank you to Joe Horn, without whom testing the code against the HP48's version of SVD would have been a major pain, even with the Jazz library.
Bug reports much appreciated!
Updated 14-FEB-2017:
Fixed a bug in the svdmmult() routine;
Fixed m x 1 and 1 x n case for pinv()
Updated 07-FEB-2017:
Fixed a bug in which SVD of 2x2 was incorrectly computed (sine and cosines swapped).
Updated 04-FEB-2017:
This update implemented matrix multiplication that will not turn small values into 0. The built-in matrix multiplication routine turns small values into 0 which defeats the point of an accurate SVD algorithm. Also fixed a sign bug.
Code:
#pragma mode( separator(.,;) integer(h32) )
#cas
// Singular Value Decomposition
// by Han Duong (14-FEB-2017)
//
// A = U*S*V
//
// Golub-Kahan bidiagonalization via
// Householder reflectors followed
// by Demmel-Kahan QR using Givens
// rotations
//
// Demmel and Kahan, Accurate singular values
// of bidiagonal matrices
// SIAM J. Sci. Stat. Comput. Vol. 11, No. 5
// pp 873-912, Sept. 1990
//
// Global Variables
// t_svdA - A (m x n)
// t_svdU - U (m x m)
// t_svdV - V (n x n)
// t_svdS - S (2 x r); r=min(m,n)
//
// several tolerances defined within
// search for: svdtol, eps
//
// Thanks to Joe Horn, whose help
// enabled me to quickly single-step
// through the HP48's SVD algorithm
// via the Jazz library as a means of
// checking the code implemented
// herein
svd2(A):=
BEGIN
local r,c,m,d;
local k;
local S;
t_svdA:=approx(A);
// r=rows; c=columns
d:=dim(A);
r:=d(1); c:=d(2);
m:=min(r,c);
t_svdU:=idenmat(r);
t_svdV:=idenmat(c);
t_svdS:=makemat(0,2,m);
// bidiagonalization
if (r>=c) then
for k from 1 to m-1 do
t_svdS[1,k]:=svdUQHQA(k,k,r,c);
t_svdS[2,k]:=svdAPPHV(k,k+1,r,c);
end;
t_svdS[1,m]:=svdUQHQA(m,m,r,c);
else
for k from 1 to m-1 do
t_svdS[1,k]:=svdAPPHV(k,k,r,c);
t_svdS[2,k]:=svdUQHQA(k+1,k,r,c);
end;
t_svdS[1,m]:=svdAPPHV(m,m,r,c);
if (m>1) then svdtoUBD(m); end;
end;
svdBDQR(m,r,c);
svdPerm(m,r,c);
// S:=makemat(0,r,c);
// for k from 1 to m do
// S[k,k]:=t_svdS[1,k];
// end;
S:=t_svdS[1];
return(t_svdU,S,t_svdV);
END;
// Householder reflection to annihilate
// subcolumns; in column j, rows i+1 to r
// are transformed to 0
// U * A = U * Q^H * Q * A
// **********************************
svdUQHQA(j,k,r,c):=
BEGIN
local a,x1,t,d;
local b,v;
local m,n;
x1:=t_svdA[j,k];
// compute norm of rest of column
if (j<r) then
v:=subMat(t_svdA,j+1,k,r,k);
// a:=trn(v)*v;
// a:=a[1,1];
a:=svdvTb(v,v);
else
a:=0;
end;
// already householder reflector
if (a==0) and (im(x1)==0) then
return(re(x1));
end;
// norm of column
a:=sqrt(a+x1*conj(x1));
if (re(x1)>=0) then a:=-a; end;
x1:=x1-a;
t_svdA[j,k]:=x1;
d:=x1*a;
v:=subMat(t_svdA,j,k,r,k);
// Q*A where Q = I-t*v*vT
for n from k+1 to c do
b:=subMat(t_svdA,j,n,r,n);
// t:=trn(v)*b;
// t:=t[1,1]/d;
t:=svdvTb(v,b)/d;
for m from j to r do
// Q*b = b+t*dot(b,v)*v
t_svdA[m,n]:=t_svdA[m,n]+t*t_svdA[m,k];
end;
end;
// U*QH
d:=conj(d);
for m from 1 to r do
b:=subMat(t_svdU,m,j,m,r);
// t:=b*v;
// t:=t[1,1]/d;
t:=svddotbv(b,v)/d;
for n from j to r do
t_svdU[m,n]:=t_svdU[m,n]+t*conj(t_svdA[n,k]);
end;
end;
return(a);
END;
// Householder reflections to annihilate
// subrows; in row j, columns j+1 to c
// are transformed to 0
// A * V = A * P * P^H * V
// **********************************
svdAPPHV(j,k,r,c):=
BEGIN
local a,x1,t,d;
local b,v;
local m,n;
x1:=t_svdA[j,k];
if (k<c) then
v:=subMat(t_svdA,j,k+1,j,c);
// a:=v*trn(v);
// a:=a[1,1];
a:=svdbvT(v,v);
else
a:=0;
end;
if (a==0) and (im(x1)==0) then
return(re(x1));
end;
a:=sqrt(a+x1*conj(x1));
if (re(x1)>=0) then a:=-a; end;
x1:=x1-a;
t_svdA[j,k]:=x1;
d:=x1*a;
v:=subMat(t_svdA,j,k,j,c);
for m from j+1 to r do
b:=subMat(t_svdA,m,k,m,c);
// t:=b*trn(v);
// t:=t[1,1]/d;
t:=svdbvT(b,v)/d;
for n from k to c do
t_svdA[m,n]:=t_svdA[m,n]+t*t_svdA[j,n];
end;
end;
d:=conj(d);
for n from 1 to c do
b:=subMat(t_svdV,k,n,c,n);
// t:=v*b;
// t:=t[1,1]/d;
t:=svddotbv(v,b)/d;
for m from k to c do
t_svdV[m,n]:=t_svdV[m,n]+t*conj(t_svdA[j,m]);
end;
end;
return(a);
END;
// convert lower bidiagonal matrix to
// upper bidiagonal using series of
// Givens rotations
// **********************************
svdtoUBD(m):=
BEGIN
local a,b,t;
local v;
local k;
for k from 1 to m-1 do
a:=t_svdS[1,k]; b:=t_svdS[2,k];
v:=svdGetRotM(a,b);
t:=t_svdS[1,k+1];
t_svdS[1,k]:=v[3];
t_svdS[2,k]:=v[2]*t;
t_svdS[1,k+1]:=v[1]*t;
svdRotUGT(v[1],v[2],k,m);
end;
END;
// determine Givens rotation matrix
// [[ cs sn ] [[ a ] [[ r ]
// [ -sn cs ]] [ b ]] = [ 0 ]]
// returns [cs,sn,r]
// **********************************
svdGetRotM(a,b):=
BEGIN
local s,t,cs,sn;
if (a==0) then
if (b==0) then
return([1,0,b]);
else
return([0,1,b]);
end;
else
if (abs(a) > abs(b)) then
t:=b/a;
s:=sqrt(1+t*t);
cs:=1/s;
sn:=t*cs;
return([cs,sn,a*s]);
else
t:=a/b;
s:=sqrt(1+t*t);
sn:=1/s;
cs:=t*sn;
return([cs,sn,b*s]);
end;
end;
END;
// U*GT where GT=G^T is a Givens rotation
// applied to submatrix at (k,k)
// **********************************
svdRotUGT(cs,sn,k,r):=
BEGIN
local m;
local a,b;
for m from 1 to r do
a:=t_svdU[m,k];
b:=t_svdU[m,k+1];
t_svdU[m,k]:=cs*a+sn*b;
t_svdU[m,k+1]:=cs*b-sn*a;
end;
END;
// G*V where G is a Givens rotation
// applied to submatrix at (k,k)
// **********************************
svdRotGV(cs,sn,k,c):=
BEGIN
local n;
local a,b;
for n from 1 to c do
a:=t_svdV[k,n];
b:=t_svdV[k+1,n];
t_svdV[k,n]:=cs*a+sn*b;
t_svdV[k+1,n]:=cs*b-sn*a;
end;
END;
// v and b are both column matrices
// compute inner product
// **********************************
svddotbv(b,v):=
BEGIN
local d,k;
local s;
d:=dim(v);
s:=0;
d:=d[1];
for k from 1 to d do
s:=s+b[1,k]*v[k,1];
end;
return(s);
END;
// v and b are both column matrices
// compute v^H * b
// **********************************
svdvTb(v,b):=
BEGIN
local m,r,s;
r:=dim(v);
r:=r[1];
s:=0;
for m from 1 to r do
s:=s+b[m,1]*conj(v[m,1]);
end;
return(s);
END;
// v and b are both row matrices
// compute b * v^H
// **********************************
svdbvT(b,v):=
BEGIN
local n,c,s;
c:=dim(v);
c:=c[2];
s:=0;
for n from 1 to c do
s:=s+b[1,n]*conj(v[1,n]);
end;
return(s);
END;
// SVD of bidiagonal matrix
// **********************************
svdBDQR(r,m,n):=
BEGIN
local thresh, svdtol;
local k,t,ks,ke;
local smax, smin;
local ksold,keold,fdir;
local qrs,abss,abse;
if (r==1) then return(0); end;
svdtol:=3.0E-13;
t:=abs(t_svdS[1,1]);
thresh:=t;
for k from 1 to r-1 do
if (t>0) then
t:=abs(t_svdS[1,k+1])*(t/(t+abs(t_svdS[2,k])));
thresh:=min(thresh,t);
end;
end;
thresh:=svdtol*thresh/sqrt(r);
smax:=0; smin:=0;
ke:=r; fdir:=1;
ksold:=r+1; keold:=0;
WHILE (ke>1) DO
smax:=abs(t_svdS[1,ke]);
// find a block to start chasing buldge
// convergence criterion 1 (Demmel-Kahan)
ks:=ke-1;
REPEAT
// FOR ks FROM ke-1 DOWNTO 1 DO
abss:=abs(t_svdS[1,ks]);
abse:=abs(t_svdS[2,ks]);
if (abse<=thresh) then
t_svdS[2,ks]:=0;
break;
end;
smax:=max(smax,max(abss,abse));
// END; // for ks;
ks:=ks-1;
UNTIL (ks==0);
// convergence from satisfying threshold?
if (ks==(ke-1)) then
ke:=ks; continue;
end;
ks:=ks+1;
// is our block a 2x2 matrix?
if (ks==(ke-1)) then
svdBDQR2X2(ks,m,n);
ke:=ke-2;
continue;
end;
// block is larger than 2x2
// find chase direction (from big to small)
// if block disjoint from previous pass
if ((ks>keold) or (ke<ksold)) then
fdir:=abs(t_svdS[1,ks]) >= abs(t_svdS[1,ke]);
end;
// check convergence
if (fdir) then
t:=svdBDFC(ks,ke);
else
t:=svdBDBC(ks,ke);
end;
smin:=t[1];
if (t[2]) then continue; end;
ksold:=ks; keold:=ke;
// do QR, with shift if applicable
// if fudge*tol*(smin/smax)<=machine precision
// then use 0-shift; else compute shift
// tol usually larger than machine precision
// and fudge=min(m,n)
if ((r*(smin/smax))<=0.01) then
qrs:=0.0;
else
qrs:=svdQRShift(ks,ke,fdir);
end;
if (qrs) then
if (fdir) then
svdQRSF(ks,ke,qrs,m,n);
k:=ke-1;
else
svdQRSB(ks,ke,qrs,m,n);
k:=ks;
end;
else
if (fdir) then
t:=svdQRF(ks,ke,m,n);
else
t:=svdQRB(ks,ke,m,n);
end;
// qrs=0; use it as temp var
k:=t[4];
qrs:=t[1]*t_svdS[1,k];
t_svdS[1,k]:=t[2]*qrs;
k:=t[5];
t_svdS[2,k]:=t[3]*qrs;
end; // if qrs
if (abs(t_svdS[2,k])<=thresh) then
t_svdS[2,k]:=0;
end;
END; // while (ke>1)
END;
// convergence test (top left to bottom right)
// **********************************
svdBDFC(ks,ke):=
BEGIN
local k,t;
local svdtol1,svdtol2,smin;
svdtol1:=1.0e-13;
svdtol2:=3.0e-13;
if (abs(t_svdS[2,ke-1]) <= svdtol1*abs(t_svdS[1,ke])) then
t_svdS[2,ke-1]:=0;
return([0,1]);
end;
t:=abs(t_svdS[1,ks]);
smin:=t;
for k from ks to ke-1 do
if (abs(t_svdS[2,k]) <= svdtol2*t) then
t_svdS[2,k]:=0;
return([0,1]);
end;
t:=abs(t_svdS[1,k+1])*(t/(t+abs(t_svdS[2,k])));
smin:=min(smin,t);
end;
return([smin,0]);
END;
// convergence test (bottom right to top left)
// **********************************
svdBDBC(ks,ke):=
BEGIN
local k,t;
local svdtol1,svdtol2,smin;
svdtol1:=1.0e-13;
svdtol2:=3.0e-13;
if (abs(t_svdS[2,ks]) <= svdtol1*abs(t_svdS[1,ks])) then
t_svdS[2,ks]:=0;
return([0,1]);
end;
t:=abs(t_svdS[1,ke]);
smin:=t;
k:=ke-1;
REPEAT
// for k from ke-1 downto ks do
if (abs(t_svdS[2,k]) <= svdtol2*t) then
t_svdS[2,k]:=0;
return([0,1]);
end;
t:=abs(t_svdS[1,k])*(t/(t+abs(t_svdS[2,k])));
smin:=min(smin,t);
// end;
k:=k-1;
UNTIL (k==(ks-1));
return([smin,0]);
END;
// compute shift for QR iteration to
// speed up convergence
// **********************************
svdQRShift(ks,ke,fdir):=
BEGIN
local f,g,h,qrs,s;
local t,t1,t2,t3;
local fa,ga,ha,fhmin,fhmax;
local eps;
eps:=1.0e-14;
if (fdir) then
s:=abs(t_svdS[1,ks]);
f:=t_svdS[1,ke-1];
g:=t_svdS[2,ke-1];
h:=t_svdS[1,ke];
else
s:=abs(t_svdS[1,ke]);
f:=t_svdS[1,ks];
g:=t_svdS[2,ks];
h:=t_svdS[1,ks+1];
end;
fa:=abs(f); ga:=abs(g); ha:=abs(h);
fhmin:=min(fa,ha);
if (fhmin==0) then
qrs:=0.0;
else
fhmax:=max(fa,ha);
t1:=1+fhmin/fhmax;
t2:=(fhmax-fhmin)/fhmax;
t3:=ga/fhmax;
t3:= t3*t3;
t:=2/(sqrt(t1*t1+t3)+sqrt(t2*t2+t3));
qrs:=t*fhmin;
if (s > 0) then
t1:=qrs/s;
t1:=t1*t1;
if (t1<eps) then qrs:=0.0; end;
end;
end;
return(qrs);
END;
// 0-shift QR iteration (forward)
// **********************************
svdQRF(ks,ke,m,n):=
BEGIN
local csR, snR, csL, snL, r;
local k,t;
csR:=1.0; csL:=1.0;
snL:=0.0;
for k from ks to ke-1 do
t:=svdGetRotM(csR*t_svdS[1,k],t_svdS[2,k]);
csR:=t[1]; snR:=t[2]; r:=t[3];
if (k>ks) then
t_svdS[2,k-1]:=snL*r;
end;
t:=svdGetRotM(csL*r,snR*t_svdS[1,k+1]);
csL:=t[1]; snL:=t[2]; r:=t[3];
t_svdS[1,k]:=r;
svdRotUGT(csL,snL,k,m);
svdRotGV(csR,snR,k,n);
end;
return([csR,csL,snL,ke,ke-1]);
END;
// 0-shift QR iteration (backward)
// **********************************
svdQRB(ks,ke,m,n):=
BEGIN
local csR, snR, csL, snL, r;
local k,t;
csR:=1.0; csL:=1.0;
snL:=0.0;
k:=ke;
REPEAT
// for k from ke downto ks+1 do
t:=svdGetRotM(csR*t_svdS[1,k],t_svdS[2,k-1]);
csR:=t[1]; snR:=t[2]; r:=t[3];
if (k<ke) then
t_svdS[2,k]:=snL*r;
end;
t:=svdGetRotM(csL*r,snR*t_svdS[1,k-1]);
csL:=t[1]; snL:=t[2]; r:=t[3];
t_svdS[1,k]:=r;
svdRotUGT(csR,-snR,k-1,m);
svdRotGV(csL,-snL,k-1,n);
// end;
k:=k-1;
UNTIL (k==ks);
return([csR,csL,snL,ks,ks]);
END;
// non-zero shift QR iteration (forward)
// **********************************
svdQRSF(ks,ke,qrs,m,n):=
BEGIN
local csR, snR, csL, snL, r;
local k,t;
local f,g;
f:=t_svdS[1,ks];
f:=(abs(f)-qrs)*(svdsgn(f)+qrs/f);
g:=t_svdS[2,ks];
for k from ks to ke-1 do
t:=svdGetRotM(f,g);
csR:=t[1]; snR:=t[2]; r:=t[3];
if (k>ks) then
t_svdS[2,k-1]:=r;
end;
f:=csR*t_svdS[1,k]+snR*t_svdS[2,k];
t_svdS[2,k]:=csR*t_svdS[2,k]-snR*t_svdS[1,k];
g:=snR*t_svdS[1,k+1];
t_svdS[1,k+1]:=csR*t_svdS[1,k+1];
svdRotGV(csR,snR,k,n);
t:=svdGetRotM(f,g);
csL:=t[1]; snL:=t[2]; r:=t[3];
t_svdS[1,k]:=r;
f:=csL*t_svdS[2,k]+snL*t_svdS[1,k+1];
t_svdS[1,k+1]:=csL*t_svdS[1,k+1]-snL*t_svdS[2,k];
if (k==(ke-1)) then
t_svdS[2,k]:=f;
else
g:=snL*t_svdS[2,k+1];
t_svdS[2,k+1]:=csL*t_svdS[2,k+1];
end;
svdRotUGT(csL,snL,k,m);
end;
END;
// non-zero shift QR iteration (backward)
// **********************************
svdQRSB(ks,ke,qrs,m,n):=
BEGIN
local csR, snR, csL, snL, r;
local k,t;
local f,g;
f:=t_svdS[1,ke];
f:=(abs(f)-qrs)*(svdsgn(f)+qrs/f);
g:=t_svdS[2,ke-1];
k:=ke;
REPEAT
// for k from ke downto ks+1 do
t:=svdGetRotM(f,g);
csL:=t[1]; snL:=t[2]; r:=t[3];
if (k<ke) then
t_svdS[2,k]:=r;
end;
f:=csL*t_svdS[1,k]+snL*t_svdS[2,k-1];
t_svdS[2,k-1]:=csL*t_svdS[2,k-1]-snL*t_svdS[1,k];
g:=snL*t_svdS[1,k-1];
t_svdS[1,k-1]:=csL*t_svdS[1,k-1];
svdRotUGT(csL,-snL,k-1,m);
t:=svdGetRotM(f,g);
csR:=t[1]; snR:=t[2]; r:=t[3];
t_svdS[1,k]:=r;
f:=csR*t_svdS[2,k-1]+snR*t_svdS[1,k-1];
t_svdS[1,k-1]:=csR*t_svdS[1,k-1]-snR*t_svdS[2,k-1];
if (k==(ks+1)) then
t_svdS[2,k-1]:=f;
else
g:=snR*t_svdS[2,k-2];
t_svdS[2,k-2]:=csR*t_svdS[2,k-2];
end;
svdRotGV(csR,-snR,k-1,n);
// end;
k:=k-1;
UNTIL (k==ks);
END;
// SVD of 2x2 upper triangular matrix
//
// Bai and Demmel, Computing the generalized
// singular value decomposition
// SIAM J. Sci. Comput. Vol. 14, No. 6
// pp 1464-1486, Nov. 1993
// **********************************
svdBDQR2X2(k,r,c):=
BEGIN
local f,g,h,fa,ga,ha;
local pmax,swap,glarge,tsign;
local d,dd,q,qq,s,ss,tt;
local spq,dpq,a;
local tmp;
local csL,snL,csR,snR,smax,smin;
local csLt,snLt,csRt,snRt;
local eps;
eps:=2.1e-16;
pmax:=1; swap:=0; glarge:=0;
f:=t_svdS[1,k];
g:=t_svdS[2,k];
h:=t_svdS[1,k+1];
fa:=abs(f); ha:=abs(h); ga:=abs(g);
if (fa<ha) then
pmax:=3;
tmp:=fa; fa:=ha; ha:=tmp;
tmp:=f; f:=h; h:=tmp;
swap:=1;
end;
if (ga==0) then
smin:=ha;
smax:=fa;
csLt:=1.0; csRt:=1.0;
snLt:=0.0; snRt:=0.0;
else
if (ga>fa) then
pmax:=2;
tmp:=fa/ga;
if (tmp<eps) then
glarge:=1;
smax:=ga;
// compute (fa*ha)/ga
if (ha>1) then
tmp:=ga/ha;
smin:=fa/tmp;
else
smin:=tmp*ha;
end;
csLt:=1.0; snLt:=h/g;
snRt:=1.0; csRt:=f/g;
end; // (fa/ga < eps)
end; // if-else (ga<=fa)
if (glarge==0) then
tmp:=fa-ha;
if (tmp==fa) then
d:=1.0;
else
d:=tmp/fa;
end;
q:=g/f;
s:=2.0-d;
qq:=q*q; ss:=s*s;
spq:=sqrt(ss+qq);
if (d==0) then
dpq:=abs(q);
else
dpq:=sqrt(d*d+qq);
end;
a:=0.5*(spq+dpq);
smin:=ha/a;
smax:=fa*a;
if (qq==0) then
if (d==0) then
tmp:=svdsgn(f)*2*svdsgn(g);
else
tmp:=g/(svdsgn(f)*tmp)+q/s;
end;
else
tmp:=(q/(spq+s) + q/(dpq+d))*(1.0+a);
end;
tt:=sqrt(tmp*tmp+4.0);
csRt:=2.0/tt; snRt:=tmp/tt;
csLt:=(csRt + snRt*q)/a;
snLt:=(h/f)*snRt/a;
end; // (glarge==0)
end; // else (ga<>0)
if swap then
tmp:=h; h:=f; f:=tmp;
csL:=snRt; snL:=csRt;
csR:=snLt; snR:=csLt;
else
csL:=csLt; snL:=snLt;
csR:=csRt; snR:=snRt;
end; // if swap
// adjust sign of smax and smin
if (pmax==1) then
tsign:=svdsgn(csR)*svdsgn(csL)*svdsgn(f);
else
if (pmax==2) then
tsign:=svdsgn(snR)*svdsgn(csL)*svdsgn(g);
else
tsign:=svdsgn(snR)*svdsgn(snL)*svdsgn(h);
end;
end;
smax:=tsign*smax;
smin:=tsign*svdsgn(f)*svdsgn(h)*smin;
t_svdS[1,k+1]:=smin;
t_svdS[1,k]:=smax;
t_svdS[2,k]:=0;
svdRotUGT(csL,snL,k,r);
svdRotGV(csR,snR,k,c);
END;
// sort singular values from largest to
// to smallest; also make them all positive
// **********************************
svdPerm(r,m,n):=
BEGIN
local j,k,rs,s;
// make all singular values positive
// and fix corresponding vector in V
for k from 1 to r do
if (t_svdS[1,k]<0) then
t_svdS[1,k]:=-t_svdS[1,k];
for j from 1 to n do
t_svdV[k,j]:=-t_svdV[k,j];
end;
end;
end;
// now sort from largest to smallest
if (r==1) then return(0); end;
k:=r;
REPEAT
// for k from r downto 2 do
rs:=1;
s:=t_svdS[1,1];
for j from 2 to k do
if (t_svdS[1,j]<=s) then
rs:=j; s:=t_svdS[1,j];
end;
end;
// need new method (?); this may be
// costly on memory
if (rs<k) then
t_svdS[1,rs]:=t_svdS[1,k];
t_svdS[1,k]:=s;
t_svdU:=swapcol(t_svdU,rs,k);
t_svdV:=swaprow(t_svdV,rs,k);
end;
// end;
k:=k-1;
UNTIL (k==1);
END;
// compute the sign of a real number
// **********************************
svdsgn(r):=
BEGIN
if (r<0) then
return(-1);
else
return(1);
end;
END;
// matrix multiplication that does not
// convert small elements to 0
svdmmult(m1,m2):=
begin
local d;
local r1,c1,r2,c2;
local j,k;
local v1, v2;
local p;
d:=dim(m1);
r1:=d[1]; c1:=d[2];
d:=dim(m2);
r2:=d[1]; c2:=d[2];
if (c1 <> r2) then
return("Invalid dimension");
end;
p:=makemat(0,r1,c2);
for j from 1 to r1 do
for k from 1 to c2 do
v1:=m1[j];
v2:=subMat(m2,1,k,r2,k);
v2:=transpose(v2);
v2:=v2(1);
p[j,k]:=v1*v2;
end;
end;
return(p);
end;
// create matrix from diagonal
// v = [ v1, v2, ..., vn ]
// r = num of rows
// c = num of columns
// **********************************
diag2(v,r,c):=
BEGIN
local tS,d,k,n;
d:=min(r,c);
tS:=makemat(0,r,c);
n:=size(v);
for k from 1 to d do
if (k>n) then
tS[k,k]:=0;
else
tS[k,k]:=v[k];
end;
end;
return(tS);
END;
// compute pseudo-inverse
pinv(M,eps):=
BEGIN
local k,tmp;
local u,s,vt;
local p;
tmp:=svd2(M);
u:=tmp(1); s:=tmp(2); vt:=tmp(3);
tmp:=size(s);
for k from 1 to tmp do
if (s[k]<eps) then
s[k]:=0;
else
s[k]:=1/s[k];
end;
end;
tmp:=dim(M);
p:=svdmmult(trn(vt),diag2(s,tmp[2],tmp[1]));
p:=svdmmult(p,trn(u));
return(p);
END;
// QR factorization with pivoting
// A*P = Q*R
qrp(M):=
BEGIN
local d,m,n,norms,v;
local a,j,k,c,r;
local matP;
local matR;
t_svdA:=approx(M);
d:=dim(M);
m:=d(1); n:=d(2);
r:=min(m,n);
t_svdU:=idenmat(m);
matP:=idenmat(n);
norms:=makelist(0,1,n);
for k from 1 to n do
v:=subMat(t_svdA,1,k,m,k);
norms[k]:=dot(v,v);
end;
for k from 1 to r do
// find largest column
a:=norms[k]; c:=k;
for j from k+1 to n do
if (norms[j]>a) then
a:=norms[j];
c:=j;
end;
end;
// swap columns if necessary
if (c>k) then
t_svdA:=swapcol(t_svdA,k,c);
matP:=swapcol(matP,k,c);
a:=norms[k];
norms[k]:=norms[c];
norms[c]:=a;
end;
// apply Householder reflection
t_svdA[k,k]:=svdUQHQA(k,k,m,n);
// adjust norms for next loop
for j from k+1 to n do
a:=t_svdA[k,j]; a:=a*conj(a);
norms[j]:=norms[j]-a;
end;
end;
matR:=makemat(0,m,n);
for k from 1 to r do
for j from k to n do
matR[k,j]:=t_svdA[k,j];
end;
end;
return(t_svdU,matR,matP);
END;
// LQ factorization with pivoting
// P*A = L*Q
lqp(M):=
BEGIN
local d,m,n,norms,v;
local a,j,k,c,r;
local matP;
local matL;
t_svdA:=approx(M);
d:=dim(M);
m:=d(1); n:=d(2);
r:=min(m,n);
t_svdV:=idenmat(n);
matP:=idenmat(m);
norms:=makelist(0,1,m);
for k from 1 to m do
v:=subMat(t_svdA,k,1,k,n);
norms[k]:=dot(v,v);
end;
for k from 1 to r do
// find largest row
a:=norms[k]; c:=k;
for j from k+1 to m do
if (norms[j]>a) then
a:=norms[j];
c:=j;
end;
end;
// swap columns if necessary
if (c>k) then
t_svdA:=swaprow(t_svdA,k,c);
matP:=swaprow(matP,k,c);
a:=norms[k];
norms[k]:=norms[c];
norms[c]:=a;
end;
// apply Householder reflection
t_svdA[k,k]:=svdAPPHV(k,k,m,n);
// adjust norms for next loop
for j from k+1 to m do
a:=t_svdA[j,k]; a:=a*conj(a);
norms[j]:=norms[j]-a;
end;
end;
matL:=makemat(0,m,n);
for k from 1 to r do
for j from k to m do
matL[j,k]:=t_svdA[j,k];
end;
end;
return(matL,t_svdV,matP);
END;
#end
Graph 3D | QPI | SolveSys
|
Search
Member List
Calendar
Help
/
SVD - Singular Value Decomposition 14-FEB-2017
