Thread Subject: how to deal with the inversion problem of a huge sparse covariance matrix

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!

Thanks in advance!

Hua

The inverse of a very large and sparse matrix is going to be just as
large as the original matrix. However, it is extremely likely to be
very dense and thus impractical to store.

You'll need to restructure your computations so that you don't
explicitly construct P=inv(C). One strategy would be to make use of
the Cholesky factorization of C, which could be (especially after
reordering rows and columns of C) reasonably sparse. Another approach
would be to adopt an iterative method for the solution of the least
squares problem. Without knowing more about exactly what you're
trying to accomplish, it's hard to provide more explicit advice.

"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).

On 26 Sep, 20:16, "Hua Wang" <ehw...@163.com> wrote:
> 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!

I'm a bit curious what kinds of situations occur where you
need to handle a huge sparse covariance matrix.

First of all, most processes of practical interest are
stationary. You don't need the long-term covariances,
which means that you can get away with a smaller covariance
matrix.

Second, the cyclostationary processes that do have some
significant long-term covariance can be handled by specialized
higher-order techniques, so you don't need the *covariance*
matrix.

The remaining processes, that are not stationary or cyclo-
stationary, would require all the covariance data to be
available, so you would not use a *sparse* covariance matrix.

All in all, I think you might want to look over the problem
statement, to see if there are other approaches than the naive
go-straight-ahead-and-implement-the-literal-expressions.

That particular approach is almost always wrong.

Rune

Rune Allnor <allnor@tele.ntnu.no> wrote in message <ace8cc81-2737-4630-83b6-76c09f427a97@p23g2000vbl.googlegroups.com>...

>
> I'm a bit curious what kinds of situations occur where you
> need to handle a huge sparse covariance matrix.

many situations:

- data are organized in clusters and variables are only correlated within a clusters.
- two data are correlated with they distance (from a certain metric) is smaller than some threshold, and zero otherwise.
...

>
> First of all, most processes of practical interest are
> stationary. You don't need the long-term covariances,
> which means that you can get away with a smaller covariance
> matrix.

You seem to infer *temporal* dimension in the OP problem. Aren't probably mistaken between correlation (time is no needed) and cross/auto-correlation ?

Bruno

On 27 Sep, 14:40, "Bruno Luong" <b.lu...@fogale.findmycountry> wrote:
> Rune Allnor <all...@tele.ntnu.no> wrote in message <ace8cc81-2737-4630-83b6-76c09f427...@p23g2000vbl.googlegroups.com>...

> > First of all, most processes of practical interest are
> > stationary. You don't need the long-term covariances,
> > which means that you can get away with a smaller covariance
> > matrix.
>
> You seem to infer *temporal* dimension in the OP problem.

No. 'Stationarity' means that you can draw any sub-sample
from a larger sample and the statistical properties will
be the same.

> Aren't probably mistaken between correlation (time is no needed) and cross/auto-correlation ?

Time has nothing to do with the statistical properties.
Assign any physical dimension you want to the data. I still
can't see why one would need large sparse correlation matrices.

Rune

Rune Allnor <allnor@tele.ntnu.no> wrote in message <524bfbe2-ff4e-4cc6-a2d0-aa05d821fbd1@l35g2000vba.googlegroups.com>...

>
> Time has nothing to do with the statistical properties.
> Assign any physical dimension you want to the data. I still
> can't see why one would need large sparse correlation matrices.
>

A very simple example is all the data are all independent. The correlation matrix is an identity matrix. May be you don't consider an identity matrix as sparse? Or you think there is never two data independent?

I would even argue in large scale problem the correlation matrix is always sparse!

Bruno

On 27 Sep, 15:03, "Bruno Luong" <b.lu...@fogale.findmycountry> wrote:
> Rune Allnor <all...@tele.ntnu.no> wrote in message <524bfbe2-ff4e-4cc6-a2d0-aa05d821f...@l35g2000vba.googlegroups.com>...
>
> > Time has nothing to do with the statistical properties.
> > Assign any physical dimension you want to the data. I still
> > can't see why one would need large sparse correlation matrices.
>
> A very simple example is all the data are all independent. The correlation matrix is an identity matrix. May be you don't consider an identity matrix as sparse? Or you think there is never two data independent?

Assume you have, say, a 1,000,000x1 data vector with independent
samples. If one takes your argument ilterally, you seem to say that
one needs a 1,000,000x1,000,000 identity matrix to represent
the covariance.

My argument is that you only need a small portion of that, say,
a 10x10 matrix. You have the same information but in a far more
convenient form.

Rune

Thank you all so much for the prompt reply!!!

For Rune Allnor & Bruno Luong,

I am dealing with a huge dataset from satellite radar images. I will combine the data from a couple of images. In each image, the data are correlated in spatial domain with a wavelength of about 20 km. That's to see, there should be a full covariance matrix for each 20km, and they are zeros for the covariance of the pixels far away from 20 km. The pixels are uncorrelated in different images, so that the covariance is also zero. In my study, each image covers about 800 km by 100 km. That's why there is a huge sparse matrix.

It is also not as easy as copy covariance of the pixels for every 20 km by 20 km. You know the covaraiance of pixel x,y can be expressed as C(x,y)=sigma(x)*sigma(y)*corr_coef(x,y). In my case, if the distances for any two pixels are the same, the correlation coefficients are the same. But the variance of each pixel is different. Therefore, covariance is different for each 20km by 20 km patch.

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

For Tim Davis & Brian Borchers,

I am using matlab R2007a. When I use backslash for A\b, it is very slow and almost impossible for my data. Is it possible for do
F=factorize(C);
x=(A'*(F\A))\(A'*(F\b)) ;

By the way, I can not find fractorize function with my matlab. I think the signal processing toolbox was not installed by the computer administrator. But I can use qr function. Could you give me more suggestions? It is better if you could tell me some functions about solve the system function: A*x=b, with giving sparse covariance matrix C?

Thank you all again!!

Cheers,
Hua

"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

On 27 Sep, 19:35, "Bruno Luong" <b.lu...@fogale.findmycountry> wrote:
> "Hua Wang" <ehw...@163.com> wrote in message <h9o54a$pn...@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.

I haven't seen anything that indicates neither the size of
these images in terms of pixels, or the dimensions of the
covarinance matrices. Nor have I seen anything to indicate
what the purpose of the analysis is.

Again, any desire to invert a large matrix is almost always
wrong. In the few cases where inversion is warranted, the
naive inversion is almost never the correct algorithm.

Of course, it could be that the OP's application is among the
handful that remains, but it might be worth the effort to
investigate alternative algorithms.

Rune

Rune Allnor <allnor@tele.ntnu.no> wrote in message <ef773a28-76ca-4e60-b6ee-dc944f4e6ba0@h30g2000vbr.googlegroups.com>...

>
> Again, any desire to invert a large matrix is almost always
> wrong.

Finite element methods, solid mechanics, thermal calculation, fluid dynamic, seismic inversion, statistics, data mining, machine learning, data assimilation, weather forecast, quantum mechanics simulation, image reconstruction in medical, photogrammetry, etc... all often call for inversion of huge system of linear equations. As far as I know they all seem to work alright.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h9oai6$o47$1@fred.mathworks.com>...
> Rune Allnor <allnor@tele.ntnu.no> wrote in message <ef773a28-76ca-4e60-b6ee-dc944f4e6ba0@h30g2000vbr.googlegroups.com>...
>
> >
> > Again, any desire to invert a large matrix is almost always
> > wrong.
>
> Finite element methods, solid mechanics, thermal calculation, fluid dynamic, seismic inversion, statistics, data mining, machine learning, data assimilation, weather forecast, quantum mechanics simulation, image reconstruction in medical, photogrammetry, etc... all often call for inversion of huge system of linear equations. As far as I know they all seem to work alright.
>
> Bruno

Do these problems call for the inversion of a large matrix? Or do they simply call for the solution of a set of large (and almost always sparse) set of linear equations?

The latter does NOT imply a need for the former. Thus the statement "any desire to invert a large matrix is almost always wrong" is true.

One might write an analytic solution using the inverse of a matrix, but that does not imply that the inverse is the only or even preferred way to numerically compute a solution. The way equations are written on paper or in a text book are almost never the way they are solved on a computer if the problem size is sufficiently large.

FWIW

Duane

Duane

On 27 Sep, 20:23, "Bruno Luong" <b.lu...@fogale.findmycountry> wrote:
> Rune Allnor <all...@tele.ntnu.no> wrote in message <ef773a28-76ca-4e60-b6ee-dc944f4e6...@h30g2000vbr.googlegroups.com>...
>
> > Again, any desire to invert a large matrix is almost always
> > wrong.
>
> Finite element methods, solid mechanics, thermal calculation, fluid dynamic,

These are forward models...

> seismic inversion,

Manual interactive process...

> statistics, data mining, machine learning, data assimilation, weather forecast, quantum mechanics simulation,

Forward models or assembling statistics...

> image reconstruction in medical, photogrammetry, etc... all often call for inversion of huge system of linear equations.

