philsupertramp/game-math
Loading...
Searching...
No Matches
svd.h
Go to the documentation of this file.
1#pragma once
2
3#include "../../Matrix.h"
4#include "gaussJordan.h"
5#include "qr.h"
6
14std::vector<Matrix<double>> svd(const Matrix<double>& A, const size_t& k, const double epsilon = 0.1e-4) {
15 // works with squared and not squared matrizes
16 size_t n_orig = A.rows();
17 size_t m_orig = A.columns();
18 size_t chosen_k = k;
19 if(k == 0) { chosen_k = n_orig < m_orig ? n_orig : m_orig; }
20 auto A_copy = A;
21 auto A_orig = A;
22 size_t n, m;
23 if(n_orig > m_orig) {
24 A_copy = A_copy.Transpose() * A_copy;
25 n = A_copy.rows();
26 m = A_copy.columns();
27 } else if(n_orig < m_orig) {
28 A_copy = A_copy * A_copy.Transpose();
29 n = A_copy.rows();
30 m = A_copy.columns();
31 } else {
32 n = n_orig;
33 m = m_orig;
34 }
35
36 auto Q = Matrix<double>::Random(n, chosen_k);
37 auto R = Matrix<double>::Random(n, chosen_k);
38 auto res = qr(Q);
39 Q = res.first;
40 auto Q_prev = Q;
41
42 for(size_t i = 0; i < 1000; i++) {
43 auto Z = A_copy * Q;
44 auto _res = qr(Z);
45 Q = _res.first;
46 R = _res.second;
47 auto err = ((Q - Q_prev)).Apply([](float val) { return val * val; }).sumElements();
48 Q_prev = Q;
49 if(err < epsilon) { break; }
50 }
51 auto singular_values = diag_elements(R).Apply([](double val) { return sqrt(val); });
52 Matrix<double> left_vecs, right_vecs;
53 if(n_orig < m_orig) {
54 left_vecs = Q.Transpose();
55 right_vecs = (gaussJordan(diag(singular_values)) * left_vecs.Transpose()) * A;
56 } else if(n_orig == m_orig) {
57 left_vecs = Q.Transpose();
58 right_vecs = left_vecs;
59 singular_values = singular_values.Apply([](double val) { return val * val; });
60 } else {
61 right_vecs = Q.Transpose();
62 left_vecs = A * (right_vecs.Transpose() * gaussJordan(diag(singular_values)));
63 }
64
65 return { left_vecs, singular_values, right_vecs };
66}
67
Definition: Matrix.h:42
constexpr Matrix< T > Transpose() const
Definition: Matrix.h:256
size_t rows() const
Definition: Matrix.h:193
size_t columns() const
Definition: Matrix.h:198
Matrix< T > Apply(const std::function< T(T)> &fun) const
Definition: Matrix.h:375
static Matrix Random(size_t rows, size_t columns, size_t element_size=1, double minValue=0.0, double maxValue=1.0)
Definition: Matrix.h:158
Matrix< double > gaussJordan(const Matrix< double > &A)
Definition: gaussJordan.h:27
Matrix< T > diag_elements(const Matrix< T > &in)
Definition: matrix_utils.h:451
Matrix< T > diag(const Matrix< T > &in)
Definition: utils.h:135
std::pair< Matrix< double >, Matrix< double > > qr(const Matrix< double > &A)
Definition: qr.h:12
std::vector< Matrix< double > > svd(const Matrix< double > &A, const size_t &k, const double epsilon=0.1e-4)
Definition: svd.h:14