Re: [eigen] Instability in LLT and LDLT methods.

[Keir Mierle <mierle@xxxxxxxxx>]

Hi,

I am working on a pivoting variant of the LDLt, but it doesn't seem
that the algorithm matches the formula on the wikipedia page

http://en.wikipedia.org/wiki/Cholesky_decomposition (see the section
avoiding square roots).

In particular, it doesn't appear the Dj is multiplied into the
product. Where does this algorithm come from?

Keir

On Wed, Jan 28, 2009 at 8:45 AM, Gael Guennebaud
<gael.guennebaud@xxxxxxxxx> wrote:
> see this thread:
> http://forum.kde.org/solved-lu-only-half-as-fast-as-lapack-t-28265-2.html#pid38539
> there is a link to source code under BSD license.
>
> that's all I have :(
>
> good luck !
>
> On Wed, Jan 28, 2009 at 5:34 PM, Keir Mierle <mierle@xxxxxxxxx> wrote:
>> I'm working on this now (pivoting in LDLt). Do you have a reference
>> for the algorithm used? It's pretty simple but a reference would help.
>>
>> Keir
>>
>> On Tue, Jan 27, 2009 at 11:42 PM, Gael Guennebaud
>> <gael.guennebaud@xxxxxxxxx> wrote:
>>> yes, indeed with LDLt we could do full pivoting and be as stable as LU
>>> for selfadjoint matrices while being faster. My initial motivation
>>> with LDLt, however, was its performance because it avoids the square
>>> roots... On the other hand, I remember my benchmark was not really in
>>> favor of the current LDLt,. I have to check again, but if so, then
>>> there is no reason not to keep the current LDLt version which could be
>>> replaced by a more complex one with full pivoting.
>>>
>>> On Wed, Jan 28, 2009 at 12:43 AM, Keir Mierle <mierle@xxxxxxxxx> wrote:
>>>> Probably it's better to do full pivoting. Apparently cholesky is
>>>> stable for semidefinite matrices when full pivoting is used:
>>>>
>>>> http://eprints.ma.man.ac.uk/1101/01/covered/MIMS_ep2008_56.pdf
>>>>
>>>> Keir
>>>>
>>>> On Tue, Jan 27, 2009 at 3:05 PM, Gael Guennebaud
>>>> <gael.guennebaud@xxxxxxxxx> wrote:
>>>>> Hi,
>>>>>
>>>>> yes, it seems the test to check whether the matrix is positive
>>>>> definite was too strict. I changed the absolute tolerance a bit, but
>>>>> we still need something better. Basically, in Cholesky we compute at
>>>>> each iteration 1/sqrt(x), and so x must be >0 with some epsilon...
>>>>>
>>>>> Gael.
>>>>>
>>>>> On Tue, Jan 27, 2009 at 7:35 PM, Keir Mierle <mierle@xxxxxxxxx> wrote:
>>>>>> Here is a testcase that fails with LLT and LDLT but works fine with
>>>>>> all of LU, SVD, and QR solving. Depends on my previous patch for QR
>>>>>> solver (or comment out the qr().solve line).
>>>>>>
>>>>>> Keir
>>>>>>
>>>>>
>>>>>
>>>>>
>>>>
>>>>
>>>>
>>>
>>>
>>>
>>
>>
>>
>
>
>