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

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

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] Instability in LLT and LDLT methods.*From*: Benoit Jacob <jacob.benoit.1@xxxxxxxxx>*Date*: Wed, 28 Jan 2009 17:00:11 +0100*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; bh=54BE52bOV2yMt+AefoQisKhvs497bhOpPdXKT1G02WY=; b=sbYwkjrDPdNo+YHdDKrbAiGd5/pT3KZNNCBkBrL5GOBN8utu4Ow3EgqH78Ov/S+Lm6 pcYZ8g21UcAi+E6xHTMu9mriTtP8z7FPdA1epGjJ//wyYqFD3Cpzx4KwiliLB0RaXlUO /GKxeJlI9HcFNxQCPq6WvHX6efg897H15agRY=*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; b=MBhIZ4UzfUjFF7fZMPLDVbgzejr08XxaFVs6RGYayV0zooj0frOU59cgkirQDpvluk DYIBJglxTkVKGtmfS86cEEWHTSoNa+Rbf5jS4bmdma3u4+awe6pjuH2dpu+oYcEK7JVS vMUZ9JXzlc6YTcV3axXImSQYwgGFyZxv59zdY=

OK so now I have a patch (attached) making LU use only partial pivoting. Results: * for an invertible matrix, partial pivoting seems to be OK. So partial pivoting can be considered for implementing inverse(), where we assume the matrix to be invertible. * for a non-invertible matrix, when one wants to compute the rank or things that depend on it, such as a basis of the kernel, partial pivoting is not enough. E.g. apply my patch, and try the test_lu unit test with 100 repetitions, it fails. So yes it would be potentially interesting to add back the template parameter "int Pivoting" to class LU and let inverse() use a partial pivoting LU. All rank-based methods (such as kernel() and image()) should only be available in the full pivoting version. Cheers, Benoit

Index: Eigen/src/LU/LU.h =================================================================== --- Eigen/src/LU/LU.h (révision 917831) +++ Eigen/src/LU/LU.h (copie de travail) @@ -340,34 +340,22 @@ IntRowVectorType cols_transpositions(matrix.cols()); int number_of_transpositions = 0; - RealScalar biggest = RealScalar(0); + RealScalar biggest = m_lu.cwise().abs().maxCoeff(); m_rank = size; for(int k = 0; k < size; ++k) { int row_of_biggest_in_corner, col_of_biggest_in_corner; RealScalar biggest_in_corner; - biggest_in_corner = m_lu.corner(Eigen::BottomRight, rows-k, cols-k) + biggest_in_corner = m_lu.block(k,k, rows-k, 1) .cwise().abs() - .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner); + .maxCoeff(&row_of_biggest_in_corner); row_of_biggest_in_corner += k; - col_of_biggest_in_corner += k; - if(k==0) biggest = biggest_in_corner; + col_of_biggest_in_corner = k; - // if the corner is negligible, then we have less than full rank, and we can finish early - if(ei_isMuchSmallerThan(biggest_in_corner, biggest)) - { - m_rank = k; - for(int i = k; i < size; i++) - { - rows_transpositions.coeffRef(i) = i; - cols_transpositions.coeffRef(i) = i; - } - break; - } - rows_transpositions.coeffRef(k) = row_of_biggest_in_corner; cols_transpositions.coeffRef(k) = col_of_biggest_in_corner; + if(k != row_of_biggest_in_corner) { m_lu.row(k).swap(m_lu.row(row_of_biggest_in_corner)); ++number_of_transpositions; @@ -376,6 +364,13 @@ m_lu.col(k).swap(m_lu.col(col_of_biggest_in_corner)); ++number_of_transpositions; } + + // if the corner is negligible, then we have less than full rank, and we can finish early + if(ei_isMuchSmallerThan(biggest_in_corner, biggest)) + { + --m_rank; continue; + } + if(k<rows-1) m_lu.col(k).end(rows-k-1) /= m_lu.coeff(k,k); if(k<size-1)

**References**:**[eigen] Instability in LLT and LDLT methods.***From:*Keir Mierle

**Re: [eigen] Instability in LLT and LDLT methods.***From:*Gael Guennebaud

**Re: [eigen] Instability in LLT and LDLT methods.***From:*Keir Mierle

**Re: [eigen] Instability in LLT and LDLT methods.***From:*Jitse Niesen

**Re: [eigen] Instability in LLT and LDLT methods.***From:*Benoit Jacob

**Re: [eigen] Instability in LLT and LDLT methods.***From:*Jitse Niesen

**Re: [eigen] Instability in LLT and LDLT methods.***From:*Benoit Jacob

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] Instability in LLT and LDLT methods.** - Next by Date:
**Re: [eigen] Instability in LLT and LDLT methods.** - Previous by thread:
**Re: [eigen] Instability in LLT and LDLT methods.** - Next by thread:
**Re: [eigen] Instability in LLT and LDLT methods.**

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