Thread Subject: cholinc for large sparse covariance matrix

Dear All,

I am using cholesky decomposition for a large variance-covariance matrix (vcm). The objetive is to simulate some noise using a known vocariance matrix. The general solution is:
1. make random noise using randn. e.g. Z=randn(nsets,rows*cols)
2. make cholesky decoposition for the covariance matrix. e.g. V=chol(vcm).
3. The simulated noise is X=Z*V.

For my study, the mathematical function of vcm is: vcm=exp(-d/alpha), where d is the distance between any two pixels in an image, and alpha is wavelength of vcm. Usually, alpla is about 20km, and the pixel size is about 100m. So it is almost impossible for me to store the full covariance matrix for a large image. My solution is set a threshold for vcm, i.e. setting vcm=0 if d>n*alpha. I can set n as 5 or so. Alternatively, I can set vcm=0 if vcm<1e-3 etc. Then there are a lot of zeros, and I can use sparse matrix to store vcm.

My problem is that if I use the above method, cholinc will be used for cholesky decomposition. But it failed if I set droptol as 0, and the error message if "vcm is not positive definite". If I use another drop tolerance, e.g. 1.0e-8, cholinc passed. For small vcm, e.g. 3000*3000, the difference between V'*V and the original vcm is smaller than 0.01. But the large matrix (9000*9000), the difference between V'*V and the original vcm is very large. I guess it is a problem of numerical errors. For my case study, the size of vcm matrix will be 20,000 by 20,000.

Could anybody give me some suggestions about setting drop tolerance and simulation? Thanks in advance!

Cheers,
Hua

"Hua Wang" <ehwang@163.com> wrote in message <hacl66$2cm$1@fred.mathworks.com>...
> Dear All,
>
> I am using cholesky decomposition for a large variance-covariance matrix (vcm). The objetive is to simulate some noise using a known vocariance matrix. The general solution is:
> 1. make random noise using randn. e.g. Z=randn(nsets,rows*cols)
> 2. make cholesky decoposition for the covariance matrix. e.g. V=chol(vcm).
> 3. The simulated noise is X=Z*V.
>
> For my study, the mathematical function of vcm is: vcm=exp(-d/alpha), where d is the distance between any two pixels in an image, and alpha is wavelength of vcm. Usually, alpla is about 20km, and the pixel size is about 100m. So it is almost impossible for me to store the full covariance matrix for a large image. My solution is set a threshold for vcm, i.e. setting vcm=0 if d>n*alpha. I can set n as 5 or so. Alternatively, I can set vcm=0 if vcm<1e-3 etc. Then there are a lot of zeros, and I can use sparse matrix to store vcm.
>
> My problem is that if I use the above method, cholinc will be used for cholesky decomposition. But it failed if I set droptol as 0, and the error message if "vcm is not positive definite". If I use another drop tolerance, e.g. 1.0e-8, cholinc passed. For small vcm, e.g. 3000*3000, the difference between V'*V and the original vcm is smaller than 0.01. But the large matrix (9000*9000), the difference between V'*V and the original vcm is very large. I guess it is a problem of numerical errors. For my case study, the size of vcm matrix will be 20,000 by 20,000.
>
> Could anybody give me some suggestions about setting drop tolerance and simulation? Thanks in advance!
>
> Cheers,
> Hua

Compute the smallest eigenvalue of your covariance matrix C (it will be negative due to numerical error), using EIGS function, then add

C = C + 2*abs(lambda)*speye(size(C))

then do the Cholesky on the modified matrix.

Bruno

>
> Compute the smallest eigenvalue of your covariance matrix C (it will be negative due to numerical error), using EIGS function, then add
>
> C = C + 2*abs(lambda)*speye(size(C))
>
> then do the Cholesky on the modified matrix.
>
> Bruno

Hi Bruno,

Thank you so much again for your prompt help!

