Assignments Assignments

SVD.inl
Go to the documentation of this file.
1
2//
3// SVD.inl
4//
5// Abstract: Singular-Value Decomposition. Numerous informative comments
6// througout code.
7//
8// Revision History Abstract:
9// xxJul1989 Ronen Barzel Initial code
10// xxFeb1990 Ronen Barzel Bug fixes
11// xxMay1993 Ronen Barzel ANSI-C backsubst
12// 21Jan1997 ChuckR Move into namespace
13// 21Jul1997 ChuckR Fix array bound error in underdetermined case
14// 13Apr1999 ChuckR Templatize
16
17
18/*************************************************************
19** Singular-Value Decomposition.
20**
21** Had to write this by hand, since the version in Numerical
22** Recipes has a bug, carried over from a bug in Golub's
23** original algorithm:
24
25 Gene H. Golub & Charles F. Van Loan
26 MATRIX COMPUTATIONS
27 Johns Hopkins University Press, 1983
28 Fifth printing, 1987
29
30 Page 293
31 Algorithm 8.3-2: The SVD Algorithm
32
33 ...Find the largest q and the smallest p such that if
34
35 +- -+
36 | A11 0 0 | p
37 | 0 A22 0 | n-p-q
38 A = | 0 0 A33 | q
39 | 0 0 0 | m-n
40 +- -+
41
42 then A33 is diagonal and A22 has a nonzero superdiagonal.
43 If q = n then quit.
44 If any diagonal entry in A22 is zero, then zero the superdiagonal
45 element in the same row and go to Repeat...
46
47 The last sentence above is the cause of the trouble, in the following case:
48
49 +- -+
50 |0 1 0|
51 A = |0 1 1|
52 |0 0 0|
53 +- -+
54
55 In this case, q is 0, A33 is null, and A22 is the same as A. The instruction
56
57 "if any diagonal entry in A22 is zero, then
58 zero the superdiagonal element in the same row"
59
60 cannot be applied -- A22 has a diagonal entry that is zero, but there
61 is no superdiagonal element in that same row. The proper thing to
62 do seems to be to special-case: If A22 has 0 in its lower-right
63 diagonal element, zero the superdiagonal element above it.
64
65
66**
67**
68**
69** Since the Num. Rec. code was cribbed from EISPACK or LINPACK
70** or some such, which in turn cribbed from other places, I would
71** consider an SVD routine suspect unless it works on the following
72** test data:
73** a = 0 1 0
74** 0 1 1
75** 0 0 0
76**
77** Anyway, Al Barr & I went through the references to figure out
78** how it all worked. I reimplimented it from scratch, first in
79** lisp, and now in C. Martin L. Livesey used the code and reported
80** various bugs, which have been fixed.
81**
82**
83** public routines:
84** num_svd()
85** num_svd_backsubst()
86**
87** Ronen Barzel July 1989
88** (bug fixes: Feb 1990)
89** (ansi-C, backsubst: May 1993)
90**
91****************************************************************/
92
93#include <stdio.h>
94#include <memory.h>
95#ifdef __APPLE__
96#include <malloc/malloc.h>
97#else // !__APPLE__
98#include <malloc.h>
99#endif // __APPLE__
100#include <math.h>
101
102#define MIN(A,B) (((A)<(B)) ? (A):(B))
103#define MAX(A,B) (((A)>(B)) ? (A):(B))
104#define ALLOC2D(m,n) (REAL *) malloc((unsigned)(m)*(n)*sizeof(REAL))
105#define ALLOC1D(n) (REAL *) malloc((unsigned)(n)*sizeof(REAL))
106#define FREE(p) free((char*)(p))
107#define CLEAR2D(a,m,n) (void) memset((char*)(a),0,(m)*(n)*sizeof(REAL))
108#define COPY2D(d,m,n,s) (void) memcpy((char*)(d),(char*)(s),(int)((m)*(n)*sizeof(REAL)))
109#define REF2D(a,m,n,i,j) (a[(n)*(i)+(j)])
110
111 template <class REAL>
112REAL fhypot (REAL a, REAL b)
113{
114 return ((REAL) sqrt (a*a + b*b));
115}
116
117/* householder transformations
118**
119** looks at the submatrix of A below and to the right of the [i,j]th element
120** (inclusive) (where i and j go from 0 to n-1). Performs a householder
121** transformation on the submatrix, to zero out all but the first
122** elements of the first column.
123**
124** Matrix u (a rotation matrix) is transformed by the inverse of the householder
125** transformation, so that (old u)(old a) = (new u)(new a)
126**
127** a is m x n, u is m x m
128**
129** If the submatrix is (X1, X2, ...) the householder transformation is
130** P_X1, X1 transforms to (|X1|,0,0,0,...), and the resulting matrix
131** is ((|X1|,0,0,...),P_X1 X2, P_X1 X3, ...)
132*/
133
134
135 template <class REAL>
136 void householder_zero_col(REAL *a,REAL *u,int i,int j,int m,int n,REAL *hv)
137 /* hv is a work vector, length m */
138{
139 REAL scale, /* a scale factor for X1 */
140 sigma, /* the norm^2 of X1 */
141 snorm, /* +/- the norm of X1 */
142 f, h, s, dot;
143 int k,l;
144
145 /* we will scale X1 by its l1-norm. squawk! */
146 for (scale=0, k=i; k<m; k++) scale += ((REAL) fabs(REF2D(a,m,n,k,j)));
147
148 /* if X1 is 0, no point doing anything else */
149 if (!scale) return;
150
151 /* divide out the l1-norm, and calculate sigma = |X|^2 */
152 for (k=i; k<m; k++) hv[k] = REF2D(a,m,n,k,j) / scale;
153 for (sigma=0, k=i; k<m; k++) sigma += hv[k]*hv[k];
154
155 /* The householder vector is X1 +/- (|X1|,0,0,...). We will
156 ** contruct this vector in hv and use it to perform the householder
157 ** transformation on the matrix. The plus or minus on the norm
158 ** is to reduce roundoff error. */
159 f = hv[i];
160 snorm = ((REAL) ((f > 0) ? -sqrt(sigma) : sqrt(sigma)));
161 h = f*snorm - sigma; /* eqn 11.2.4 in Press et.al. */
162 hv[i] = f - snorm;
163
164 /* set the first column of a to be the (|X1|,0,...) -- this is
165 ** what we would get by performing the transformation, but this
166 ** way's faster */
167 REF2D(a,m,n,i,j) = scale * snorm;
168 for (k=i+1; k<m; k++) REF2D(a,m,n,k,j) = 0;
169
170 /* Now perform the householder transformation on the rest of
171 ** the columns. If the householder vector is X, and -half its norm^2
172 ** is h, the householder transform is P_ij : delta_ij + x_i x_j / h.
173 ** We don't actually create the P matrix, we just do the calcuation.
174 */
175 for (k=j+1; k<n; k++) {
176 for (dot=0, l=i; l<m; l++) dot += REF2D(a,m,n,l,k)*hv[l];
177 s = dot/h;
178 for (l=i; l<m; l++) REF2D(a,m,n,l,k) += s*hv[l];
179 }
180
181 /* Similarly, perform the householder transformation on (all)
182 ** the rows of u. Note that it's the rows of u rather than
183 ** the columns, because we want U to invert what we've done to a.
184 */
185 for (k=0; k<m; k++) {
186 for (dot=0, l=i; l<m; l++) dot += REF2D(u,m,m,k,l)*hv[l];
187 s = dot/h;
188 for (l=i; l<m; l++) REF2D(u,m,m,k,l) += s*hv[l];
189 }
190}
191
192/* this is the same as householder_zero_col, but with rows
193** and cols swapped.
194*/
195
196 template <class REAL>
197 void householder_zero_row(REAL *a,REAL *v,int i,int j,int m,int n,REAL *hv)
198{
199 REAL scale, sigma,snorm, f, h, s, dot;
200 int k,l;
201
202 for (scale=0, k=j; k<n; k++) scale += ((REAL) fabs(REF2D(a,m,n,i,k)));
203 if (!scale) return;
204
205 for (k=j; k<n; k++) hv[k] = REF2D(a,m,n,i,k) / scale;
206 for (sigma=0, k=j; k<n; k++) sigma += hv[k]*hv[k];
207
208 f = hv[j];
209 snorm = ((REAL) ((f > 0) ? -sqrt(sigma) : sqrt(sigma)));
210 h = f*snorm - sigma;
211 hv[j] = f - snorm;
212
213 REF2D(a,m,n,i,j) = scale * snorm;
214 for (k=j+1; k<n; k++) REF2D(a,m,n,i,k) = 0;
215
216 for (k=i+1; k<m; k++) {
217 for (dot=0, l=j; l<n; l++) dot += REF2D(a,m,n,k,l)*hv[l];
218 s = dot/h;
219 for (l=j; l<n; l++) REF2D(a,m,n,k,l) += s*hv[l];
220 }
221
222 for (k=0; k<n; k++) {
223 for (dot=0, l=j; l<n; l++) dot += REF2D(v,n,n,l,k)*hv[l];
224 s = dot/h;
225 for (l=j; l<n; l++) REF2D(v,n,n,l,k) += s*hv[l];
226 }
227}
228
229/*
230** performs a Givens rotation on the ith and jth columns of a, looking
231** at the n elements starting at start. a is mm x nn.
232*/
233 template <class REAL>
234 void rotate_cols(int i,int j,REAL cos,REAL sin,REAL* a,int start,int n,int mm,int nn)
235{
236 int end = start+n, k;
237 REAL x,y;
238 for (k=start; k<end; k++) {
239 x = REF2D(a,mm,nn,k,i);
240 y = REF2D(a,mm,nn,k,j);
241 REF2D(a,mm,nn,k,i) = cos*x + sin*y;
242 REF2D(a,mm,nn,k,j) = -sin*x + cos*y;
243 }
244}
245
246
247/*
248** performs a Givens rotation on the ith and jth rows of a, looking
249** at the n elements starting at start. a is mm x nn.
250*/
251 template <class REAL>
252 void rotate_rows(int i,int j, REAL cos,REAL sin,REAL* a,int start,int n,int mm,int nn)
253{
254 int end = start+n, k;
255 REAL x,y;
256 for (k=start; k<end; k++) {
257 x = REF2D(a,mm,nn,i,k);
258 y = REF2D(a,mm,nn,j,k);
259 REF2D(a,mm,nn,i,k) = cos*x + sin*y;
260 REF2D(a,mm,nn,j,k) = -sin*x + cos*y;
261 }
262}
263
264/*
265** This takes a submatrix of b (from p through z inclusive). The submatrix
266** must be bidiagonal, with a zero in the upper-left corner element. Mucks
267** with the sumatrix so that the top row is completely 0, accumulating
268** the rotations into u. b is m x n. u is m x min.
269**
270** Suppose the matrix looks like
271** 0 R 0 0 0...
272** 0 X X 0 0...
273** 0 0 X X 0...
274** ...
275** Where R & X's are different & non-zero. We can rotate the first
276** and second rows, giving
277** 0 0 R 0 0...
278** 0 X X 0 0...
279** 0 0 X X 0...
280** ...
281** with new R's and X's. We rotate the first and third rows, moving
282** R over again, etc. till R falls off the end.
283*/
284
285 template <class REAL>
286 void clr_top_supdiag_elt(REAL * b,int p,int z,REAL * u,int m,int n,int min)
287{
288 int i;
289 REAL r, x, hypot, cos, sin;
290
291 for (i=p+1; i<=z; i++) {
292 r = REF2D(b,m,n,p,i);
293 x = REF2D(b,m,n,i,i);
294 if (r==0) break;
295 hypot = ((REAL) sqrt (r*r + x*x)); //hypot = pythag(r,x);
296 cos = x/hypot;
297 sin = r/hypot;
298 /* update the diagonal and superdiagonal elements */
299 REF2D(b,m,n,p,i) = 0;
300 REF2D(b,m,n,i,i) = hypot;
301 /* Rotate the remainder of rows p and i (only need to
302 ** rotate one more element, since the rest are zero) */
303 if (i != z) rotate_rows(p,i,cos,sin,b,i+1,1,m,n);
304 /* accumulate the rotation into u */
305 rotate_cols(i,p,cos,sin,u,0,m,m,min);
306 }
307}
308
309/*
310** This takes a submatrix of b (from p through z inclusive). The submatrix
311** must be bidiagonal, with a zero in the lower-right corner element. Mucks
312** with the sumatrix so that the right column is completely 0, accumulating
313** the rotations into v. b is m x n. v is min x n
314**
315** Suppose the matrix looks like
316** X X 0 0
317** 0 X X 0
318** 0 0 X R
319** 0 0 0 0
320** Where R & X's are different & non-zero. We can rotate the last
321** and second-to-last columns, yielding
322** X X 0 0
323** 0 X X R
324** 0 0 X 0
325** 0 0 0 0
326** with new R's and X's. We rotate the last and third-to-last cols, moving
327** R over again, etc. till R falls off the end.
328*/
329
330 template <class REAL>
331 void clr_bot_supdiag_elt(REAL * b,int p,int z,REAL * v,int m,int n,int min)
332{
333 int i;
334 REAL r, x, hypot, cos, sin;
335
336 for (i=z-1; i>=p; i--) {
337 r = REF2D(b,m,n,i,z);
338 x = REF2D(b,m,n,i,i);
339 if (r==0) break;
340 hypot = ((REAL) sqrt (r*r + x*x)); //hypot = pythag(r,x);
341 cos = x/hypot;
342 sin = r/hypot;
343 /* update the diagonal and superdiagonal elements */
344 REF2D(b,m,n,i,z) = 0;
345 REF2D(b,m,n,i,i) = hypot;
346 /* Rotate the remainder of cols z and i (only need to
347 ** rotate one more element, since the rest are zero) */
348 if (i != p) rotate_cols(i,z,cos,sin,b,i-1,1,m,n);
349 /* accumulate the rotation into v */
350 rotate_rows(i,z,cos,sin,v,0,n,min,n);
351 }
352}
353
354/*
355** This takes a submatrix of b (from p through z inclusive). The submatrix
356** must be bidiagonal except that the topmost subdiagonal element is non-zero.
357** Mucks with the submatrix to make it bidiagonal, accumulating the rotations
358** into u and v. b is m x n u is m x min v is min x n
359**
360** Suppose the matrix looks like
361** X X 0 0 0...
362** R X X 0 0...
363** 0 0 X X 0...
364** ...
365** Where R & X's are different & non-zero. We can rotate the first and
366** second rows, giving
367** X X R 0 0...
368** 0 X X 0 0...
369** 0 0 X X 0...
370** ...
371** with new R &X's. Now rotate the second and third columns, getting
372** X X 0 0 0...
373** 0 X X 0 0...
374** 0 R X X 0...
375** ...
376** which is the same as the initial problem, but decreased in
377** size by 1. Eventually, we'll reach the situation where we have
378** ...
379** ... X X
380** ... R X
381** and the row rotation will eliminate R.
382*/
383
384 template <class REAL>
385 void clr_top_subdiag_elt(REAL * b,int p,int z,REAL * u,REAL * v,int m,int n,int min)
386{
387 int i;
388 REAL x, r, hypot, cos, sin;
389
390 for (i=p; ; i++) {
391 /* figure out row rotation to zero out the subdiagonal element */
392 x = REF2D(b,m,n,i,i);
393 r = REF2D(b,m,n,i+1,i);
394 hypot = fhypot(x,r);
395 cos = x/hypot;
396 sin = r/hypot;
397 /* rotate the leading elements of the row */
398 REF2D(b,m,n,i,i) = hypot;
399 REF2D(b,m,n,i+1,i) = 0;
400 /* rotate out the rest of the row */
401 rotate_rows(i,i+1,cos,sin,b,i+1,(i+1==z)?1:2,m,n);
402 /* accumulate transformation into columns of u */
403 rotate_cols(i,i+1,cos,sin,u,0,m,m,min);
404
405 /* end with a row rotation */
406 if (i+1==z) break;
407
408 /* figure out column rotation */
409 x = REF2D(b,m,n,i,i+1);
410 r = REF2D(b,m,n,i,i+2);
411 hypot = fhypot(x,r);
412 cos = x/hypot;
413 sin = r/hypot;
414 /* rotate the leading elements of the column */
415 REF2D(b,m,n,i,i+1) = hypot;
416 REF2D(b,m,n,i,i+2) = 0;
417 /* rotate the rest of the column */
418 rotate_cols(i+1,i+2,cos,sin,b,i+1,2,m,n);
419 /* accumulate transformation into columns of v */
420 rotate_rows(i+1,i+2,cos,sin,v,0,n,min,n);
421 }
422}
423
424/*
425** This is the first part of an implicit-shift QR step. We do some
426** magic eigenvalue calculation, to calculate a Givens rotation for
427** the top-left corner.
428**
429** This rotation is described as part of Golub's Algorithm 3.3-1
430**
431** This is also described as the implicit QL algorithm for symmetric
432** tridiagonal matrices, in Numerical Recipes. Here's what's going
433** on, as far as I can tell:
434**
435** Hypothetically, one could do a series of QR decompositions of a symmetric
436** tridiagonal matrix A.
437** - Ignoring the R's, one serially computes An+1 = Qn^t An Q
438** - Eventually, An will approach diagonal
439** - The rate of convergence goes like A[i,j] = u_i/u_j, where
440** u_i and u_j are the eigenvalues of A
441** - To make it converge faster, we shift the eigenvalues by some amount
442** ("uu" in the code) by decomposing the matrix A - uu I.
443** - uu is computed so as to attempt to minimize (u_i-uu)/(u_j-uu), which
444** is the convergence rate of the shifted matrix. Ideally uu would
445** be equal to the smallest eigenvalue. Since we don't know the
446** smallest eigenvalue, we look at the eigenvalues of the
447** trailing 2x2 submatrix.
448** - Rather than actually computing the matrix A - uu I, we just keep
449** track of the shift in the calculation of the coefficients for the
450** rotations that make up the Q's. Hence the "implicit" in the
451** name of the algorithm.
452**
453** For SVD, we are looking at a bidiagonal matrix, rather than a tridiagonal
454** one. Thus we will do one more level of implicitness, by computing the
455** coefficients we WOULD get were we to consider the tridiagonal matrix
456** T = B^t B.
457**
458** This particular routine just performs the initial rotation on the
459** bidiagonal matrix, based on the eigenvalue-shift stuff. This
460** makes the matrix no loger bidiagonal. Other routines will have to
461** clean up the non-bidiagonal terms. The net rotation that is performed
462** by this routine and the cleanup routines is the Q at one step of
463** the above iteration. We don't ever explicitly compute Q, we
464** just keep updating the U and V matrices.
465**
466** b is m x n, v is min x n
467*/
468
469 template <class REAL>
470 void golub_kahn_svd_rot(REAL * b,int p,int q,REAL * v,int m,int n,int min)
471{
472 REAL t11, t12, t22;
473 REAL uu;
474 REAL b11, b12;
475 REAL dm,fm,dn,fn;
476 REAL hypot,cos,sin;
477 REAL d, s, y, z;
478
479 /* grab the last diagonal and superdiagonal elements of b22 */
480 fm = (q-2) < 0 ? 0 : REF2D(b,m,n,q-2,q-1);
481 dm = REF2D(b,m,n,q-1,q-1); fn = REF2D(b,m,n,q-1,q);
482 dn = REF2D(b,m,n,q,q);
483
484 /* create the trailing 2x2 submatrix of T = b22^t b22 */
485 t11 = dm*dm + fm * fm;
486 t12 = dm * fn; t22 = dn * dn + fn * fn;
487
488 /* find the eigenvalues of the T matrix.
489 **
490 ** The quadratic equation for the eigenvalues has coefficients:
491 ** a = 1
492 ** b = -(t11+t22)
493 ** c = t11 t22 - t12 t12
494 **
495 ** b^2 - 4ac = (t11^2 + t22^2 + 2 t11 t12) - 4 t11 t22 + 4 t12 t12
496 ** = (t11 - t22)^2 + 4 t12 t12
497 **
498 ** sqrt(b^2-4ac) = sqrt((t11 - t22)^2 + 4 t12 t12) -- use "pythag()"
499 **
500 ** using quadratic formula, we have the eigenvalues:
501 ** (u1,u2) = .5 (t11 + t22 +/- pythag(...))
502 ** = t22 + .5 (t11 - t22 +/- pythag(...))
503 **
504 ** We propogate the .5 into the pythag term, yielding golub's equation 8.2-2
505 ** for the "wilkinson shift".
506 ** [Note: Golub says to use (t22 + d - signum(d) s) to find the eigenvalue
507 ** closer to t22. He's right. ]
508 **/
509 d = ((REAL) 0.5)*(t11 - t22);
510 s = fhypot (d,t12);
511 uu = t22 + d + ((d > 0) ? -s : s);
512
513 /* grab the leading 2x1 of b */
514 b11 = REF2D(b,m,n,p,p); b12 = REF2D(b,m,n,p,p+1);
515
516 /* make the leading 2x1 submatrix of T */
517 t11 = b11 * b11;
518 t12 = b11 * b12;
519
520 /* calculate the rotation that would zero the off-diagonal term of the
521 ** shifted matrix T - uu I
522 */
523 y = t11 - uu;
524 z = t12;
525 hypot = fhypot (y,z);
526 cos = y / hypot;
527 sin = z / hypot;
528
529 /* perform the rotation on B. This sprays some glop into the upper-left
530 ** subdiagonal element.
531 */
532 rotate_cols(p,p+1,cos,sin,b,p,2,m,n);
533 /* accumulate the rotation into the rows of v */
534 rotate_rows(p,p+1,cos,sin,v,0,n,min,n);
535}
536
537
538/* bidiagonalize
539**
540** Given a (m x n)
541** computes u (m x min) orthogonal cols
542** b (m x n) bidiagonal
543** v (min x n) orthogonal rows
544** such that a = u b v (looking at only the min x min leading part of b)
545**
546** Works by starting with B = A, then performing a series of Householder
547** transformations to zero-out the off-diagonal elements, while
548** accumulating the transformations into U and V.
549**
550** b may point to the same array as a, in which case the result
551** overwrites the original data.
552*/
553
554 template <class REAL>
555 void bidiagonalize(const REAL *a,int m,int n,REAL *u,REAL *b,REAL *v)
556{
557 int i,j, min=MIN(m,n);
558 REAL *usave = u, *vsave = v, *h = ALLOC1D(MAX(m,n));
559
560 /* we need square matrices to accumulate the transformations */
561 if (min != m) u = ALLOC2D(m,m);
562 if (min != n) v = ALLOC2D(n,n);
563
564 /* start off with u and v equal to the identity, and b equal to a */
565 CLEAR2D(u,m,m);
566 CLEAR2D(v,n,n);
567 for (i=0; i<m; i++) REF2D(u,m,m,i,i) = 1;
568 for (i=0; i<n; i++) REF2D(v,n,n,i,i) = 1;
569 if (b != a) COPY2D(b,m,n,a);
570
571 /* walk down the diagonal */
572 for (i = 0; i<min; i++) {
573 /* zero the entries below b[i,i] */
574 householder_zero_col(b,u,i,i,m,n,h);
575 /* zero the entries to the right of b[i,i+1] */
576 householder_zero_row(b,v,i,i+1,m,n,h);
577 }
578 /* when m < n (matrix wider than tall), the above
579 leaves a non-0 element in the [m-1,m]th spot.
580 This is the bottom superdiagonal element of an (m+1)x(m+1);
581 use the clear routine to get rid of it. */
582 if (min != n)
583 clr_bot_supdiag_elt(b,0,min,v,m,n,n);
584
585 /* For non-square arrays, lop off the parts we don't need */
586 if (min!=m) {
587 for (i=0; i<m; i++) for (j=0; j<min; j++)
588 REF2D(usave,m,min,i,j)=REF2D(u,m,m,i,j);
589 FREE(u);
590 }
591 if (min!=n) {
592 for (i=0; i<n; i++) for (j=0; j<min; j++)
593 REF2D(vsave,min,n,j,i)=REF2D(v,n,n,j,i);
594 FREE(v);
595 }
596 FREE(h);
597}
598
599/*
600** Finds the SVD of a bidiagonal matrix.
601**
602** Zeroes the superdiagonal of a bidiagonal matrix, accumulating
603** left- and right-hand transforms into u and v. The matrix
604** is modified in place.
605**
606** That is, given u, b, v where b is bidiagonal and u & v are
607** rotations, modify the matrices so that (old) u b v = (new) u b v,
608** where the new b is diagonal, and u & v are still rotations.
609**
610** b is m x n, u is m x min , and v is min x n
611**
612** This is Golub's SVD algorithm (Algorithm 8.3-2)
613*/
614
615 template <class REAL>
616 void bidiagonal_svd(REAL *b,int m,int n,REAL *u,REAL *v)
617{
618 REAL anorm, t;
619 int p, q, i, j, z, iter, min=MIN(m,n);
620
621 /* use the max of the sum of each diagonal element and the
622 ** superdiagonal element above it as a norm for the array.
623 ** The use of this norm comes from the Numerical Recipes (Press et al)
624 ** code. squawk!
625 */
626
627 // CHUCKR: this is where the array access if fouling up
628 //for (anorm=REF2D(b,m,n,0,0), i=1; i<n; i++) {
629 for (anorm = REF2D(b,m,n,0,0), i=1 ; i < min ; i++)
630 {
631 t = ((REAL) (fabs(REF2D(b,m,n,i,i)) + fabs(REF2D(b,m,n,i-1,i))));
632
633 if (t > anorm)
634 anorm = t;
635 }
636
637 /* Each iteration of this loop will zero out the superdiagonal
638 ** element above b[q][q] -- i.e. b[q][q] will be a singular value.
639 */
640 for (q = min-1; q>=0; q--) {
641
642 /* Each iteration will make the superdiagonal elements smaller,
643 ** till they drop below numerical noise. It ought to converge
644 ** by 30 tries. squawk!
645 */
646 for (iter=0; iter<30; iter++) {
647 /* Find a block that has a zero on the diagonal or superdiagonal.
648 ** That is, we are breaking up b into three submatrices,
649 ** b11, b22, b33, such that b33 is diagonal and b22 has no zeros
650 ** on the superdiagonal.
651 ** b33 goes from q+1 to n-1 -- it's the part we already did
652 ** b22 goes from p through q -- it's a bidiagonal block
653 */
654
655 /* sweep back till we reach a numerically 0 superdiagonal element */
656 for (p = q; p>0; p--) if (REF2D(b,m,n,p-1,p) + anorm == anorm) {
657 REF2D(b,m,n,p-1,p) = 0;
658 break;
659 }
660
661 /* if b22 is 1x1, i.e. there is a 0 above b[q,q], then
662 ** b[q,q] is the singular value. Go on to the next
663 ** singular value. (But first, make sure it's
664 ** positive. If it's negative, negate it and the
665 ** corresponding row of v)
666 */
667
668 if (p==q) {
669 if (REF2D(b,m,n,q,q) < 0) {
670 REF2D(b,m,n,q,q) *= -1;
671 for (j=0; j<n; j++) REF2D(v,min,n,q,j) *= -1;
672 }
673 break;
674 }
675
676 /* check for zero on the diagonal */
677 for (z= -1, i=q; i>=p; i--) if (REF2D(b,m,n,i,i)+anorm==anorm) {
678 z = i;
679 break;
680 }
681 if (z >= 0) {
682 /* get rid of zero on the diagonl. where is it? */
683 if (z == q)
684 /* lower-right corner. clr element above it */
685 clr_bot_supdiag_elt(b,0,z,v,m,n,min);
686 else
687 /* in the middle. clr elt to its right */
688 clr_top_supdiag_elt(b,z,q,u,m,n,min);
689 }
690 else {
691 /* b22 has no zeroes on the diagonal or superdiagonal.
692 ** Do magic rotation, leaving glop in the uppermost
693 ** subdiagonal element
694 */
695 golub_kahn_svd_rot(b,p,q,v,m,n,min);
696 /* get rid of the glop in the uppermost subdiagonal */
697 clr_top_subdiag_elt(b,p,q,u,v,m,n,min);
698 }
699 }
701// (void) fprintf(stderr,"svd: didn't converge after 30 iterations\n");
702 }
703}
704
705
706/**************************************************************************
707**
708** SVD
709**
710** Given a (m x n)
711** computes u (m x min) column-orthonormal
712** w (min x min) diagonal
713** vt (min x n) column-orthonormal
714** where min=min(m,n),
715** such that a = u w vt
716**
717** w is returned as a vector of length min
718**
719** t
720** NOTE: the SVD is commonly represented as U W V
721** where U is m x min and V is n x min. This routine
722** computes the min x n transpose of V, placing the result
723** in (the memory pointed to by) parameter "vt"
724**
725****************************************************************************/
726
727 template <class REAL>
728 void num_svd(const REAL *a,int m,int n,REAL *u,REAL *w,REAL *vt)
729{
730 int i, j, k, min = MIN(m,n);
731 REAL *g, *p, wi;
732
733 g = ALLOC2D(m,n);
734 bidiagonalize(a,m,n,u,g,vt);
735
736
737 bidiagonal_svd(g,m,n,u,vt);
738 for (i=0; i<min; i++) w[i] = REF2D(g,m,n,i,i);
739
740 /* insertion sort singular values. Note that since
741 ** the svd algorithm produces values sorted or nearly sorted,
742 ** an insertion sort is efficient.
743 */
744 for (i=1; i<min; i++) {
745 if (w[i] > w[i-1]) {
746 /* save w[i], and col i of u and row i of vt. use "g" as buffer */
747 wi = w[i];
748 p=g;
749 for (k=0; k<m; k++) *p++ = REF2D(u,m,min,k,i);
750 for (k=0; k<n; k++) *p++ = REF2D(vt,min,n,i,k);
751 /* slide columns over */
752 for (j=i; j>0 && wi > w[j-1]; j--) {
753 w[j] = w[j-1];
754 for (k=0; k<m; k++) REF2D(u,m,min,k,j)=REF2D(u,m,min,k,j-1);
755 for (k=0; k<n; k++) REF2D(vt,min,n,j,k)=REF2D(vt,min,n,j-1,k);
756 }
757 /* put original i stuff where we ended up */
758 w[j] = wi;
759 p=g;
760 for (k=0; k<m; k++) REF2D(u,m,min,k,j) = *p++;
761 for (k=0; k<n; k++) REF2D(vt,min,n,j,k) = *p++;
762 }
763 }
764
765 FREE(g);
766}
767
768/**************************************************************************
769**
770** SVD_BACKSUBST
771**
772** Given the svd of a matrix A, and a vector B, solves A X = B for
773** a vector X. Takes a condition-threshold factor "eps", singular
774** values less than "eps" times the largest value are considered to be
775** zero.
776**
777** input:
778** u (m x min) column-orthonormal
779** w (min-element vector) diagonal elements
780** vt (min x n) column-orthonormal
781** where min=min(m,n),
782** such that a = u w vt, as decomposed by num_svd()
783** b (m-element vector)
784**
785** output:
786** x (n-element vector)
787**
788****************************************************************************/
789
790 template <class REAL>
791 void num_svd_backsubst(int m, int n, const REAL *u, const REAL *w, const REAL *vt,
792 const REAL b[], REAL x[], REAL eps)
793{
794 const int min = MIN(m,n);
795 const REAL thresh = eps * w[0]; /* singular vals are sorted, w[0] is max */
796 REAL *tmp = ALLOC1D(min);
797 int i,j,k;
798
799 for (j=0; j<min; j++) {
800 REAL s = 0;
801 if (w[j] >= thresh) {
802 for (i=0; i<m; i++) s += REF2D(u,m,min,i,j) * b[i];
803 s /= w[j];
804 }
805 tmp[j] = s;
806 }
807 for (k=0; k<n; k++) {
808 REAL s=0;
809 for (j=0; j<min; j++) s += REF2D(vt,min,n,j,k) * tmp[j];
810 x[k] = s;
811 }
812 FREE(tmp);
813}
814
815#undef MIN
816#undef MAX
817#undef ALLOC2D
818#undef ALLOC1D
819#undef FREE
820#undef CLEAR2D
821#undef COPY2D
822#undef REF2D
FREE
#define FREE(p)
Definition SVD.inl:106
bidiagonal_svd
void bidiagonal_svd(REAL *b, int m, int n, REAL *u, REAL *v)
Definition SVD.inl:616
bidiagonalize
void bidiagonalize(const REAL *a, int m, int n, REAL *u, REAL *b, REAL *v)
Definition SVD.inl:555
ALLOC2D
#define ALLOC2D(m, n)
Definition SVD.inl:104
num_svd
void num_svd(const REAL *a, int m, int n, REAL *u, REAL *w, REAL *vt)
Definition SVD.inl:728
golub_kahn_svd_rot
void golub_kahn_svd_rot(REAL *b, int p, int q, REAL *v, int m, int n, int min)
Definition SVD.inl:470
rotate_cols
void rotate_cols(int i, int j, REAL cos, REAL sin, REAL *a, int start, int n, int mm, int nn)
Definition SVD.inl:234
fhypot
REAL fhypot(REAL a, REAL b)
Definition SVD.inl:112
ALLOC1D
#define ALLOC1D(n)
Definition SVD.inl:105
clr_top_subdiag_elt
void clr_top_subdiag_elt(REAL *b, int p, int z, REAL *u, REAL *v, int m, int n, int min)
Definition SVD.inl:385
CLEAR2D
#define CLEAR2D(a, m, n)
Definition SVD.inl:107
COPY2D
#define COPY2D(d, m, n, s)
Definition SVD.inl:108
clr_top_supdiag_elt
void clr_top_supdiag_elt(REAL *b, int p, int z, REAL *u, int m, int n, int min)
Definition SVD.inl:286
householder_zero_row
void householder_zero_row(REAL *a, REAL *v, int i, int j, int m, int n, REAL *hv)
Definition SVD.inl:197
num_svd_backsubst
void num_svd_backsubst(int m, int n, const REAL *u, const REAL *w, const REAL *vt, const REAL b[], REAL x[], REAL eps)
Definition SVD.inl:791
REF2D
#define REF2D(a, m, n, i, j)
Definition SVD.inl:109
householder_zero_col
void householder_zero_col(REAL *a, REAL *u, int i, int j, int m, int n, REAL *hv)
Definition SVD.inl:136
MAX
#define MAX(A, B)
Definition SVD.inl:103
MIN
#define MIN(A, B)
Definition SVD.inl:102
rotate_rows
void rotate_rows(int i, int j, REAL cos, REAL sin, REAL *a, int start, int n, int mm, int nn)
Definition SVD.inl:252
clr_bot_supdiag_elt
void clr_bot_supdiag_elt(REAL *b, int p, int z, REAL *v, int m, int n, int min)
Definition SVD.inl:331
b
long b
Definition jpegint.h:371