R help - Jul 2008 - spectral decomposition for near-singular pd matrices

Hi,

I have a matrix M, quite a few (< 1/4th) of its eigen values are of
O(10^-16). Analytically I know that M is positive definite, but
numerically of course it is not. Some of the small (O(10^-16)) eigen
values (as obtained from eigen()) are negative. It is the
near-singularity that is causing the numerical errors. I could use
svd(), but then the left ($u) and right ($v) vectors are not
identical, again numerical errors. Is there any function that imposes
pd property while calculating the eigen decomposition.

Thanks,
PK

How large is your matrix?

Are the very small eigenvalues well separated?

If your matrix is not very small and the lower eigenvalues are clustered, this
may be a really hard problem! You may need a special purpose algorithm and/or
higher precision arithmetic.
If your matrix is A and there exists a sequence of matrices A1,A2,...An = A such
that A1 is "good", A2 is a bit worse (and is not very different from
A1), etc., you may be able to compute the eigensystem for A1 and then track it
down to A2, A3,...,An. I have worked on such problems some 15 years ago but I
believe that by now there should be packages doing this (I do not know whether
they exist in R).


--- On Thu, 17/7/08, Prasenjit Kapat <kapatp at gmail.com> wrote:
> From: Prasenjit Kapat <kapatp at gmail.com>
> Subject: [R] spectral decomposition for near-singular pd matrices
> To: r-help at stat.math.ethz.ch
> Received: Thursday, 17 July, 2008, 7:27 AM
> Hi,
> 
> I have a matrix M, quite a few (< 1/4th) of its eigen
> values are of
> O(10^-16). Analytically I know that M is positive definite,
> but
> numerically of course it is not. Some of the small
> (O(10^-16)) eigen
> values (as obtained from eigen()) are negative. It is the
> near-singularity that is causing the numerical errors. I
> could use
> svd(), but then the left ($u) and right ($v) vectors are
> not
> identical, again numerical errors. Is there any function
> that imposes
> pd property while calculating the eigen decomposition.
> 
> Thanks,
> PK
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained,
> reproducible code.

Moshe Olshansky <m_olshansky <at> yahoo.com> writes:
> How large is your matrix?
Right now I am looking at sizes between 30x30 to 150x150, though it will 
increase in future.
> Are the very small eigenvalues well separated?
> 
> If your matrix is not very small and the lower eigenvalues are clustered, 
this may be a really hard problem!

Here are the deciles of the eigenvalues (using eigen()) from a typical random 
generation (30x30):
          Min          10th         20th         30th         40th
-1.132398e-16 -6.132983e-17 3.262002e-18 1.972702e-17 8.429709e-17
         50th         60th          70th         80th      90th     Max
 2.065645e-13  1.624553e-09 4.730935e-06   0.00443353 0.9171549 16.5156

I am guessing they are *not* "well-separated." 
> You may need a special purpose algorithm and/or higher precision
arithmetic.
> If your matrix is A and there exists a sequence of matrices A1,A2,...An = A
such that A1 is "good", A2 is a bit
> worse (and is not very different from A1), etc., you may be able to compute
the eigensystem for A1 and then
> track it down to A2, A3,...,An. I have worked on such problems some 15
years
ago but I believe that by now
> there should be packages doing this (I do not know whether they exist in
R).
I will have to think on the possibility of such a product-decomposition. But 
even if it were possible, it would be an infinite product (just thinking out 
loud).

Thanks again,
PK

PS: Kindly CC me when replying.