Thread Subject: calculating the inverse efficiently (not for solving equations :-) )

Hi,

well most of you probably think "Oh no, not another one
trying to solve a system of linear equations using the
matrix inverse", but luckily for those guys I don't need it
for that :-)

So my problem is that after some calculations I got stuck
with an equation of the following kind:

(A*B^-1*C)*x=y

As the matrix B is multiplied with A and C there's no way
(or at least I don't know of any) to get around the problem
of knowing the full inverse of B.

Some important properties for B are, that it's symmetric and
sparse. Yes, I know having the inverse will destroy the
sparisty :-)


So far the possibilities for calculating the inverse of B, I
have come up with are:
*LU decomposition: tends to fail because B is rather large
and therefore my computer runs out of memory
*Cholesky decomposition: not much more efficient than the LU
decomposition (n^3/3 compared to 2n^3/3)

I also thought of using eigenvalues, but that wouldn't help
much since you need the inverse of the matrix with the
eigenvectors.

I hope that I could give anyone enough insight into my
problem...

Does anyone know another efficient way for calculating the
inverse for sparse symmetric matrices, or maybe to simplify
or transform the upper equation to make it less costly?

Thanks in advance!!

Albert Buehrli

"M E" <gitsnedbutzi@hotmail.com> wrote in message <fsu1r6$ib3
$1@fred.mathworks.com>...
> ........
> I also thought of using eigenvalues, but that wouldn't help
> much since you need the inverse of the matrix with the
> eigenvectors.
> ........
> Albert Buehrli
---------
  If B is symmetric, real, and non-singular, you CAN use eigenvectors to find its
inverse. Just take the reciprocals of its eigenvalues:

 [v,d] = eig(B);
 Binv = v*diag(1./diag(d))*v';

  However, I am not claiming this is necessarily the most efficient way to solve
your problem.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid>
wrote in message <fsu4bu$l24$1@fred.mathworks.com>...
> If B is symmetric, real, and non-singular, you CAN use
eigenvectors to find its
> inverse. Just take the reciprocals of its eigenvalues:
>
> [v,d] = eig(B);
> Binv = v*diag(1./diag(d))*v';
>
> However, I am not claiming this is necessarily the most
efficient way to solve
> your problem.
>
> Roger Stafford
>

Indeed, yes you're right. Forgot that for symmetric matrices
you don't need the inverse of V, but only the transpose.

Maybe still not the most efficient way, but at least I'm one
step closer to an efficient solution.

Thanks!!

> Indeed, yes you're right. Forgot that for symmetric matrices
> you don't need the inverse of V, but only the transpose.
>
> Maybe still not the most efficient way, but at least I'm one
> step closer to an efficient solution.
>
> Thanks!!

Btw. has anyone experience with how long it takes to
calculate all eigenvalues of the large sparse matrix? Is
this fast?

"M E" <gitsnedbutzi@hotmail.com> wrote in message
<fsu6on$qaa$1@fred.mathworks.com>...
> > Indeed, yes you're right. Forgot that for symmetric matrices
> > you don't need the inverse of V, but only the transpose.
> >
> > Maybe still not the most efficient way, but at least I'm one
> > step closer to an efficient solution.
> >
> > Thanks!!
>
> Btw. has anyone experience with how long it takes to
> calculate all eigenvalues of the large sparse matrix? Is
> this fast?

No, I don't think this is the right way to go.

You don't want to use eig to do this if your
matrix is large and sparse, as the eigenvector
set will NOT be sparse, and you need them
too.

Best is probably to use chol, perhaps with a
permutation to minimize the fill-in. I think
that the current best method is amd if I
remember right. A quick test shows that for
some matrices, the difference can be
significant.

n = 1000;
A = sprand(n,n,.001);
A = (A + A') + 3*speye(n,n);

timeit(@() chol(A))
ans =
      0.56226

p = amd(A);
timeit(@() chol(A(p,p)))
ans =
    0.0038918

timeit(@() amd(A))
ans =
    0.0014697

R = chol(A);
nnz(R)
ans =
       31133

r(:,p) = chol(A(p,p));
nnz(r)
ans =
        3940

Most of the time, the difference will not be
that dramatic. You will still need to undo the
permutation. This should work:

R(:,p) = chol(A(p,p));

HTH,
John

"M E" <gitsnedbutzi@hotmail.com> wrote in message
<fsu1r6$ib3$1@fred.mathworks.com>...
> Hi,
>
> well most of you probably think "Oh no, not another one
> trying to solve a system of linear equations using the
> matrix inverse", but luckily for those guys I don't need it
> for that :-)
>
> So my problem is that after some calculations I got stuck
> with an equation of the following kind:
>
> (A*B^-1*C)*x=y
>
> As the matrix B is multiplied with A and C there's no way
> (or at least I don't know of any) to get around the problem
> of knowing the full inverse of B.
>
> Some important properties for B are, that it's symmetric and
> sparse. Yes, I know having the inverse will destroy the
> sparisty :-)
>
>
> So far the possibilities for calculating the inverse of B, I
> have come up with are:
> *LU decomposition: tends to fail because B is rather large
> and therefore my computer runs out of memory
> *Cholesky decomposition: not much more efficient than the LU
> decomposition (n^3/3 compared to 2n^3/3)
>
> I also thought of using eigenvalues, but that wouldn't help
> much since you need the inverse of the matrix with the
> eigenvectors.
>
> I hope that I could give anyone enough insight into my
> problem...
>
> Does anyone know another efficient way for calculating the
> inverse for sparse symmetric matrices, or maybe to simplify
> or transform the upper equation to make it less costly?
>
> Thanks in advance!!
>
> Albert Buehrli
>
>

If you ever find yourself multiplying by the inverse, then
you know one thing for certain. You know, for sure, that
you don't know what you're doing. One needs never, ever,
ever multiply by the inverse. Multiplying by the inverse is
mathematically equivalent to solving a system of linear
equations, but it is an awful computational replacement for
the latter.

You really *are* solving a system of equations, either B*z=C
for z, or you are solving A=w*B for w. Both are systems of
linear equations.

Instead of inv(B)*C, you should use B\C instead. Or,
replace A*inv(B) by A/B.

You gave as the equation:

(A*B^-1*C)*x=y

does this mean you know A, B, C, and x, and you want to
compute y? Or are you solving for x? If it's the latter,
then what you are should do is this:

x = ((A/B)*C)\y

or

x = (A*(B\C))\y

If you have x and want to compute y, then do

y = ((A/B)*C)*x

or

y = (A*(B\C))*x

"Tim Davis" <davis@cise.ufl.edu> wrote in message <fsukql$df3
$1@fred.mathworks.com>...
> If you ever find yourself multiplying by the inverse, then
> you know one thing for certain. You know, for sure, that
> you don't know what you're doing. One needs never, ever,
> ever multiply by the inverse. Multiplying by the inverse is
> mathematically equivalent to solving a system of linear
> equations, but it is an awful computational replacement for
> the latter.
> ...........
---------
  Tim, it is best not to make claims that can be misunderstood, as in: "If you
ever find yourself multiplying by the inverse, then you know one thing for
certain. You know, for sure, that you don't know what you're doing. One
needs never, ever, ever multiply by the inverse." There is a Murphy's Law-like
principle that states that someone, sometime, somewhere is surely going to
come up with an example for which multiplying by the inverse of a square
matrix is the very, very best way of carrying out a given task.

  How about this one for starters. We are given a non-singular 3 x 3 square
matrix M that transforms 3 x 1 vectors, x, to 3 x 1 vectors, y, in 3D space.
We are also given a point cloud of 1,000,000 different points of y-values and
are to determine the x-point cloud they came from by way of multiplication
by M. I claim the very best way to do this is to determine, once and for all, by
whatever means are necessary, fair or foul, the explicit, full inverse of M, and
then, having accomplished that despicable operation, to then proceed with
the 1,000,000 matrix multiplications of it by y-points necessary to determine
all the x-points.

  At this juncture, I can envision you claiming, "Foul, foul, that isn't what I
meant in my statement!" But that is exactly the point of my comment here.
Your extreme statement could easily be interpreted in just such a manner. It
is always best to carefully quantify such statements so as not to allow of such
misinterpretations.

Roger Stafford

@John:

Thanks for your hints on using amd. I'll have to try that
one.

@Tim:

Unfortunately, the mldivide also runs out of memory, so I
can't use that one.


Since memory seems to be the main issue, I also thought of
using blockwise inversion. There I would only have to
compute the inverse of A and (D-C*A^-1*B), where I could
reduce the size from n to n/2.

For those who don't know it:
[A, B; C, D]^-1 = [A^-1+ A^-1*B*(D-C*A^-1*B)^-1*C*A^-
1, ...; ..., ...]

Has anyone any experience on using blockwise inversion?

Thx in advance