HP Forums

Full Version: SVD - Singular Value Decomposition 14-FEB-2017
You're currently viewing a stripped down version of our content. View the full version with proper formatting.
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
Code updated to fix major bug (bad indexing) which caused backward sweeps to fail.
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
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)
thank you Han, very interesting!

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

Salvo
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)
(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
Hopefully squashed the last bug (f and h were not swapped within 2X2 routine; this lead to sign problems)... updates in first post.
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!
(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.)
Reference URL's