Matrix Cosine is Slooooooow!

[rdbyk]

To Reproduce

--> A=rand(500,500);
--> tic;expm(A);toc
Elapsed time is 0.382739 seconds.
--> tic;cosm(A);toc
Elapsed time is 11.684408 seconds.

Desktop (please complete the following information):

• OS: Ubuntu 14.04 LTS

[Nelson-numerical-software]

cosm is equivalent to 0.5*(expm(iA)+expm(-iA))
currently expm(real) ok but expm(complex) is sloooooow :)

[rdbyk]

Aah I see, so a little improvement could be achieved by

--> A=rand(500,500);
--> tic;B=expm(i*A);COSM=0.5*(B+inv(B));toc
Elapsed time is 7.523827 seconds.

[Nelson-numerical-software]

other way is to use expoKit (Cpp version) but there is some licence limitation (ACM)

[rdbyk]

For real arguments, it could be done as

Nelson:

--> A=rand(500,500);
--> tic;COSM=real(expm(i*A));toc
Elapsed time is 6.079377 seconds.

Regarding the implementation of expm itself, why not looking at Scilab, Octave, or Julia ...

Scilab:

--> A=rand(500,500);
--> tic;COSM=real(expm(%i*A));toc
 ans  =
   0.702032 

Octave:

octave:1> A=rand(500,500);
octave:2> tic;COSM=real(expm(i*A));toc
Elapsed time is 1.49243 seconds.

Julia:

julia> A=rand(500,500);

julia> tic();COSM=real(expm(I*A));toc();
elapsed time: 2.118478498 seconds

julia> tic();COSM=real(expm(I*A));toc();
elapsed time: 0.231352753 seconds

julia> B=rand(500,500);

julia> tic();COSM=real(expm(I*B));toc();
elapsed time: 0.305516143 seconds

[Nelson-numerical-software]

from my memory, scilab or nsp uses expoKit.

for real, I will continue to use current code.
for complex, I will try to do something as:
[V, D] = eig(A);
expmcomplex = V * diag(exp(diag(D))) / V;
but I need to implement diag and exp ...

[rdbyk]

It seems, that in Scilab you did it by yourself a long time ago ;-), cf. here

expoKit is not mentioned in the Copyright/License header !?

[Nelson-numerical-software]

speed optimization about expm using blas, lapack (mkl on windows).

previous version on Windows:

-->  A = rand(500,500) + i;
--> tic();expm(A);toc();
Elapsed time is 22.779935 seconds.
-->  A = rand(500,500);
--> tic();expm(A);toc();
Elapsed time is 0.564563 seconds.

New version using MKL on windows

--> A = rand(500,500) + i;
--> tic();expm(A);toc();
Elapsed time is 1.798484 seconds.
--> A = rand(500,500);
--> tic();expm(A);toc();
Elapsed time is 0.242285 seconds.

for real numbers, expm is very fast

[Nelson-numerical-software]

expm (complex numbers) will be more optimized in a next version.

[Nelson-numerical-software]

About cosm, on overloading branch with cosm updated:
on my notebook pc:
--> A = rand(500, 500) + i;
tic();cosm(A);toc

Elapsed time is 1.482141 seconds.

previous version build 1068:
A = rand(500, 500) + i;
tic();cosm(A);toc

Elapsed time is 2.245938 seconds