|Re: [eigen] Another LDLt issue|
[ Thread Index |
| More lists.tuxfamily.org/eigen Archives
- To: eigen@xxxxxxxxxxxxxxxxxxx
- Subject: Re: [eigen] Another LDLt issue
- From: Benoit Jacob <jacob.benoit.1@xxxxxxxxx>
- Date: Mon, 30 Mar 2009 15:21:36 +0200
- Dkim-signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:mime-version:received:in-reply-to:references :date:message-id:subject:from:to:content-type :content-transfer-encoding; bh=XDJ6N35lB+3NSJPU/2Rd0kTOB8RiPE3fBkyS1Wh4OAs=; b=fXNhW1Feu+E3kNPoUcZFH3cAMIQSU0+21p6Y8jyFmyroTySGkgKfZGegB1snqNR/jp X6fFgcZITz6cEeHIlEIAfX2F1uVOQ/3EhWOnEYLKVRj+1vVf2yy1Ha41BKhtZliZn7+r OLJFyc8VeNGJV895UiJPjDcQHi00IC9rIi71k=
- Domainkey-signature: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:in-reply-to:references:date:message-id:subject:from:to :content-type:content-transfer-encoding; b=Xt9phDG6NOPmpQ70hURhptq61keZySZcxjoCpR+VAEW6luyuJvM6ZaWQjH/M6Aiia2 yG5GfHXMMlRiRcRjyVnCr59VcjlFwfENWZtNPxK7lD/4PEOCq4XVeRravqRgfapzHozo Ms+8Zspp05cMOlRXZ8L1+Hcp/ZHwcttnr9C3k=
Oops. In fact, invertibility won't help, and I'm not sure at all that
LDLt (even pivoting) exists for all self-adjoint matrices. I mean, I
"this factorization can be used for any square, symmetrical matrix"
but I don't see anything there backing this assertion, it seems to me
that it's wrong, and at least the formulas that they give fail for the
>So a generalization of LDLt handling non-positive matrices would have
>to assume invertibility. It would work for self-adjoint invertible
>matrices. Why not. Are you still interested in that?
>The next question would then be: implement as separate dec or as a
>variant of LDLt controlled by a template parameter?
2009/3/30 Benoit Jacob <jacob.benoit.1@xxxxxxxxx>:
> I had a closer look.
> What happens is that the "positive" assumption allows to restrict the
> pivot lookup to the diagonal. So we get full pivoting at the price of
> partial pivoting.
> What allows us to do that is that the norm of a positive matrix is
> bounded by a multiple of the norm of its diagonal. In particular, if
> the diagonal is zero, then the whole matrix is zero.
> This is really characteristic of positive matrices, it doesn't work
> for general self-adjoint matrices, e.g.
> 0 1
> 1 0
> So a generalization of LDLt handling non-positive matrices would have
> to assume invertibility. It would work for self-adjoint invertible
> matrices. Why not. Are you still interested in that?
> The next question would then be: implement as separate dec or as a
> variant of LDLt controlled by a template parameter?
> By the way, Keir, at line 162:
> m_matrix.block(j, j, size-j, size-j).fill(0); // Zero unreliable data.
> Unless this is backed by experiment, I kind of disagree with this
> line, I don't see how zeroing this noise can be useful, and I can
> think of disadvantages of zeroing: it makes the reconstruction of the
> original matrix more inaccurate, not to mention that it generates more
> 2009/3/30 Benoit Jacob <jacob.benoit.1@xxxxxxxxx>:
>> This is interesting. Indeed it seems that this works for all
>> self-adjoint matrices. I'll try to have a look at it.
>> 2009/3/29 <w.h.greene@xxxxxxxxx>:
>>> I just finished reading a long recent thread on some issues with LDLt.
>>> In some numerical experiments I was doing today, I've come across a
>>> related issue. It appears that the current implementation of LDLt requires
>>> that the matrix be positive definite. This is not strictly necessary for
>>> the LDLt factorization to succeed. It is necessary only that the diagonal
>>> term not be zero during factorization. A negative diagonal term is not,
>>> by itself, a problem.
>>> From reading the previous posts, it appears that a main reason for the
>>> LDLt implementation compared with LLt was to avoid the performance
>>> penalty of a square root. In fact, I think the main benefit of LDLt is that
>>> it can factor both negative- and positive-definite symmetric matrices.
>>> (As an aside, these routinely occur in dynamic analysis of mechanical
>>> I suggest the test in the LDLt factorization routine be changed to
>>> check that abs(diag_term) > eps.
>>> Bill Greene