Calculate larger exp(A)!?
Magnus Ögren
I'm working (presently in Matlab) on a quantum-physics problem where I
need access to the matrix exp(A), i.e. I do not have any c-number 'v'
vector to multiply it with.
From scenario 1 on your help page with F.A.Q. I read "Simply use the
built-in function expm. Expokit provides another function if your
matrix is sparse (see scenario 3 below)." But in scenario 3 I only see
the case with a 'v' vector i.e. exp(A)*v
The matrix A=[A11 A12; A21 A22] have a structure such that A11 and A22
is diagonal and A12 and A21 is sparse i.e. >50% zeros.
So I have a sparse A matrix but no vector 'v', is there any algorithm
that would be better for me than expm(A)?
I guess simply expm(sparse(A)) does not make any improvements.
Presently I'm working in matlab and have calculated matrices in the
size of n=6000.
Would going to fortran or C++ substantially allow me to use a larger n?
It is the size of the matrix (rather than the computational time) that
is limiting my accuracy at the moment. Is it simply the RAM memory of
the PC that sets the maximum size n?
Best regards!
Magnus (PhD student in physics)
Arjan van Dijk
> So I have a sparse A matrix but no vector 'v', is there any algorithm
> that would be better for me than expm(A)?
To my opinion it is no problem that you don't have a vector "v".
Make a loop over the dimension of "v" and calculate expm(v)
for vectors "v" where all elements are zero except for 1, which you
set to 1. The resulting vectors are the columns of the exponential
matrix you want.
Regards,
Arjan
Roger B. Sidje
>
> Dear exp(A) experts!
>
>
> I'm working (presently in Matlab) on a quantum-physics problem where I
> need access to the matrix exp(A), i.e. I do not have any c-number 'v'
> vector to multiply it with.
What do you do with E=exp(A) in full? Do you use ALL the entries of E
once you get it?
> From scenario 1 on your help page with F.A.Q. I read "Simply use the
> built-in function expm. Expokit provides another function if your
> matrix is sparse (see scenario 3 below)." But in scenario 3 I only see
> the case with a 'v' vector i.e. exp(A)*v
>
> The matrix A=[A11 A12; A21 A22] have a structure such that A11 and A22
> is diagonal and A12 and A21 is sparse i.e. >50% zeros.
> So I have a sparse A matrix but no vector 'v', is there any algorithm
> that would be better for me than expm(A)?
%A==sparse(A); % this is assumed!
E = zeros(size(A));
I = eye(size(A));
for j = 1:n
v = I(:, j)
E(:, j) = expv(1, A, v, 15, 1e-5); % for example.
end
But I still wonder what you do with E. It is quite unusual to need
exp(A) in full for large sparse problems.
>
> I guess simply expm(sparse(A)) does not make any improvements.
You mean expm(full(A)). expm doesn't run with a sparse matrix.
>
> Presently I'm working in matlab and have calculated matrices in the
> size of n=6000.
> Would going to fortran or C++ substantially allow me to use a larger n?
>
> It is the size of the matrix (rather than the computational time) that
> is limiting my accuracy at the moment. Is it simply the RAM memory of
> the PC that sets the maximum size n?
Matlab is noticeably slower. But its convenience tends to compensate
for its slow speed if you only use very large n occasionally. You
should consider switching to the compiled languages if you are
consistently going to be using large n, especially in your context
where you need exp(A) in full.
The way to decide is to try expv and see if you can put up its run
time. If you start getting impatient of waiting, then it is time to
try the compiled languages.
Roger B. Sidje
Roger B. Sidje
>
> %A==sparse(A); % this is assumed!
> E = zeros(size(A));
> I = eye(size(A));
> for j = 1:n
> v = I(:, j)
> E(:, j) = expv(1, A, v, 15, 1e-5); % for example.
> end
>
Alternatively, a slightly more efficient way without the intermediate I
%A==sparse(A); % this is assumed!
E = zeros(size(A));
v = zeros(n, 1);
for j = 1:n
v(j) = 1;
E(:, j) = expv(1, A, v, 15, 1e-5); % for example.
v(j) = 0;
end
once you get it
h m ogren
Thanks for your ideas and comments!
Why do I need the full exp(A) matrix?
I'm solving linear ODE's for apples, pears and bananas... i.e. symbols
instead of c-numbers.
That is I express a certain symbol (creation and annihilation
operators for atoms) at time t (assumed to be one below) as a linear
combination of such symbols at time zero. It is therefore the rows of
the exp(A) matrix that are interesting for me, since the scalar
product of certain rows will give me my observables, i.e. certain
combinations of pairs of symbols.
If you (R. B. S.) are still working at U.Q. I can come in for an
informal presentation of my problem at a group meeting or so in a
couple of weeks if you like, since I'm going to be down in Brisbane
for a couple of months now.
The method of calculating exp(A)*v, v being 'all Cartesian base
vectors' indeed gives the full exp(A) matrix. In my case I'm using the
rows, so it is better to calculate exp( transpose(A) )*v. It takes
substantially longer time then expm(full(A)), but can handle larger
matrices :-)
I still think may be ways to simplify the problem due to the
symmetries of the A matrix, such that I do not need to calculate
exp(A) of the full A matrix?
If I write A in terms of four block matrices
A22=conj(A11)=-A11 since A11 is pure imaginary (and also diagonal)
A21=q*conj(A12), where q=+-1, since A12 is also real for most
interesting cases (*): A21=q*A12, i.e.:
A=[A11 A12; q*A12 -A11]
(*) Then A is antihermitian, and "The matrix exponential of an
antihermitian matrix is a unitary matrix".
Now my numerical calculations of the matrix M=exp(A)=[M11 M12; M21
M22] (that is neither real or sparse) suggests that:
M22=conj(M11) and M21=q*conj(M12) (**)
It is always the case that M^\dagger=M^{-1}, i.e. M is unitary as is
often the case in quantum physics (though my approach involves some
approximations, so it was not clear for me a priori)
In my view (**) indicates that a simpler procedure (e.g. to only
calculate exp of block matrices) may be possible? (Besides it have
already reduced the calculation of M=exp(A) to half the number of
rows.)
Another thing that could lead to improvements are the way the (sparse)
matrix is represented in Matlab.
In general a matrix with some (e.g. only one) complex or pure
imaginary elements uses two times the memory of a pure real matrix.
Since at the moment it is the size of the memory for the sparse A
matrix that is my limitation (i.e. the RAM memory of the computer) and
my matrix is real (in most cases of interest) except for the diagonal
which is always pure imaginary, I then guess one could maybe modify
expv.m such that it takes as input one real sparse matrix and one
complex vector (the imaginary diagonal), in this way the memory
required for sparse(A) would be a factor of two less.
I have also tried to perform my calculations with 'single precision'
i.e. 'A=single(sparse(A))' with good result, meaning that I did not
see any differences in the elements of M=exp(A) (checked only a few
cases...) within 'format short'. My calculations will result in
'qualitative images' rather than precise numbers, maybe it is then OK
to use single precision which reduce the memory limitations by a
factor of two. Do you have any experiences of how good/bad single
precision is for exp(A)?
Another way to reduce the memory limitation is off course to make the
A matrix more sparse (it is about 25% nnz now) by truncation of the
smallest elements in A12, which is say 1/100 of the largest. But
clearly this has to be done with care, since it is actually ignoring
some (small) physical effects present.
Best regards!
Magnus
Roger B. Sidje
>
> Dear exp(A) experts!
>
> Thanks for your ideas and comments!
>
> Why do I need the full exp(A) matrix?
> I'm solving linear ODE's for apples, pears and bananas... i.e. symbols
> instead of c-numbers.
> That is I express a certain symbol (creation and annihilation
> operators for atoms) at time t (assumed to be one below) as a linear
> combination of such symbols at time zero. It is therefore the rows of
> the exp(A) matrix that are interesting for me, since the scalar
> product of certain rows will give me my observables, i.e. certain
> combinations of pairs of symbols.
> If you (R. B. S.) are still working at U.Q. I can come in for an
> informal presentation of my problem at a group meeting or so in a
> couple of weeks if you like, since I'm going to be down in Brisbane
> for a couple of months now.
Too bad, I am not at UQ anymore. I am at the University of Minnesota
in Minneapolis for the time being.
>
> The method of calculating exp(A)*v, v being 'all Cartesian base
> vectors' indeed gives the full exp(A) matrix. In my case I'm using the
> rows, so it is better to calculate exp( transpose(A) )*v. It takes
> substantially longer time then expm(full(A)), but can handle larger
> matrices :-)
>
> I still think may be ways to simplify the problem due to the
> symmetries of the A matrix, such that I do not need to calculate
> exp(A) of the full A matrix?
I get that feeling too. But as always, the devil is in the details, so
let's try to work them out to see if anything comes up. You said
above:
>That is I express a certain symbol (creation and annihilation
> operators for atoms) at time t (assumed to be one below) as a linear
> combination of such symbols at time zero. It is therefore the rows of
> the exp(A) matrix that are interesting for me,
That is, doing w_i' = e_i' * exp(A), to get the i-th row of exp(A).
Then,
> since the scalar
> product of certain rows will give me my observables, i.e. certain
> combinations of pairs of symbols.
That is, doing w_i' * w_j, for some selected i, j, ... right?
And putting all together, you are doing
scalar = (e_i' * exp(A)) * (e_j' * exp(A))' = (exp(A')*e_i)' * (exp(A')*e_j),
for selected pairs of (known) indices i and j.
So the first optimization that comes to mind is whether you really
need exp(A) in full. That would be the case only if i and j run
through 1:n, meaning with E=exp(A), you need half of ALL the entries
of E' * E. That would indeed be surprising (in this context of large
sparse problems). If i and j are restricted to a few combinations,
there is definitely some savings by not computing exp(A) in full.
[Another thing is to avoid the undue A'. That is, reframe the problem
to label the symbols column-wise. There is not particular reason to
drag transpose(A) in the computations. Yeah, it means a different
mindset from what you started, but if it saves a few cycles, I am
taker... especially considering that it is only a one-off deal -- when
generating the actual data for computational purposes -- or the
transpose can simply be embedded in your tailor-made matrix-vector
product, see below.]
> If I write A in terms of four block matrices
> A=[A11 A12; A21 A22]
> There are additional symmetries
> A22=conj(A11)=-A11 since A11 is pure imaginary (and also diagonal)
> A21=q*conj(A12), where q=+-1, since A12 is also real for most
> interesting cases (*): A21=q*A12, i.e.:
> A=[A11 A12; q*A12 -A11]
>
> (*) Then A is antihermitian, and "The matrix exponential of an
> antihermitian matrix is a unitary matrix".
>
> Now my numerical calculations of the matrix M=exp(A)=[M11 M12; M21
> M22] (that is neither real or sparse) suggests that:
> M22=conj(M11) and M21=q*conj(M12) (**)
> It is always the case that M^\dagger=M^{-1}, i.e. M is unitary as is
> often the case in quantum physics (though my approach involves some
> approximations, so it was not clear for me a priori)
> In my view (**) indicates that a simpler procedure (e.g. to only
> calculate exp of block matrices) may be possible? (Besides it have
> already reduced the calculation of M=exp(A) to half the number of
> rows.)
Your A matrix suggests that the Lanczos approach could be apdated to
your case (it is only available in the Fortran version). The Matlab
version only has the Arnoldi approach. But the other optimization of
not computing exp(A) in full should already be significant.
>
> Another thing that could lead to improvements are the way the (sparse)
> matrix is represented in Matlab.
> In general a matrix with some (e.g. only one) complex or pure
> imaginary elements uses two times the memory of a pure real matrix.
> Since at the moment it is the size of the memory for the sparse A
> matrix that is my limitation (i.e. the RAM memory of the computer) and
> my matrix is real (in most cases of interest) except for the diagonal
> which is always pure imaginary, I then guess one could maybe modify
> expv.m such that it takes as input one real sparse matrix and one
> complex vector (the imaginary diagonal), in this way the memory
> required for sparse(A) would be a factor of two less.
Feel free to do so. That's the point of having the (open) source. You
can change expv to myexpv(t, diagA, offA, ...), and define your
matrix-vector product as needed. (Only a couple of places to tune in
expv to suit your new A*v, or rather transpose(A)*v, with only
offA=sparse(A12) and diagA the vector that defines A11.)
>
> I have also tried to perform my calculations with 'single precision'
> i.e. 'A=single(sparse(A))' with good result, meaning that I did not
> see any differences in the elements of M=exp(A) (checked only a few
> cases...) within 'format short'. My calculations will result in
> 'qualitative images' rather than precise numbers, maybe it is then OK
> to use single precision which reduce the memory limitations by a
> factor of two. Do you have any experiences of how good/bad single
> precision is for exp(A)?
Nope. I haven't seen any good or bad stories about single precision
usage of the package. However, it is a tad slower because the hardware
is 32bit (double precision), and trying to do single precision (16bit)
means that padding is needed to actually fit the data in the physical
registers when comes the crunch time. Such a padding brings some
overhead. In the old days with 16bit in the hardware, it used to be
that double precision was dearer because two registers will be needed
for computations. But the situtation is now reversed with 32bit in the
hardware, although those dazzingling fast CPUs that we now have make
the cost of padding to be only noise. So you can do that if it saves
memory. But I would rather play safe (in double precision) and tune
the matrix-vector product per above.
---
RBS
Roger B. Sidje
(real*4) is 32-bit, while double precision (real*8) is 64-bit. The
later has been pretty common on supercomputers, but has just started
to make its way into some of today's personal systems. [I wrote
something about this a while back somewhere else w.r.t the PCarray
(PowerChallenge array --
http://www.maths.uq.edu.au/~rbs/sgi/node5.html#SECTION00032000000000000000
]
---
RBS