Pivoted Cholesky

Start by reading the existing ndarray_linalg::cholesky entry points, especially factorizec and CholeskyFactorized, then compare their behavior with LAPACK's dpstrf example. Define how the factor, pivot, and rank should be returned and how a truncated factor could initialize CholeskyFactorized; done means the pivoted PSD case and the reduced-system reuse path are specified and tested.

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.