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 <numeric>
      36             : #include <openGPMP/linalg/eigen.hpp>
      37             : #include <vector>
      38             : 
      39           7 : gpmp::linalg::Eigen::Eigen(const std::vector<std::vector<double>> &mat)
      40           7 :     : matrix(mat), size(mat.size()) {
      41             :     // Check if the matrix is square
      42             :     /*if (size != mat[0].size()) {
      43             :         throw std::invalid_argument(
      44             :             "Error: Eigenvalue computation requires a square matrix.");
      45             :     }*/
      46           7 : }
      47             : 
      48           1 : double gpmp::linalg::Eigen::power_iter(int max_iters, double tolerance) const {
      49             :     // initial guess
      50           1 :     std::vector<double> x(size, 1.0);
      51           1 :     std::vector<double> y(size);
      52             : 
      53           7 :     for (int iter = 0; iter < max_iters; ++iter) {
      54             :         // perform matrix-vector multiplication: y = A * x
      55          28 :         for (int i = 0; i < size; ++i) {
      56          21 :             y[i] = 0.0;
      57          84 :             for (int j = 0; j < size; ++j) {
      58          63 :                 y[i] += matrix[i][j] * x[j];
      59             :             }
      60             :         }
      61             : 
      62             :         // find the maximum element in y
      63           7 :         double max_element = 0.0;
      64          28 :         for (int i = 0; i < size; ++i) {
      65          21 :             max_element = std::max(max_element, std::abs(y[i]));
      66             :         }
      67             : 
      68             :         // normalize y
      69          28 :         for (int i = 0; i < size; ++i) {
      70          21 :             y[i] /= max_element;
      71             :         }
      72             : 
      73             :         // check for convergence
      74           7 :         double error = 0.0;
      75          28 :         for (int i = 0; i < size; ++i) {
      76          21 :             error += std::abs(y[i] - x[i]);
      77             :         }
      78             : 
      79           7 :         if (error < tolerance) {
      80             :             // eigenvalue found
      81           1 :             return max_element;
      82             :         }
      83             : 
      84           6 :         x = y;
      85             :     }
      86             : 
      87             :     std::cerr << "Warning: Power iteration did not converge within the"
      88           0 :                  "specified iterations\n";
      89             :     // not converged
      90           0 :     return 0.0;
      91           1 : }
      92             : 
      93           6 : std::vector<double> gpmp::linalg::Eigen::qr_algorithm(int max_iters,
      94             :                                                       double tolerance) const {
      95             :     // initialize a copy of the matrix to avoid modifying the original
      96           6 :     std::vector<std::vector<double>> hessenberg_mtx = matrix;
      97             : 
      98             :     // QR Iteration
      99           6 :     for (int iter = 0; iter < max_iters; ++iter) {
     100             :         // QR decomposition of the Hessenberg matrix
     101          18 :         for (int i = 0; i < size - 1; ++i) {
     102             :             // perform Givens rotation to introduce zeros below the subdiagonal
     103          12 :             double a = hessenberg_mtx[i][i];
     104          12 :             double b = hessenberg_mtx[i + 1][i];
     105          12 :             double r = std::hypot(a, b);
     106          12 :             double c = a / r;
     107          12 :             double s = b / r;
     108             : 
     109             :             // apply Givens rotation to the submatrix
     110          48 :             for (int j = 0; j < size; ++j) {
     111             :                 double temp1 =
     112          36 :                     c * hessenberg_mtx[i][j] + s * hessenberg_mtx[i + 1][j];
     113             :                 double temp2 =
     114          36 :                     -s * hessenberg_mtx[i][j] + c * hessenberg_mtx[i + 1][j];
     115          36 :                 hessenberg_mtx[i][j] = temp1;
     116          36 :                 hessenberg_mtx[i + 1][j] = temp2;
     117             :             }
     118             :         }
     119             : 
     120             :         // check for convergence (off-diagonal elements close to zero)
     121           6 :         double offDiagonalSum = 0.0;
     122          12 :         for (int i = 2; i < size; ++i) {
     123           6 :             offDiagonalSum += std::abs(hessenberg_mtx[i][i - 1]);
     124             :         }
     125             : 
     126           6 :         if (offDiagonalSum < tolerance) {
     127             :             // extract eigenvalues from the diagonal of the Hessenberg matrix
     128           6 :             std::vector<double> eigenvalues;
     129          24 :             for (int i = 0; i < size; ++i) {
     130          18 :                 eigenvalues.push_back(hessenberg_mtx[i][i]);
     131             :             }
     132             : 
     133           6 :             return eigenvalues;
     134           6 :         }
     135             :     }
     136             : 
     137             :     std::cerr << "Warning: QR iteration did not converge within the "
     138           0 :                  "specified iterations\n";
     139             :     // empty vector if not converged
     140           0 :     return std::vector<double>();
     141           6 : }
     142             : 
     143           1 : double gpmp::linalg::Eigen::determinant() const {
     144             :     // compute the determinant as the product of eigenvalues
     145           1 :     std::vector<double> eigenvalues = qr_algorithm();
     146           1 :     double determinant = std::accumulate(eigenvalues.begin(),
     147             :                                          eigenvalues.end(),
     148             :                                          1.0,
     149             :                                          std::multiplies<double>());
     150           1 :     return determinant;
     151           1 : }
     152             : 
     153             : std::vector<double>
     154           1 : gpmp::linalg::Eigen::sensitivity(double perturbation) const {
     155           1 :     std::vector<double> originalEigenvalues = qr_algorithm();
     156             : 
     157           1 :     std::vector<double> perturbedEigenvalues;
     158           4 :     for (int i = 0; i < size; ++i) {
     159             :         // perturb the matrix
     160           3 :         std::vector<std::vector<double>> perturbedMatrix = matrix;
     161           3 :         perturbedMatrix[i][i] += perturbation;
     162             : 
     163             :         // calculate eigenvalues of perturbed matrix
     164             :         std::vector<double> perturbedVals =
     165           3 :             Eigen(perturbedMatrix).qr_algorithm();
     166           3 :         perturbedEigenvalues.push_back(perturbedVals[i]);
     167           3 :     }
     168             : 
     169             :     // calculate sensitivity: (perturbed eigenvalue - original eigenvalue) /
     170             :     // perturbation
     171           1 :     std::vector<double> sensitivity;
     172           4 :     for (int i = 0; i < size; ++i) {
     173             :         double sens =
     174           3 :             (perturbedEigenvalues[i] - originalEigenvalues[i]) / perturbation;
     175           3 :         sensitivity.push_back(sens);
     176             :     }
     177             : 
     178           2 :     return sensitivity;
     179           1 : }
     180             : 
     181           0 : void gpmp::linalg::Eigen::schur_decomp(
     182             :     std::vector<std::vector<double>> &schurMatrix,
     183             :     double tolerance) const {
     184             :     // Hessenberg matrix is already calculated in QR iteration
     185           0 :     schurMatrix = matrix;
     186             : 
     187           0 :     for (int i = size - 1; i > 0; --i) {
     188             :         // check for small subdiagonal elements to introduce a zero
     189           0 :         if (std::abs(schurMatrix[i][i - 1]) < tolerance) {
     190           0 :             schurMatrix[i][i - 1] = 0.0;
     191             :         } else {
     192             :             // Apply Givens rotation to introduce a zero below the subdiagonal
     193           0 :             double a = schurMatrix[i - 1][i - 1];
     194           0 :             double b = schurMatrix[i][i - 1];
     195           0 :             double r = std::hypot(a, b);
     196           0 :             double c = a / r;
     197           0 :             double s = b / r;
     198             : 
     199           0 :             for (int j = 0; j < size; ++j) {
     200             :                 double temp1 =
     201           0 :                     c * schurMatrix[i - 1][j] + s * schurMatrix[i][j];
     202             :                 double temp2 =
     203           0 :                     -s * schurMatrix[i - 1][j] + c * schurMatrix[i][j];
     204           0 :                 schurMatrix[i - 1][j] = temp1;
     205           0 :                 schurMatrix[i][j] = temp2;
     206             :             }
     207             :         }
     208             :     }
     209           0 : }
     210             : 
     211             : std::vector<std::vector<double>>
     212           0 : gpmp::linalg::Eigen::jordan_normal_form(double tolerance) const {
     213             :     std::vector<std::vector<double>> jordanMatrix(
     214           0 :         size,
     215           0 :         std::vector<double>(size, 0.0));
     216           0 :     std::vector<double> eigenvalues = qr_algorithm();
     217             : 
     218           0 :     for (int i = 0; i < size; ++i) {
     219           0 :         jordanMatrix[i][i] = eigenvalues[i];
     220             : 
     221             :         // find the size of the Jordan block
     222           0 :         int jordanBlockSize = 1;
     223           0 :         while (i + jordanBlockSize < size &&
     224           0 :                std::abs(matrix[i + jordanBlockSize][i]) < tolerance) {
     225           0 :             ++jordanBlockSize;
     226             :         }
     227             : 
     228             :         // fill the Jordan block
     229           0 :         for (int j = 1; j < jordanBlockSize; ++j) {
     230           0 :             jordanMatrix[i + j][i + j - 1] = 1.0;
     231             :         }
     232             :     }
     233             : 
     234           0 :     return jordanMatrix;
     235           0 : }
     236             : 
     237           0 : double gpmp::linalg::Eigen::rayleigh_iter(int maxIterations,
     238             :                                           double tolerance) const {
     239             :     // initial guess
     240           0 :     std::vector<double> x(size, 1.0);
     241           0 :     double lambdaPrev = 0.0;
     242             : 
     243           0 :     for (int iter = 0; iter < maxIterations; ++iter) {
     244             :         // Rayleigh quotient iteration
     245           0 :         std::vector<double> Ax(size, 0.0);
     246             : 
     247             :         // perform matrix-vector multiplication: Ax = A * x
     248           0 :         for (int i = 0; i < size; ++i) {
     249           0 :             for (int j = 0; j < size; ++j) {
     250           0 :                 Ax[i] += matrix[i][j] * x[j];
     251             :             }
     252             :         }
     253             : 
     254             :         // compute the Rayleigh quotient
     255             :         double lambda =
     256           0 :             std::inner_product(x.begin(), x.end(), Ax.begin(), 0.0) /
     257           0 :             std::inner_product(x.begin(), x.end(), x.begin(), 0.0);
     258             : 
     259             :         // check for convergence
     260           0 :         if (std::abs(lambda - lambdaPrev) < tolerance) {
     261           0 :             return lambda;
     262             :         }
     263             : 
     264             :         // update the guess vector
     265           0 :         double normAx = std::sqrt(
     266             :             std::inner_product(Ax.begin(), Ax.end(), Ax.begin(), 0.0));
     267           0 :         for (int i = 0; i < size; ++i) {
     268             :             // normalize each element of x by dividing by its magnitude
     269           0 :             x[i] = Ax[i] / normAx;
     270             :         }
     271             : 
     272           0 :         lambdaPrev = lambda;
     273           0 :     }
     274             : 
     275             :     std::cerr << "Warning: Rayleigh quotient iteration did not converge within "
     276           0 :                  "the specified iterations."
     277           0 :               << std::endl;
     278             :     // if not converged
     279           0 :     return 0.0;
     280           0 : }