SVD-Impl.h
Go to the documentation of this file.
1 /*************************************************************************
2 > File Name: SVD-Impl.h
3 > Project Name: CubbyFlow
4 > Author: Chan-Ho Chris Ohk
5 > Purpose: Singular value decomposition (SVD).
6 > Created Time: 2017/11/19
7 > Copyright (c) 2018, Chan-Ho Chris Ohk
8 *************************************************************************/
9 #ifndef CUBBYFLOW_SVD_IMPL_H
10 #define CUBBYFLOW_SVD_IMPL_H
11 
12 namespace CubbyFlow
13 {
14  namespace Internal
15  {
16  template <typename T>
17  inline T Sign(T a, T b)
18  {
19  return (b >= 0.0) ? std::fabs(a) : -std::fabs(a);
20  }
21 
22  template <typename T>
23  inline T Pythag(T a, T b)
24  {
25  T at = std::fabs(a);
26  T bt = std::fabs(b);
27  T ct;
28  T result;
29 
30  if (at > bt)
31  {
32  ct = bt / at;
33  result = at * std::sqrt(1 + ct * ct);
34  }
35  else if (bt > 0)
36  {
37  ct = at / bt;
38  result = bt * std::sqrt(1 + ct * ct);
39  }
40  else
41  {
42  result = 0;
43  }
44 
45  return result;
46  }
47  }
48 
49  template <typename T>
50  void SVD(const MatrixMxN<T>& a, MatrixMxN<T>& u, VectorN<T>& w, MatrixMxN<T>& v)
51  {
52  const int m = static_cast<int>(a.Rows());
53  const int n = static_cast<int>(a.Cols());
54 
55  int i, j = 0, jj = 0, k = 0, l = 0, nm = 0;
56  T c = 0, f = 0, h = 0, s = 0, x = 0, y = 0, z = 0;
57  T anorm = 0, g = 0, scale = 0;
58 
59  if (m < n)
60  {
61  throw std::invalid_argument("Number of rows of input matrix must greater than or equal to columns.");
62  }
63 
64  // Prepare workspace
65  VectorN<T> rv1(n, 0);
66  u = a;
67  w.Resize(n, 0);
68  v.Resize(n, n, 0);
69 
70  // Householder reduction to bidiagonal form
71  for (i = 0; i < n; i++)
72  {
73  // left-hand reduction
74  l = i + 1;
75  rv1[i] = scale * g;
76  g = s = scale = 0;
77 
78  if (i < m)
79  {
80  for (k = i; k < m; k++)
81  {
82  scale += std::fabs(u(k, i));
83  }
84 
85  if (scale)
86  {
87  for (k = i; k < m; k++)
88  {
89  u(k, i) /= scale;
90  s += u(k, i) * u(k, i);
91  }
92 
93  f = u(i, i);
94  g = -Internal::Sign(std::sqrt(s), f);
95  h = f * g - s;
96  u(i, i) = f - g;
97 
98  if (i != n - 1)
99  {
100  for (j = l; j < n; j++)
101  {
102  for (s = 0, k = i; k < m; k++)
103  {
104  s += u(k, i) * u(k, j);
105  }
106 
107  f = s / h;
108 
109  for (k = i; k < m; k++)
110  {
111  u(k, j) += f * u(k, i);
112  }
113  }
114  }
115 
116  for (k = i; k < m; k++)
117  {
118  u(k, i) *= scale;
119  }
120  }
121  }
122 
123  w[i] = scale * g;
124 
125  // right-hand reduction
126  g = s = scale = 0;
127 
128  if (i < m && i != n - 1)
129  {
130  for (k = l; k < n; k++)
131  {
132  scale += std::fabs(u(i, k));
133  }
134 
135  if (scale)
136  {
137  for (k = l; k < n; k++)
138  {
139  u(i, k) /= scale;
140  s += u(i, k) * u(i, k);
141  }
142 
143  f = u(i, l);
144  g = -Internal::Sign(std::sqrt(s), f);
145  h = f * g - s;
146  u(i, l) = f - g;
147 
148  for (k = l; k < n; k++)
149  {
150  rv1[k] = static_cast<T>(u(i, k)) / h;
151  }
152 
153  if (i != m - 1)
154  {
155  for (j = l; j < m; j++)
156  {
157  for (s = 0, k = l; k < n; k++)
158  {
159  s += u(j, k) * u(i, k);
160  }
161 
162  for (k = l; k < n; k++)
163  {
164  u(j, k) += s * rv1[k];
165  }
166  }
167  }
168 
169  for (k = l; k < n; k++)
170  {
171  u(i, k) *= scale;
172  }
173  }
174  }
175 
176  anorm = std::max(anorm, (std::fabs(static_cast<T>(w[i])) + std::fabs(rv1[i])));
177  }
178 
179  // accumulate the right-hand transformation
180  for (i = n - 1; i >= 0; i--)
181  {
182  if (i < n - 1)
183  {
184  if (g)
185  {
186  for (j = l; j < n; j++)
187  {
188  v(j, i) = ((u(i, j) / u(i, l)) / g);
189  }
190 
191  // T division to avoid underflow
192  for (j = l; j < n; j++)
193  {
194  for (s = 0, k = l; k < n; k++)
195  {
196  s += u(i, k) * v(k, j);
197  }
198 
199  for (k = l; k < n; k++)
200  {
201  v(k, j) += s * v(k, i);
202  }
203  }
204  }
205 
206  for (j = l; j < n; j++)
207  {
208  v(i, j) = v(j, i) = 0;
209  }
210  }
211 
212  v(i, i) = 1;
213  g = rv1[i];
214  l = i;
215  }
216 
217  // accumulate the left-hand transformation
218  for (i = n - 1; i >= 0; i--)
219  {
220  l = i + 1;
221  g = w[i];
222 
223  if (i < n - 1)
224  {
225  for (j = l; j < n; j++)
226  {
227  u(i, j) = 0;
228  }
229  }
230 
231  if (g)
232  {
233  g = 1 / g;
234 
235  if (i != n - 1)
236  {
237  for (j = l; j < n; j++)
238  {
239  for (s = 0, k = l; k < m; k++)
240  {
241  s += u(k, i) * u(k, j);
242  }
243 
244  f = (s / u(i, i)) * g;
245 
246  for (k = i; k < m; k++)
247  {
248  u(k, j) += f * u(k, i);
249  }
250  }
251  }
252 
253  for (j = i; j < m; j++)
254  {
255  u(j, i) = u(j, i) * g;
256  }
257  }
258  else
259  {
260  for (j = i; j < m; j++)
261  {
262  u(j, i) = 0;
263  }
264  }
265 
266  ++u(i, i);
267  }
268 
269  // diagonalize the bidiagonal form
270  for (k = n - 1; k >= 0; k--)
271  {
272  // loop over singular values
273  for (int its = 0; its < 30; its++)
274  {
275  // loop over allowed iterations
276  int flag = 1;
277 
278  for (l = k; l >= 0; l--)
279  {
280  // test for splitting
281  nm = l - 1;
282 
283  if (std::fabs(rv1[l]) + anorm == anorm)
284  {
285  flag = 0;
286  break;
287  }
288 
289  if (std::fabs(static_cast<T>(w[nm])) + anorm == anorm)
290  {
291  break;
292  }
293  }
294 
295  if (flag)
296  {
297  c = 0;
298  s = 1;
299 
300  for (i = l; i <= k; i++)
301  {
302  f = s * rv1[i];
303 
304  if (std::fabs(f) + anorm != anorm)
305  {
306  g = w[i];
307  h = Internal::Pythag(f, g);
308  w[i] = static_cast<T>(h);
309  h = 1 / h;
310  c = g * h;
311  s = -f * h;
312 
313  for (j = 0; j < m; j++)
314  {
315  y = u(j, nm);
316  z = u(j, i);
317  u(j, nm) = y * c + z * s;
318  u(j, i) = z * c - y * s;
319  }
320  }
321  }
322  }
323 
324  z = w[k];
325 
326  if (l == k)
327  {
328  // convergence
329  if (z < 0)
330  {
331  // make singular value nonnegative
332  w[k] = -z;
333 
334  for (j = 0; j < n; j++)
335  {
336  v(j, k) = -v(j, k);
337  }
338  }
339 
340  break;
341  }
342 
343  if (its >= 30)
344  {
345  throw "No convergence after 30 iterations";
346  }
347 
348  // shift from bottom 2 x 2 minor
349  x = w[l];
350  nm = k - 1;
351  y = w[nm];
352  g = rv1[nm];
353  h = rv1[k];
354  f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2 * h * y);
355  g = Internal::Pythag(f, static_cast<T>(1));
356  f = ((x - z) * (x + z) + h * ((y / (f + Internal::Sign(g, f))) - h)) / x;
357 
358  // next QR transformation
359  c = s = 1;
360 
361  for (j = l; j <= nm; j++)
362  {
363  i = j + 1;
364  g = rv1[i];
365  y = w[i];
366  h = s * g;
367  g = c * g;
368  z = Internal::Pythag(f, h);
369  rv1[j] = z;
370  c = f / z;
371  s = h / z;
372  f = x * c + g * s;
373  g = g * c - x * s;
374  h = y * s;
375  y = y * c;
376 
377  for (jj = 0; jj < n; jj++)
378  {
379  x = v(jj, j);
380  z = v(jj, i);
381  v(jj, j) = x * c + z * s;
382  v(jj, i) = z * c - x * s;
383  }
384 
385  z = Internal::Pythag(f, h);
386  w[j] = z;
387 
388  if (z)
389  {
390  z = 1 / z;
391  c = f * z;
392  s = h * z;
393  }
394 
395  f = (c * g) + (s * y);
396  x = (c * y) - (s * g);
397 
398  for (jj = 0; jj < m; jj++)
399  {
400  y = u(jj, j);
401  z = u(jj, i);
402  u(jj, j) = y * c + z * s;
403  u(jj, i) = z * c - y * s;
404  }
405  }
406 
407  rv1[l] = 0;
408  rv1[k] = f;
409  w[k] = x;
410  }
411  }
412  }
413 
414  template <typename T, size_t M, size_t N>
415  void SVD(const Matrix<T, M, N>& a, Matrix<T, M, N>& u, Vector<T, N>& w, Matrix<T, N, N>& v)
416  {
417  const int m = static_cast<int>(M);
418  const int n = static_cast<int>(N);
419 
420  int i, its, j = 0, jj = 0, k = 0, l = 0, nm = 0;
421  T c = 0, f = 0, h = 0, s = 0, x = 0, y = 0, z = 0;
422  T anorm = 0, g = 0, scale = 0;
423 
424  static_assert(m >= n, "Number of rows of input matrix must greater than or equal to columns.");
425 
426  // Prepare workspace
427  Vector<T, N> rv1;
428  u = a;
429  w = Vector<T, N>();
430  v = Matrix<T, N, N>();
431 
432  // Householder reduction to bidiagonal form
433  for (i = 0; i < n; i++)
434  {
435  // left-hand reduction
436  l = i + 1;
437  rv1[i] = scale * g;
438  g = s = scale = 0;
439 
440  if (i < m)
441  {
442  for (k = i; k < m; k++)
443  {
444  scale += std::fabs(u(k, i));
445  }
446 
447  if (scale)
448  {
449  for (k = i; k < m; k++)
450  {
451  u(k, i) /= scale;
452  s += u(k, i) * u(k, i);
453  }
454 
455  f = u(i, i);
456  g = -Internal::Sign(std::sqrt(s), f);
457  h = f * g - s;
458  u(i, i) = f - g;
459 
460  if (i != n - 1)
461  {
462  for (j = l; j < n; j++)
463  {
464  for (s = 0, k = i; k < m; k++)
465  {
466  s += u(k, i) * u(k, j);
467  }
468 
469  f = s / h;
470 
471  for (k = i; k < m; k++)
472  {
473  u(k, j) += f * u(k, i);
474  }
475  }
476  }
477 
478  for (k = i; k < m; k++)
479  {
480  u(k, i) *= scale;
481  }
482  }
483  }
484 
485  w[i] = scale * g;
486 
487  // right-hand reduction
488  g = s = scale = 0;
489 
490  if (i < m && i != n - 1)
491  {
492  for (k = l; k < n; k++)
493  {
494  scale += std::fabs(u(i, k));
495  }
496 
497  if (scale)
498  {
499  for (k = l; k < n; k++)
500  {
501  u(i, k) /= scale;
502  s += u(i, k) * u(i, k);
503  }
504 
505  f = u(i, l);
506  g = -Internal::Sign(std::sqrt(s), f);
507  h = f * g - s;
508  u(i, l) = f - g;
509 
510  for (k = l; k < n; k++)
511  {
512  rv1[k] = static_cast<T>(u(i, k)) / h;
513  }
514 
515  if (i != m - 1)
516  {
517  for (j = l; j < m; j++)
518  {
519  for (s = 0, k = l; k < n; k++)
520  {
521  s += u(j, k) * u(i, k);
522  }
523 
524  for (k = l; k < n; k++)
525  {
526  u(j, k) += s * rv1[k];
527  }
528  }
529  }
530 
531  for (k = l; k < n; k++)
532  {
533  u(i, k) *= scale;
534  }
535  }
536  }
537  anorm = std::max(anorm, (std::fabs(static_cast<T>(w[i])) + std::fabs(rv1[i])));
538  }
539 
540  // accumulate the right-hand transformation
541  for (i = n - 1; i >= 0; i--)
542  {
543  if (i < n - 1)
544  {
545  if (g)
546  {
547  for (j = l; j < n; j++)
548  {
549  v(j, i) = ((u(i, j) / u(i, l)) / g);
550  }
551 
552  // T division to avoid underflow
553  for (j = l; j < n; j++)
554  {
555  for (s = 0, k = l; k < n; k++)
556  {
557  s += u(i, k) * v(k, j);
558  }
559 
560  for (k = l; k < n; k++)
561  {
562  v(k, j) += s * v(k, i);
563  }
564  }
565  }
566 
567  for (j = l; j < n; j++)
568  {
569  v(i, j) = v(j, i) = 0;
570  }
571  }
572 
573  v(i, i) = 1;
574  g = rv1[i];
575  l = i;
576  }
577 
578  // accumulate the left-hand transformation
579  for (i = n - 1; i >= 0; i--)
580  {
581  l = i + 1;
582  g = w[i];
583 
584  if (i < n - 1)
585  {
586  for (j = l; j < n; j++)
587  {
588  u(i, j) = 0;
589  }
590  }
591 
592  if (g)
593  {
594  g = 1 / g;
595 
596  if (i != n - 1)
597  {
598  for (j = l; j < n; j++)
599  {
600  for (s = 0, k = l; k < m; k++)
601  {
602  s += u(k, i) * u(k, j);
603  }
604 
605  f = (s / u(i, i)) * g;
606 
607  for (k = i; k < m; k++)
608  {
609  u(k, j) += f * u(k, i);
610  }
611  }
612  }
613 
614  for (j = i; j < m; j++)
615  {
616  u(j, i) = u(j, i) * g;
617  }
618  }
619  else
620  {
621  for (j = i; j < m; j++)
622  {
623  u(j, i) = 0;
624  }
625  }
626 
627  ++u(i, i);
628  }
629 
630  // diagonalize the bidiagonal form
631  for (k = n - 1; k >= 0; k--)
632  {
633  // loop over singular values
634  for (its = 0; its < 30; its++)
635  {
636  // loop over allowed iterations
637  int flag = 1;
638 
639  for (l = k; l >= 0; l--)
640  {
641  // test for splitting
642  nm = l - 1;
643 
644  if (std::fabs(rv1[l]) + anorm == anorm)
645  {
646  flag = 0;
647  break;
648  }
649 
650  if (std::fabs(static_cast<T>(w[nm])) + anorm == anorm)
651  {
652  break;
653  }
654  }
655 
656  if (flag)
657  {
658  c = 0;
659  s = 1;
660 
661  for (i = l; i <= k; i++)
662  {
663  f = s * rv1[i];
664 
665  if (std::fabs(f) + anorm != anorm)
666  {
667  g = w[i];
668  h = Internal::Pythag(f, g);
669  w[i] = static_cast<T>(h);
670  h = 1 / h;
671  c = g * h;
672  s = -f * h;
673 
674  for (j = 0; j < m; j++)
675  {
676  y = u(j, nm);
677  z = u(j, i);
678  u(j, nm) = y * c + z * s;
679  u(j, i) = z * c - y * s;
680  }
681  }
682  }
683  }
684 
685  z = w[k];
686 
687  if (l == k)
688  {
689  // convergence
690  if (z < 0)
691  {
692  // make singular value nonnegative
693  w[k] = -z;
694 
695  for (j = 0; j < n; j++)
696  {
697  v(j, k) = -v(j, k);
698  }
699  }
700 
701  break;
702  }
703 
704  if (its >= 30)
705  {
706  throw "No convergence after 30 iterations";
707  }
708 
709  // shift from bottom 2 x 2 minor
710  x = w[l];
711  nm = k - 1;
712  y = w[nm];
713  g = rv1[nm];
714  h = rv1[k];
715  f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2 * h * y);
716  g = Internal::Pythag(f, static_cast<T>(1));
717  f = ((x - z) * (x + z) + h * ((y / (f + Internal::Sign(g, f))) - h)) / x;
718 
719  // next QR transformation
720  c = s = 1;
721 
722  for (j = l; j <= nm; j++)
723  {
724  i = j + 1;
725  g = rv1[i];
726  y = w[i];
727  h = s * g;
728  g = c * g;
729  z = Internal::Pythag(f, h);
730  rv1[j] = z;
731  c = f / z;
732  s = h / z;
733  f = x * c + g * s;
734  g = g * c - x * s;
735  h = y * s;
736  y = y * c;
737 
738  for (jj = 0; jj < n; jj++)
739  {
740  x = v(jj, j);
741  z = v(jj, i);
742  v(jj, j) = x * c + z * s;
743  v(jj, i) = z * c - x * s;
744  }
745 
746  z = Internal::Pythag(f, h);
747  w[j] = z;
748 
749  if (z) {
750  z = 1 / z;
751  c = f * z;
752  s = h * z;
753  }
754 
755  f = (c * g) + (s * y);
756  x = (c * y) - (s * g);
757 
758  for (jj = 0; jj < m; jj++)
759  {
760  y = u(jj, j);
761  z = u(jj, i);
762  u(jj, j) = y * c + z * s;
763  u(jj, i) = z * c - y * s;
764  }
765  }
766 
767  rv1[l] = 0;
768  rv1[k] = f;
769  w[k] = x;
770  }
771  }
772  }
773 }
774 
775 #endif
CubbyFlow::Internal::Sign
T Sign(T a, T b)
Definition: SVD-Impl.h:17
CubbyFlow::Vector
Generic statically-sized N-D vector class.
Definition: Vector.h:33
CubbyFlow::MatrixMxN::Cols
size_t Cols() const
Returns number of columns of this matrix.
Definition: MatrixMxN-Impl.h:209
CubbyFlow::VectorN
General purpose dynamically-sizedN-D vector class.
Definition: VectorN.h:27
CubbyFlow::Matrix
Static-sized M x N matrix class.
Definition: Matrix.h:30
CubbyFlow::MatrixMxN::Resize
void Resize(size_t m, size_t n, const T &s=T(0))
Resizes to m x n matrix with initial value s.
Definition: MatrixMxN-Impl.h:59
CubbyFlow
Definition: pybind11Utils.h:24
CubbyFlow::VectorN::Resize
void Resize(size_t n, const T &val=0)
Resizes to n dimensional vector with initial value val.
Definition: VectorN-Impl.h:57
CubbyFlow::Internal::Pythag
T Pythag(T a, T b)
Definition: SVD-Impl.h:23
CubbyFlow::SVD
void SVD(const MatrixMxN< T > &a, MatrixMxN< T > &u, VectorN< T > &w, MatrixMxN< T > &v)
Singular value decomposition (SVD).
Definition: SVD-Impl.h:50
CubbyFlow::MatrixMxN
M x N matrix class.
Definition: MatrixMxN.h:26
CubbyFlow::MatrixMxN::Rows
size_t Rows() const
Returns number of rows of this matrix.
Definition: MatrixMxN-Impl.h:203
Generated on Sun Apr 22 2018 23:43:20 for CubbyFlow by   doxygen 1.8.14