Skip to content

Confusing documentation for ldlt #725

Description

@stephen-huan

However, since pivoting is on by default, the factorization is internally
represented as `A == P'*L*D*L'*P` with a permutation matrix `P`;
using just `L` without accounting for `P` will give incorrect answers.
To include the effects of permutation, it is typically preferable to extract
"combined" factors like `PtL = F.PtL` (the equivalent of
`P'*L`) and `LtP = F.UP` (the equivalent of `L'*P`).
The complete list of supported factors is `:L, :PtL, :D, :UP, :U, :LD, :DU, :PtLD, :DUP`.
The permutation vector is available as `F.p`, defined such that `L*D*L' == A[p, p]`,
The `LD` component can be materialized as a sparse matrix using `sparse(F.LD)`,
Other components cannot be materialized directly, but can be reconstructed from
`sparse(F.LD)` and `F.p` if needed.

Reading the documentation and seeing F.PtL = P'*L, F.UP = L'*P, etc. naturally seems to imply that F.LD = L * D. It is not; it is L (which with its unit diagonal replaced by D as opposed to L * D which is each row of L scaled by the corresponding diagonal entry of D (which happens to also have diagonal D). If LD was actually L * D, then L would be unrecoverable from L * D when D has diagonal entry 0 (of course, that row doesn't affect the factorization). Since LD is just L's diagonal replaced with D, it suffices to read D off the diagonal of LD and replace the diagonal with 1's for L.

Maybe this convention is obvious to experts, but it tripped me up. It's common to reduce the L D L' factorization to the Cholesky factorization by L D L' = (L D^{1 / 2}) (L D^{1 / 2})', so it's not entirely unreasonable to multiply L by D.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions