[eigen] Solving positive semidefinite systems |

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

*To*: eigen@xxxxxxxxxxxxxxxxxxx
*Subject*: [eigen] Solving positive semidefinite systems
*From*: Edward Lam <edward@xxxxxxxxxx>
*Date*: Mon, 19 Jun 2017 16:38:59 -0400
*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=sidefx.com; s=google; h=to:from:subject:message-id:date:user-agent:mime-version :content-language:content-transfer-encoding; bh=lDai9WKSVt7f65VyuX9j/elLN5W7p5vfalPsa5nRVhc=; b=q3ogEktn132q2FaX2fg7LsHYRsZ/Ju0bWi+MRBUO2tDQPH6pLWvHq79US1iLXK097T kGBOnna7KlEqJahQcX3yC5vuODyCCmocQ2E3neUac9Q5IFURNTbruXPqUYKw6oUk35vE XciZV7aGtFnaLvXgiMIUOpB87fQyRP0DNoCY+ui2pXEobcmg/WVlUJVh4KSsqbNWTcdC ONXdYTlAb//RXyCnZSbPoOUNmlZ7soG772h/Zn7MHi+M9cqCqz/RayPmgqYFdriKFtfX Jz484wRqSILIalaxUTzK9LBdSxZ5n3VQtW+FmLau35KQsxomO5M+pw2HNWziwOc7zKJR +A1Q==

Hi,
First of all, I'd like to thank you all for the wonderful library!

`In trying to convert some particular old code to use Eigen, I've run into a case
``where Eigen's Cholesky solvers don't have a away to ignore A(k,k) if it's 0 (or
``sufficiently close to 0). Instead, they either silently give me bad results or
``give up altogether (depending on if it's the sparse or dense variants).
`

`Golub and Van Loan's Matrix Computations, chapter 4.2.8 (The semdidefinite case)
``gives a very easy modification to the regular Cholesky decomposition that allows
``it to be used for positive semidefinite matrices. All it does is skip the update
``iteration if A(k,k) is found to be 0. In practice, the implementation that I'm
``used to checks if A(k,k) is less than some epsilon, and if it is, sets the
``current diagonal entry and the rest of the row (or column) to be exactly 0 and
``then skip the rest of the loop. Then during back substitution, it checks whether
``the diagonal value is 0, and if so, then it sets the corresponding component of
``the solution to 0 instead of dividing by it. This trick is similarly used in the
``LU matrix solvers as well with the exception that we just set the 0 diagonal
``entries to 0 during factorization and simply skip the divide by 0's during back
``substitution.
`

`Can such a modification be considered for Eigen's dense/sparse Cholesky/LU
``solvers? I have managed to work around the problem by using other Eigen solvers
``like FullPivLU (for dense matrices) or SparseQR (for sparse matrices). However,
``they're much slower.
`
Thanks,
-Edward