So by a dozen or so applications you mention only a handful
actually are inversion problems, and of those at moset one
or two rely on the inversion of large matrices.

> As far as I know they all seem to work alright.

Most of them do. But the problems you list have very little
in common with the still non-specified problem the OP works
with.

Rune

>
> I haven't seen anything that indicates neither the size of
> these images in terms of pixels, or the dimensions of the
> covarinance matrices. Nor have I seen anything to indicate
> what the purpose of the analysis is.

Sorry that I forgot the pixel size. In my dataset, the original pixel size is about 90 m. It is to slow to process such high resolution data. After multilooking (i.e. average for the neighbor pixels), the pixel size is about 360 m. So there are about 300*2400 pixels for each image. After removing some pixels with NaN values, the dimension of the final covariance matrix is about 40,000 by 40,000 with 15 million non-zeros.

> Again, any desire to invert a large matrix is almost always
> wrong. In the few cases where inversion is warranted, the
> naive inversion is almost never the correct algorithm.
>
> Of course, it could be that the OP's application is among the
> handful that remains, but it might be worth the effort to
> investigate alternative algorithms.

In my study, I got deformation from different radar images. Because of different geometry of the radar images, the deformation field is not continuous from one image to the other. My idea is to divide my study area into triangular mesh, and invert the velocities on the nodes of the mesh using all radar data. If I know the radar geometry, it is reasonable to solve it.

I am also thinking to use cleverer algorithm, e.g. dividing the image into different patches. But I some parameters, e.g. orbital errors of the radar image, should be the same for the patches within a image. I am not clear how to add such constraints if I use patch strategy. Alternatively, could you give me some suggestions?

Thanks a lot!

Hua

On 28 Sep, 11:19, "Hua Wang" <ehw...@163.com> wrote:
> > I haven't seen anything that indicates neither the size of
> > these images in terms of pixels, or the dimensions of the
> > covarinance matrices. Nor have I seen anything to indicate
> > what the purpose of the analysis is.
>
> Sorry that I forgot the pixel size. In my dataset, the original pixel size is about 90 m. It is to slow to process such high resolution data. After multilooking (i.e. average for the neighbor pixels), the pixel size is about 360 m. So there are about 300*2400 pixels for each image. After removing some pixels with NaN values, the dimension of the final covariance matrix is about 40,000 by 40,000 with 15 million non-zeros.

The numbers alone go a very long way to indicate that the
problem statement is almost certainly wrong.

Rune

To Hua:

You probably consider the step of adding a preprocessing of the radar images data and compress them to a space of fewer dimensions before carry out the regression.

If your data is correlated spatially you might consider to project on few smoother modes. Take a look at references on EOF (Emperical Orthogonal Functions) where people do this kind of reduction/simplification without too much affected by the final accuracy.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h9q13q$pnf$1@fred.mathworks.com>...
> To Hua:
>
> You probably consider the step of adding a preprocessing of the radar images data and compress them to a space of fewer dimensions before carry out the regression.
>
> If your data is correlated spatially you might consider to project on few smoother modes. Take a look at references on EOF (Emperical Orthogonal Functions) where people do this kind of reduction/simplification without too much affected by the final accuracy.
>
> Bruno

Thanks Bruno!

I've done such redunction already. The original resolution of radar data should be 20 by 20 m. I can not reduce more data. The high resolution is one of the advantage that we use radar data. And the objective of my study is to get a high resolution velocity field using radar data.

I will try your suggestions about pcg. But I am not sure whether it is possilbe to get d=A'*(C\b) before pcg. Anyway, I will try it soon. Many thanks!

Cheers,
Hua

> > Sorry that I forgot the pixel size. In my dataset, the original pixel size is about 90 m. It is to slow to process such high resolution data. After multilooking (i.e. average for the neighbor pixels), the pixel size is about 360 m. So there are about 300*2400 pixels for each image. After removing some pixels with NaN values, the dimension of the final covariance matrix is about 40,000 by 40,000 with 15 million non-zeros.
>
> The numbers alone go a very long way to indicate that the
> problem statement is almost certainly wrong.
>
Sorry, what do you mean?

Hua,

another technique that you might consider is multigrid-like regression. I'm not expert on this but I have seen some people using with big success.

Of course there is multi-grid algorithm for matrix/linear inversion. I have no idea if such tool exists in Matlab.

Bruno

