Computation of Eigenvalues and Eigenvectors in C++ by Using Eigen C++ Matrix Library – Fusion of Engineering, Control, Coding, Machine Learning, and Science

In this scientific computing, numerical analysis, and control engineering tutorial, we will learn how to compute eigenvalues and eigenvectors of matrices in C++ by using the Eigen matrix library. The YouTube video accompanying this post is given below.

Eigenvector and Eigenvalues in Eigen C++ Matrix Library

The first step is to install the Eigen matrix library. If you are using the VS Code editor, you can follow this tutorial to do that. On the other hand, if you are using Microsoft Visual Studio, then follow this tutorial. The tutorial presented here is based on the VS Code editor. The GitHub page with all the code files is given here. The first step is to include the necessary libraries and namespaces:

#include<iostream>

#include<vector>
#include <iterator>

#include<Eigen/Dense>
#include <Eigen/Eigenvalues> 

using namespace std;
using namespace Eigen;

The include statement

#include<Eigen/Dense>

is used to include a set of header files that are used for dense matrix computations. On the other hand, the include statement

#include <Eigen/Eigenvalues>

is used to include the header files for algorithms that are used for eigenvalue computations. The code presented below should be placed inside the “main()” C++ function.

First, define a random matrix whose eigenvalues and eigenvectors we will compute. Also, we define complex matrices for storing the eigenvalues and eigenvectors.

        // create a random matrix
    MatrixXd randomMatrix;
    randomMatrix.setRandom(3,3);
    // this vector/matrix will be used to store the eigenvalues
    
    // The matrix has to be complex
    // This is the typedef for MatrixXcd
    // typedef Matrix< std::complex<double> , Dynamic , Dynamic > MatrixXcd
    // That is, this matrix can store std::complex<double> data type
    MatrixXcd eigenvalueMatrix;
    // this matrix will be used to store the eigenvectors
    // the matrix needs to be complex since the eigenvalues can be complex
    MatrixXcd eigenvectorMatrix;
    // print the random matrix
    cout <<endl<< randomMatrix << endl;

The data type “MatrixXd” is a type definition for a matrix of double real entries whose dimensions can be changed during the compile time. On the other hand, “MatrixXcd” is a type definition for a matrix of double complex entries whose dimensions can be changed during the compile time. The complex matrix “eigenvalueMatrix” is used to store the complex eigenvalues. The matrix “eigenvectorMatrix” is used to store the eigenvectors of the matrix “randomMatrix” column-wise.

Next, we compute the eigenvalues and eigenvectors

 EigenSolver<MatrixXd> eigenValueSolver(randomMatrix);
    
    eigenvalueMatrix=eigenValueSolver.eigenvalues();
    // the eigenvectors are stored column wise
    eigenvectorMatrix=eigenValueSolver.eigenvectors();
    // print the eigenvalues
    cout <<endl<<"Computed eigenvalues:"<<endl<<eigenvalueMatrix<<endl;
    // print the eigenvectors
    cout <<endl<<"Computed eigenvectors:"<<endl<<eigenvectorMatrix<<endl;

   // extract real and imaginary parts of the eigenvalues
    MatrixXd realPartEigenValues;
    MatrixXd imagPartEigenValues;

    realPartEigenValues=eigenvalueMatrix.real();
    imagPartEigenValues=eigenvalueMatrix.imag();

    // print the real parts of the eigenvalues
    cout<<endl<<realPartEigenValues<<endl;
    // print the complex parts of the eigenvalues
    cout<<endl<<imagPartEigenValues<<endl;

The following code line

 EigenSolver<MatrixXd> eigenValueSolver(randomMatrix);

is used to create the EigenSolver object. We pass to the constructor the matrix “randomMatrix”. The eigenvalues and eigenvectors are computed as follows

    eigenvalueMatrix=eigenValueSolver.eigenvalues();
    eigenvectorMatrix=eigenValueSolver.eigenvectors();

The eigenvectors are stored column-wise in the matrix “eigenvectorMatrix”. We extract the real and imaginary parts of the computed eigenvalues as follows

    realPartEigenValues=eigenvalueMatrix.real();
    imagPartEigenValues=eigenvalueMatrix.imag();

The following piece of code will double-check the result. That is, it will pick one computed eigenvalue and eigenvector and it will check if the defining equation is satisfied.

  // let us double check one eigenvalue and eigenvector
    MatrixXcd testLeft;
    MatrixXcd testRight;

    // definition of the eigenvalue and the eigenvector
    // test1Right = A*v 
    // test2Right = lambda*v
    // A*v=lambda*v 
    // if everything is computed correctly
    // testLeft = testRight
    testLeft=randomMatrix*eigenvectorMatrix.col(0);
    testRight=eigenvalueMatrix(0,0)*eigenvectorMatrix.col(0);
    cout<<endl<<"Left side vector:"<<endl<<testLeft<<endl;
    cout<<endl<<"Right side vector:"<<endl<<testRight<<endl;
    cout<<endl<<"Left minus right side vector:"<<endl<<testLeft-testRight<<endl;

This code will test if this equation

(1) $\begin{align*}A\mathbf{v}=\lambda \mathbf{v}\end{align*}$

is satisfied, where $A$ is the random matrix, $\mathbf{v}$ is the computed eigenvector, and $\lambda$ is the computed eigenvalue. We test is the left-hand side of this equation is equal to the right-hand side. That is if $A\mathbf{v}$ equals $\lambda \mathbf{v}$ . We also test if the left-hand side of this expression:

(2) $\begin{align*}A\mathbf{v}-\lambda \mathbf{v}=0\end{align*}$

is equal to zero. By running this piece of code, if we can see that the eigenvalues and eigenvectors are correctly computed. Next, we explain how to perform basic operations on computed eigenvalues. The following piece of code will double the computed eigenvalues and it will store them in C++ vector STL object.

// let us do some manipulations with eigenvalues
    // for example let us define a vector structure of eigenvalues
    // we need std:: complex<double>

    std::vector<std::complex<double>> doubledEigenValueVector;
    // eigenvalues are of the type complex <double>
    complex<double> eigenValue;
   
    for (int i=0; i<eigenvalueMatrix.rows(); i++)
    {
        eigenValue=eigenvalueMatrix(i,0);
        cout<<endl<<" Real part of the eigenvalue "<<i<<" is "<<eigenValue.real()<<endl;
        cout<<" Imaginary part of the eigenvalue "<<i<<" is "<<eigenValue.imag()<<endl;
        // double the eigenvalues
        doubledEigenValueVector.push_back(2.0*eigenValue);
    }
    // print the results
    for (auto it=doubledEigenValueVector.begin(); 
    it !=doubledEigenValueVector.end();
    it++)
    {
        cout<<endl<<"Original eigenvalue: "<<eigenvalueMatrix(distance(doubledEigenValueVector.begin(), it),0)<<endl;
        cout<<"Doubled eigenvalue: "<< *it <<endl;
    }

The above code is self explanatory. Here, we should just say that eigenvalues as complex numbers are stored by using the STL “complex <double>” data type.

Eigenvector and Eigenvalues in Eigen C++ Matrix Library

admin

You might also like

Install Arduino in Linux (Ubuntu)

Wiring diagram for the Differential Wheeled Robot – Ver 1.0

Simulating State-Space Models in Eigen C++ Matrix Library – Object-Oriented Program