NumCpp  2.9.0
A Templatized Header Only C++ Implementation of the Python NumPy Library
SVDClass.hpp
Go to the documentation of this file.
1 #pragma once
30 
31 #include <cmath>
32 #include <limits>
33 #include <string>
34 
36 #include "NumCpp/Core/Types.hpp"
37 #include "NumCpp/NdArray.hpp"
38 
39 namespace nc
40 {
41  namespace linalg
42  {
43  // =============================================================================
44  // Class Description:
47  class SVD
48  {
49  public:
50  // =============================================================================
51  // Description:
56  explicit SVD(const NdArray<double>& inMatrix) :
57  m_(inMatrix.shape().rows),
58  n_(inMatrix.shape().cols),
59  u_(inMatrix),
60  v_(n_, n_),
61  s_(1, n_),
62  eps_(std::numeric_limits<double>::epsilon())
63  {
64  decompose();
65  reorder();
66  tsh_ = 0.5 * std::sqrt(m_ + n_ + 1.) * s_.front() * eps_;
67  }
68 
69  // =============================================================================
70  // Description:
75  const NdArray<double>& u() noexcept
76  {
77  return u_;
78  }
79 
80  // =============================================================================
81  // Description:
86  const NdArray<double>& v() noexcept
87  {
88  return v_;
89  }
90 
91  // =============================================================================
92  // Description:
97  const NdArray<double>& s() noexcept
98  {
99  return s_;
100  }
101 
102  // =============================================================================
103  // Description:
111  NdArray<double> solve(const NdArray<double>& inInput, double inThresh = -1.)
112  {
113  double ss = 0.;
114 
115  if (inInput.size() != m_)
116  {
117  THROW_INVALID_ARGUMENT_ERROR("bad sizes.");
118  }
119 
120  NdArray<double> returnArray(1, n_);
121 
122  NdArray<double> tmp(1, n_);
123 
124  tsh_ = (inThresh >= 0. ? inThresh : 0.5 * sqrt(m_ + n_ + 1.) * s_.front() * eps_);
125 
126  for (uint32 j = 0; j < n_; j++)
127  {
128  ss = 0.;
129  if (s_[j] > tsh_)
130  {
131  for (uint32 i = 0; i < m_; i++)
132  {
133  ss += u_(i, j) * inInput[i];
134  }
135  ss /= s_[j];
136  }
137  tmp[j] = ss;
138  }
139 
140  for (uint32 j = 0; j < n_; j++)
141  {
142  ss = 0.;
143  for (uint32 jj = 0; jj < n_; jj++)
144  {
145  ss += v_(j, jj) * tmp[jj];
146  }
147 
148  returnArray[j] = ss;
149  }
150 
151  return returnArray;
152  }
153 
154  private:
155  // =============================================================================
156  // Description:
164  static double SIGN(double inA, double inB) noexcept
165  {
166  return inB >= 0 ? (inA >= 0 ? inA : -inA) : (inA >= 0 ? -inA : inA);
167  }
168 
169  // =============================================================================
170  // Description:
173  void decompose()
174  {
175  bool flag = true;
176  uint32 i = 0;
177  uint32 its = 0;
178  uint32 j = 0;
179  uint32 jj = 0;
180  uint32 k = 0;
181  uint32 l = 0;
182  uint32 nm = 0;
183 
184  double anorm = 0.;
185  double c = 0.;
186  double f = 0.;
187  double g = 0.;
188  double h = 0.;
189  double ss = 0.;
190  double scale = 0.;
191  double x = 0.;
192  double y = 0.;
193  double z = 0.;
194 
195  NdArray<double> rv1(n_, 1);
196 
197  for (i = 0; i < n_; ++i)
198  {
199  l = i + 2;
200  rv1[i] = scale * g;
201  g = ss = scale = 0.;
202 
203  if (i < m_)
204  {
205  for (k = i; k < m_; ++k)
206  {
207  scale += std::abs(u_(k, i));
208  }
209 
210  if (scale != 0.)
211  {
212  for (k = i; k < m_; ++k)
213  {
214  u_(k, i) /= scale;
215  ss += u_(k, i) * u_(k, i);
216  }
217 
218  f = u_(i, i);
219  g = -SIGN(std::sqrt(ss), f);
220  h = f * g - ss;
221  u_(i, i) = f - g;
222 
223  for (j = l - 1; j < n_; ++j)
224  {
225  for (ss = 0., k = i; k < m_; ++k)
226  {
227  ss += u_(k, i) * u_(k, j);
228  }
229 
230  f = ss / h;
231 
232  for (k = i; k < m_; ++k)
233  {
234  u_(k, j) += f * u_(k, i);
235  }
236  }
237 
238  for (k = i; k < m_; ++k)
239  {
240  u_(k, i) *= scale;
241  }
242  }
243  }
244 
245  s_[i] = scale * g;
246  g = ss = scale = 0.;
247 
248  if (i + 1 <= m_ && i + 1 != n_)
249  {
250  for (k = l - 1; k < n_; ++k)
251  {
252  scale += std::abs(u_(i, k));
253  }
254 
255  if (scale != 0.)
256  {
257  for (k = l - 1; k < n_; ++k)
258  {
259  u_(i, k) /= scale;
260  ss += u_(i, k) * u_(i, k);
261  }
262 
263  f = u_(i, l - 1);
264  g = -SIGN(std::sqrt(ss), f);
265  h = f * g - ss;
266  u_(i, l - 1) = f - g;
267 
268  for (k = l - 1; k < n_; ++k)
269  {
270  rv1[k] = u_(i, k) / h;
271  }
272 
273  for (j = l - 1; j < m_; ++j)
274  {
275  for (ss = 0., k = l - 1; k < n_; ++k)
276  {
277  ss += u_(j, k) * u_(i, k);
278  }
279 
280  for (k = l - 1; k < n_; ++k)
281  {
282  u_(j, k) += ss * rv1[k];
283  }
284  }
285 
286  for (k = l - 1; k < n_; ++k)
287  {
288  u_(i, k) *= scale;
289  }
290  }
291  }
292 
293  anorm = std::max(anorm, (std::abs(s_[i]) + std::abs(rv1[i])));
294  }
295 
296  for (i = n_ - 1; i != static_cast<uint32>(-1); --i)
297  {
298  if (i < n_ - 1)
299  {
300  if (g != 0.)
301  {
302  for (j = l; j < n_; ++j)
303  {
304  v_(j, i) = (u_(i, j) / u_(i, l)) / g;
305  }
306 
307  for (j = l; j < n_; ++j)
308  {
309  for (ss = 0., k = l; k < n_; ++k)
310  {
311  ss += u_(i, k) * v_(k, j);
312  }
313 
314  for (k = l; k < n_; ++k)
315  {
316  v_(k, j) += ss * v_(k, i);
317  }
318  }
319  }
320 
321  for (j = l; j < n_; ++j)
322  {
323  v_(i, j) = v_(j, i) = 0.;
324  }
325  }
326 
327  v_(i, i) = 1.;
328  g = rv1[i];
329  l = i;
330  }
331 
332  for (i = std::min(m_, n_) - 1; i != static_cast<uint32>(-1); --i)
333  {
334  l = i + 1;
335  g = s_[i];
336 
337  for (j = l; j < n_; ++j)
338  {
339  u_(i, j) = 0.;
340  }
341 
342  if (g != 0.)
343  {
344  g = 1. / g;
345 
346  for (j = l; j < n_; ++j)
347  {
348  for (ss = 0., k = l; k < m_; ++k)
349  {
350  ss += u_(k, i) * u_(k, j);
351  }
352 
353  f = (ss / u_(i, i)) * g;
354 
355  for (k = i; k < m_; ++k)
356  {
357  u_(k, j) += f * u_(k, i);
358  }
359  }
360 
361  for (j = i; j < m_; ++j)
362  {
363  u_(j, i) *= g;
364  }
365  }
366  else
367  {
368  for (j = i; j < m_; ++j)
369  {
370  u_(j, i) = 0.;
371  }
372  }
373 
374  ++u_(i, i);
375  }
376 
377  for (k = n_ - 1; k != static_cast<uint32>(-1); --k)
378  {
379  for (its = 0; its < 30; ++its)
380  {
381  flag = true;
382  for (l = k; l != static_cast<uint32>(-1); --l)
383  {
384  nm = l - 1;
385  if (l == 0 || std::abs(rv1[l]) <= eps_ * anorm)
386  {
387  flag = false;
388  break;
389  }
390 
391  if (std::abs(s_[nm]) <= eps_ * anorm)
392  {
393  break;
394  }
395  }
396 
397  if (flag)
398  {
399  c = 0.;
400  ss = 1.;
401  for (i = l; i < k + 1; ++i)
402  {
403  f = ss * rv1[i];
404  rv1[i] = c * rv1[i];
405 
406  if (std::abs(f) <= eps_ * anorm)
407  {
408  break;
409  }
410 
411  g = s_[i];
412  h = pythag(f, g);
413  s_[i] = h;
414  h = 1. / h;
415  c = g * h;
416  ss = -f * h;
417 
418  for (j = 0; j < m_; ++j)
419  {
420  y = u_(j, nm);
421  z = u_(j, i);
422  u_(j, nm) = y * c + z * ss;
423  u_(j, i) = z * c - y * ss;
424  }
425  }
426  }
427 
428  z = s_[k];
429  if (l == k)
430  {
431  if (z < 0.)
432  {
433  s_[k] = -z;
434  for (j = 0; j < n_; ++j)
435  {
436  v_(j, k) = -v_(j, k);
437  }
438  }
439  break;
440  }
441 
442  if (its == 29)
443  {
444  THROW_INVALID_ARGUMENT_ERROR("no convergence in 30 svdcmp iterations");
445  }
446 
447  x = s_[l];
448  nm = k - 1;
449  y = s_[nm];
450  g = rv1[nm];
451  h = rv1[k];
452  f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2. * h * y);
453  g = pythag(f, 1.);
454  f = ((x - z) * (x + z) + h * ((y / (f + SIGN(g, f))) - h)) / x;
455  c = ss = 1.;
456 
457  for (j = l; j <= nm; j++)
458  {
459  i = j + 1;
460  g = rv1[i];
461  y = s_[i];
462  h = ss * g;
463  g = c * g;
464  z = pythag(f, h);
465  rv1[j] = z;
466  c = f / z;
467  ss = h / z;
468  f = x * c + g * ss;
469  g = g * c - x * ss;
470  h = y * ss;
471  y *= c;
472 
473  for (jj = 0; jj < n_; ++jj)
474  {
475  x = v_(jj, j);
476  z = v_(jj, i);
477  v_(jj, j) = x * c + z * ss;
478  v_(jj, i) = z * c - x * ss;
479  }
480 
481  z = pythag(f, h);
482  s_[j] = z;
483 
484  if (z != 0.)
485  {
486  z = 1. / z;
487  c = f * z;
488  ss = h * z;
489  }
490 
491  f = c * g + ss * y;
492  x = c * y - ss * g;
493 
494  for (jj = 0; jj < m_; ++jj)
495  {
496  y = u_(jj, j);
497  z = u_(jj, i);
498  u_(jj, j) = y * c + z * ss;
499  u_(jj, i) = z * c - y * ss;
500  }
501  }
502  rv1[l] = 0.;
503  rv1[k] = f;
504  s_[k] = x;
505  }
506  }
507  }
508 
509  // =============================================================================
510  // Description:
513  void reorder()
514  {
515  uint32 i = 0;
516  uint32 j = 0;
517  uint32 k = 0;
518  uint32 ss = 0;
519  uint32 inc = 1;
520 
521  double sw = 0.;
522  NdArray<double> su(m_, 1);
523  NdArray<double> sv(n_, 1);
524 
525  do
526  {
527  inc *= 3;
528  ++inc;
529  } while (inc <= n_);
530 
531  do
532  {
533  inc /= 3;
534  for (i = inc; i < n_; ++i)
535  {
536  sw = s_[i];
537  for (k = 0; k < m_; ++k)
538  {
539  su[k] = u_(k, i);
540  }
541 
542  for (k = 0; k < n_; ++k)
543  {
544  sv[k] = v_(k, i);
545  }
546 
547  j = i;
548  while (s_[j - inc] < sw)
549  {
550  s_[j] = s_[j - inc];
551 
552  for (k = 0; k < m_; ++k)
553  {
554  u_(k, j) = u_(k, j - inc);
555  }
556 
557  for (k = 0; k < n_; ++k)
558  {
559  v_(k, j) = v_(k, j - inc);
560  }
561 
562  j -= inc;
563 
564  if (j < inc)
565  {
566  break;
567  }
568  }
569 
570  s_[j] = sw;
571 
572  for (k = 0; k < m_; ++k)
573  {
574  u_(k, j) = su[k];
575  }
576 
577  for (k = 0; k < n_; ++k)
578  {
579  v_(k, j) = sv[k];
580  }
581  }
582  } while (inc > 1);
583 
584  for (k = 0; k < n_; ++k)
585  {
586  ss = 0;
587 
588  for (i = 0; i < m_; i++)
589  {
590  if (u_(i, k) < 0.)
591  {
592  ss++;
593  }
594  }
595 
596  for (j = 0; j < n_; ++j)
597  {
598  if (v_(j, k) < 0.)
599  {
600  ss++;
601  }
602  }
603 
604  if (ss > (m_ + n_) / 2)
605  {
606  for (i = 0; i < m_; ++i)
607  {
608  u_(i, k) = -u_(i, k);
609  }
610 
611  for (j = 0; j < n_; ++j)
612  {
613  v_(j, k) = -v_(j, k);
614  }
615  }
616  }
617  }
618 
619  // =============================================================================
620  // Description:
628  static double pythag(double inA, double inB) noexcept
629  {
630  const double absa = std::abs(inA);
631  const double absb = std::abs(inB);
632  return (absa > absb ? absa * std::sqrt(1. + utils::sqr(absb / absa))
633  : (absb == 0. ? 0. : absb * std::sqrt(1. + utils::sqr(absa / absb))));
634  }
635 
636  private:
637  // ===============================Attributes====================================
638  const uint32 m_;
639  const uint32 n_;
640  NdArray<double> u_;
641  NdArray<double> v_;
642  NdArray<double> s_;
643  double eps_;
644  double tsh_;
645  };
646  } // namespace linalg
647 } // namespace nc
#define THROW_INVALID_ARGUMENT_ERROR(msg)
Definition: Error.hpp:36
size_type size() const noexcept
Definition: NdArrayCore.hpp:4105
const_reference front() const noexcept
Definition: NdArrayCore.hpp:2661
Definition: SVDClass.hpp:48
const NdArray< double > & u() noexcept
Definition: SVDClass.hpp:75
const NdArray< double > & s() noexcept
Definition: SVDClass.hpp:97
NdArray< double > solve(const NdArray< double > &inInput, double inThresh=-1.)
Definition: SVDClass.hpp:111
const NdArray< double > & v() noexcept
Definition: SVDClass.hpp:86
SVD(const NdArray< double > &inMatrix)
Definition: SVDClass.hpp:56
constexpr auto j
Definition: Constants.hpp:45
constexpr double c
speed of light
Definition: Constants.hpp:40
dtype f(GeneratorType &generator, dtype inDofN, dtype inDofD)
Definition: f.hpp:58
constexpr dtype sqr(dtype inValue) noexcept
Definition: sqr.hpp:44
Definition: Coordinate.hpp:45
NdArray< dtype > max(const NdArray< dtype > &inArray, Axis inAxis=Axis::NONE)
Definition: max.hpp:44
auto abs(dtype inValue) noexcept
Definition: abs.hpp:49
NdArray< dtype > min(const NdArray< dtype > &inArray, Axis inAxis=Axis::NONE)
Definition: min.hpp:44
auto sqrt(dtype inValue) noexcept
Definition: sqrt.hpp:48
Shape shape(const NdArray< dtype > &inArray) noexcept
Definition: Functions/Shape.hpp:42
std::uint32_t uint32
Definition: Types.hpp:40