From e96431a577356fdd51e5664c38a1bfbc7f1fc956 Mon Sep 17 00:00:00 2001 From: Relintai Date: Tue, 7 Feb 2023 23:23:48 +0100 Subject: [PATCH] Reworked svd in LinAlg. Also fixed the usage order of the parameters in a few helper methods. --- mlpp/lin_alg/lin_alg.cpp | 296 +++++++++++++++++++++++++++++++++++++-- mlpp/lin_alg/lin_alg.h | 12 +- mlpp/pca/pca.cpp | 9 +- mlpp/pca/pca.h | 5 +- 4 files changed, 295 insertions(+), 27 deletions(-) diff --git a/mlpp/lin_alg/lin_alg.cpp b/mlpp/lin_alg/lin_alg.cpp index d579b6c..55cf20d 100644 --- a/mlpp/lin_alg/lin_alg.cpp +++ b/mlpp/lin_alg/lin_alg.cpp @@ -680,6 +680,43 @@ real_t MLPPLinAlg::det(std::vector> A, int d) { return deter; } +real_t MLPPLinAlg::detm(const Ref &A, int d) { + ERR_FAIL_COND_V(!A.is_valid(), 0); + + real_t deter = 0; + Ref B; + B.instance(); + B->resize(Size2i(d, d)); + B->fill(0); + + /* This is the base case in which the input is a 2x2 square matrix. + Recursion is performed unless and until we reach this base case, + such that we recieve a scalar as the result. */ + if (d == 2) { + return A->get_element(0, 0) * A->get_element(1, 1) - A->get_element(0, 1) * A->get_element(1, 0); + } else { + for (int i = 0; i < d; i++) { + int sub_i = 0; + for (int j = 1; j < d; j++) { + int sub_j = 0; + for (int k = 0; k < d; k++) { + if (k == i) { + continue; + } + + B->set_element(sub_i, sub_j, A->get_element(j, k)); + sub_j++; + } + sub_i++; + } + + deter += Math::pow(static_cast(-1), static_cast(i)) * A->get_element(0, i) * detm(B, d - 1); + } + } + + return deter; +} + real_t MLPPLinAlg::trace(std::vector> A) { real_t trace = 0; for (int i = 0; i < A.size(); i++) { @@ -755,6 +792,76 @@ std::vector> MLPPLinAlg::pinverse(std::vector MLPPLinAlg::cofactorm(const Ref &A, int n, int i, int j) { + Ref cof; + cof.instance(); + cof->resize(A->size()); + + int sub_i = 0; + int sub_j = 0; + + for (int row = 0; row < n; row++) { + for (int col = 0; col < n; col++) { + if (row != i && col != j) { + cof->set_element(sub_i, sub_j++, A->get_element(row, col)); + + if (sub_j == n - 1) { + sub_j = 0; + sub_i++; + } + } + } + } + + return cof; +} +Ref MLPPLinAlg::adjointm(const Ref &A) { + Ref adj; + + ERR_FAIL_COND_V(!A.is_valid(), adj); + + Size2i a_size = A->size(); + + ERR_FAIL_COND_V(a_size.x != a_size.y, adj); + + //Resizing the initial adjoint matrix + + adj.instance(); + adj->resize(a_size); + + // Checking for the case where the given N x N matrix is a scalar + if (a_size.y == 1) { + adj->set_element(0, 0, 1); + return adj; + } + + if (a_size.y == 2) { + adj->set_element(0, 0, A->get_element(1, 1)); + adj->set_element(1, 1, A->get_element(0, 0)); + + adj->set_element(0, 1, -A->get_element(0, 1)); + adj->set_element(1, 0, -A->get_element(1, 0)); + + return adj; + } + + for (int i = 0; i < a_size.y; i++) { + for (int j = 0; j < a_size.x; j++) { + Ref cof = cofactorm(A, a_size.y, i, j); + // 1 if even, -1 if odd + int sign = (i + j) % 2 == 0 ? 1 : -1; + adj->set_element(j, i, sign * detm(cof, int(a_size.y) - 1)); + } + } + return adj; +} +Ref MLPPLinAlg::inversem(const Ref &A) { + return scalar_multiplym(1 / detm(A, int(A->size().y)), adjointm(A)); +} +Ref MLPPLinAlg::pinversem(const Ref &A) { + return matmultm(inversem(matmultm(transposem(A), A)), transposem(A)); +} + std::vector> MLPPLinAlg::zeromat(int n, int m) { std::vector> zeromat; zeromat.resize(n); @@ -772,7 +879,7 @@ Ref MLPPLinAlg::zeromatm(int n, int m) { Ref mat; mat.instance(); - mat->resize(Size2i(n, m)); + mat->resize(Size2i(m, n)); mat->fill(0); return mat; @@ -781,7 +888,7 @@ Ref MLPPLinAlg::onematm(int n, int m) { Ref mat; mat.instance(); - mat->resize(Size2i(n, m)); + mat->resize(Size2i(m, n)); mat->fill(1); return mat; @@ -790,7 +897,7 @@ Ref MLPPLinAlg::fullm(int n, int m, int k) { Ref mat; mat.instance(); - mat->resize(Size2i(n, m)); + mat->resize(Size2i(m, n)); mat->fill(k); return mat; @@ -938,6 +1045,21 @@ std::vector> MLPPLinAlg::identity(real_t d) { return identityMat; } +Ref MLPPLinAlg::identitym(int d) { + Ref identity_mat; + identity_mat.instance(); + identity_mat->resize(Size2i(d, d)); + identity_mat->fill(0); + + real_t *im_ptr = identity_mat->ptrw(); + + for (int i = 0; i < d; i++) { + im_ptr[identity_mat->calculate_index(i, i)] = 1; + } + + return identity_mat; +} + std::vector> MLPPLinAlg::cov(std::vector> A) { MLPPStat stat; std::vector> covMat; @@ -1198,6 +1320,146 @@ MLPPLinAlg::EigenResultOld MLPPLinAlg::eigen_old(std::vector return res; } +MLPPLinAlg::EigenResult MLPPLinAlg::eigen(Ref A) { + EigenResult res; + + ERR_FAIL_COND_V(!A.is_valid(), res); + + /* + A (the entered parameter) in most use cases will be X'X, XX', etc. and must be symmetric. + That simply means that 1) X' = X and 2) X is a square matrix. This function that computes the + eigenvalues of a matrix is utilizing Jacobi's method. + */ + + real_t diagonal = true; // Perform the iterative Jacobi algorithm unless and until we reach a diagonal matrix which yields us the eigenvals. + + HashMap val_to_vec; + Ref a_new; + Ref eigenvectors = identitym(A->size().y); + Size2i a_size = A->size(); + + do { + real_t a_ij = A->get_element(0, 1); + real_t sub_i = 0; + real_t sub_j = 1; + for (int i = 0; i < a_size.y; ++i) { + for (int j = 0; j < a_size.x; ++j) { + real_t ca_ij = A->get_element(i, j); + real_t abs_ca_ij = ABS(ca_ij); + + if (i != j && abs_ca_ij > a_ij) { + a_ij = ca_ij; + sub_i = i; + sub_j = j; + } else if (i != j && abs_ca_ij == a_ij) { + if (i < sub_i) { + a_ij = ca_ij; + sub_i = i; + sub_j = j; + } + } + } + } + + real_t a_ii = A->get_element(sub_i, sub_i); + real_t a_jj = A->get_element(sub_j, sub_j); + real_t a_ji = A->get_element(sub_j, sub_i); + real_t theta; + + if (a_ii == a_jj) { + theta = M_PI / 4; + } else { + theta = 0.5 * atan(2 * a_ij / (a_ii - a_jj)); + } + + Ref P = identitym(A->size().y); + P->set_element(sub_i, sub_j, -Math::sin(theta)); + P->set_element(sub_i, sub_i, Math::cos(theta)); + P->set_element(sub_j, sub_j, Math::cos(theta)); + P->set_element(sub_j, sub_i, Math::sin(theta)); + + a_new = matmultm(matmultm(inversem(P), A), P); + + Size2i a_new_size = a_new->size(); + + for (int i = 0; i < a_new_size.y; ++i) { + for (int j = 0; j < a_new_size.x; ++j) { + if (i != j && Math::round(a_new->get_element(i, j)) == 0) { + a_new->set_element(i, j, 0); + } + } + } + + bool non_zero = false; + for (int i = 0; i < a_new_size.y; ++i) { + for (int j = 0; j < a_new_size.x; ++j) { + if (i != j && Math::round(a_new->get_element(i, j)) == 0) { + non_zero = true; + } + } + } + + if (non_zero) { + diagonal = false; + } else { + diagonal = true; + } + + if (a_new == A) { + diagonal = true; + for (int i = 0; i < a_new_size.y; ++i) { + for (int j = 0; j < a_new_size.x; ++j) { + if (i != j) { + a_new->set_element(i, j, 0); + } + } + } + } + + eigenvectors = matmultm(eigenvectors, P); + A = a_new; + + } while (!diagonal); + + Ref a_new_prior = a_new; + + Size2i a_new_size = a_new->size(); + + // Bubble Sort. Should change this later. + for (int i = 0; i < a_new_size.y - 1; ++i) { + for (int j = 0; j < a_new_size.x - 1 - i; ++j) { + if (a_new->get_element(j, j) < a_new->get_element(j + 1, j + 1)) { + real_t temp = a_new->get_element(j + 1, j + 1); + a_new->set_element(j + 1, j + 1, a_new->get_element(j, j)); + a_new->set_element(j, j, temp); + } + } + } + + for (int i = 0; i < a_new_size.y; ++i) { + for (int j = 0; j < a_new_size.x; ++j) { + if (a_new->get_element(i, i) == a_new_prior->get_element(j, j)) { + val_to_vec[i] = j; + } + } + } + + Ref eigen_temp = eigenvectors->duplicate(); + + Size2i eigenvectors_size = eigenvectors->size(); + + for (int i = 0; i < eigenvectors_size.y; ++i) { + for (int j = 0; j < eigenvectors_size.x; ++j) { + eigenvectors->set_element(i, j, eigen_temp->get_element(i, val_to_vec[j])); + } + } + + res.eigen_vectors = eigenvectors; + res.eigen_values = a_new; + + return res; +} + MLPPLinAlg::SDVResultOld MLPPLinAlg::SVD(std::vector> A) { EigenResultOld left_eigen = eigen_old(matmult(A, transpose(A))); EigenResultOld right_eigen = eigen_old(matmult(transpose(A), A)); @@ -1219,27 +1481,31 @@ MLPPLinAlg::SDVResultOld MLPPLinAlg::SVD(std::vector> A) { } MLPPLinAlg::SDVResult MLPPLinAlg::svd(const Ref &A) { - /* - EigenResultOld left_eigen = eigen(matmult(A, transpose(A))); - EigenResultOld right_eigen = eigen(matmult(transpose(A), A)); + SDVResult res; - std::vector> singularvals = sqrt(left_eigen.eigen_values); - std::vector> sigma = zeromat(A.size(), A[0].size()); - for (int i = 0; i < singularvals.size(); i++) { - for (int j = 0; j < singularvals[i].size(); j++) { - sigma[i][j] = singularvals[i][j]; + ERR_FAIL_COND_V(!A.is_valid(), res); + + Size2i a_size = A->size(); + + EigenResult left_eigen = eigen(matmultm(A, transposem(A))); + EigenResult right_eigen = eigen(matmultm(transposem(A), A)); + + Ref singularvals = sqrtm(left_eigen.eigen_values); + Ref sigma = zeromatm(a_size.y, a_size.x); + + Size2i singularvals_size = singularvals->size(); + + for (int i = 0; i < singularvals_size.y; ++i) { + for (int j = 0; j < singularvals_size.x; ++j) { + sigma->set_element(i, j, singularvals->get_element(i, j)); } } - SDVResult res; res.U = left_eigen.eigen_vectors; res.S = sigma; res.Vt = right_eigen.eigen_vectors; return res; - */ - - return SDVResult(); } std::vector MLPPLinAlg::vectorProjection(std::vector a, std::vector b) { diff --git a/mlpp/lin_alg/lin_alg.h b/mlpp/lin_alg/lin_alg.h index 6bbc466..8782d7b 100644 --- a/mlpp/lin_alg/lin_alg.h +++ b/mlpp/lin_alg/lin_alg.h @@ -76,6 +76,7 @@ public: Ref absm(const Ref &A); real_t det(std::vector> A, int d); + real_t detm(const Ref &A, int d); real_t trace(std::vector> A); @@ -84,6 +85,11 @@ public: std::vector> inverse(std::vector> A); std::vector> pinverse(std::vector> A); + Ref cofactorm(const Ref &A, int n, int i, int j); + Ref adjointm(const Ref &A); + Ref inversem(const Ref &A); + Ref pinversem(const Ref &A); + std::vector> zeromat(int n, int m); std::vector> onemat(int n, int m); std::vector> full(int n, int m, int k); @@ -109,6 +115,7 @@ public: real_t norm_2(std::vector> A); std::vector> identity(real_t d); + Ref identitym(int d); std::vector> cov(std::vector> A); @@ -121,14 +128,13 @@ public: EigenResultOld eigen_old(std::vector> A); -/* struct EigenResult { Ref eigen_vectors; Ref eigen_values; }; - EigenResult eigen(const Ref &A); -*/ + EigenResult eigen(Ref A); + struct SDVResultOld { std::vector> U; std::vector> S; diff --git a/mlpp/pca/pca.cpp b/mlpp/pca/pca.cpp index 039bda9..ce8914d 100644 --- a/mlpp/pca/pca.cpp +++ b/mlpp/pca/pca.cpp @@ -11,12 +11,6 @@ #include #include - - -MLPPPCA::MLPPPCA(std::vector> inputSet, int k) : - inputSet(inputSet), k(k) { -} - std::vector> MLPPPCA::principalComponents() { MLPPLinAlg alg; MLPPData data; @@ -52,3 +46,6 @@ real_t MLPPPCA::score() { return 1 - num / den; } +MLPPPCA::MLPPPCA(std::vector> inputSet, int k) : + inputSet(inputSet), k(k) { +} diff --git a/mlpp/pca/pca.h b/mlpp/pca/pca.h index 36e9d5d..f5e43c7 100644 --- a/mlpp/pca/pca.h +++ b/mlpp/pca/pca.h @@ -12,13 +12,13 @@ #include - class MLPPPCA { public: - MLPPPCA(std::vector> inputSet, int k); std::vector> principalComponents(); real_t score(); + MLPPPCA(std::vector> inputSet, int k); + private: std::vector> inputSet; std::vector> X_normalized; @@ -27,5 +27,4 @@ private: int k; }; - #endif /* PCA_hpp */