# Matrix Exponential

I am quite new to ROOT and I mainly use it for finance quantitative analysis. I find that there is no function to calculate matrix exponential, so I just wrote a simple piece of code based on TMatrixD. Hope we can have a formal one in the future

``````/*
\file: MatrixExponential.cpp

Calculation of matrix exponential

author: Rucheng YANG(Brandon),GATECH

website:http://vschon.com/2013/05/12/matrix-exponential/

date: 2013/5/11

Dependency: ROOT
*/
#include "TMath.h"
#include "TMatrixD.h"

#include <iostream>

using namespace std;

TMatrixD MatrixExponential(TMatrixD A,Int_t m)
{
//Function to Calculate the exp(A) using Pade Approximation
/*Reference: http://www.maths.manchester.ac.uk/~higham/talks/exp09.pdf
How and How Not to Compute the Exponential of a Matrix by Nick Higham
*/

//Input:
//A: Square Matrix
//m: Expansion degree

//Get size of A
int NRow=A.GetNrows();
int NCol=A.GetNcols();

//P,Q: Expansion matrix
TMatrixD P(NRow,NCol);
TMatrixD Q(NRow,NCol);
P.Zero();Q.Zero();
TMatrixD A_j(NRow,NCol);
TMatrixD A_minus_j(NRow,NCol);
TMatrixD zeroMat(NRow,NCol);
TMatrixD A_minus(NRow,NCol);
zeroMat.Zero();
A_j.UnitMatrix();
A_minus_j.UnitMatrix();
P.UnitMatrix();
Q.UnitMatrix();

//C: Coefficient
Double_t C=0;

for(int j=1;j<=m;j++)
{
//A_j:calculating A to the power j
//A_minus_j:Calculatig (-A) to the power of j
A_j=A_j*A;
A_minus_j=((j%2)==0?A_j:(zeroMat-A_j));

//Calculating coefficient C
C=TMath::Factorial(2*m-j)*TMath::Factorial(m)/(TMath::Factorial(2*m)*TMath::Factorial(j)*TMath::Factorial(m-j));

P=P+C*A_j;
Q=Q+C*A_minus_j;
}
return P*Q.Invert();
}

int main()
{
//Test the function, input 0 to quit
Int_t dim;
cout<<"Matrix dim=: "<<endl;
cin>>dim;
while(1)
{
if(dim==0){break;}
TMatrixD TestMat(dim,dim);
for(int i=0;i<dim;i++){
for(int j=0;j<dim;j++){
TestMat[i][j]=i+j+1;
}
}
TestMat.Print();

TestMat=MatrixExponential(TestMat,50);
TestMat.Print();
cout<<"Matrix dim=: "<<endl;
cin>>dim;
}
return 0;
}``````

Hi,

thanks for sharing! Note that the GSL library (which is practically a dependency of ROOT already) includes an undocumented function to calculate exponentials of matrices with different methods. That function will pick an optimal parameter set for PadÃ© or Taylor series methods based on the input matrix.

``````/* Calculate the matrix exponential by the scaling and
* squaring method described in Moler + Van Loan,
* SIAM Rev 20, 801 (1978). The mode argument allows
* choosing an optimal strategy, from the table
* given in the paper, for a given precision.
*
*/
int gsl_linalg_exponential_ss(
const gsl_matrix * A,
gsl_matrix * eA,
gsl_mode_t mode
);``````

mode is a GSL precision constant.
So in the spirit of reducing custom implementations (aka prefer standing on the shoulders of giants)

``````#include <gsl/gsl_linalg.h>
#include <TMatrixD.h>

TMatrixD GSL_exp(const TMatrixD& a) {
const size_t Nx = a.GetNrows();
gsl_matrix_const_view va =
gsl_matrix_const_view_array(a.GetMatrixArray(), Nx, Nx);

TMatrixD b(Nx, Nx);
gsl_matrix_view vb = gsl_matrix_view_array(b.GetMatrixArray(), Nx, Nx);

gsl_linalg_exponential_ss((const gsl_matrix*)&va, (gsl_matrix*)&vb, 0);

return b;
}``````

I am not sure this method seems only faster because of better memory usage patterns.

Hi, honk, thanks for your feedback. I did not notice that there is such a function in GSL library. It is much convenient than my little toy code:)