On 28 Sep, 12:50, "Hua Wang" <ehw...@163.com> wrote:
> > > Sorry that I forgot the pixel size. In my dataset, the original pixel size is about 90 m. It is to slow to process such high resolution data. After multilooking (i.e. average for the neighbor pixels), the pixel size is about 360 m. So there are about 300*2400 pixels for each image. After removing some pixels with NaN values, the dimension of the final covariance matrix is about 40,000 by 40,000 with 15 million non-zeros.
>
> > The numbers alone go a very long way to indicate that the
> > problem statement is almost certainly wrong.
>
> Sorry, what do you mean?

A 40k x 40k matrix will contain 1.6 billion elements.
Just storing those elements takes >>10 GBytes of memory.

I know, you have only a few nonzero elements in all this,
but just about any non-trivial operation on a sparse
matrix will destroy sparsity. If you try to invert the
matrix (or even just reformulate it) the above amounts of
memory will be required to hold the result. And a lot of
algorithms produce several similar-sized marices, either
as part of the end result or as intermediate variables.

Even if you had the memory required to handle all those
data, there is no reason to trust the computed results,
due to accumulated errors from numerical inaccuracies.

Rune

> > > > Sorry that I forgot the pixel size. In my dataset, the original pixel size is about 90 m. It is to slow to process such high resolution data. After multilooking (i.e. average for the neighbor pixels), the pixel size is about 360 m. So there are about 300*2400 pixels for each image. After removing some pixels with NaN values, the dimension of the final covariance matrix is about 40,000 by 40,000 with 15 million non-zeros.
> >
>
> A 40k x 40k matrix will contain 1.6 billion elements.
> Just storing those elements takes >>10 GBytes of memory.
>

Sorry that you misunderstood my explanation!
If the covariance matrix is 40k x 40k, but there is only 15 million non-zeros in this sparse matrix. Why does it need 10G memory to store? It should be 15M * sizeof(double).

Hua

On 28 Sep, 14:41, "Hua Wang" <ehw...@163.com> wrote:
> > > > > Sorry that I forgot the pixel size. In my dataset, the original pixel size is about 90 m. It is to slow to process such high resolution data. After multilooking (i.e. average for the neighbor pixels), the pixel size is about 360 m. So there are about 300*2400 pixels for each image. After removing some pixels with NaN values, the dimension of the final covariance matrix is about 40,000 by 40,000 with 15 million non-zeros.
>
> > A 40k x 40k matrix will contain 1.6 billion elements.
> > Just storing those elements takes >>10 GBytes of memory.
>
> Sorry that you misunderstood my explanation!
> If the covariance matrix is 40k x 40k, but there is only 15 million non-zeros in this sparse matrix. Why does it need 10G memory to store? It should be 15M * sizeof(double).

You need some overhead for some indexing stuff as well,
but yes. As long as the matrix remains sparse, you don't
need a lot of memory.

What I tried to explain in the rest of the post, is that
if you try to invert this sparse matrix (or do just about
anything else to it), the result is no longer sparse. And
this dense result will require such huge amounts of memory.

Rune

> 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, ...)
>

For Bruno,

Huge thanks to Bruno! I think your method can solve my problem. It tooks less than one minute and 2G memory to solve my system function when I use your above three lines. I am still not clear why it is so fast to solve d=A'*(C\b), but it can not solve W=C\I, where I is a unit matrix. Is it because C\b will not be stored for the calculation of d, but C\I has to be stored?

I've not get the final result yet. PCG gave an warning message that it could be convergent. But I do not think it is due to Bruno's method. In my plan, there is still something to do to improve my system model. I am so happy that it seems the above three lines can solve it at last.

Thank you all very much! I will let you know if I can get the final result.

Cheers,
Hua

> 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, ...)
>

I compared the above method with the following solutions:

w=cholinc(C,'0');
x=(w\A)\(w\b);

I found that Bruno's pcg result is the same with x=(A'*(C\A))\d, but it is different from the result by cholinc method. Which one is better?

Hua

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

>
> I found that Bruno's pcg result is the same with x=(A'*(C\A))\d, but it is different from the result by cholinc method. Which one is better?

If your matrices are ill-conditioned then the variation is expected. The ill-conditioning can come from A or C. I'm still surprised that your C matrix is singular. That means your data are redundant (in the radar framework it means the data resolution is coarser than the images).

If you don't put any other a priori knowledge in the solution of your inversion problem, then mathematically a better solution is the one that gives smaller 2-norm of the residue:
|x' A'*inv(C)*A*x - inv(C)*b|.

However the better solution (in the sense above) might be non-physical and not preferable. It is up to modeler to put the constraint of the norm, to regularize, etc... in order to obtain what he expected. There is no generic answer I'm afraid.

