SmartInv

Large sparse matrix inversion. Returns block diagonal, tridiagonal or pentadiagonal elements.

Description

Returns block mono, tri or penta diagonal elements of the inverse of a symetric square matrix.

Useful for large sparse matrix which LU decomposition is easy to compute for. Not the fastest way to compute inverse matrix, but avoid memory problem of a full matrix storage.

Optional progressive diagonal computation display. Useful to quickly observe an important modification on the diagonal.

% Q = smartinv(N) returns the N^-1. N is a square symetric matrix n x n.
%
% Q = smartinv(N,blocksize) returns the block (blocksize x blocksize)
% diagonal elements of N^-1.
%
% Q = smartinv(N,blocksize,type) returns the block type-diagonal
% elements of N^-1 type can be mono, tri or penta. The type
% determines the inclusion of diagonal block computation,
% i.e there is redundancy for tri (25%) and for penta (~45%).
%
% Q = smartinv(N,blocksize,position) returns the block diagonal elements
% of N^-1 at position. Q is blocksize x blocksize if position
% is a single value or n x n if posistion is a vector. If
% position is not set, the block values will be computed in an
% order so the curve shape will appear progressively (if graph
% displaymode is used).
%
% Q = smartinv(N,blocksize,...,'displaymode',displaymodevalue) specifies
% process graphical display: none, waitbar (default), graph.
% If position exists and is scalar, no display is set up at all.
% Graph mode will draw a curve for each diagonal element of block
% (i.e. a block represents an epoch, and each element of a block
% a variable to display).
%
% Q = smartinv(N,blocksize,...,'linelegend',linelegendvalue) where
% linelegendvalue is a cell array of string specifying legend
% strings to display for each block element (used for graph
% display mode).
%
% Q = smartinv(N,blocksize,...,'linecolor',linecolorvalue) where
% linecolorvalue is a cell array of color string or RGB value
% to display for each block element (used for graph display mode).
%
% Q = smartinv(N,blocksize,...,'linefunction',linefunctionname) where
% linefunctionname is a function called to display value in
% graph displaymode. linefunctionname can be a string (the same
% function will be used for each element in the block) or a cell
% array of strings of blocksize length.

• MATLAB release: MATLAB 7.5 (R2007b)
• Other Files: license.txt, ; smartinv.m

[Tim Davis]

This function attempts to do the right thing when computing selected entries of the inverse, by solving a sequence of linear equations after doing an LU factorization. However, the implementation is flawed.

One suggestion: the statement [L,U,P]=lu(A) does not allow LU to permute A to improve sparsity. It would be better to use the 4-output LU, [L,U,P,Q]=lu(A), except in your case you'd need to use another name for Q since you have a Q variable already.

You would then need to use the column permutation Q in the forward/backsolve steps.

Including a "waitbar" in a function like this is not a good idea.

There seems to be a bug. The documentation says that s=smartinv(A) returns the diagonal of the inverse, but in fact I get the whole inverse back.

If I do s=smartinv(A,1), I get a diagonal matrix, but not the diagonal entries of the inverse, which is what the documentatioon said I would get. (I tried with this matrix, which is included in MATLAB):

load west0479
A=west0479;

With smartinv(A,1) or smartinv(A,1,'mono'), I should get the diagonal of inv(A) as the result, but I don't. SOmething else is returned.

Finally, the method is very slow. Computing the entire inverse takes 0.13 seconds using S=inv(A). Using S=smartinv(A) I get the same result, but it takes 15 seconds.

In the profiler, all the time is spent in this one line of code:

 Q(q0+1:q0+w,q0+1:q0+w) = insblock;

that's a very bad idea. You should never use sparse submatrix assignment, if at all possible. A better approach is to build the matrix all at once, from the blocks, using the "sparse" command.

So there appears to be at least three serious bugs in the code:
(1) does not perform according to the documentation (smartinv(A) returns inv(A), not diag(inv(A))).
(2) very slow as compared with just "inv"
(3) cannot compute the diagonal of the inverse, in spite of the fact that the documentation says it can do so. Instead it returns something else altogether.

