SVD - Singular Value Decomposition 14-FEB-2017
10-20-2015, 05:26 PM (This post was last modified: 02-15-2017 04:56 AM by Han.)
Post: #1
 Han Senior Member Posts: 1,843 Joined: Dec 2013
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.

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
10-21-2015, 05:09 PM
Post: #2
 Han Senior Member Posts: 1,843 Joined: Dec 2013
RE: SVD - Singular Value Decomposition
Code updated to fix major bug (bad indexing) which caused backward sweeps to fail.

Graph 3D | QPI | SolveSys
10-21-2015, 06:20 PM
Post: #3
 Werner Senior Member Posts: 561 Joined: Dec 2013
RE: SVD - Singular Value Decomposition
Looks like a transcript of the LAPACK GESVD to me, nice job.
The 48/49 use the older LINPACK version (but without the shift calculation bug) that does not have the zero-shift QR nor bulge chasing in two directions, and different (simpler - too simple in some cases) convergence criteria.
I have the sources - no need for Jazz.

Cheers, Werner
10-27-2015, 07:07 PM (This post was last modified: 10-28-2015 04:41 AM by Han.)
Post: #4
 Han Senior Member Posts: 1,843 Joined: Dec 2013
RE: SVD - Singular Value Decomposition
Updated to include:

pinv(A,eps) -- calculates the pseudo-inverse of a matrix using the SVD
qrp(A) -- computes the QR factorization of A with pivoting, so that AP = QR, where P is a permutation

Edit: fixed a bug in the dimension of the P matrix for qrp(A)

Graph 3D | QPI | SolveSys
10-28-2015, 09:45 AM
Post: #5
 salvomic Senior Member Posts: 1,392 Joined: Jan 2015
RE: SVD - Singular Value Decomposition
thank you Han, very interesting!

I'm trying some examples and with them it seems to be very good

Salvo

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
11-09-2015, 03:36 PM
Post: #6
 Han Senior Member Posts: 1,843 Joined: Dec 2013
RE: SVD - Singular Value Decomposition
updated 9-nov-2015 to add lqp, which gives the LQ factorization of a matrix with pivoting (see first thread for info on the commands)

Graph 3D | QPI | SolveSys
02-08-2016, 03:13 PM
Post: #7
 Werner Senior Member Posts: 561 Joined: Dec 2013
RE: SVD - Singular Value Decomposition
(10-21-2015 06:20 PM)Werner Wrote:  The 48/49 use the older LINPACK version (but without the shift calculation bug) that does not have the zero-shift QR nor bulge chasing in two directions, and different (simpler - too simple in some cases) convergence criteria.

Well there goes my memory.
The 48/49/50 use the LAPACK version, as well.
How can you know?
1. look at the source code, or:
2. calculate the SVD of
Code:
   [ a^2 1  0  0  ]   [  0  1  a  0  ]   [  0  0  1  1  ]   [  0  0  0 a^2 ]
with a very small a so that 1+a = 1. (I used a=10^-20)
When the smallest singular value is about a^3, the LAPACK code is used.
When it's a^2/SQRT(2), the LINPACK code is used.
(the LINPACK code would set the offdiagonal a to zero since the convergence test used is 1 + 1 + a = 1 + 1 (adding the neighbouring diagonal elements)

Cheers, Werner
02-03-2017, 07:21 AM (This post was last modified: 02-03-2017 02:37 PM by Han.)
Post: #8
 Han Senior Member Posts: 1,843 Joined: Dec 2013
RE: SVD - Singular Value Decomposition
Hopefully squashed the last bug (f and h were not swapped within 2X2 routine; this lead to sign problems)... updates in first post.

Graph 3D | QPI | SolveSys
02-24-2020, 04:52 PM
Post: #9
 mark4flies Member Posts: 112 Joined: Jan 2015
RE: SVD - Singular Value Decomposition 14-FEB-2017
I love this program. I am a statistician and use the SVD all the time. I am no longer able to get this program to run since I installed recent firmware updates. Is there a new version that is compatible? Is there another way to make it work on the new firmware?

Thanks!
07-24-2020, 04:05 AM
Post: #10
 Han Senior Member Posts: 1,843 Joined: Dec 2013
RE: SVD - Singular Value Decomposition 14-FEB-2017
(02-24-2020 04:52 PM)mark4flies Wrote:  I love this program. I am a statistician and use the SVD all the time. I am no longer able to get this program to run since I installed recent firmware updates. Is there a new version that is compatible? Is there another way to make it work on the new firmware?

Thanks!

Does opening and closing the program's source file make any difference? Which firmware are you running? (Sorry for the late reply.)

Graph 3D | QPI | SolveSys
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)