Why is "expm" computation for symbolic matrix slow?

I am trying to do matrix exponential computation symbolically because I want a high-precision result. But it is taking forever to run:

>> digits(32)

>> a = sym(1.3651817760E+10)*sym(log(2));
>> b = sym(1.0E-04);
>> c = a*b;

>> A = [c c 0 0 0 0 0 0 0 0 0 0;
        0 c a c 0 0 0 0 0 0 0 a;
        0 0 c a c 0 0 0 0 0 0 0;
        0 0 0 c c c c 0 0 0 0 c;
        0 0 0 0 0 c b c 0 0 0 0;
        0 0 0 0 0 0 c c c 0 0 0;
        0 0 0 0 0 0 0 0 c c c 0;
        0 0 0 0 a 0 0 0 0 c c c;
        0 0 0 b 0 0 0 0 0 0 c c;
        0 0 0 0 0 0 0 0 0 0 c a;
        0 0 0 0 0 0 0 0 0 0 c c;
        0 0 0 0 0 0 0 0 0 0 0 c];

>> expm(A)

Answer 1

[MathWorks Support Team]

1. Why is computing the matrix exponential of A slow?

In the script, you defined matrix A in term of sym (each sym being applied to a MATLAB double). The only 'non-trivial' (i.e. exact irrational) input value is log(2).

Each double is approximated by a rational number. Thus, the matrix ending up in the symbolic engine consists of exact symbolic data (rational numbers and the occasional log(2)).

Note that the sparsity of matrix A does not help: for exact input data, the matrix has to be diagonalized symbolically. The resulting matrix of eigenvectors is not sparse. It contains many rational numbers with huge numerators and denominators.

The runtime is mainly determined by a final inversion of the (dense) eigenvector matrix (this is needed algorithmically for "expm"), and therefore causing the slow performance.

1. About using symbolic variables to get high-precision result:

Numerical exponential of a matrix is a very stable process. Notice that the input data consist of symbolic rational approximations of MATLAB doubles,and therefore contains numerical input errors even with the symbolic computation.

2. The workaround:

Thus, you should get a much faster (and same precision) computation when you use "vpa" to the input of "expm" when all input values are of purely numerical nature. (The resulting numerical vpa algorithm in "expm" is completely different from the "sym" algorithm and numerically stable, so you do not have to worry about numerical errors)

>> A = vpa(A);
>> E = expm(A);