PositiveFactorizations is a package for computing a positive definite
matrix decomposition (factorization) from an arbitrary symmetric
input. The motivating application is optimization (Newton or
quasi-Newton methods), in which the canonical search direction `-H\g`

(`H`

being the Hessian and `g`

the gradient) may not be a descent
direction if `H`

is not positive definite. This package provides an
efficient approach to computing `-Htilde\g`

, where `Htilde`

is equal
to `H`

if `H`

is positive definite, and otherwise is a
positive definite matrix that is "spiritually like `H`

."

The approach favored here is different from the well-known
Gill-Murray-Wright approach of computing the Cholesky factorization of
`H+E`

, where `E`

is a minimal correction needed to make `H+E`

positive-definite (sometimes known as GMW81). See the discussion
starting
here
for justification; briefly, the idea of a small correction conflates
large negative eigenvalues with numerical roundoff error, which (when
stated that way) seems like a puzzling choice. In contrast, this
package provides methods that are largely equivalent to taking the
absolute value of the diagonals D in an LDLT factorization, and setting
any "tiny" diagonals (those consistent with roundoff error) to 1. For
a diagonal matrix with some entries negative, this results in
approximately twice the correction used in GMW81.

Given a symmetric matrix `H`

, compute a positive factorization `F`

like this:

```
F = cholesky(Positive, H, [pivot=Val{false}])
```

Pivoting (turned on with `Val{true}`

) can make the correction smaller
and increase accuracy, but is not necessary for existence or stability.

For a little more information, call `ldlt`

instead:

```
F, d = ldlt(Positive, H, [pivot=Val{false}])
```

`F`

will be the same as for `cholesky`

, but this also returns `d`

, a
vector of `Int8`

with values +1, 0, or -1 indicating the sign of the
diagonal as encountered during processing (so in order of rows/columns
if not using pivoting, in order of pivot if using pivoting). This
output can be useful for determining whether the original matrix was
already positive (semi)definite.

Note that `cholesky`

/`ldlt`

can be used with any matrix, even
those which lack a conventional LDLT factorization. For example, the
matrix `[0 1; 1 0]`

is factored as `L = [1 0; 0 1]`

(the identity matrix),
with all entries of `d`

being 0. Symmetry is assumed but not checked;
only the lower-triangle of the input is used.

`cholesky`

is recommended because it is very efficient. A slower alternative is to use `eigen`

:

```
F = eigen(Positive, H)
```

which may be easier to reason about from the standpoint of fundamental linear algebra.

01/26/2016

about 1 month ago

30 commits