Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices |

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices*From*: Ben Goodrich <bgokgm@xxxxxxxxxxxxxx>*Date*: Sat, 27 Feb 2010 19:04:38 -0500*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=googlemail.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=jXhG9+ndeNIyrM59oeqcJdgGBOSlMmzDB6xa5YMjh8g=; b=cIjRkxznjyXilS/kHYouX1cx0QhywY4OJ6NhOHOQsFXrE3OYUMQuLjiK2Ztvv21evy 7VeJJp57oeeycqas+yLJ6oBSBhoHPdtLB0ZA3fqxxsGQX84ycxmnR4PZBgMZcTTLNvW3 9A6SHRw6hAnLVyJbYnVTP08MvnmlqKjBTBYSk=*Domainkey-signature*: a=rsa-sha1; c=nofws; d=googlemail.com; s=gamma; h=mime-version:in-reply-to:references:date:message-id:subject:from:to :content-type:content-transfer-encoding; b=fTA+VbKXWLfa4tTjAiNnXdKsK8aU4p6R8r21U2gHRG0L4LxZcIsnBivXtG1aYSp5be tFE/6RAdgEXF3AqJ0jWrw97pwo58eU8THaVxSL9dbdEZgFCFXVZ3vjh8DbYGSYtbef0O 833ClUZt6CZQdWVPwSNdoBv/Zj92awima26dQ=

Hi Benoit, On Sat, Feb 27, 2010 at 5:27 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote: > 2010/2/27 Ben Goodrich <bgokgm@xxxxxxxxxxxxxx>: >> On Wed, Feb 24, 2010 at 12:45 PM, Ben Goodrich <bgokgm@xxxxxxxxxxxxxx> wrote: >>> This is a great example, thank you. In my opinion, it would be good to >>> have a way in eigen to define D uniquely by utilizing the unit lower >>> triangular restriction on L. > > Why do you care about D? What is the use case? Definitely not a standard use case :) Basically, it is not solving AX = B systems, but rather numerical optimization / statistics where the A matrix is a function of a parameter vector and the objective function being maximized, f(theta), is defined with respect to the D in "the" LDLt decomposition of A(theta). [Or at least one could define the objective function that way. I have a paper under review that operationalizes the objective function on the eigenvalues of A(theta), and I am trying to see if the LDLt decomposition would work better or worse for the revision. It seems to have some pros and cons in my testing so far. For what it is worth, although I am interested in this for my own research, the paper by Benson and Vanderbei I referred to earlier on this thread http://www.pages.drexel.edu/~hvb22/dimacsppr3.pdf talks in more general terms about how D is useful for operationalizing the constraint in semidefinite programming problems (mine is just one specific instance of that) because the diagonal cells of D are convex functions of the cells of A.] One big problem was the situation I outlined earlier in test_noise.cc VectorXd D(5); D << // feed data MatrixXd L(5,5); L << // feed data MatrixXd A = L * D.asDiagonal() * L.transpose(); A_factored = A.ldlt(); The numerical noise causes ldlt() to pivot differently and thus A_factored.vectorD() is very different from the D that "actually" generated A. However if I do A.diagonal().setOnes(); A_factored = A.ldlt(); then A_factored.vectorD() is essentially the same as the D that "actually" generated A (which in symbolic math would have had 1 in all its diagonal cells). Thus, I am worried that numerical noise can cause (more) discontinuities in my objective function due to different ways of pivoting during the optimization and would like to have a "consistent" (maybe textbook was a poor choice of word) way to define D so that f(theta) is smooth and always approximately the same as f(theta + epsilon). > The vectorD and matrixL methods are zero-cost accessors that just > return stuff that has previously been computed during the > decomposition. So, if you wanted to implement something like that, > that would rather be > - either as an option to the decomposition itself > - or as new separate methods. > > I still don't understand what is "textbook form" and why it's useful :/ Right. And I am basically just trying to figure out the fastest way to get a consistent D that has decent numerical precision even when A(theta) is nearly singular (but not indefinite). If I am understanding correctly, I need to unpivot in basically the same way that reconstructMatrix() does. But unpivoting is going to be computationally expensive when it has to be done repeatedly in an optimization context, so I was hoping someone had a short-cut relative to the steps that first came to my mind :) . So, it sounds like I just need to try to implement it every way that I can think of and see how well / fast they work out. Thanks, Ben

**Follow-Ups**:**Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices***From:*Ben Goodrich

**Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices***From:*Benoit Jacob

**Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices***From:*Jitse Niesen

**References**:**[eigen] [patch] LDLt decomposition with rank-deficient matrices***From:*Ben Goodrich

**Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices***From:*Benoit Jacob

**Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices***From:*Ben Goodrich

**Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices***From:*Gael Guennebaud

**Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices***From:*Jitse Niesen

**Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices***From:*Ben Goodrich

**Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices***From:*Jitse Niesen

**Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices***From:*Ben Goodrich

**Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices***From:*Ben Goodrich

**Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices***From:*Benoit Jacob

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] Eigen3 ->Eigen2 performance regression: patch.** - Next by Date:
**Re: [eigen] merged the eigen-strides fork** - Previous by thread:
**Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices** - Next by thread:
**Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices**

Mail converted by MHonArc 2.6.19+ | http://listengine.tuxfamily.org/ |