Main Page   Class Hierarchy   Compound List   File List   Compound Members   File Members

# Factor.cc

Go to the documentation of this file.
```00001 /*
00002     File:       Factor.cc
00003
00004     Function:   Implements some matrix factorizations: SVD, QR
00005
00006     Author:     Andrew Willmott
00007
00008     Notes:      Hopefully this is a bit more understandable
00009                 than most of the other SVD code out there.
00010
00011                 This should also make it a bit easier to optimize
00012                 for modern compilers, but that's a task for another
00013                 day.
00014 */
00015
00016 #include "vl/Factor.h"
00017
00018 #define DBG_COUT if (0) cerr
00019
00020 /*
00021     NOTE
00022
00023     Householder transformations zero all elements of a vector v apart
00024     from the first element.
00025
00026         H = I - 2uut
00027             where ||u||^2 = 1
00028             HtH = HHt = HH = I
00029             - it's a reflection, so we would expect HH = I
00030
00031         Hv = [c 0 0 0 0]
00032             where c = +/- ||v||
00033             and u = f(v - [c 0 0..])
00034             where f = 1/sqrt(2c(c - v[0]))
00035
00036         For convenience, we can write
00037             Hv = v - 2u(u.v)
00038         or
00039             Hv = v - 2w(w.v)/g
00040             where g = c(c - v[0])
00041             w = v - [c 0 0..]
00042 */
00043
00044 //#define DO_VL_SUB
00045
00046 #ifdef DO_VL_SUB
00047 // experiment to see just how much slower using the generic sub-matrix
00048 // calls makes things.
00049
00050 Double LeftHouseholder(Matd &A, Matd &U, const Int i)
00051 // Zeroes out those elements below the diagonal in column i by use of a
00052 // Householder transformation matrix H: A' = HA. U is replaced with UH,
00053 // so that U'A' = UH HA = UA
00054 {
00055     Assert(i < A.Rows(), "bad i");
00056
00057     Int         j, k;
00058     Double      scale, vSqrLen, oldDiag, newDiag, g, dotProd;
00059     Int         clen = A.Rows() - i;
00060     SubVecd     Ai = sub(col(A, i), i, clen);
00061
00062     // scale factor is the sum of abs of elements below and including
00063     // the diagonal in column i
00064     scale = 0.0;
00065     for (k = 0; k < Ai.Elts(); k++)
00066         scale += abs(Ai[k]);
00067
00068     // return if nothing to eliminate...
00069     if (scale == 0.0)
00070         return(0.0);
00071
00072     // scale the column, find sum of squares
00073     Ai = Ai / scale;
00074     vSqrLen = sqrlen(Ai);
00075
00076     oldDiag = Ai[0];            // v[0]
00077     newDiag = sqrt(vSqrLen);    // c
00078     if (oldDiag > 0.0)
00079         newDiag = -newDiag;
00080
00081     g = vSqrLen - oldDiag * newDiag;    // c(c - v0)
00082
00083     // replace column i of A with w
00084     Ai[0] = oldDiag - newDiag;      // v[0] -= c
00085
00086     // Apply H to A by transforming remaining columns of A
00087     for (j = i + 1; j < A.Cols(); j++)
00088     {
00089         SubVecd     Aj = sub(col(A, j), i, clen);
00090
00091         // Hv = v - w(v.w) / g
00092         Aj = Aj - Ai * (dot(Ai, Aj) / g);
00093     }
00094
00095     // Apply H to rows of U
00096     for (j = 0; j < A.Rows(); j++)
00097     {
00098         SubVecd     Uj = sub(row(U, j), i, clen);
00099
00100         // vH = v - (v.w)w / g
00101         Uj = Uj - Ai * (dot(Ai, Uj) / g);
00102     }
00103
00104     return(newDiag * scale);
00105 }
00106
00107 #else
00108
00109 Double LeftHouseholder(Matd &A, Matd &U, const Int i)
00110 // Zeroes out those elements below the diagonal in column i by use of a
00111 // Householder transformation matrix H: A' = HA. U is replaced with UH,
00112 // so that U'A' = UH HA = UA
00113 {
00114     Int     j, k;
00115     Double  scale, vSqrLen, oldDiag, newDiag, g, dotProd;
00116
00117     Assert(i < A.Rows(), "bad i");
00118
00119     // scale factor is the sum of abs of elements below and including
00120     // the diagonal in column i
00121     scale = 0.0;
00122     for (k = i; k < A.Rows(); k++)
00123         scale += abs(A[k][i]);
00124
00125     // return if nothing to eliminate...
00126     if (scale == 0.0)
00127         return(0.0);
00128
00129     // scale the column, find sum of squares
00130     vSqrLen = 0.0;
00131     for (k = i; k < A.Rows(); k++)
00132     {
00133         A[k][i] /= scale;
00134         vSqrLen += sqr(A[k][i]);
00135     }
00136
00137     oldDiag = A[i][i];          // v[0]
00138     newDiag = sqrt(vSqrLen);    // c
00139     if (oldDiag > 0.0)
00140         newDiag = -newDiag;
00141
00142     g = vSqrLen - oldDiag * newDiag;    // c(c - v0)
00143
00144     // replace column i of A with w
00145     A[i][i] = oldDiag - newDiag;        // v[0] -= c
00146
00147     // Apply H to A by transforming remaining columns of A
00148     for (j = i + 1; j < A.Cols(); j++)
00149     {
00150         // Hv = v - w(v.w) / g
00151
00152         // find dot product of the columns [i:i..m] (w) and [j:i..m]
00153         dotProd = 0.0;
00154         for (k = i; k < A.Rows(); k++)
00155             dotProd += A[k][j] * A[k][i];
00156
00157         dotProd /= g;
00158
00159         // calculate A's new column j
00160         for (k = i; k < A.Rows(); k++)
00161             A[k][j] -= A[k][i] * dotProd;
00162     }
00163
00164     // Apply H to rows of U
00165     for (j = 0; j < A.Rows(); j++)
00166     {
00167         // vH = v - (v.w)w / g
00168
00169         dotProd = 0.0;
00170         for (k = i; k < A.Rows(); k++)
00171             dotProd += A[k][i] * U[j][k];
00172
00173         dotProd /= g;
00174
00175         for (k = i; k < A.Rows(); k++)
00176             U[j][k] -=  A[k][i] * dotProd;
00177     }
00178
00179     return(newDiag * scale);
00180 }
00181
00182 #endif
00183
00184 Double RightHouseholder(Matd &A, Matd &V, const Int i)
00185 // Zeroes out those elements to the right of the super-diagonal in row i
00186 // by use of a Householder transformation matrix H: A' = AH. V is
00187 // replaced with VH, so that A'V't = AH (HV)t = AVt
00188 {
00189     Int     j, k;
00190     Double  scale, vSqrLen, oldSuperDiag, newSuperDiag, g, dotProd;
00191
00192     Assert(i < A.Cols() - 1, "bad i");
00193
00194     // scale factor is the sum of abs of elements to the right of the
00195     // diagonal in row i
00196     scale = 0.0;
00197     for (k = i + 1; k < A.Cols(); k++)
00198         scale += abs(A[i][k]);
00199
00200     // return if nothing to eliminate...
00201     if (scale == 0.0)
00202         return(0.0);
00203
00204     vSqrLen = 0.0;
00205     for (k = i + 1; k < A.Cols(); k++)
00206     {
00207         A[i][k] /= scale;
00208         vSqrLen += sqr(A[i][k]);
00209     }
00210
00211     oldSuperDiag = A[i][i + 1];
00212     newSuperDiag = sqrt(vSqrLen);
00213     if (oldSuperDiag > 0.0)
00214         newSuperDiag = -newSuperDiag;
00215
00216     g = oldSuperDiag * newSuperDiag - vSqrLen;
00217     A[i][i + 1] = oldSuperDiag - newSuperDiag;
00218
00219     // transform the remaining rows below i
00220     for (j = i + 1; j < A.Rows(); j++)
00221     {
00222         dotProd = 0.0;
00223         for (k = i + 1; k < A.Cols(); k++)
00224             dotProd += A[i][k] * A[j][k];
00225
00226         dotProd /= g;
00227
00228         for (k = i + 1; k < A.Cols(); k++)
00229             A[j][k] += A[i][k] * dotProd;
00230     }
00231
00232     // Accumulate the transform in V
00233     for (j = 1; j < A.Cols(); j++)
00234     {
00235         dotProd = 0.0;
00236         for (k = i + 1; k < A.Cols(); k++)
00237             dotProd += A[i][k] * V[j][k];
00238
00239         dotProd /= g;
00240
00241         for (k = i + 1; k < A.Cols(); k++)
00242             V[j][k] += A[i][k] * dotProd;
00243     }
00244
00245     // return new super-diagonal element A[i][i+1]
00246     return(newSuperDiag * scale);
00247 }
00248
00249 Void Bidiagonalize(Matd &A, Matd &U, Matd &V, Vecd &diagonal, Vecd &superDiag)
00250 // bidiagonalize matrix A by using householder transformations to eliminate
00251 // columns below the diagonal and rows to the right of the super-diagonal.
00252 // A is modified, and the diagonal and superDiag set.
00253 {
00254     Int     i;
00255     Matd    Us;
00256
00257     Assert(A.Rows() >= A.Cols(), "matrix must have rows >= cols");
00258
00259     diagonal.SetSize(A.Cols());
00260     superDiag.SetSize(A.Cols() - 1);
00261     U.SetSize(A.Rows(), A.Cols());
00262     Us.SetSize(A.Rows(), A.Rows());
00263     V.SetSize(A.Cols(), A.Cols());
00264     Us = vl_I;
00265     V = vl_I;
00266
00267     for (i = 0; i < A.Cols(); i++)
00268     {
00269         diagonal[i] = LeftHouseholder(A, Us, i);
00270
00271         if (i < A.Cols() - 1)
00272             superDiag[i] = RightHouseholder(A, V, i);
00273     }
00274
00275     U = sub(Us, 0, 0, A.Rows(), A.Cols());
00276 }
00277
00278 Double QRFactorization(Matd &A, Matd &Q, Matd &R)
00279 // Factor A into an orthogonal matrix Q and an upper-triangular matrix R.
00280 // Destroys A.
00281 {
00282     Double  normAcc = 0.0, diagElt;
00283     Int     i, j;
00284     Matd    Qs;
00285
00286     Assert(A.Rows() >= A.Cols(), "matrix must have rows >= cols");
00287
00288     Qs.SetSize(A.Rows(), A.Rows());
00289     Q.SetSize(A.Rows(), A.Cols());
00290     R.SetSize(A.Cols(), A.Cols());
00291     Qs = vl_I;
00292     R = vl_0;
00293
00294     // for each column
00295     for (i = 0; i < A.Cols(); i++)
00296     {
00297         // apply householder
00298         diagElt = LeftHouseholder(A, Qs, i);
00299         // copy over row i of A to R
00300         R[i][i] = diagElt;
00301         j = A.Cols() - i - 1;
00302         if (j)
00303             last(R[i], j) = last(A[i], j);
00304         normAcc = Max(normAcc, abs(diagElt));
00305     }
00306
00307     Q = sub(Qs, 0, 0, A.Rows(), A.Cols());
00308
00309     return(normAcc);
00310 }
00311
00312 /*
00313     NOTE
00314
00315     Givens rotations perform a rotation in a 2D plane. Can be
00316     used to zero an entry of a matrix.
00317
00318     We pick axes i and j, then form G such that it is basically
00319     I with entries:
00320
00321         i    j
00322     i:  c   -s
00323          ...
00324     j:  s    c
00325
00326
00327     s^2 + c^2 = 1
00328     inv(G) = Gt
00329
00330     A' = GA
00331         modifies rows i and j only:
00332             row'[i] = c * row[i] - s * row[j]
00333             row'[j] = c * row[j] + s * row[i]
00334         can force:
00335         A'[i][j] = 0
00336             if c * A[i][j] - s * A[j][j] = 0
00337             if c = A[j][j], s = A[i][j]
00338             must normalise to retain identity:
00339                 norm = sqrt(c^2 + s^2), c /= norm, s /= norm
00340             A'[j][j] = c * (norm * c) + s * (norm * s) = norm
00341         A'[i][i] = 0
00342             if c * A[i][i] - s * A[j][i] = 0
00343             if c = A[j][i], s = A[i][i]
00344
00345     A' = AGt
00346
00347
00348 */
00349
00350 Void RotateRight(Matd& U, const Int i, const Int j, const Double c, const Double s)
00351 // rotate U by the given Givens rotation: U' = UGt
00352 // where G is defined as above
00353 {
00354     Vecd    temp;
00355
00356     temp = col(U, i);
00357     col(U, i) =  c * col(U, i) - s * col(U, j);
00358     col(U, j) =  c * col(U, j) + s * temp;
00359 }
00360
00361 Void ClearSuperEntry(Vecd &diagonal, Vecd &superDiag, Matd &U,
00362             const Int k, const Int l, const Double prec)
00363 // We have a zero diagonal element at diag[l]. This means we can zero
00364 // out the super-diagonal entry to its right with a series of rotations.
00365 // Each rotation clears an entry (starting with the original
00366 // super-diagonal entry) but creates another entry immediately to its
00367 // right which must in turn be zeroed. Eventually we run off the end of
00368 // the matrix, and the process terminates.
00369 {
00370     Double  c = 0.0, s = 1.0;
00371     Double  f;
00372     Int     i;
00373     Double  norm;
00374
00375     // We zero the superdiagonal element l by using the row immediately
00376     // below it in a Givens rotation. Unfortunately, the superdiagonal
00377     // in this row in turn creates another entry A[l][l+2]. We must keep
00378     // applying rotations in the plane of the axes l and l + n to keep
00379     // zeroing each new entry created until we reach the right hand side
00380     // of the matrix.
00381
00382     // initialise: start with f being the original super diagonal entry
00383     // we're eliminating
00384     f = superDiag[l];
00385     superDiag[l] = 0.0;
00386
00387     // at each stage, f = A[l][i], rotate in the l/i plane
00388     for (i = l + 1; true; i++)
00389     {
00390         if (abs(f) + prec == prec)  // if f is zero, we don't have to eliminate further
00391             break;
00392
00393         // to eliminate f at each stage we pick s = -f, c = di
00394         s = -f;
00395         c = diagonal[i];
00396
00397         // normalise
00398         norm = sqrt(sqr(c) + sqr(s));
00399         c /= norm;
00400         s /= norm;
00401
00402         // apply inverse rotation to U
00403         RotateRight(U, i, l, c, s);
00404
00405         // do the rotation, zeroing f and creating f' immediately to its right
00406         diagonal[i] = norm;     // di' = c * di - s * f = norm
00407         if (i == k)             // done?
00408             break;
00409         f = s * superDiag[i];   // f'  = c * 0  + s * superdiag[i]
00410         superDiag[i] *= c;      // ei' = c * ei - s * 0
00411     }
00412 }
00413
00414 Int FindSplit(Vecd &diagonal, Vecd &superDiag, Matd &U, const Int k, const Double prec)
00415 // Check for a zero entry along the superdiagonal; if we find one, we can
00416 // QR iterate on two separate matrices above and below it.
00417 // If there is a zero on the *diagonal*, we can call ClearSuperEntry to
00418 // zero out the corresponding superdiagonal entry to its right.
00419 {
00420     Int     l;
00421
00422     for (l = k - 1; l >= 0; l--)
00423         if (superDiag[l] + prec == prec)
00424             // can split here
00425             return(l + 1);
00426         else if (diagonal[l] + prec == prec)
00427         {
00428             // can create a split here
00429             DBG_COUT << "creating split at " << l << endl;
00430             DBG_COUT << "diagonal " << diagonal << endl;
00431             DBG_COUT << "superDiag " << superDiag << endl;
00432             ClearSuperEntry(diagonal, superDiag, U, k, l, prec);
00433             DBG_COUT << "done: " << l << endl;
00434             DBG_COUT << "diagonal " << diagonal << endl;
00435             DBG_COUT << "superDiag " << superDiag << endl;
00436             DBG_COUT << endl;
00437             return(l + 1);
00438         }
00439
00440     return(0);
00441 }
00442
00443 Void Diagonalize(Vecd &diagonal, Vecd &superDiag, Matd &U, Matd &V)
00444 // Diagonalise the bidiagonal matrix A = (diagonal, superDiag), accumulating
00445 // transforms into U on the left and Vt on the right.
00446 //
00447 // diag(A) = diagonal and diag(A, 1) = superDiag, that is to say, diagonal[0]
00448 // is A[0][0], and superDiag[0] is A[0][1].
00449 {
00450     Int     i, j, k, l;
00451     Double  prec;
00452
00453     // work out a good precision value
00454     prec = abs(diagonal[0]);
00455     for (i = 1; i < diagonal.Elts(); i++)
00456         prec = Max(prec, abs(diagonal[i]) + abs(superDiag[i - 1]));
00457
00458     // work our way up from the bottom of the matrix
00459     for (k = diagonal.Elts() - 1; k >= 0; k--)
00460     {
00461         while (1)
00462         {
00463             // if there's a split, start from there rather than A[0][0]
00464             l = FindSplit(diagonal, superDiag, U, k, prec);
00465
00466             DBG_COUT << "QR-shift A " << l << ":" << k << endl;
00467
00468             DBG_COUT << "diagonal " << diagonal << endl;
00469             DBG_COUT << "superDiag " << superDiag << endl;
00470
00471             // are we done?
00472             if (l == k)     // last super-diag entry must be zero -- what we wanted.
00473                 break;
00474
00475             // QR iterate over the sub-matrix of A, A[l..k][l..k], until
00476             // we've nuked super diagonal entry A[k - 1][k]
00477
00478             // For extra stability, we shift the QR computation.  We
00479             // want the shift to be as close as possible to the largest
00480             // eigenvalue, which of course we don't know yet... so we
00481             // take the eigenvalues of the 2x2 matrix at the bottom of
00482             // our window of A, and use those.
00483
00484             Double  shift;
00485             Double  e0;
00486             Double  e1 = superDiag[k - 1];
00487             Double  d1 = diagonal[k - 1];
00488             Double  d2 = diagonal[k];
00489             Double  dl = diagonal[l];
00490
00491             if (k > 1)
00492                 e0 = superDiag[k - 2];
00493             else
00494                 e0 = 0.0;
00495
00496             //      d0  e0
00497             //          d1  e1
00498             //   k:         d2
00499
00500             shift = (d1 - d2) * (d1 + d2) + (e0 - e1) * (e0 + e1);
00501             shift /= 2.0 * e1 * d1;
00502             shift += sign(shift) * sqrt(sqr(shift) + 1.0);
00503             shift = ((dl - d2) * (dl + d2) + e1 * (d1 / shift - e1)) / dl;
00504
00505
00506             Double      cos_th, sin_th, cos_ph, sin_ph;
00507             Double      d1n, norm, e2, f1, g0, g1, e1n;
00508
00509             d1 = dl;    // d1 == A[i - 1][i - 1]: initialise it to that
00510             e1 = superDiag[l];
00511             // the first rotation is the QR-shifted one. After that, we
00512             // are continually eliminating below-diagonal and
00513             // above-super diagonal elements created by the previous
00514             // rotation, until we hit the bottom-right of the matrix and
00515             // we have a bidiagonal matrix again. The QR-shift will
00516             // ensure that the new matrix has smaller super-diagonal
00517             // elements, however.
00518             g0 = e1;
00519             e0 = shift;
00520
00521             for (i = l + 1; true; i++)
00522             {
00523                 // rotate in plane of axes i - 1, i:
00524                 //    d0  e0  g0      d0  e0'  0
00525                 //        d1  e1   ->     d1'  e1'
00526                 //  i:        d2          f1   d2'
00527                 // we are rotating on the right, and so affecting columns
00528                 d2 = diagonal[i];
00529
00530                 // eliminate g0 using e0 (except for the first iteration,
00531                 // where e0 and g0 would be off the top of the matrix, and
00532                 // we're performing the QR-shift.)
00533                 sin_th = g0;
00534                 cos_th = e0;
00535
00536                 norm = sqrt(sqr(cos_th) + sqr(sin_th));
00537                 cos_th /= norm;
00538                 sin_th /= norm;
00539
00540                 // perform the rotation
00541                 if (i > 1)
00542                     superDiag[i - 2] = norm;    // e0' = cos_th * e0 + sin_th * g0
00543                 d1n = cos_th * d1 + sin_th * e1;
00544                 e1n = -sin_th * d1 + cos_th * e1;
00545
00546                 f1 = d2 * sin_th;           // f1  = c *  0 + s * d2
00547                 d2 *= cos_th;               // d2' = c * d2 + s * 0
00548
00549                 // Accumulate rotation in V
00550                 RotateRight(V, i, i - 1, cos_th, sin_th);
00551
00552                 // in eliminating g0, we've created f1: eet must be
00553                 // destroyed
00554
00555                 // rotate in plane of axes i - 1, i:
00556                 //      d0  e0             d0   e0
00557                 //      d1  e1      ->     d1'  e1'  g1'
00558                 //  i:  f1  d2  e2         0    d2'  e2'
00559                 // we are rotating on the left, and so affecting rows
00560
00561                 // standard rotation to eliminate the f1 we've just created
00562                 cos_ph = d1n;
00563                 sin_ph = f1;
00564
00565                 norm = sqrt(sqr(cos_ph) + sqr(sin_ph));
00566
00567                 if (norm == 0.0)
00568                 {
00569                     // rotation angle is arbitrary
00570                     cos_ph = cos_th;
00571                     sin_ph = sin_th;
00572                 }
00573                 else
00574                 {
00575                     cos_ph /= norm;
00576                     sin_ph /= norm;
00577                 }
00578
00579                 // as usual, for the elimination pair, f1 = 0, d1 = norm
00580                 diagonal[i - 1] = norm;     // d1'
00581                 // rotate e1 and d2
00582                 e1 =  cos_ph * e1n + sin_ph * d2;
00583                 d2 = -sin_ph * e1n + cos_ph * d2;
00584
00585                 // Accumulate rotation into U
00586                 RotateRight(U, i, i - 1, cos_ph, sin_ph);
00587
00588                 if (i == k)
00589                     break;
00590
00591                 // rotate g1 and e2
00592                 e2 = superDiag[i];
00593                 g1 = sin_ph * e2;
00594                 e2 = cos_ph * e2;
00595
00596                 d1 = d2;        // the new d1 will be the old d2
00597                 e0 = e1;
00598                 e1 = e2;
00599                 g0 = g1;
00600             }
00601
00602             if (l > 0)
00603                 superDiag[l - 1] = 0.0;     // Supposed to be eliminated by now
00604             superDiag[k - 1] = e1;          // write in the last superdiagonal
00605             diagonal[k] = d2;               // and diagonal
00606         }
00607
00608         // Force singular value to be +ve if needs be
00609         if (diagonal[k] < 0)
00610         {
00611             DBG_COUT << "flipping " << k << endl;
00612             diagonal[k] = -diagonal[k];
00613             // flip the corresponding row of Vt to balance out
00614             col(V, k) = -col(V, k);
00615         }
00616         DBG_COUT << "done: " << endl;
00617         DBG_COUT << "diagonal " << diagonal << endl;
00618         DBG_COUT << "superDiag " << superDiag << endl;
00619         DBG_COUT << endl;
00620     }
00621 }
00622
00623 Void SVDFactorization(Matd &A, Matd &U, Matd &V, Vecd &diagonal)
00624 // pretty easy...
00625 {
00626     Vecd    superDiag;
00627
00628     Bidiagonalize(A, U, V, diagonal, superDiag);
00629     Diagonalize(diagonal, superDiag, U, V);
00630 }
```

Generated at Sat Aug 5 00:16:46 2000 for Class Library by 1.1.0 written by Dimitri van Heesch, © 1997-2000