Line data    Source code

       1             : /*************************************************************************
       2             :  *
       3             :  *  Project
       4             :  *                         _____ _____  __  __ _____
       5             :  *                        / ____|  __ \|  \/  |  __ \
       6             :  *  ___  _ __   ___ _ __ | |  __| |__) | \  / | |__) |
       7             :  * / _ \| '_ \ / _ \ '_ \| | |_ |  ___/| |\/| |  ___/
       8             :  *| (_) | |_) |  __/ | | | |__| | |    | |  | | |
       9             :  * \___/| .__/ \___|_| |_|\_____|_|    |_|  |_|_|
      10             :  *      | |
      11             :  *      |_|
      12             :  *
      13             :  * Copyright (C) Akiel Aries, <akiel@akiel.org>, et al.
      14             :  *
      15             :  * This software is licensed as described in the file LICENSE, which
      16             :  * you should have received as part of this distribution. The terms
      17             :  * among other details are referenced in the official documentation
      18             :  * seen here : https://akielaries.github.io/openGPMP/ along with
      19             :  * important files seen in this project.
      20             :  *
      21             :  * You may opt to use, copy, modify, merge, publish, distribute
      22             :  * and/or sell copies of the Software, and permit persons to whom
      23             :  * the Software is furnished to do so, under the terms of the
      24             :  * LICENSE file. As this is an Open Source effort, all implementations
      25             :  * must be of the same methodology.
      26             :  *
      27             :  *
      28             :  *
      29             :  * This software is distributed on an AS IS basis, WITHOUT
      30             :  * WARRANTY OF ANY KIND, either express or implied.
      31             :  *
      32             :  ************************************************************************/
      33             : #include <cmath>
      34             : #include <iostream>
      35             : #include <openGPMP/linalg/svd.hpp>
      36             : #include <stdexcept>
      37             : 
      38           0 : gpmp::linalg::SVD::SVD(std::vector<std::vector<double>> &matrix) {
      39           0 :     if (matrix.empty() || matrix[0].empty()) {
      40           0 :         throw std::invalid_argument("Error: Matrix for SVD is empty.");
      41             :     }
      42             : 
      43           0 :     computeSVD(matrix);
      44           0 : }
      45             : 
      46           0 : std::vector<double> gpmp::linalg::SVD::getSingularValues() const {
      47           0 :     return singularValues_;
      48             : }
      49             : 
      50             : std::vector<std::vector<double>>
      51           0 : gpmp::linalg::SVD::getLeftSingularVectors() const {
      52           0 :     return leftSingularVectors_;
      53             : }
      54             : 
      55             : std::vector<std::vector<double>>
      56           0 : gpmp::linalg::SVD::getRightSingularVectors() const {
      57           0 :     return rightSingularVectors_;
      58             : }
      59             : 
      60           0 : void gpmp::linalg::SVD::computeSVD(std::vector<std::vector<double>> &matrix) {
      61           0 :     size_t m = matrix.size();
      62           0 :     size_t n = matrix[0].size();
      63             : 
      64           0 :     singularValues_.resize(std::min(m, n));
      65             :     leftSingularVectors_ =
      66           0 :         std::vector<std::vector<double>>(m, std::vector<double>(m, 0.0));
      67             :     rightSingularVectors_ =
      68           0 :         std::vector<std::vector<double>>(n, std::vector<double>(n, 0.0));
      69             : 
      70           0 :     std::vector<double> superdiagonal(n - 1, 0.0);
      71             : 
      72           0 :     bidiagonalize(const_cast<std::vector<std::vector<double>> &>(matrix),
      73           0 :                   singularValues_,
      74             :                   superdiagonal);
      75             : 
      76             :     // left/right singular vectors as identity matrices
      77           0 :     for (size_t i = 0; i < m; ++i) {
      78           0 :         leftSingularVectors_[i][i] = 1.0;
      79             :     }
      80             : 
      81           0 :     for (size_t i = 0; i < n; ++i) {
      82           0 :         rightSingularVectors_[i][i] = 1.0;
      83             :     }
      84           0 : }
      85             : 
      86           0 : void gpmp::linalg::SVD::bidiagonalize(std::vector<std::vector<double>> &matrix,
      87             :                                       std::vector<double> &diagonal,
      88             :                                       std::vector<double> &superdiagonal) {
      89           0 :     size_t m = matrix.size();
      90           0 :     size_t n = matrix[0].size();
      91             : 
      92           0 :     diagonal.resize(std::min(m, n));
      93           0 :     superdiagonal.resize(std::min(m, n) - 1);
      94             : 
      95           0 :     for (size_t k = 0; k < std::min(m, n); ++k) {
      96           0 :         double alpha = 0.0;
      97           0 :         for (size_t i = k; i < m; ++i) {
      98           0 :             alpha += matrix[i][k] * matrix[i][k];
      99             :         }
     100           0 :         alpha = (matrix[k][k] > 0) ? -std::sqrt(alpha) : std::sqrt(alpha);
     101             : 
     102           0 :         double beta = alpha * (alpha - matrix[k][k]);
     103           0 :         diagonal[k] = -alpha;
     104             : 
     105             :         // make a copy of the current column of the matrix
     106           0 :         std::vector<double> columnCopy(m - k);
     107           0 :         for (size_t i = k; i < m; ++i) {
     108           0 :             columnCopy[i - k] = matrix[i][k];
     109             :         }
     110             : 
     111           0 :         computeHouseholderReflection(columnCopy, superdiagonal, beta);
     112           0 :         applyHouseholder(matrix, superdiagonal, beta, k, k + 1);
     113           0 :     }
     114           0 : }
     115             : 
     116           0 : void gpmp::linalg::SVD::computeHouseholderReflection(
     117             :     const std::vector<double> &x,
     118             :     std::vector<double> &v,
     119             :     double &beta) {
     120           0 :     size_t n = x.size();
     121           0 :     double sigma = 0.0;
     122             : 
     123           0 :     for (size_t i = 1; i < n; ++i) {
     124           0 :         sigma += x[i] * x[i];
     125             :     }
     126             : 
     127           0 :     v.resize(n);
     128           0 :     v[0] = 1.0;
     129           0 :     for (size_t i = 1; i < n; ++i) {
     130           0 :         v[i] = x[i] / sigma;
     131             :     }
     132             : 
     133           0 :     beta = 2.0 / (sigma + 1.0);
     134           0 : }
     135             : 
     136           0 : void gpmp::linalg::SVD::applyHouseholder(
     137             :     std::vector<std::vector<double>> &matrix,
     138             :     const std::vector<double> &v,
     139             :     double beta,
     140             :     size_t fromRow,
     141             :     size_t fromCol) {
     142           0 :     size_t m = matrix.size();
     143           0 :     size_t n = matrix[0].size();
     144             : 
     145           0 :     for (size_t i = fromRow; i < m; ++i) {
     146           0 :         double dotProduct = 0.0;
     147           0 :         for (size_t j = fromCol; j < n; ++j) {
     148           0 :             dotProduct += matrix[i][j] * v[j - fromCol];
     149             :         }
     150             : 
     151           0 :         dotProduct *= beta;
     152             : 
     153           0 :         for (size_t j = fromCol; j < n; ++j) {
     154           0 :             matrix[i][j] -= dotProduct * v[j - fromCol];
     155             :         }
     156             :     }
     157           0 : }