Class for computing eigenvalues of a matrix. More...

#include <eigen.hpp>

Public Member Functions
	Eigen (const std::vector< std::vector< double >> &mat)
	Constructor for the Eigen class. More...

double	power_iter (int max_iters=100, double tolerance=1e-10) const
	Perform the power iteration method to find the dominant eigenvalue. More...

std::vector< double >	qr_algorithm (int max_iters=100, double tolerance=1e-10) const
	Calculate all eigenvalues of the matrix using the QR algorithm. More...

std::vector< double >	sensitivity (double perturbation) const
	Calculates the sensitivity of the eigenvalues of a matrix to small perturbations. More...

double	determinant () const

void	schur_decomp (std::vector< std::vector< double >> &schurMatrix, double tolerance=1e-10) const
	Computes the Schur decomposition of a matrix. More...

std::vector< std::vector< double > >	jordan_normal_form (double tolerance=1e-10) const
	Computes the Jordan normal form of a matrix. More...

double	rayleigh_iter (int maxIterations, double tolerance) const
	Perform Rayleigh Quotient Iteration to approximate the dominant eigenvalue. More...

Private Attributes
std::vector< std::vector< double > >	matrix

int	size

Detailed Description

Class for computing eigenvalues of a matrix.

Definition at line 45 of file eigen.hpp.

Constructor & Destructor Documentation

◆ Eigen()

gpmp::linalg::Eigen::Eigen ( const std::vector< std::vector< double >> & mat )

Constructor for the Eigen class.

Parameters

mat	The square matrix for eigenvalue computation

Definition at line 39 of file eigen.cpp.

     : matrix(mat), size(mat.size()) {
     // Check if the matrix is square
     /*if (size != mat[0].size()) {
         throw std::invalid_argument(
             "Error: Eigenvalue computation requires a square matrix.");
     }*/
 }

Member Function Documentation

◆ determinant()

double gpmp::linalg::Eigen::determinant ( ) const

Definition at line 143 of file eigen.cpp.

                                           {
     // compute the determinant as the product of eigenvalues
     std::vector<double> eigenvalues = qr_algorithm();
     double determinant = std::accumulate(eigenvalues.begin(),
                                          eigenvalues.end(),
                                          1.0,
                                          std::multiplies<double>());
     return determinant;
 }

◆ jordan_normal_form()

std::vector< std::vector< double > > gpmp::linalg::Eigen::jordan_normal_form ( double tolerance = 1e-10 ) const

Computes the Jordan normal form of a matrix.

This method computes the Jordan normal form of a matrix The Jordan normal form is a block diagonal matrix, where each block is either a diagonal matrix or a Jordan block A Jordan block is a square matrix with ones on the superdiagonal and the eigenvalue on the diagonal

Parameters

tolerance The tolerance for determining if a subdiagonal element is zero

Returns: The Jordan normal form of the matrix

The returned matrix will have the same dimensions as the original matrix

Definition at line 212 of file eigen.cpp.

                                                           {
     std::vector<std::vector<double>> jordanMatrix(
         size,
         std::vector<double>(size, 0.0));
     std::vector<double> eigenvalues = qr_algorithm();
  
     for (int i = 0; i < size; ++i) {
         jordanMatrix[i][i] = eigenvalues[i];
  
         // find the size of the Jordan block
         int jordanBlockSize = 1;
         while (i + jordanBlockSize < size &&
                std::abs(matrix[i + jordanBlockSize][i]) < tolerance) {
             ++jordanBlockSize;
         }
  
         // fill the Jordan block
         for (int j = 1; j < jordanBlockSize; ++j) {
             jordanMatrix[i + j][i + j - 1] = 1.0;
         }
     }
  
     return jordanMatrix;
 }

◆ power_iter()

double gpmp::linalg::Eigen::power_iter	(	int	max_iters = `100`,
		double	tolerance = `1e-10`
	)		const

Perform the power iteration method to find the dominant eigenvalue.

Parameters

max_iters	Maximum number of iterations for power iteration
tolerance	Tolerance for convergence

Returns: The dominant eigenvalue of the matrix

