Matrix Exponential

ROOT
1

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 :smiley:

/*
\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;
}
2

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.
 *
 * exceptions: GSL_ENOTSQR, GSL_EBADLEN
 */
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.

3

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:)