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.

$\endgroup$
  • 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 '13 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 '13 at 3:53
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
$\endgroup$
  • 1
    $\begingroup$ An answer from The Man himself. Welcome to Scicomp, Gael! $\endgroup$ – Damien May 17 '13 at 11:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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