Essentia  2.1-beta6-dev
jama_svd.h
Go to the documentation of this file.
1 #ifndef JAMA_SVD_H
2 #define JAMA_SVD_H
3 
4 
5 #include "tnt_array1d.h"
6 #include "tnt_array1d_utils.h"
7 #include "tnt_array2d.h"
8 #include "tnt_array2d_utils.h"
9 #include "tnt_math_utils.h"
10 
11 #include <algorithm>
12 // for min(), max() below
13 #include <cmath>
14 // for abs() below
15 
16 using namespace TNT;
17 using namespace std;
18 
19 namespace JAMA
20 {
38 template <class Real>
39 class SVD
40 {
41 
42 
45  int m, n;
46 
47  public:
48 
49 
50  SVD (const Array2D<Real> &Arg) {
51 
52 
53  m = Arg.dim1();
54  n = Arg.dim2();
55  int nu = min(m,n);
56  s = Array1D<Real>(min(m+1,n));
57  U = Array2D<Real>(m, nu, Real(0));
58  V = Array2D<Real>(n,n);
59  Array1D<Real> e(n);
60  Array1D<Real> work(m);
61  Array2D<Real> A(Arg.copy());
62  int wantu = 1; /* boolean */
63  int wantv = 1; /* boolean */
64  int i=0, j=0, k=0;
65 
66  // Reduce A to bidiagonal form, storing the diagonal elements
67  // in s and the super-diagonal elements in e.
68 
69  int nct = min(m-1,n);
70  int nrt = max(0,min(n-2,m));
71  for (k = 0; k < max(nct,nrt); k++) {
72  if (k < nct) {
73 
74  // Compute the transformation for the k-th column and
75  // place the k-th diagonal in s[k].
76  // Compute 2-norm of k-th column without under/overflow.
77  s[k] = 0;
78  for (i = k; i < m; i++) {
79  s[k] = hypot(s[k],A[i][k]);
80  }
81  if (s[k] != 0.0) {
82  if (A[k][k] < 0.0) {
83  s[k] = -s[k];
84  }
85  for (i = k; i < m; i++) {
86  A[i][k] /= s[k];
87  }
88  A[k][k] += 1.0;
89  }
90  s[k] = -s[k];
91  }
92  for (j = k+1; j < n; j++) {
93  if ((k < nct) && (s[k] != 0.0)) {
94 
95  // Apply the transformation.
96 
97  double t = 0;
98  for (i = k; i < m; i++) {
99  t += A[i][k]*A[i][j];
100  }
101  t = -t/A[k][k];
102  for (i = k; i < m; i++) {
103  A[i][j] += t*A[i][k];
104  }
105  }
106 
107  // Place the k-th row of A into e for the
108  // subsequent calculation of the row transformation.
109 
110  e[j] = A[k][j];
111  }
112  if (wantu & (k < nct)) {
113 
114  // Place the transformation in U for subsequent back
115  // multiplication.
116 
117  for (i = k; i < m; i++) {
118  U[i][k] = A[i][k];
119  }
120  }
121  if (k < nrt) {
122 
123  // Compute the k-th row transformation and place the
124  // k-th super-diagonal in e[k].
125  // Compute 2-norm without under/overflow.
126  e[k] = 0;
127  for (i = k+1; i < n; i++) {
128  e[k] = hypot(e[k],e[i]);
129  }
130  if (e[k] != 0.0) {
131  if (e[k+1] < 0.0) {
132  e[k] = -e[k];
133  }
134  for (i = k+1; i < n; i++) {
135  e[i] /= e[k];
136  }
137  e[k+1] += 1.0;
138  }
139  e[k] = -e[k];
140  if ((k+1 < m) & (e[k] != 0.0)) {
141 
142  // Apply the transformation.
143 
144  for (i = k+1; i < m; i++) {
145  work[i] = 0.0;
146  }
147  for (j = k+1; j < n; j++) {
148  for (i = k+1; i < m; i++) {
149  work[i] += e[j]*A[i][j];
150  }
151  }
152  for (j = k+1; j < n; j++) {
153  double t = -e[j]/e[k+1];
154  for (i = k+1; i < m; i++) {
155  A[i][j] += t*work[i];
156  }
157  }
158  }
159  if (wantv) {
160 
161  // Place the transformation in V for subsequent
162  // back multiplication.
163 
164  for (i = k+1; i < n; i++) {
165  V[i][k] = e[i];
166  }
167  }
168  }
169  }
170 
171  // Set up the final bidiagonal matrix or order p.
172 
173  int p = min(n,m+1);
174  if (nct < n) {
175  s[nct] = A[nct][nct];
176  }
177  if (m < p) {
178  s[p-1] = 0.0;
179  }
180  if (nrt+1 < p) {
181  e[nrt] = A[nrt][p-1];
182  }
183  e[p-1] = 0.0;
184 
185  // If required, generate U.
186 
187  if (wantu) {
188  for (j = nct; j < nu; j++) {
189  for (i = 0; i < m; i++) {
190  U[i][j] = 0.0;
191  }
192  U[j][j] = 1.0;
193  }
194  for (k = nct-1; k >= 0; k--) {
195  if (s[k] != 0.0) {
196  for (j = k+1; j < nu; j++) {
197  double t = 0;
198  for (i = k; i < m; i++) {
199  t += U[i][k]*U[i][j];
200  }
201  t = -t/U[k][k];
202  for (i = k; i < m; i++) {
203  U[i][j] += t*U[i][k];
204  }
205  }
206  for (i = k; i < m; i++ ) {
207  U[i][k] = -U[i][k];
208  }
209  U[k][k] = 1.0 + U[k][k];
210  for (i = 0; i < k-1; i++) {
211  U[i][k] = 0.0;
212  }
213  } else {
214  for (i = 0; i < m; i++) {
215  U[i][k] = 0.0;
216  }
217  U[k][k] = 1.0;
218  }
219  }
220  }
221 
222  // If required, generate V.
223 
224  if (wantv) {
225  for (k = n-1; k >= 0; k--) {
226  if ((k < nrt) & (e[k] != 0.0)) {
227  for (j = k+1; j < nu; j++) {
228  double t = 0;
229  for (i = k+1; i < n; i++) {
230  t += V[i][k]*V[i][j];
231  }
232  t = -t/V[k+1][k];
233  for (i = k+1; i < n; i++) {
234  V[i][j] += t*V[i][k];
235  }
236  }
237  }
238  for (i = 0; i < n; i++) {
239  V[i][k] = 0.0;
240  }
241  V[k][k] = 1.0;
242  }
243  }
244 
245  // Main iteration loop for the singular values.
246 
247  int pp = p-1;
248  int iter = 0;
249  double eps = pow(2.0,-52.0);
250  while (p > 0) {
251  int k=0;
252  int kase=0;
253 
254  // Here is where a test for too many iterations would go.
255 
256  // This section of the program inspects for
257  // negligible elements in the s and e arrays. On
258  // completion the variables kase and k are set as follows.
259 
260  // kase = 1 if s(p) and e[k-1] are negligible and k<p
261  // kase = 2 if s(k) is negligible and k<p
262  // kase = 3 if e[k-1] is negligible, k<p, and
263  // s(k), ..., s(p) are not negligible (qr step).
264  // kase = 4 if e(p-1) is negligible (convergence).
265 
266  for (k = p-2; k >= -1; k--) {
267  if (k == -1) {
268  break;
269  }
270  if (abs(e[k]) <= eps*(abs(s[k]) + abs(s[k+1]))) {
271  e[k] = 0.0;
272  break;
273  }
274  }
275  if (k == p-2) {
276  kase = 4;
277  } else {
278  int ks;
279  for (ks = p-1; ks >= k; ks--) {
280  if (ks == k) {
281  break;
282  }
283  double t = (ks != p ? abs(e[ks]) : 0.) +
284  (ks != k+1 ? abs(e[ks-1]) : 0.);
285  if (abs(s[ks]) <= eps*t) {
286  s[ks] = 0.0;
287  break;
288  }
289  }
290  if (ks == k) {
291  kase = 3;
292  } else if (ks == p-1) {
293  kase = 1;
294  } else {
295  kase = 2;
296  k = ks;
297  }
298  }
299  k++;
300 
301  // Perform the task indicated by kase.
302 
303  switch (kase) {
304 
305  // Deflate negligible s(p).
306 
307  case 1: {
308  double f = e[p-2];
309  e[p-2] = 0.0;
310  for (j = p-2; j >= k; j--) {
311  double t = hypot(s[j],f);
312  double cs = s[j]/t;
313  double sn = f/t;
314  s[j] = t;
315  if (j != k) {
316  f = -sn*e[j-1];
317  e[j-1] = cs*e[j-1];
318  }
319  if (wantv) {
320  for (i = 0; i < n; i++) {
321  t = cs*V[i][j] + sn*V[i][p-1];
322  V[i][p-1] = -sn*V[i][j] + cs*V[i][p-1];
323  V[i][j] = t;
324  }
325  }
326  }
327  }
328  break;
329 
330  // Split at negligible s(k).
331 
332  case 2: {
333  double f = e[k-1];
334  e[k-1] = 0.0;
335  for (j = k; j < p; j++) {
336  double t = hypot(s[j],f);
337  double cs = s[j]/t;
338  double sn = f/t;
339  s[j] = t;
340  f = -sn*e[j];
341  e[j] = cs*e[j];
342  if (wantu) {
343  for (i = 0; i < m; i++) {
344  t = cs*U[i][j] + sn*U[i][k-1];
345  U[i][k-1] = -sn*U[i][j] + cs*U[i][k-1];
346  U[i][j] = t;
347  }
348  }
349  }
350  }
351  break;
352 
353  // Perform one qr step.
354 
355  case 3: {
356 
357  // Calculate the shift.
358 
359  double scale = max(max(max(max(
360  abs(s[p-1]),abs(s[p-2])),abs(e[p-2])),
361  abs(s[k])),abs(e[k]));
362  double sp = s[p-1]/scale;
363  double spm1 = s[p-2]/scale;
364  double epm1 = e[p-2]/scale;
365  double sk = s[k]/scale;
366  double ek = e[k]/scale;
367  double b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/2.0;
368  double c = (sp*epm1)*(sp*epm1);
369  double shift = 0.0;
370  if ((b != 0.0) || (c != 0.0)) {
371  shift = sqrt(b*b + c);
372  if (b < 0.0) {
373  shift = -shift;
374  }
375  shift = c/(b + shift);
376  }
377  double f = (sk + sp)*(sk - sp) + shift;
378  double g = sk*ek;
379 
380  // Chase zeros.
381 
382  for (j = k; j < p-1; j++) {
383  double t = hypot(f,g);
384  double cs = f/t;
385  double sn = g/t;
386  if (j != k) {
387  e[j-1] = t;
388  }
389  f = cs*s[j] + sn*e[j];
390  e[j] = cs*e[j] - sn*s[j];
391  g = sn*s[j+1];
392  s[j+1] = cs*s[j+1];
393  if (wantv) {
394  for (i = 0; i < n; i++) {
395  t = cs*V[i][j] + sn*V[i][j+1];
396  V[i][j+1] = -sn*V[i][j] + cs*V[i][j+1];
397  V[i][j] = t;
398  }
399  }
400  t = hypot(f,g);
401  cs = f/t;
402  sn = g/t;
403  s[j] = t;
404  f = cs*e[j] + sn*s[j+1];
405  s[j+1] = -sn*e[j] + cs*s[j+1];
406  g = sn*e[j+1];
407  e[j+1] = cs*e[j+1];
408  if (wantu && (j < m-1)) {
409  for (i = 0; i < m; i++) {
410  t = cs*U[i][j] + sn*U[i][j+1];
411  U[i][j+1] = -sn*U[i][j] + cs*U[i][j+1];
412  U[i][j] = t;
413  }
414  }
415  }
416  e[p-2] = f;
417  iter = iter + 1;
418  }
419  break;
420 
421  // Convergence.
422 
423  case 4: {
424 
425  // Make the singular values positive.
426 
427  if (s[k] <= 0.0) {
428  s[k] = (s[k] < 0.0 ? -s[k] : 0.0);
429  if (wantv) {
430  for (i = 0; i <= pp; i++) {
431  V[i][k] = -V[i][k];
432  }
433  }
434  }
435 
436  // Order the singular values.
437 
438  while (k < pp) {
439  if (s[k] >= s[k+1]) {
440  break;
441  }
442  double t = s[k];
443  s[k] = s[k+1];
444  s[k+1] = t;
445  if (wantv && (k < n-1)) {
446  for (i = 0; i < n; i++) {
447  t = V[i][k+1]; V[i][k+1] = V[i][k]; V[i][k] = t;
448  }
449  }
450  if (wantu && (k < m-1)) {
451  for (i = 0; i < m; i++) {
452  t = U[i][k+1]; U[i][k+1] = U[i][k]; U[i][k] = t;
453  }
454  }
455  k++;
456  }
457  iter = 0;
458  p--;
459  }
460  break;
461  }
462  }
463  }
464 
465 
466  void getU (Array2D<Real> &A)
467  {
468  int minm = min(m+1,n);
469 
470  A = Array2D<Real>(m, minm);
471 
472  for (int i=0; i<m; i++)
473  for (int j=0; j<minm; j++)
474  A[i][j] = U[i][j];
475 
476  }
477 
478  /* Return the right singular vectors */
479 
480  void getV (Array2D<Real> &A)
481  {
482  A = V;
483  }
484 
488  {
489  x = s;
490  }
491 
496  void getS (Array2D<Real> &A) {
497  A = Array2D<Real>(n,n);
498  for (int i = 0; i < n; i++) {
499  for (int j = 0; j < n; j++) {
500  A[i][j] = 0.0;
501  }
502  A[i][i] = s[i];
503  }
504  }
505 
508  double norm2 () {
509  return s[0];
510  }
511 
514  double cond () {
515  return s[0]/s[min(m,n)-1];
516  }
517 
522  int rank ()
523  {
524  double eps = pow(2.0,-52.0);
525  double tol = max(m,n)*s[0]*eps;
526  int r = 0;
527  for (int i = 0; i < s.dim(); i++) {
528  if (s[i] > tol) {
529  r++;
530  }
531  }
532  return r;
533  }
534 };
535 
536 }
537 #endif
538 // JAMA_SVD_H
Definition: jama_svd.h:40
int rank()
Definition: jama_svd.h:522
SVD(const Array2D< Real > &Arg)
Definition: jama_svd.h:50
double cond()
Definition: jama_svd.h:514
int m
Definition: jama_svd.h:45
void getU(Array2D< Real > &A)
Definition: jama_svd.h:466
void getSingularValues(Array1D< Real > &x)
Definition: jama_svd.h:487
void getS(Array2D< Real > &A)
Definition: jama_svd.h:496
void getV(Array2D< Real > &A)
Definition: jama_svd.h:480
double norm2()
Definition: jama_svd.h:508
Array2D< Real > U
Definition: jama_svd.h:43
Array1D< Real > s
Definition: jama_svd.h:44
int dim() const
Definition: tnt_array1d.h:221
int dim1() const
Definition: tnt_array2d.h:231
int dim2() const
Definition: tnt_array2d.h:234
Array2D copy() const
Definition: tnt_array2d.h:180
Definition: jama_cholesky.h:9
Definition: tnt_array1d.h:36
Real hypot(const Real &a, const Real &b)
Definition: tnt_math_utils.h:18
float Real
Definition: types.h:69