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