# Computing Matrix Exponentials using C++ with boost library

## Introduction

Evaluating the exponential of a matrix is an important problem that arises in physics, mathematics and engineering. For example, in quantum theory, a central problem consists in solving the ODE of Schrödinger equation,

where is a hermitian matrix and is a complex vector.

If the Hamiltonian, is time independent, the solution for Eq. (??) is given by

where is the final state, and is the initial state.

There are several methods to compute matrix exponential. Methods involving approximation theory, differential equations, the matrix eigenvalues, and the matrix characteristic polynomial have been proposed.

## Computed by Taylor Series Expansion

A direct way to define the matrix exponential <p class="ql-center-displayed-equation" style="line-height: 18px;">     </p>

where t is a real number is undoubtedly  through the exponential power series expansion,

whose convergence is guaranteed for any square matrix . If we momentarily ignore the efficiency, we can simply sum up the series until adding another term which does not alter the numbers stored in computer. That is, if

, then we can find such that

where is some prescribed error tolerance. Concern over where to truncate the series is of importance if efficiency is being considered. Unfortunately, such an algorithm is know to be inefficiency even in the scalar case, and there are lots of papers\cite{Moler2003, Liou1966} concerning the truncation error of Taylor series. As a result, directly summing up all Taylor series seems not work very well.

## Computed by Diagonalization of the Matrix

For most physicists, another instinctive method is to diagonalize the matrix , such that

then using the power series definition of <p class="ql-center-displayed-equation" style="line-height: 17px;">     </p>

implies

where is the matrix whose columns are eigenvectors of A, , and are the eigenvalues of . The exponential of is very easy to compute and we have a satisfactory build-in routine in lots of computer languages for computing the exponential of a scalar,

The difficulty with this approach is not find the eigenvalues of the matrix in itself, but occurs when does not have a complete set of linearly independent eigenvectors. In this case, is defective and there is no invertible matrix of <p class="ql-center-displayed-equation" style="line-height: 12px;">     </p>

, therefore, the algorithm will break down. In real computing world, the difficulties occur even when is nearly defective. Define the condition number as . While is nearly defective, then will be very large. The errors of computing <p class="ql-center-displayed-equation" style="line-height: 17px;">     </p>

, including the roundoff errors from the eigenvalues computation, may be magnified by . As a result, the computed exponential of a matrix will most likely be inaccurate when is very large. An example had been demonstrated in this paper\cite{Moler2003}.

The most easy method for computing numerically might be Taylor series, but we have discussed that this approach is inefficiency and not accurate. A Padé approximation often gives better result of a given function – here than truncating its Taylor series, and it may still work where the Taylor series does not converge well. Since it is defined as a rational function which is a ratio of polynomial series, it can be calculated numerically easily.  Therefore, Padé approximations are used extensively in computer calculation. A Padé approximation approximates a function in only one variable, an approximation of a function in two variables is called a Chisholm approximation, and in multiple variables is called a Canterbury approximation.

A Padé rational approximation to is the quotient of two polynomials and of degrees and respectively. This is so-called (p,q)-degree type Padé approximation. We use the notation to denote this quotient:

where x is a scalar. If x is a matrix – , the quotient is defined as

The polynomials used in Eq. (??) are

and

In the case , the approximation will reduce to the Taylor (Maclaurin) expansion for . There are unknown coefficients in , and unknown coefficients in , hence the rational function has unknown coefficients. Assume that <p class="ql-center-displayed-equation" style="line-height: 18px;">     </p>

is analytic and has the Maclaurin expansion,

then is said to be a Padé approximation to the series .
Since the highest possible oder of nonzero derivative of is , that is

as a result, the first derivatives of and are to agree at ,

Eq. (??) implies that,

Multiply on Eq. (??) giving

and we obtain

When the left side of Eq. (??) is multiplied out out and the coefficients of the power of are set equal to zero for , the result is a system of linear equations:

and

The equations in Eq. (??) involve only the  unknowns , and have to be solved first. Then the equations in Eq. (??) are used successively to find .

Go back to our original problem, and setting as gives us,

Solving Eq. (??) and Eq. (??) together with Eq. (??)  gives us,

and

It has been discussed that there are several reasons \cite{Moler2003} why the diagonal approximants are preferred over the off diagonal approximants for stability and economy of computation. For , we have , , and . As noted in Sidje’s thesis\cite{Sidje1994}

, Eq. (??) can be written as the following irreducible form for economical computing reason,

Since Padé approximation is only accurate near the origin so that the approximation of is not valid, when is too large. Fortunately, this problem will be solved when we introduce so-called scaling and squaring technology\cite{Moler2003 ,Sidje1998}. We make use of the exponential property

where is chosen such that . The idea is to choose to be a nature number for which can be reliably and efficiently computed, and then to compute the result by repeated squaring.

Because the result is evaluated by repeated squaring the exponential of scaled matrix, the drawback of this algorithm may come from the fact that if , then the roundoff errors may be large. The error analysis discussed in these papers\cite{Moler2003, Sidje1998} has shown that if <p class="ql-center-displayed-equation" style="line-height: 19px;">     </p>

then

where

Therefore, a value of is generally satisfactory for computer computing while using double precision.

## Matrix Exponential Source Code

