Problems with inversion using chptrf/chptri

[thombos]

Hello,

i have a hermitian matrix in packed storage that I'd like to invert. I have implemented some code using chptrf/chptri, which seems to work fine for a number of test matrices (compared with matlab). In my application however the inverse obtained with these routines is often wrong (very wrong). I can unpack the matrix in to full storage and use cgesvd, no problem. Is there a restriction, which hermitian matrices this routine can invert and which not ? I have not done the work to check whether the factorization is bad or whether the inverse is bad.

Thanks, Thomas

Edit: I used a matrix from my application and factored it in matlab using ldl(). The inverse I got using that is slightly different from my svd based pseudo-inverse, but not way off, like the result from the lapack routines I get. I compared the Ds from matlab ldl() and chptrf and they look similar (well they are permuted differently), but the chptrf produced distinctively more negative values.

Edit 2: The rcond value estimated with chpcon for a problem that "fails" is 7.79277e-11

Edit 3: I added noise to improve the conditioning and now it works reasonably well, so this might be my fault after all. If I check the conditioning after chptrf, is there anything I can do to improve the conditioning/regularize the inverse ? I'm not very familiar with LDL'/UDU' decompositions, and its also baffling that the matlab ldl() routine seems to deal with this a bit better.

[Julien Langou]

Hello Thomas,

* I am surprised that Matlab ldl() and LAPACK CHPTRF/CHPTRI gives different results. I would have assume Matlab ldl() is using CHETRF/CHETRI (which is the unpacked version of CHPTRF/CHPTRI).

* It seems like your problem is ill-conditioned. (See the RCOND value you are giving.) This means that two different (correct) backward stable algorithms can give two answers that are very far from each other while however correct in the backward stable sense.

* Adding noise removes the small singular values and so improves the condition number and so should solve your problem to some extent, yes.

Julien.

[thombos]

Thanks for your reply. For LDL inversion, is there a way of regularizing similar to say truncated SVD ? I'm not sure how the factorization works, and whether "tampering" with D would achieve any sort of regularizing effect...

[Julien Langou]

This can turn out to be an excellent idea. I am not aware of anything related to this. I will need to check the bibliography.

L is supposed/assumed to be well conditionned independently of the condition number of A. (This is not always the case but that's a standard assumption.) Therefore all the ill-conditioning of A is in D (if it's not in L). Therefore, yes, you can "filter" the small singular values of A by using its LDLT factorization. This needs some checks but that should work almost all the time. As long as L is well conditioned.

Julien.