Thread Subject: ilu vs luinc

Dear All,

I'm solving a large linear system of equations as
B*x=y;

I got the solution about two years ago by the suggestions of some experts in the forum by the bicg method as following (where vcm is the variance-covariance matrix of y):

tol=1.0e-6;
maxit=500;
d=B'*(vcm\y);
nbb=B'*(vcm\B);
[l1,u1]=luinc(nbb,tol);
x=bicg(nbb,d,tol,maxit,l1,u1);

I always works well in the past two years. I got to know that ilu works faster than luinc, and there will be a warning message while using matlab 2011 because luinc will not be supported in the future.

I tried icu by [l2,u2]=icu((nbb,struct('droptol',1.0e-6)); I had thought it should equal to the above luinc function. It is indeed much faster, but the solution of x is obvious wrong. Bicg can not be convergent after using icu, but it converge after a couple of iteration while using luinc. The maximum difference of L and U (e.g. max(max(l1-l2))) are about 1000.

Something must be wrong for my usage in icu. Could you please shed light on it? Many thanks!

Another question is about the inverse of nbb, a large sparse matrix. I know it is an old question for inversion of large sparse matrix. And I did find a solution to invert nbb by Tim Davis. http://www.mathworks.es/matlabcentral/newsreader/view_thread/158111
But the solution can only get the diag element of nbb, e.g. the variance of x. Is there a way to get the off-diag elements (covariance of x)?

Hua

> I tried icu by [l2,u2]=icu((nbb,struct('droptol',1.0e-6));

I tried icu again by [l2,u2]=ilu(nbb,struct('type','ilutp','droptol',1.0e-6)). It gave the same results as luinc. They both used about 8 seconds. I got to know that 'ilutp' use the same method with luinc. Is there a faster way to get the same (proper) answer using ilu?

Comments for inverse of nbb is still very welcome. I'm thinking if it is possible to get the inverse if nbb based on l2/u2?

Thanks!
Hua

"Hua Wang" <ehwang@163.com> wrote in message <j42lrc$1mm$1@newscl01ah.mathworks.com>...
>
>
> Another question is about the inverse of nbb, a large sparse matrix. I know it is an old question for inversion of large sparse matrix. And I did find a solution to invert nbb by Tim Davis. http://www.mathworks.es/matlabcentral/newsreader/view_thread/158111
> But the solution can only get the diag element of nbb, e.g. the variance of x. Is there a way to get the off-diag elements (covariance of x)?

No, for the simple reason that because if the matrix is large and its inverse is likely to be full therefore method will run into the storage of such matrix.

However under some assumption, the inverse could approximate with sparse matrix. A starting reference is: http://dl.acm.org/citation.cfm?id=317422

Bruno

To inverse the sparse normal matrix, I find that it is much faster to use cholinc.

It takes 15s for the following codes:

tc=cholinc(nbb,1e-6);
lp=inv(tc');
vcm1=lp'*lp;

It takes about the same time to replace cholinc by ichol.
tc=ichol(nbb,struct('type','ict','droptol',1e-6)).
lp=inv(tc');
vcm2=lp'*lp;

It takes about 30s for vcm3=inv(nbb).

Although the same droptol is used in cholinc and ichol, the peak-to-peak value of vcm1-vcm3 is about 0.04, but that of vcm2-vcm3 is about 0.4. As you know, cholinc will not be supported in the future version. How can I get the same accuracy while suing ichol?

Thanks!

> It takes about the same time to replace cholinc by ichol.
> tc=ichol(nbb,struct('type','ict','droptol',1e-6)).
> lp=inv(tc');
> vcm2=lp'*lp;

I am processing a large sparse matrix for nbb, the normal function nbb (8953 by 8953). I need its inverse to get the uncertainty of the unknown parameters. It's reasonable for the first two steps (ichol/inv). But it is very slow for the last multiplication. Is there a fast way to do so? It's urgent! Thanks!

"Hua Wang" <ehwang@163.com> wrote in message <j469kc$j19$1@newscl01ah.mathworks.com>...
> To inverse the sparse normal matrix, I find that it is much faster to use cholinc.
>
> It takes 15s for the following codes:
>
> tc=cholinc(nbb,1e-6);
> lp=inv(tc');
> vcm1=lp'*lp;
>
> It takes about the same time to replace cholinc by ichol.
> tc=ichol(nbb,struct('type','ict','droptol',1e-6)).
> lp=inv(tc');
> vcm2=lp'*lp;
>
> It takes about 30s for vcm3=inv(nbb).
>
> Although the same droptol is used in cholinc and ichol, the peak-to-peak value of vcm1-vcm3 is about 0.04, but that of vcm2-vcm3 is about 0.4. As you know, cholinc will not be supported in the future version. How can I get the same accuracy while suing ichol?
>
> Thanks!

Hi Hua.

Note that by default ichol returns lower triangular factors while cholinc returns upper triangular factors. To get upper triangular factors out of ichol, add the 'shape' parameter to your call:

tc=ichol(nbb,struct('type','ict','droptol',1e-6,'shape','upper'));

On the other hand, you could (and probably should) use the lower triangular factor natively and use
lp = inv(tc); % Note the absence of the '
instead of
lp = inv(tc');
in your computation. Of course, you should be careful considering the inverses of any of these factors, especially when the matrix you are considering may be poorly conditioned. Also, as others have pointed out, the complete inverse is very likely to be dense (hence the time required by inv).

Take care.
Pat.

> On the other hand, you could (and probably should) use the lower triangular factor natively and use
> lp = inv(tc); % Note the absence of the '

Thanks Pat!

I noticed this problem, and used the correct codes like you suggested. But it is still very slow to do inversion for the lower triangular matrix inv(lp). Now I am using the method in Factorize package to solve it.

[l,g,p]=chol(nbb,'lower');
chold=p*(l'\(l\(p'*speye(size(nbb)))));

It's faster than explicitly inverse the lower triangular matrix, although it is still not satisfactory as speye used.