1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119

// Copyright 2018-2020 argmin developers
//
// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
// http://apache.org/licenses/LICENSE-2.0> or the MIT license <LICENSE-MIT or
// http://opensource.org/licenses/MIT>, at your option. This file may not be
// copied, modified, or distributed except according to those terms.

//! Modified Cholesky decompositions
//!
//! Given a symmetric matrix A which is potentially not positive definite, a modified Cholesky
//! algorithm obtains the Cholesky decomposition `LL^T` of the positive definite matrix
//! `P(A + E)P^T` where `E` is symmetric and `>= 0`, `P` is a permutation matrix and `L` is lower
//! triangular.
//! If `A` is already positive definite, then `E = 0`.
//! The perturbation `E` should be as small as possible for `A + E` to be "sufficiently positive
//! definite".
//! This is used in optimization methods where indefinite Hessians can be problematic.
//!
//! This crate implements the algorithms by Gill, Murray and Wright
//! ([GMW81](trait.ModCholeskyGMW81.html)) and Schnabel and Eskow
//! ([SE90](trait.ModCholeskySE90.html) and [SE99](trait.ModCholeskySE99.html)).
//! All algorithms are currently based on `ndarray` but will also be implemented for `nalgebra` in
//! the future.
//!
//! # Example
//!
//! ```rust
//! # use modcholesky::utils::{diag_mat_from_arr, index_to_permutation_mat};
//! # use approx::AbsDiffEq;
//! use modcholesky::ModCholeskySE99;
//!
//! let a = ndarray::arr2(&[[1.0, 1.0, 2.0],
//!                         [1.0, 1.0, 3.0],
//!                         [2.0, 3.0, 1.0]]);
//!
//! // Perform modified Cholesky decomposition
//! // The `Decomposition` struct holds L, E and P
//! let decomp = a.mod_cholesky_se99();
//!
//! println!("L:\n{:?}", decomp.l);
//! println!("E:\n{:?}", decomp.e);
//! println!("P:\n{:?}", decomp.p);
//!
//! # let l = decomp.l;
//! #
//! # let res_l: ndarray::Array2<f64> = ndarray::arr2(&[
//! #     [1.732050807568877, 0.0, 0.0],
//! #     [0.5773502691896257, 1.698920954907997, 0.0],
//! #     [1.154700538379251, 1.37342077428181, 0.006912871809428971],
//! # ]);
//! # let res = res_l.dot(&res_l.t());
//! #
//! # assert!(res_l.abs_diff_eq(&l, 1e-12));
//! #
//! # let e = diag_mat_from_arr(decomp.e.as_slice().unwrap());
//! # let p = index_to_permutation_mat(decomp.p.as_slice().unwrap());
//! # let paptpept = p.dot(&a.dot(&p.t())) + p.dot(&e.dot(&p.t()));
//! # assert!(paptpept.abs_diff_eq(&l.dot(&l.t()), 1e-12));
//! # assert!(l.dot(&(l.t())).abs_diff_eq(&res, 1e-12));
//! ```
//!
//! # TODOs
//!
//! * Be generic over element type (f64/f32/complex32/complex64)
//! * Implement algorithms for nalgebra types
//!
//! # References
//!
//! * Philip E. Gill, Walter Murray and Margaret H. Wright.
//!   Practical Optimization.
//!   Emerald Group Publishing Limited. ISBN 978-0122839528. 1982
//! * Semyon Aranovich Gershgorin.
//!   Über die Abgrenzung der Eigenwerte einer Matrix.
//!   Izv. Akad. Nauk. USSR Otd. Fiz.-Mat. Nauk, 6: 749–754, 1931.
//! * Robert B. Schnabel and Elizabeth Eskow.
//!   A new modified Cholesky factorization.
//!   SIAM J. Sci. Stat. Comput. Vol. 11, No. 6, pp. 1136-1158, November 1990
//! * Elizabeth Eskow and Robert B. Schnabel.
//!   Algorithm 695: Software for a new modified Cholesky factorization.
//!   ACM Trans. Math. Softw. Vol. 17, p. 306-312, 1991
//! * Sheung Hun Cheng and Nicholas J. Higham.
//!   A modified Cholesky algorithm based on a symmetric indefinite factorization.
//!   SIAM J. Matrix Anal. Appl. Vol. 19, No. 4, pp. 1097-1110, October 1998
//! * Robert B. Schnabel and Elizabeth Eskow.
//!   A revised modified Cholesky factorization.
//!   SIAM J. Optim. Vol. 9, No. 4, pp. 1135-1148, 1999
//! * Jorge Nocedal and Stephen J. Wright.
//!   Numerical Optimization.
//!   Springer. ISBN 0-387-30303-0. 2006.

// necessary to get clippy to shut up about the s! macro
#![allow(clippy::deref_addrof)]
// I really do not like the a..=b syntax
#![allow(clippy::range_plus_one)]
// allow single-character names
#![allow(clippy::many_single_char_names)]

mod gershgorin;
mod gmw81;
mod se90;
mod se99;
pub mod utils;

pub use crate::gershgorin::*;
pub use crate::gmw81::ModCholeskyGMW81;
pub use crate::se90::ModCholeskySE90;
pub use crate::se99::ModCholeskySE99;

pub struct Decomposition<L, E, P> {
    pub l: L,
    pub e: E,
    pub p: P,
}

impl<L, E, P> Decomposition<L, E, P> {
    pub fn new(l: L, e: E, p: P) -> Self {
        Decomposition { l, e, p }
    }
}