How can I get the smallest eigenvalue? When I use the function d=eigs(C,5,'sm') and d=eigs(C,1,'sm'), they give different smallest value. For the first, it gives five eigenvaues, and the smallest is -0.2*10^(-3). For the second, it gives -7*10^(-5).

What does lambda mean in your formula? Is it the smallest eigenvalue? Thanks!

Hua

"Hua Wang" <ehwang@163.com> wrote in message <haco7n$1m6$1@fred.mathworks.com>...
> >
> > Compute the smallest eigenvalue of your covariance matrix C (it will be negative due to numerical error), using EIGS function, then add
> >
> > C = C + 2*abs(lambda)*speye(size(C))
> >
> > then do the Cholesky on the modified matrix.
> >
> > Bruno
>
> Hi Bruno,
>
> Thank you so much again for your prompt help!
>
> How can I get the smallest eigenvalue? When I use the function d=eigs(C,5,'sm') and d=eigs(C,1,'sm'), they give different smallest value. For the first, it gives five eigenvaues, and the smallest is -0.2*10^(-3). For the second, it gives -7*10^(-5).

Try with eigs(C,1'sr') (left-most eigenvalue). Don't worry too much about the difference in values, which is probably the magnitude of the error of the numerical method due to round-off (also depending also how the Arnoldi iteration strated from a random vector).

>
> What does lambda mean in your formula? Is it the smallest eigenvalue? Thanks!

Yes, *I* would pick the left-most.

Bruno

> > >
> > > C = C + 2*abs(lambda)*speye(size(C))

Great!
The difference between V'*V and C is only about 10^(-15) now even if without using droptol for cholinc, e.g. V=cholinc(C,0). I have a couple of more questions about this.

1. cholinc gives a warning message "Warning: Incomplete Cholesky with a zero or unspecified drop tolerance calls complete sparse Cholesky". I guess it does not matter because I set 0 and cholinc uses complete sparse choleskty instead. Do you think so? Or any suggestions about droptol?

2. I use a function to calculate the covariance matrix. I wonder whether I can add "C = C + 2*abs(lambda)*speye(size(C))" in that function. Then, I will use the modified C wherever covariance matrix is used. Is it correct?

Many thanks!

Hua

>
> 2. I use a function to calculate the covariance matrix. I wonder whether I can add "C = C + 2*abs(lambda)*speye(size(C))" in that function. Then, I will use the modified C wherever covariance matrix is used. Is it correct?

Yes. The difference should be the order of numerical errors, and there shouldn't be any side issue.

% This code is probably more robust to warranty positiveness:
s1=eigs(C,1,'sm');
s2=eigs(C,1,'sa');
lambda = -min([s1 s2]);
lambda = max(lambda, eps(max(diag(C)))*length(C));

C = C+2*lambda*speye(size(C))

Bruno

> Yes. The difference should be the order of numerical errors, and there shouldn't be any side issue.
>
> % This code is probably more robust to warranty positiveness:
> s1=eigs(C,1,'sm');
> s2=eigs(C,1,'sa');
> lambda = -min([s1 s2]);
> lambda = max(lambda, eps(max(diag(C)))*length(C));
>
> C = C+2*lambda*speye(size(C))
>

Thank you so much!

Cheers,
Hua

>
> % This code is probably more robust to warranty positiveness:
> s1=eigs(C,1,'sm');
> s2=eigs(C,1,'sa');
> lambda = -min([s1 s2]);
> lambda = max(lambda, eps(max(diag(C)))*length(C));
>
> C = C+2*lambda*speye(size(C))
>
> Bruno

Hi Bruno,

I find it is quite slow to find the eigs when C is a very big matrix (9000 by 9000). It takes 10 minutes to finish the above two EIGS function. Is there a faster way to to use this function?

Thanks!
Hua