Definition at line 48 of file eigen.cpp.

                                                                         {
     // initial guess
     std::vector<double> x(size, 1.0);
     std::vector<double> y(size);
  
     for (int iter = 0; iter < max_iters; ++iter) {
         // perform matrix-vector multiplication: y = A * x
         for (int i = 0; i < size; ++i) {
             y[i] = 0.0;
             for (int j = 0; j < size; ++j) {
                 y[i] += matrix[i][j] * x[j];
             }
         }
  
         // find the maximum element in y
         double max_element = 0.0;
         for (int i = 0; i < size; ++i) {
             max_element = std::max(max_element, std::abs(y[i]));
         }
  
         // normalize y
         for (int i = 0; i < size; ++i) {
             y[i] /= max_element;
         }
  
         // check for convergence
         double error = 0.0;
         for (int i = 0; i < size; ++i) {
             error += std::abs(y[i] - x[i]);
         }
  
         if (error < tolerance) {
             // eigenvalue found
             return max_element;
         }
  
         x = y;
     }
  
     std::cerr << "Warning: Power iteration did not converge within the"
                  "specified iterations\n";
     // not converged
     return 0.0;
 }

Referenced by main().

◆ qr_algorithm()

std::vector< double > gpmp::linalg::Eigen::qr_algorithm	(	int	max_iters = `100`,
		double	tolerance = `1e-10`
	)		const

Calculate all eigenvalues of the matrix using the QR algorithm.

Parameters

max_iters	Maximum number of iterations for QR algorithm
tolerance	Tolerance for convergence

Returns: A vector containing all eigenvalues of the matrix

Definition at line 93 of file eigen.cpp.

                                                                               {
     // initialize a copy of the matrix to avoid modifying the original
     std::vector<std::vector<double>> hessenberg_mtx = matrix;
  
     // QR Iteration
     for (int iter = 0; iter < max_iters; ++iter) {
         // QR decomposition of the Hessenberg matrix
         for (int i = 0; i < size - 1; ++i) {
             // perform Givens rotation to introduce zeros below the subdiagonal
             double a = hessenberg_mtx[i][i];
             double b = hessenberg_mtx[i + 1][i];
             double r = std::hypot(a, b);
             double c = a / r;
             double s = b / r;
  
             // apply Givens rotation to the submatrix
             for (int j = 0; j < size; ++j) {
                 double temp1 =
                     c * hessenberg_mtx[i][j] + s * hessenberg_mtx[i + 1][j];
                 double temp2 =
                     -s * hessenberg_mtx[i][j] + c * hessenberg_mtx[i + 1][j];
                 hessenberg_mtx[i][j] = temp1;
                 hessenberg_mtx[i + 1][j] = temp2;
             }
         }
  
         // check for convergence (off-diagonal elements close to zero)
         double offDiagonalSum = 0.0;
         for (int i = 2; i < size; ++i) {
             offDiagonalSum += std::abs(hessenberg_mtx[i][i - 1]);
         }
  
         if (offDiagonalSum < tolerance) {
             // extract eigenvalues from the diagonal of the Hessenberg matrix
             std::vector<double> eigenvalues;
             for (int i = 0; i < size; ++i) {
                 eigenvalues.push_back(hessenberg_mtx[i][i]);
             }
  
             return eigenvalues;
         }
     }
  
     std::cerr << "Warning: QR iteration did not converge within the "
                  "specified iterations\n";
     // empty vector if not converged
     return std::vector<double>();
 }

Referenced by sensitivity().

◆ rayleigh_iter()

double gpmp::linalg::Eigen::rayleigh_iter	(	int	maxIterations,
		double	tolerance
	)		const

Perform Rayleigh Quotient Iteration to approximate the dominant eigenvalue.

This method iteratively refines an estimate of the dominant eigenvalue of the matrix using the Rayleigh quotient iteration. The method returns the final estimate of the dominant eigenvalue

Parameters

maxIterations	Maximum number of iterations
tolerance	Convergence tolerance

Returns: The dominant eigenvalue estimate

Note: The iteration may not converge for all matrices. A warning is issued if convergence is not achieved

Definition at line 237 of file eigen.cpp.

                                                                   {
     // initial guess
     std::vector<double> x(size, 1.0);
     double lambdaPrev = 0.0;
  
     for (int iter = 0; iter < maxIterations; ++iter) {
         // Rayleigh quotient iteration
         std::vector<double> Ax(size, 0.0);
  
         // perform matrix-vector multiplication: Ax = A * x
         for (int i = 0; i < size; ++i) {
             for (int j = 0; j < size; ++j) {
                 Ax[i] += matrix[i][j] * x[j];
             }
         }
  
         // compute the Rayleigh quotient
         double lambda =
             std::inner_product(x.begin(), x.end(), Ax.begin(), 0.0) /
             std::inner_product(x.begin(), x.end(), x.begin(), 0.0);
  
         // check for convergence
         if (std::abs(lambda - lambdaPrev) < tolerance) {
             return lambda;
         }
  
         // update the guess vector
         double normAx = std::sqrt(
             std::inner_product(Ax.begin(), Ax.end(), Ax.begin(), 0.0));
         for (int i = 0; i < size; ++i) {
             // normalize each element of x by dividing by its magnitude
             x[i] = Ax[i] / normAx;
         }
  
         lambdaPrev = lambda;
     }
  
     std::cerr << "Warning: Rayleigh quotient iteration did not converge within "
                  "the specified iterations."
               << std::endl;
     // if not converged
     return 0.0;
 }

