Pivoted Cholesky #252
Description
Currently, Cholesky will fail on PSD but not PD matrices, because it calls ?pptrf.
use ndarray::{array, Array, Axis}; use ndarray_linalg::cholesky::*; fn main() -> Result<(), Box<dyn std::error::Error>> { let x1 = array![1., 2., 3., 4., 7.]; let x2 = array![-1., 2., 3., 5., 8.]; let x3 = array![1., -2., 3., 6., 9.]; let x4 = &x1 * 1. + &x2 * 2. - &x3 * 3.; let xs = [x1, x2, x3, x4]; let xs = xs .iter() .map(|x| x.view().insert_axis(Axis(1))) .collect::<Vec<_>>(); let X = ndarray::stack(Axis(1), &xs)?; println!("{:?}", X); let XTX = X.t().dot(&X); println!("{:?}", XTX); let chol = XTX.factorizec(UPLO::Lower)?; // Error: Lapack { return_code: 4 } Ok(()) } However, if we allow pivoting, then we can return a cholesky factor U, pivot matrix P such that P U^T U P^T = A for an input matrix A that's merely PSD (this also returns the rank r).
It would be nice if ndarray-linalg could also provide this pivoted version, e.g., as shown here in python:
from scipy.linalg.lapack import dpstrf import numpy as np xt = np.array([ [1, 2, 3, 4, 7], [-1, 2, 3, 5, 8], [1, -2, 3, 6, 9]]).astype(float) x = np.insert(xt, len(xt), xt[0] + 2 * xt[1] - 3 * xt[2], axis=0).T xtx = x.T.dot(x) U, P, rank, info = dpstrf(xtx) assert info > 0 assert rank == 3 P -= 1 U = np.triu(U) invP = np.empty_like(P) invP[P] = np.arange(len(P), dtype=int) print(np.linalg.norm(U.T.dot(U)[np.ix_(invP, invP)] - xtx, ord='fro')) # row indexing inverts permutations # 3.552713678800501e-14 An interesting design question would be what the interface should be. Clearly this routine should return the factor, pivot, and rank in some form. But it'd be nice if I could take my pivoted Cholesky output, truncate U to its leading principal minor of order r, and initialize a CholeskyFactorized struct directly, so that I can just re-use existing methods for solving the reduced subsystem.