expm.h File Reference

Matrix exponential by scaling and squaring with a diagonal Pade approximant. More...

#include <cmath>
#include <cstddef>
#include <limits>
#include <vector>
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/linalg.h"
#include "line/util/matrix.h"

Include dependency graph for expm.h:

Go to the source code of this file.

Namespaces
namespace	line

Functions
template<class T>
Matrix< T >	line::expm (const Matrix< T > &A)
	Matrix exponential exp(A).
template<class T>
Matrix< T >	line::expm (const Matrix< T > &A, const T &t)
	exp(t A), the form every MAP descriptor actually needs.

Detailed Description

Matrix exponential by scaling and squaring with a diagonal Pade approximant.

This is the numerical primitive behind MATLAB's expm(), which the kpctoolbox MAP counting-process descriptors (map_cdf, map_pdf, map_acfc, map_count_var, map_varcount, map_count_moment) and the LRU(m)-MAP TTL cache approximation (cache_lrum_map_levelstats) all call. It has no MATLAB source file of its own in the tree; the reference is the algorithm of

N. J. Higham, "The scaling and squaring method for the matrix exponential revisited", SIAM J. Matrix Anal. Appl. 26(4):1179-1193, 2005,

which is what MATLAB implements. Given A, pick a scaling s so that ||2^-s A||_1 <= theta_m, evaluate the [m/m] Pade approximant

r_m(X) = D_m(X)^-1 N_m(X), N_m(X) = sum_k c_k X^k, D_m(X) = sum_k (-1)^k c_k X^k, c_0 = 1, c_k = c_{k-1} (m-k+1) / (k (2m-k+1)),

and square the result s times: exp(A) = r_m(2^-s A)^(2^s).

ARITHMETIC: the answer is an approximation controlled by a tolerance – the Pade truncation error is nonzero for a general A no matter how the arithmetic is carried out – so this is gated on num_traits<T>::has_transcendental and refuses to instantiate at exact rational arithmetic. It is nevertheless exact (to rounding) on the two cases that matter for testing: A = 0 returns the identity by a shortcut, and a nilpotent A with A^q = 0, q <= 2m+1, is reproduced exactly because the Pade error term is O(X^(2m+1)).

PRECISION AWARENESS: Higham's theta table is calibrated for double unit roundoff, and using it unchanged at Real<50> would silently deliver only double accuracy. The bound on the [m/m] Pade error,

|exp(x) - r_m(x)| ~ C_m |x|^(2m+1), C_m = (m!)^2 / ((2m)! (2m+1)!),

is inverted here against the working precision of T, so theta_m shrinks as the precision grows (theta_13 = 5.1 at double, 0.28 at 50 digits, 0.0039 at 100 digits). At double the resulting thresholds agree with Higham's table to within a few percent.

Definition in file expm.h.