"Hua Wang" <ehwang@163.com> wrote in message <had2c2$kiu$1@fred.mathworks.com>...
> >
> > % This code is probably more robust to warranty positiveness:
> > s1=eigs(C,1,'sm');
> > s2=eigs(C,1,'sa');
> > lambda = -min([s1 s2]);
> > lambda = max(lambda, eps(max(diag(C)))*length(C));
> >
> > C = C+2*lambda*speye(size(C))
> >
> > Bruno
>
> Hi Bruno,
>
> I find it is quite slow to find the eigs when C is a very big matrix (9000 by 9000). It takes 10 minutes to finish the above two EIGS function. Is there a faster way to to use this function?
>

Try this if it helps you:

% Hilbert matrix
m=1:100;
C = 1./(bsxfun(@plus,m',m)-1);
lambda = smeig(C); % function defined below
lambda = max(lambda, eps(max(diag(C)))*length(C));

C = C+2*lambda*speye(size(C));
L = chol(C)

end

% Estimate the smallest eigenvalue by Rayleigh-Ritz iteration
function lambda = smeig(C)

wn = warning('off','MATLAB:nearlySingularMatrix');
n = size(C,1);
x = randn(n,1);
maxit = 10;
lambda = Inf;
for k=1:maxit
    z = C\x;
    nz = norm(z);
    if nz==0
        lambda = 0;
        break
    end
    lambda = min(lambda, 1/nz);
    x = lambda*z;
end

warning(wn.state,'MATLAB:nearlySingularMatrix');

end % smeig

% Bruno

> lambda = max(lambda, eps(max(diag(C)))*length(C));

Thanks Bruno!

What's the meaning of eps(max(diag(C)))? Can I just use eps(double)? When I add this to my function, there is an error that "Subscription must be positive or logical values". It is quite strange. There must be something wrong before this step. I met the same problem before, but I can not remember how I fixed it.

Also, if I move this statement outside of the function. It give an error message "error using eps, class must be 'single' or 'double'".

btw, I've sent you an email about Costantini's thesis.

Cheers,
Hua

"Hua Wang" <ehwang@163.com> wrote in message <had5j6$g31$1@fred.mathworks.com>...
> > lambda = max(lambda, eps(max(diag(C)))*length(C));
>
> Thanks Bruno!
>
> What's the meaning of eps(max(diag(C)))?

This returns the smallest number such that when adding to C the largest diagonal elements must change. The purpose is to avoid adding a diagonal matrix that does nothing to the original one. The number is conservative (large), however it is small enough compared to the norm of the covariance matrix.

>
> btw, I've sent you an email about Costantini's thesis.
>

Thanks.

Bruno

> >
> > What's the meaning of eps(max(diag(C)))?
>
> This returns the smallest number such that when adding to C the largest diagonal elements must change. The purpose is to avoid adding a diagonal matrix that does nothing to the original one. The number is conservative (large), however it is small enough compared to the norm of the covariance matrix.
>
There is an error when I use this statement. The error is "eps: class must be 'single' or 'double'". When I run max(diag(C)), the anser is "(1,1), 1". You know, the maximum covariance coefficient is 1. I am not sure whether it is due to sparse matrix I used? It seems correct if C is not a sparse matrix. If I change the sparse diag(C) as full matrix, eps(full(diag(C))) is the same with eps('double') then. Shall I use eps('double') directly? Hope to hear from you!

Hua

> lambda = smeig(C); % function defined below
> lambda = max(lambda, eps(max(diag(C)))*length(C));
>
> C = C+2*lambda*speye(size(C));
> L = chol(C)
>

I tested the above codes and your SMEIG function in my program. It is indeed must faster than EIGS. For the test data, it takes 140 s but EIGS takes 392 s. But they give different smallest eigs.
For s1 = eigs(C,1, 'sm'), s1=0.056
For s2 = eigs(C,1, 'sa'), s2=0
and for your smeig, lambda = 0.0018.
Do you think the above difference is important?

Shall I use lambda = - smeig(C)? Since lambda= - min[s1 s2].

It seems that chol also works for the sparse matrix when I use matlab-2009a. And the function w=chol(C) is much faster than w=cholinc(C,'0'). What is the difference with them? The outputs are not a very different. Thanks!

Hua