I finish with a note on a small technical trick: for 'y' a vector and 'A' a large sparse matrix, instead of calculate:

A' * y

It is faster to do

(y'*A)'

Bruno

Sorry for the typo, should read:

>... a better solution is the one that gives smaller 2-norm of the residue:
> |x' *A'*inv(C)*A*x - A'*inv(C)*b|.

|x' *A'*inv(C)*A*x - d|.

% Bruno

>I'm still surprised that your C matrix is singular. That means your data are redundant (in the radar framework it means the data resolution is coarser than the images).
>

Thanks Bruno!

How do you know whether C matrix is singular? I suppose you remember my another post about cholinc decomposition of C matrix. In that post, C is calculated by exp(-d/alpha), where d is distance between two pixels and alpha is a wavelength (about 20 km). If I calculate C using this formular, C is not spare matrix although most of the values are very small if d is much bigger than alpha. To save the memory, I set a threshold of d, for instance, 5 * alpha. So that C=0 if d>5*alpha. Then C becomes a sparse matrix. This is the reason that sometimes C will fail for cholesky decomposition.

I also find that if the threshold is 5*alpha, C usually will pass cholesky decomposition. So that i need not do EIGS correction any more.

I did not do EIGS correction, and do least-squares directly. Then the difference exists between pcg and cholinc methods.

Do you think the above theshold strategy will make C singular?

Hua

"Hua Wang" <ehwang@163.com> wrote in message <haki2n$14o$1@fred.mathworks.com>...
> >I'm still surprised that your C matrix is singular. That means your data are redundant (in the radar framework it means the data resolution is coarser than the
>
> Do you think the above theshold strategy will make C singular?

Not sure. For the convolution problem, usually the condition number increases when the kernel get broader. It seems fishy to me that C becomes more "singular" when thresholding using smaller alpha. The opposite seems more natural for me. But you must know it better than me, after all it's *your* data.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message ...
> I finish with a note on a small technical trick: for 'y' a vector and 'A' a large sparse matrix, instead of calculate:
>
> A' * y
>
> It is faster to do
>
> (y'*A)'
>

Are you sure that's still true in the most recent MATLAB version? I rewrote A*x, A'*x, etc, a while back, and I think A'*x computes the result without transposing A (or at least my code does that; I can't recall if MATLAB exploits that feature yet, or not).

My SFMULT code in SuiteSparse (on my web page) is what's used inside MATLAB. It can do A'*x without computing A' first. It appears in MATLAB 7.6 and later.

There will be some differences between the two statements, perhaps, but I'd be surprised if it was very significant.

"Tim Davis" <davis@cise.ufl.edu> wrote in message <halhsr$iv2$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message ...
> > I finish with a note on a small technical trick: for 'y' a vector and 'A' a large sparse matrix, instead of calculate:
> >
> > A' * y
> >
> > It is faster to do
> >
> > (y'*A)'
> >
>
> Are you sure that's still true in the most recent MATLAB version? I rewrote A*x, A'*x, etc, a while back, and I think A'*x computes the result without transposing A (or at least my code does that; I can't recall if MATLAB exploits that feature yet, or not).
>
> My SFMULT code in SuiteSparse (on my web page) is what's used inside MATLAB. It can do A'*x without computing A' first. It appears in MATLAB 7.6 and later.
>
> There will be some differences between the two statements, perhaps, but I'd be surprised if it was very significant.

You are absolutely correct Tim. There are a significant speed difference on v7.3 (between two syntax) and very similar speed on v7.9. Thanks for pointing out the improvement.

Bruno

"Hua Wang" <ehwang@163.com> wrote in message <hailt4$evb$1@fred.mathworks.com>...
> > 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, ...)
> >
>
> I compared the above method with the following solutions:
>
> w=cholinc(C,'0');
> x=(w\A)\(w\b);
>
> I found that Bruno's pcg result is the same with x=(A'*(C\A))\d, but it is different from the result by cholinc method. Which one is better?
>

I've not solve this problem yet. And I met the save problem for another application which has smaller vcm matrix than the above. The condition number is very big (around 10^5). Also the result from least-sqaures inversion (e.g. Bruno's suggestion for pcg method) and my "x=(w\A)\(w\b) are quite different. The least-squares results seems making no sense because I expect the first parameter be a positive value (around 5-10), but it give an negative value.

I have put my data at http://homepages.see.leeds.ac.uk/~earhw/down/ . Could anybody test it and give me some suggestion? Anybody know how to set the preconditional matrix to drop the condition number? Many thanks!

Hua