In this section, we will show the source code of matrix exponential in Listing ?? which is implemented based on irreducible -degree Padé approximation described in the previous section and Fortran source code of Expokit software package\cite{Sidje1998}. Our source code is written by C++, and we take the advantage of Boost C++ Libraries which have a lots of useful Basic Linear Algebra routines (uBLAS). Our source code is released under GPLv2, its later version, or Boost Software License Version 1.0.

/*!
*   \file expm.hpp
*
*   Implement matrix exponential using pade approximation.
*
*  \author Tsai, Dung-Bang  蔡東邦
*
*/
//  Department of Physics,
//  National Taiwan University.
//
//  E-Mail : dbtsai [_at_] dbtsai [_dot_] org
//  Begine : 2007/11/20
//  Last modify : 2007/11/26
//  Version : v0.4
//
//  expm_pad computes the matrix exponential exp(H) for general matrixs,
//  including complex and real matrixs using the irreducible (p,p) degree
//  rational Pade approximation to the exponential
//  exp(z) = r(z) = (+/-)( I+2*(Q(z)/P(z))).
//
//  Usage :
//
//   U = expm_pad(H, t, p),
//
//  where t is a real number which is default set to 1.0 such that U=exp(t*H),
//  and p is internally set to 6 (recommended and gererally satisfactory).
//
//
//  Reference :
//  EXPOKIT, Software Package for Computing Matrix Exponentials.
//  ACM - Transactions On Mathematical Software, 24(1):130-156, 1998
//
// Use, modification and distribution are subject to the Boost Software License,
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
//

#ifndef _BOOST_UBLAS_EXPM_
#define _BOOST_UBLAS_EXPM_
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/traits.hpp>

namespace boost { namespace numeric { namespace ublas {

template MATRIX expm_pad(const MATRIX H, typename type_traits::real_type t = 1.0, const int p = 6){
typedef typename MATRIX::value_type value_type;
typedef typename MATRIX::size_type size_type;
typedef typename type_traits<value_type>::real_type real_value_type;
assert(H.size1() == H.size2());
assert(p >= 1);
const size_type n = H.size1();
const identity_matrix<value_type> I(n);
matrix<value_type> U(n,n),H2(n,n),P(n,n),Q(n,n);
real_value_type norm = 0.0;
vector<real_value_type> c(p+1);
c(0)=1;
for(size_type i = 0; i < p; ++i)
c(i+1) = c(i) * ((p - i)/((i + 1.0) * (2.0 * p - i)));
// Calcuate the infinty norm of H, which is defined as the largest row sum of a matrix
for(size_type i=0; i<n; ++i) {
real_value_type temp = 0.0;
for(size_type j = 0; j < n; j++)
temp += std::abs(H(i, j));
norm = t * std::max<real_value_type>(norm, temp);
}
// If norm = 0, and all H elements are not NaN or infinity but zero,
// then U should be identity.
if (norm == 0.0) {
bool all_H_are_zero = true;
for(size_type i = 0; i < n; i++)
for(size_type j = 0; j < n; j++)
if( H(i,j) != value_type(0.0) )
all_H_are_zero = false;
if( all_H_are_zero == true ) return I;
// Some error happens, H has elements which are NaN or infinity.
std::cerr<<"Null input error in the template expm_pad.\n";
std::cout << "Null INPUT : " << H <<"\n";
exit(0);
}
// Scaling, seek s such that || H*2^(-s) || < 1/2, and set scale = 2^(-s) int s = 0; real_value_type scale = 1.0; if(norm > 0.5) {
s = std::max(0, static_cast((log(norm) / log(2.0) + 2.0)));
scale /= real_value_type(std::pow(2.0, s));
U.assign((scale * t) * H); // Here U is used as temp value due to that H is const
}
else
U.assign(H);

// Horner evaluation of the irreducible fraction, see the following ref above.
// Initialise P (numerator) and Q (denominator)
H2.assign( prod(U, U) );
Q.assign( c(p)*I );
P.assign( c(p-1)*I );
size_type odd = 1;
for( size_type k = p - 1; k > 0; --k) {
( odd == 1 ) ?
( Q = ( prod(Q, H2) + c(k-1) * I ) ) :
( P = ( prod(P, H2) + c(k-1) * I ) ) ;
odd = 1 - odd;
}
( odd == 1 ) ? ( Q = prod(Q, U) ) : ( P = prod(P, U) );
Q -= P;
// In origine expokit package, they use lapack ZGESV to obtain inverse matrix,
// and in that ZGESV routine, it uses LU decomposition for obtaing inverse matrix.
// Since in ublas, there is no matrix inversion template, I simply use the build-in
// LU decompostion package in ublas, and back substitute by myself.

// Implement Matrix Inversion
permutation_matrix<size_type> pm(n);
int res = lu_factorize(Q, pm);
if( res != 0) {
std::cerr << "Matrix inversion error in the template expm_pad.\n";
exit(0);
}
// H2 is not needed anymore, so it is temporary used as identity matrix for substituting.
H2.assign(I);
lu_substitute(Q, pm, H2);
(odd == 1) ?
( U.assign( -(I + real_value_type(2.0) * prod(H2, P))) ):
( U.assign(   I + real_value_type(2.0) * prod(H2, P) ) );
// Squaring
for(size_type i = 0; i < s; ++i)
U = (prod(U,U));
return U;
}

}}}

#endif


BOOST LIBRARY · NUMERICAL ANALYSIS
English