◆ schur_decomp()

void gpmp::linalg::Eigen::schur_decomp	(	std::vector< std::vector< double >> &	schurMatrix,
		double	tolerance = `1e-10`
	)		const

Computes the Schur decomposition of a matrix.

This method computes the Schur decomposition of a matrix The Schur decomposition is a block triangular decomposition of a matrix, where the diagonal blocks are square and the off-diagonal blocks are zero or upper triangular

Parameters

schurMatrix	The output matrix to store the Schur decomposition
tolerance	The tolerance for determining if a subdiagonal element is zero

The schurMatrix parameter must be pre-allocated with the same dimensions as the matrix object This method modifies the contents of the schurMatrix parameter

Definition at line 181 of file eigen.cpp.

                             {
     // Hessenberg matrix is already calculated in QR iteration
     schurMatrix = matrix;
  
     for (int i = size - 1; i > 0; --i) {
         // check for small subdiagonal elements to introduce a zero
         if (std::abs(schurMatrix[i][i - 1]) < tolerance) {
             schurMatrix[i][i - 1] = 0.0;
         } else {
             // Apply Givens rotation to introduce a zero below the subdiagonal
             double a = schurMatrix[i - 1][i - 1];
             double b = schurMatrix[i][i - 1];
             double r = std::hypot(a, b);
             double c = a / r;
             double s = b / r;
  
             for (int j = 0; j < size; ++j) {
                 double temp1 =
                     c * schurMatrix[i - 1][j] + s * schurMatrix[i][j];
                 double temp2 =
                     -s * schurMatrix[i - 1][j] + c * schurMatrix[i][j];
                 schurMatrix[i - 1][j] = temp1;
                 schurMatrix[i][j] = temp2;
             }
         }
     }
 }

◆ sensitivity()

std::vector< double > gpmp::linalg::Eigen::sensitivity ( double perturbation ) const

Calculates the sensitivity of the eigenvalues of a matrix to small perturbations.

This method calculates the sensitivity of each eigenvalue of a matrix to small perturbations The sensitivity is defined as the change in the eigenvalue divided by the amount of perturbation

Parameters

perturbation The amount of perturbation to apply to the matrix

Returns: A vector containing the sensitivity of each eigenvalue

The size of the returned vector will be the same as the size of the matrix

Definition at line 154 of file eigen.cpp.

                                                       {
     std::vector<double> originalEigenvalues = qr_algorithm();
  
     std::vector<double> perturbedEigenvalues;
     for (int i = 0; i < size; ++i) {
         // perturb the matrix
         std::vector<std::vector<double>> perturbedMatrix = matrix;
         perturbedMatrix[i][i] += perturbation;
  
         // calculate eigenvalues of perturbed matrix
         std::vector<double> perturbedVals =
             Eigen(perturbedMatrix).qr_algorithm();
         perturbedEigenvalues.push_back(perturbedVals[i]);
     }
  
     // calculate sensitivity: (perturbed eigenvalue - original eigenvalue) /
     // perturbation
     std::vector<double> sensitivity;
     for (int i = 0; i < size; ++i) {
         double sens =
             (perturbedEigenvalues[i] - originalEigenvalues[i]) / perturbation;
         sensitivity.push_back(sens);
     }
  
     return sensitivity;
 }

References qr_algorithm().

Member Data Documentation

◆ matrix

std::vector<std::vector<double> > gpmp::linalg::Eigen::matrix

private

Definition at line 48 of file eigen.hpp.

◆ size

int gpmp::linalg::Eigen::size

private

Definition at line 50 of file eigen.hpp.

The documentation for this class was generated from the following files:

include/openGPMP/linalg/eigen.hpp
modules/linalg/eigen.cpp

Public Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ Eigen()

Member Function Documentation

◆ determinant()

◆ jordan_normal_form()

◆ power_iter()

◆ qr_algorithm()

◆ rayleigh_iter()

◆ schur_decomp()

◆ sensitivity()

Member Data Documentation

◆ matrix

◆ size