5
$\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
  • 4
    $\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$ 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

2 Answers 2

4
$\begingroup$

Short Answer: LAPACK's dsytf2 (for symmetric full) and dsptrf (for symmetric packed, which is the same layout that Bierman uses in its Kalman filter subroutines) actually computes $UDU'$ decomposition as it is used in estimation community.

Longer Version:

Cholesky decomposition routines of LAPACK (such as dpotrf) compute only $LL'$ and $U'U$ (not $UU'$). It is interesting to note until version 3.6 of LAPACK, the documentation of LAPACK wrongfully stated that when the UPLO='U', the Cholesky decomposition returns $UU'$ decomposition. But, that wrong statement was fixed after 3.6. Currently (as of 3.12), LAPACK Cholesky decomposition routines can only return $U'U$ or $LL'$ decomposition based on UPLO argument (same as matlab's chol function).

On the other hand, $LDL$ decomposition routines of LAPACK actually compute $UDU'$ when the UPLO is set to 'U' in the arguments. This fact is also explicitly stated in these functions' documentation.

However, one should also note that LAPACK LDL routines execute much more complicated algorithms as the LDL routines are designed for indefinite matrices. In estimation algorithms, we generally set the row/col to zero when the pivot is zero and continue to decompose the rest of the matrix as it is known that the input (the covariance matrix) is not indefinite. However, LDL routines of LAPACK, perform row/col rooks in such cases. Therefore, one may need to use syconv routine to convert permutated $LDL'$/$UDU'$ matrices into unpermutated forms.

$\endgroup$
7
$\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
  • 1
    $\begingroup$ I could not understand how this answer actually answers the question. The question asks if there is a routine that computes UDU', the answer says eigen can only computes U'DU (which is a completely different decomposition than UDU'), and this answer was chosen as accepted. $\endgroup$
    – tantuni
    Jan 13 at 10:59

Your Answer

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

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