Path: news.mathworks.com!not-for-mail
From: "Bruno Luong" <b.luong@fogale.findmycountry>
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 17:35:07 +0000 (UTC)
Organization: FOGALE nanotech
Lines: 35
Message-ID: <h9o7ob$j3l$1@fred.mathworks.com>
References: <h9llp6$mho$1@fred.mathworks.com> <h9o54a$pnq$1@fred.mathworks.com>
Reply-To: "Bruno Luong" <b.luong@fogale.findmycountry>
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 1254072907 19573 172.30.248.37 (27 Sep 2009 17:35:07 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Sun, 27 Sep 2009 17:35:07 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 390839
Xref: news.mathworks.com comp.soft-sys.matlab:573238


"Hua Wang" <ehwang@163.com> wrote in message <h9o54a$pnq$1@fred.mathworks.com>...

> 
> I am not sure whether it is clear enough for the explanation of the covariance matrix. Sorry that my english is not good enough.

Perfectly clear. You have actually confirmed the two cases that I assumed above where the covariance matrix is sparse - namely distance (within an image) and clusters (different images).

I have no idea by which trick Rune want to employ to reduce it to a 10 x 10 matrix.

Factorize can be found here

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

If I was you, I'll use an iterative method to solve weighted A*x = b.  It is equaivalent to minimize the quadratic function:

J(x) := 1/2 x'*A'*inv(C)*A*x - A'*inv(C)*b

Let's denote
H = A'*inv(C)*A
d = A'*inv(C)*b (computed as  A' * (C\b) )

J(x) = 1/2 x'*H*x - d

Now use PCG function to minimize J (or solve H*x=d), but you don't need to compute explicitly H. All you need is a function that is able to calculate H*x for any given x. This can be done by few simple lines:

d = A' * (C\b);
hfun = @(x) A'*(C \ (A*x)); % compute H*x
pcg(hfun, d, ...)

As Tim suggest, and if you can afford memory to store the factorization, you can use FACTORIZE or INVERSE in replace of two calls of (C \ y) by :

S = inverse(C) % <- call once
S*y % instead of (C \ y)

Bruno