4
$\begingroup$

I would like to compute the decomposition of a real symmetric positive definite matrix $\mathbf{A} = \mathbf{UDU}^\top$.

LINPACK seems to have it as DSIFA, but I cannot find an equivalent routine in LAPACK. It also doesn't appear to implemented in Eigen.

What other packages and routines support this decomposition?

Background:

In my field, the most common of Kalman Filter implemented in the "real world" is via UDU as evidenced in multiple books on the subject. I understand that there is an equivalent $\mathbf{LDL}^\top$ form, but I am trying to stick with the convention in my field.

Improve this question
$\endgroup$
2
  • 2
    $\begingroup$ The algorithm for $UDU^T$ factorization is essentially equivalent to an $LDL^T$ factorization, but executed "up" the diagonal of the original matrix rather than down it. I haven't checked the details, but it should therefore be possible to reverse the order of the entries in the original matrix, i.e., $\mathcal{R}(A)(i,j) := A((n-1)-i,(n-1)-j)$, run a standard $LDL^T$ factorization, and then reverse the memory of each of the resulting factors. $\endgroup$
    – Jack Poulson
    May 14, 2013 at 3:50
  • 1
    $\begingroup$ Now, if there were only some way to hijack the "twisted factorizations" implemented in LAPACK... BTW, DSIFA doesn't do a genuine diagonal factorization, as the $\mathbf D$ factor it returns is in fact block-diagonal, at least if the input matrix is symmetric-indefinite (Bunch-Parlett). $\endgroup$
    – J. M.
    May 14, 2013 at 3:53

1 Answer 1

Reset to default
4
$\begingroup$

Eigen's LDLT class actually performs a U^T.D.U factorization if the input matrix is column-major with the symmetric entries stored in the upper triangular part:

MatrixXd A;
// fill at least the upper triangular part of A
LDLT<MatrixXd,Upper> udu(A); // Only the upper part of A is read to form U
udu.matrixU(); // returns an expression of the triangular matrix U (column major)
udu.matrixL(); // returns an expression of the triangular matrix U^T (row major)
udu.vectorD(); // expression of the diagonal coefficients of the matrix D as a vector
udu.vectorD().asDiagonal(); // an expression of the diagonal matrix D
Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ An answer from The Man himself. Welcome to Scicomp, Gael! $\endgroup$
    – Damien
    May 17, 2013 at 11:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.