[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

singular_value_decomposition.hxx
1 /************************************************************************/
2 /* */
3 /* Copyright 2007 by Ullrich Koethe */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35 
36 
37 #ifndef VIGRA_SINGULAR_VALUE_DECOMPOSITION_HXX
38 #define VIGRA_SINGULAR_VALUE_DECOMPOSITION_HXX
39 
40 #include "matrix.hxx"
41 #include "array_vector.hxx"
42 
43 
44 namespace vigra
45 {
46 
47 namespace linalg
48 {
49 
50  /** Singular Value Decomposition.
51  \ingroup MatrixAlgebra
52 
53  For an m-by-n matrix \a A with m >= n, the singular value decomposition is
54  an m-by-n orthogonal matrix \a U, an n-by-n diagonal matrix S, and
55  an n-by-n orthogonal matrix \a V so that A = U*S*V'.
56 
57  To save memory, this functions stores the matrix \a S in a column vector of
58  appropriate length (a diagonal matrix can be obtained by <tt>diagonalMatrix(S)</tt>).
59  The singular values, sigma[k] = S(k, 0), are ordered so that
60  sigma[0] >= sigma[1] >= ... >= sigma[n-1].
61 
62  The singular value decomposition always exists, so this function will
63  never fail (except if the shapes of the argument matrices don't match).
64  The effective numerical rank of A is returned.
65 
66  (Adapted from JAMA, a Java Matrix Library, developed jointly
67  by the Mathworks and NIST; see http://math.nist.gov/javanumerics/jama).
68 
69  <b>\#include</b> <<a href="singular__value__decomposition_8hxx-source.html">vigra/singular_value_decomposition.hxx</a>> or<br>
70  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
71  Namespaces: vigra and vigra::linalg
72  */
73 template <class T, class C1, class C2, class C3, class C4>
74 unsigned int
77 {
78  typedef T Real;
79  typedef MultiArrayShape<2>::type Shape;
80 
81  const MultiArrayIndex rows = rowCount(A);
82  const MultiArrayIndex cols = columnCount(A);
83  vigra_precondition(rows >= cols,
84  "singularValueDecomposition(): Input matrix A must be rectangular with rowCount >= columnCount.");
85  vigra_precondition(rowCount(S) == cols && columnCount(S) == 1,
86  "singularValueDecomposition(): Output S must be column vector with rowCount == columnCount(A).");
87  vigra_precondition(rowCount(U) == rows && columnCount(U) == cols,
88  "singularValueDecomposition(): Output matrix U must have the same dimensions as input matrix A.");
89  vigra_precondition(rowCount(V) == cols && columnCount(V) == cols,
90  "singularValueDecomposition(): Output matrix V must be square with n = columnCount(A).");
91 
92  MultiArrayIndex m = rows;
93  MultiArrayIndex n = cols;
94  MultiArrayIndex nu = n;
95 
96  U.init(0.0);
97  S.init(0.0);
98  V.init(0.0);
99 
100  ArrayVector<Real> e((unsigned int)n);
101  ArrayVector<Real> work((unsigned int)m);
102  Matrix<Real> a(A);
104 
105  MultiArrayIndex i=0, j=0, k=0;
106 
107  // Reduce a to bidiagonal form, storing the diagonal elements
108  // in s and the super-diagonal elements in e.
109  MultiArrayIndex nct = std::min(m-1,n);
110  MultiArrayIndex nrt = std::max((MultiArrayIndex)0,n-2);
111  for (k = 0; k < std::max(nct,nrt); ++k)
112  {
113  if (k < nct)
114  {
115  // Compute the transformation for the k-th column and
116  // place the k-th diagonal in s(k).
117  // Compute 2-norm of k-th column without under/overflow.
118  s(k) = 0.0;
119  for (i = k; i < m; ++i)
120  {
121  s(k) = hypot(s(k), a(i, k));
122  }
123  if (s(k) != 0.0)
124  {
125  if (a(k, k) < 0.0)
126  {
127  s(k) = -s(k);
128  }
129  for (i = k; i < m; ++i)
130  {
131  a(i, k) /= s(k);
132  }
133  a(k, k) += 1.0;
134  }
135  s(k) = -s(k);
136  }
137  for (j = k+1; j < n; ++j)
138  {
139  if ((k < nct) && (s(k) != 0.0))
140  {
141  // Apply the transformation.
142  Real t(0.0);
143  for (i = k; i < m; ++i)
144  {
145  t += a(i, k)*a(i, j);
146  }
147  t = -t/a(k, k);
148  for (i = k; i < m; ++i)
149  {
150  a(i, j) += t*a(i, k);
151  }
152  }
153 
154  // Place the k-th row of a into e for the
155  // subsequent calculation of the row transformation.
156 
157  e[j] = a(k, j);
158  }
159  if (k < nct)
160  {
161  // Place the transformation in U for subsequent back
162  // multiplication.
163 
164  for (i = k; i < m; ++i)
165  {
166  U(i, k) = a(i, k);
167  }
168  }
169  if (k < nrt)
170  {
171  // Compute the k-th row transformation and place the
172  // k-th super-diagonal in e[k].
173  // Compute 2-norm without under/overflow.
174  e[k] = 0;
175  for (i = k+1; i < n; ++i)
176  {
177  e[k] = hypot(e[k],e[i]);
178  }
179  if (e[k] != 0.0)
180  {
181  if (e[k+1] < 0.0)
182  {
183  e[k] = -e[k];
184  }
185  for (i = k+1; i < n; ++i)
186  {
187  e[i] /= e[k];
188  }
189  e[k+1] += 1.0;
190  }
191  e[k] = -e[k];
192  if ((k+1 < m) & (e[k] != 0.0))
193  {
194  // Apply the transformation.
195  for (i = k+1; i < m; ++i)
196  {
197  work[i] = 0.0;
198  }
199  for (j = k+1; j < n; ++j)
200  {
201  for (i = k+1; i < m; ++i)
202  {
203  work[i] += e[j]*a(i, j);
204  }
205  }
206  for (j = k+1; j < n; ++j)
207  {
208  Real t(-e[j]/e[k+1]);
209  for (i = k+1; i < m; ++i)
210  {
211  a(i, j) += t*work[i];
212  }
213  }
214  }
215  // Place the transformation in V for subsequent
216  // back multiplication.
217  for (i = k+1; i < n; ++i)
218  {
219  V(i, k) = e[i];
220  }
221  }
222  }
223 
224  // Set up the final bidiagonal matrix of order p.
225 
226  MultiArrayIndex p = n;
227  if (nct < n)
228  {
229  s(nct) = a(nct, nct);
230  }
231  if (m < p)
232  {
233  s(p-1) = 0.0;
234  }
235  if (nrt+1 < p)
236  {
237  e[nrt] = a(nrt, p-1);
238  }
239  e[p-1] = 0.0;
240 
241  // Generate U.
242  for (j = nct; j < nu; ++j)
243  {
244  for (i = 0; i < m; ++i)
245  {
246  U(i, j) = 0.0;
247  }
248  U(j, j) = 1.0;
249  }
250  for (k = nct-1; k >= 0; --k)
251  {
252  if (s(k) != 0.0)
253  {
254  for (j = k+1; j < nu; ++j)
255  {
256  Real t(0.0);
257  for (i = k; i < m; ++i)
258  {
259  t += U(i, k)*U(i, j);
260  }
261  t = -t/U(k, k);
262  for (i = k; i < m; ++i)
263  {
264  U(i, j) += t*U(i, k);
265  }
266  }
267  for (i = k; i < m; ++i )
268  {
269  U(i, k) = -U(i, k);
270  }
271  U(k, k) = 1.0 + U(k, k);
272  for (i = 0; i < k-1; ++i)
273  {
274  U(i, k) = 0.0;
275  }
276  }
277  else
278  {
279  for (i = 0; i < m; ++i)
280  {
281  U(i, k) = 0.0;
282  }
283  U(k, k) = 1.0;
284  }
285  }
286 
287  // Generate V.
288  for (k = n-1; k >= 0; --k)
289  {
290  if ((k < nrt) & (e[k] != 0.0))
291  {
292  for (j = k+1; j < nu; ++j)
293  {
294  Real t(0.0);
295  for (i = k+1; i < n; ++i)
296  {
297  t += V(i, k)*V(i, j);
298  }
299  t = -t/V(k+1, k);
300  for (i = k+1; i < n; ++i)
301  {
302  V(i, j) += t*V(i, k);
303  }
304  }
305  }
306  for (i = 0; i < n; ++i)
307  {
308  V(i, k) = 0.0;
309  }
310  V(k, k) = 1.0;
311  }
312 
313  // Main iteration loop for the singular values.
314 
315  MultiArrayIndex pp = p-1;
316  int iter = 0;
317  Real eps = NumericTraits<Real>::epsilon()*2.0;
318  while (p > 0)
319  {
320  MultiArrayIndex k=0;
321  int kase=0;
322 
323  // Here is where a test for too many iterations would go.
324 
325  // This section of the program inspects for
326  // negligible elements in the s and e arrays. On
327  // completion the variables kase and k are set as follows.
328 
329  // kase = 1 if s(p) and e[k-1] are negligible and k<p
330  // kase = 2 if s(k) is negligible and k<p
331  // kase = 3 if e[k-1] is negligible, k<p, and
332  // s(k), ..., s(p) are not negligible (qr step).
333  // kase = 4 if e(p-1) is negligible (convergence).
334 
335  for (k = p-2; k >= -1; --k)
336  {
337  if (k == -1)
338  {
339  break;
340  }
341  if (abs(e[k]) <= eps*(abs(s(k)) + abs(s(k+1))))
342  {
343  e[k] = 0.0;
344  break;
345  }
346  }
347  if (k == p-2)
348  {
349  kase = 4;
350  }
351  else
352  {
353  MultiArrayIndex ks;
354  for (ks = p-1; ks >= k; --ks)
355  {
356  if (ks == k)
357  {
358  break;
359  }
360  Real t( (ks != p ? abs(e[ks]) : 0.) +
361  (ks != k+1 ? abs(e[ks-1]) : 0.));
362  if (abs(s(ks)) <= eps*t)
363  {
364  s(ks) = 0.0;
365  break;
366  }
367  }
368  if (ks == k)
369  {
370  kase = 3;
371  }
372  else if (ks == p-1)
373  {
374  kase = 1;
375  }
376  else
377  {
378  kase = 2;
379  k = ks;
380  }
381  }
382  ++k;
383 
384  // Perform the task indicated by kase.
385 
386  switch (kase)
387  {
388  case 1: // Deflate negligible s(p).
389  {
390  Real f(e[p-2]);
391  e[p-2] = 0.0;
392  for (j = p-2; j >= k; --j)
393  {
394  Real t( hypot(s(j),f));
395  Real cs(s(j)/t);
396  Real sn(f/t);
397  s(j) = t;
398  if (j != k)
399  {
400  f = -sn*e[j-1];
401  e[j-1] = cs*e[j-1];
402  }
403  for (i = 0; i < n; ++i)
404  {
405  t = cs*V(i, j) + sn*V(i, p-1);
406  V(i, p-1) = -sn*V(i, j) + cs*V(i, p-1);
407  V(i, j) = t;
408  }
409  }
410  break;
411  }
412  case 2: // Split at negligible s(k).
413  {
414  Real f(e[k-1]);
415  e[k-1] = 0.0;
416  for (j = k; j < p; ++j)
417  {
418  Real t(hypot(s(j),f));
419  Real cs( s(j)/t);
420  Real sn(f/t);
421  s(j) = t;
422  f = -sn*e[j];
423  e[j] = cs*e[j];
424  for (i = 0; i < m; ++i)
425  {
426  t = cs*U(i, j) + sn*U(i, k-1);
427  U(i, k-1) = -sn*U(i, j) + cs*U(i, k-1);
428  U(i, j) = t;
429  }
430  }
431  break;
432  }
433  case 3: // Perform one qr step.
434  {
435  // Calculate the shift.
436  Real scale = std::max(std::max(std::max(std::max(
437  abs(s(p-1)),abs(s(p-2))),abs(e[p-2])),
438  abs(s(k))),abs(e[k]));
439  Real sp = s(p-1)/scale;
440  Real spm1 = s(p-2)/scale;
441  Real epm1 = e[p-2]/scale;
442  Real sk = s(k)/scale;
443  Real ek = e[k]/scale;
444  Real b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/2.0;
445  Real c = (sp*epm1)*(sp*epm1);
446  Real shift = 0.0;
447  if ((b != 0.0) || (c != 0.0))
448  {
449  shift = VIGRA_CSTD::sqrt(b*b + c);
450  if (b < 0.0)
451  {
452  shift = -shift;
453  }
454  shift = c/(b + shift);
455  }
456  Real f = (sk + sp)*(sk - sp) + shift;
457  Real g = sk*ek;
458 
459  // Chase zeros.
460  for (j = k; j < p-1; ++j)
461  {
462  Real t = hypot(f,g);
463  Real cs = f/t;
464  Real sn = g/t;
465  if (j != k)
466  {
467  e[j-1] = t;
468  }
469  f = cs*s(j) + sn*e[j];
470  e[j] = cs*e[j] - sn*s(j);
471  g = sn*s(j+1);
472  s(j+1) = cs*s(j+1);
473  for (i = 0; i < n; ++i)
474  {
475  t = cs*V(i, j) + sn*V(i, j+1);
476  V(i, j+1) = -sn*V(i, j) + cs*V(i, j+1);
477  V(i, j) = t;
478  }
479  t = hypot(f,g);
480  cs = f/t;
481  sn = g/t;
482  s(j) = t;
483  f = cs*e[j] + sn*s(j+1);
484  s(j+1) = -sn*e[j] + cs*s(j+1);
485  g = sn*e[j+1];
486  e[j+1] = cs*e[j+1];
487  if (j < m-1)
488  {
489  for (i = 0; i < m; ++i)
490  {
491  t = cs*U(i, j) + sn*U(i, j+1);
492  U(i, j+1) = -sn*U(i, j) + cs*U(i, j+1);
493  U(i, j) = t;
494  }
495  }
496  }
497  e[p-2] = f;
498  iter = iter + 1;
499  break;
500  }
501  case 4: // Convergence.
502  {
503  // Make the singular values positive.
504  if (s(k) <= 0.0)
505  {
506  s(k) = (s(k) < 0.0 ? -s(k) : 0.0);
507  for (i = 0; i <= pp; ++i)
508  {
509  V(i, k) = -V(i, k);
510  }
511  }
512 
513  // Order the singular values.
514 
515  while (k < pp)
516  {
517  if (s(k) >= s(k+1))
518  {
519  break;
520  }
521  Real t = s(k);
522  s(k) = s(k+1);
523  s(k+1) = t;
524  if (k < n-1)
525  {
526  for (i = 0; i < n; ++i)
527  {
528  t = V(i, k+1); V(i, k+1) = V(i, k); V(i, k) = t;
529  }
530  }
531  if (k < m-1)
532  {
533  for (i = 0; i < m; ++i)
534  {
535  t = U(i, k+1); U(i, k+1) = U(i, k); U(i, k) = t;
536  }
537  }
538  ++k;
539  }
540  iter = 0;
541  --p;
542  break;
543  }
544  default:
545  vigra_fail("vigra::svd(): Internal error.");
546  }
547  }
548  Real tol = std::max(m,n)*s(0)*eps;
549  unsigned int rank = 0;
550  for (MultiArrayIndex i = 0; i < n; ++i)
551  {
552  if (s(i) > tol)
553  {
554  ++rank;
555  }
556  }
557  return rank; // effective rank
558 }
559 
560 } // namespace linalg
561 
563 
564 } // namespace vigra
565 
566 #endif // VIGRA_SINGULAR_VALUE_DECOMPOSITION_HXX

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.7.1 (Thu Jun 14 2012)