Allow fallback value for Cholesky decomposition

This is useful in case where the values are very close to zero or
possibly even slightly negative. This can quite easily happen due to
numerical errors.

A common strategy is to replace these values with a small epsilon value
that keeps the matrix SPD. Some libraries even do this by default (such
as https://github.com/STOR-i/GaussianProcesses.jl/issues/1).

We point the user to `LU` decomposition and also make it clearer that
the method is basically a hack. The public method no longer takes an
`Option` which didn't really make sense. A private method is used to not
repeat implementation in `new`.
This commit is contained in:
Mateusz Kowalczyk 2021-08-30 11:52:26 +09:00
parent a4c2ca941d
commit fc3526b58a
No known key found for this signature in database
GPG Key ID: 2CFE145ADE8C2E62
1 changed files with 44 additions and 11 deletions

View File

@ -171,7 +171,32 @@ where
/// ///
/// Returns `None` if the input matrix is not definite-positive. The input matrix is assumed /// Returns `None` if the input matrix is not definite-positive. The input matrix is assumed
/// to be symmetric and only the lower-triangular part is read. /// to be symmetric and only the lower-triangular part is read.
pub fn new(mut matrix: OMatrix<T, D, D>) -> Option<Self> { pub fn new(matrix: OMatrix<T, D, D>) -> Option<Self> {
Self::new_internal(matrix, None)
}
/// Attempts to approximate the Cholesky decomposition of `matrix` by
/// replacing non-positive values on the diagonals during the decomposition
/// with the given `substitute`.
///
/// [`try_sqrt`](ComplexField::try_sqrt) will be applied to the `substitute`
/// when it has to be used.
///
/// If your input matrix results only in positive values on the diagonals
/// during the decomposition, `substitute` is unused and the result is just
/// the same as if you used [`new`](Cholesky::new).
///
/// This method allows to compensate for matrices with very small or even
/// negative values due to numerical errors but necessarily results in only
/// an approximation: it is basically a hack. If you don't specifically need
/// Cholesky, it may be better to consider alternatives like the
/// [`LU`](crate::linalg::LU) decomposition/factorization.
pub fn new_with_substitute(matrix: OMatrix<T, D, D>, substitute: T) -> Option<Self> {
Self::new_internal(matrix, Some(substitute))
}
/// Common implementation for `new` and `new_with_substitute`.
fn new_internal(mut matrix: OMatrix<T, D, D>, substitute: Option<T>) -> Option<Self> {
assert!(matrix.is_square(), "The input matrix must be square."); assert!(matrix.is_square(), "The input matrix must be square.");
let n = matrix.nrows(); let n = matrix.nrows();
@ -187,17 +212,25 @@ where
col_j.axpy(factor.conjugate(), &col_k, T::one()); col_j.axpy(factor.conjugate(), &col_k, T::one());
} }
let diag = unsafe { matrix.get_unchecked((j, j)).clone() }; let sqrt_denom = |v: T| {
if !diag.is_zero() { if v.is_zero() {
if let Some(denom) = diag.try_sqrt() { return None;
unsafe {
*matrix.get_unchecked_mut((j, j)) = denom.clone();
}
let mut col = matrix.slice_range_mut(j + 1.., j);
col /= denom;
continue;
} }
v.try_sqrt()
};
let diag = unsafe { matrix.get_unchecked((j, j)).clone() };
if let Some(denom) =
sqrt_denom(diag).or_else(|| substitute.clone().and_then(sqrt_denom))
{
unsafe {
*matrix.get_unchecked_mut((j, j)) = denom.clone();
}
let mut col = matrix.slice_range_mut(j + 1.., j);
col /= denom;
continue;
} }
// The diagonal element is either zero or its square root could not // The diagonal element is either zero or its square root could not