And 2 serious drawbacks:
(1) Not really a bug, but a serious performance issue: doesn't allow for pivoting to maintain sparsity, and is thus unsuitable for large matrices.
(2) Remove the waitbar. If your function is used within other codes, the end user isn't interested in the "Inverting matrix..." message.

[Rouzaud Denis]

Hello,
Thanks for your message and to take the time to explain!

Well, effectively I compute inv(west0479) and it does not return the same as smartinv(west0479).
I missed an important explanation in the documentation. I made this code for blcok symetric matrices. I tried with the ones I work on, and I got differences about e-17.
Do you agree with the method for symetric matrices ?

I added the [L,U,P,Q] decomposition. It effectively seems to improve the computation time about 30% (I tried with a135000x135000 matrix). Thanks !

About the drawbacks, I removed the default option (displaymode) of the waitbar and I corrected Q = smartinv(N) definition. However smartinv(N) is really not what this function is made for. If can invert it in one single step, smartinv is really useless!

I know this code is very slow. But it allwos the computation of the block pentadiagonal elements of the inverse of a matrix 135000x135000, which would need 165 Go of RAM to be fully compute.

Update is following.
Thanks again
Denis

[Tim Davis]

There's still something wrong with this submission. I computed the exact inverse of this matrix:

http://www.cise.ufl.edu/research/sparse/matrices/HB/bcsstk28.html

and then extracted the tridiagonal entries of the inversion (entries on the diagonal, subdiagonal, and +1 diagonal).

I then computed S = smartinv(A,1,'tri'), excepting to get the same thing.

But I don't I get a matrix with zeros on the +1 and -1 diagonals, every 15 entries or so. So this revised submission is still not computing what it says it computes.

[Rouzaud Denis]

Hello again, thanks again.
 
There is an indexing used for the displaymode 'graph' that is used to display progressively the diagonal of the inverse while computation. There was an error in the use of the unique function, which is now corrected. It was leading to missing values in the computed inverse!
 
It seems that the computation works for non symetric matrices. Here are the result I get now.
 
Let me know what you think.
 
Thanks for your help.
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
load bcsstk28
N=Problem.A;
sm = smartinv(N,1,'tri');
sq = inv(N);
dm = sq-sm;
dmax = abs(dm(1,1));
P = [ 1 1];
for i = 2:size(dm)
if abs(dm(i,i))>dmax,dmax=abs(dm(i,i));p=[i i];end
if abs(dm(i-1,i))>dmax,dmax=abs(dm(i-1,i));p=[i-1 i];end
if abs(dm(i,i-1))>dmax,dmax=abs(dm(i,i-1));p=[i i-1];end
end
dm(p(1),p(2))
sq(p(1),p(2))
 
ans =
    2.664812814856532e-013
 
ans =
   0.033834948807814
 
clear
load west0479
N = west0479;
sm = smartinv(N,1,'penta');
sq = inv(N);
dm = sq-sm;
dmax = abs(dm(1,1));
P = [ 1 1];
if abs(dm(1,2))>dmax,dmax=abs(dm(1,2));p=[1 2];end
if abs(dm(2,1))>dmax,dmax=abs(dm(2,1));p=[2 1];end
for i = 3:size(dm)
if abs(dm(i,i))>dmax,dmax=abs(dm(i,i));p=[i i];end
if abs(dm(i-1,i))>dmax,dmax=abs(dm(i-1,i));p=[i-1 i];end
if abs(dm(i,i-1))>dmax,dmax=abs(dm(i,i-1));p=[i i-1];end
if abs(dm(i-2,i))>dmax,dmax=abs(dm(i-2,i));p=[i-2 i];end
if abs(dm(i,i-2))>dmax,dmax=abs(dm(i,i-2));p=[i i-2];end
end
dm(p(1),p(2))
sq(p(1),p(2))
 
ans =
    6.366462912410498e-011
 
ans =
    5.218763111339072e+003