Path: news.mathworks.com!not-for-mail
From: "Tim Davis" <davis@cise.ufl.edu>
Newsgroups: comp.soft-sys.matlab
Subject: Re: how to deal with the inversion problem of a huge sparse covariance matrix
Date: Sun, 27 Sep 2009 01:45:03 +0000 (UTC)
Organization: University of Florida
Lines: 28
Message-ID: <h9mg2v$9gk$1@fred.mathworks.com>
References: <h9llp6$mho$1@fred.mathworks.com>
Reply-To: "Tim Davis" <davis@cise.ufl.edu>
NNTP-Posting-Host: webapp-02-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1254015903 9748 172.30.248.37 (27 Sep 2009 01:45:03 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Sun, 27 Sep 2009 01:45:03 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 45902
Xref: news.mathworks.com comp.soft-sys.matlab:573146


"Hua Wang" <ehwang@163.com> wrote in message <h9llp6$mho$1@fred.mathworks.com>...
> Dear All,
> 
> I am processing some data using weighted least-squares (WLS) method. As you know, the solution of a system function A*x=b with weighting matrix P is that: x=inv(A'*P*A)*(A'*P*b);
> 
> I am solving the system function using Cholesky decomposition. The formula is:
> w=cholinc(P,'0');
> x=(w*A)\(w*b);
> 
> I am not sure whether the above method is possible for the back slash operator. But first of all, I met the problem of calculating weight matrix P from the covariance matrix C, i.e. P=inv(C). It is almost impossible to get the inverse of the huge sparse matrix C.
> 
> Could anybody give me some suggestions to solve such equation A*x=B, giving huge sparse covariance matrix C? It is urgent for me!

The short answer is never multiply by the inverse.  Ever.  If P=inv(C), then do this (first cut, not as efficient as possible, but much better than inv):

x = (A'*(C\A))\(A'*(C\b)) ;

but this is less efficient than using the Cholesky factorization of C.  Using the toolbox "Don't let that inv go past your eyes, to solve that system, FACTORIZE":

F=factorize(C);
x=(A'*(F\A))\(A'*(F\b)) ;

You can find the FACTORIZE toolbox here:

 http://www.mathworks.com/matlabcentral/fileexchange/24119

however, it would be better to reformulate this to use a sparse QR factorization
instead.  You're solving with the normal equations and Cholesky factorization, which is not as accurate as using a sparse QR.  If A is rectangular, x=A\b uses a sparse QR (be sure to use R2009b for that, for best performance; otherwise it will be very slow), whereas x=(A'*A)\(A'*b) uses sparse Cholesky for the Normal Equations (fast but not as accurate).