| Re: [eigen] [RFC] Full pivoting LDLT |
[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen Archives
]
- To: eigen@xxxxxxxxxxxxxxxxxxx
- Subject: Re: [eigen] [RFC] Full pivoting LDLT
- From: Keir Mierle <mierle@xxxxxxxxx>
- Date: Wed, 4 Feb 2009 10:19:16 -0800
- 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=+LPcyc5K9G2JvRs+7Vk6FVxKYfv5mBMcKzNQW9uwPoI=; b=RMv5OpgUxs0ghE6FdTQEgMDiQcMoOmgdiNMEEGDUTcF9uzLecdT6mPOiqOSXMUKmaD SpyTmtWTFua20ZcBniHtdOy6MOoH218WA/E4skIrW59+0nz8JlB93/egglaVXVXMUWNJ GTPrIcYElHxgUcieYR4Xr7TpKcpOGJHpYxUUs=
- 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=eWiYCEnYaRUa+RTsLe+fw3831/J7K8db9T/6HKta0xFSY+K0yoVp3ddfYz6vAqLT8d HiZK1bEtJoSVyhWgYUBTgO2C+otdXYSOWW4SB4LAeD+tSI74YVAD6r+pASOjRP1aGc3W yAWlRnvDUEJ7Hy2DRGUug7wMZj5i0u0lqWdnY=
On Wed, Feb 4, 2009 at 10:10 AM, Jitse Niesen <jitse@xxxxxxxxxxxxxxxxx> wrote:
> On Tue, 3 Feb 2009, Keir Mierle wrote:
>
>> http://codereview.appspot.com/14042/diff/1/4
>>
>> There are a couple of issues.
>>
>> 1. Pivoting is not optional. Any objections?
>> 2. If the matrix is only semidefinite, then matrixL() gives wrong results
>> (1's in diagonal where they shouldn't be).
>>
>> Minus objections, I'll commit this.
>
> I don't understand issue 2. I thought the L in the LDL^T decomposition has
> per definition 1s on the diagonal?
>
> It looks like you assume that the matrix is positive semi-definite. However,
> I believe that the LDL^T decomposition is also useful for symmetric
> indefinite matrices.
If A is symmetric indefinite, then a very different algorithm is
required because D is block-diagonal as 1x1 and 2x2 blocks. I did not
implement this because (for my use cases) A is always (semi)definite,
and it would substantially slow the code.
> Solving Ax = b with A symmetric is twice as fast with
> LDL^T than with LU (if I remember correctly).
>
> Why do you store P in two places: m_p and m_transpositions? Is using
> m_transpositions faster than m_p?
I left that in for the non-solve-in-place case. If you are solving in
place, you need to apply the transpositions via m_transpositions. If
you are solving into a new matrix, then m_p should be used instead
(see LU.h's use of m_p in solve()). I haven't implemented solve()
using m_p yet though (feel free!). m_p is a direct mapping from one
permutation to another; transpositions is the series of swaps to apply
to get to the other permutation.
> I was a bit surprised that you use both the triangles in m_matrix during the
> computation, since one half is just a rescaling of the other half. So I
> tried using only the lower triangle, but my timings (using g++ 4.1) are
> inconclusive on whether this is actually faster. I attach the diff for your
> entertainment, but I don't have any good reason for it.
Hmm, that is interesting that it's slower. I'll look at this tonight.
Keir
>
> Cheers,
> Jitse