I haven't had a chance to review the code, but I thought I'd respond to this portion of your comment:

I wrote some functions for creating various ranked matrices but I think there might be a better way to do this. I added a rank function to help create random matrices of various ranks.

Probably the simplest approach is to take advantage of the properties of the rank of matrix products and the fact that matrices generated with ndarray_linalg::generate::random are almost always full-rank. So, given two random matrices of shape m ×ばつ r and r ×ばつ n, where r <= min(m, n), their product will almost always have rank r.

use ndarray::prelude::*;
use ndarray::linalg::general_mat_mul;
use ndarray_linalg::{Scalar, generate::random};
/// Returns an array with the specified shape and rank.
///
/// # Panics
///
/// Panics if the rank is impossible to achieve for the given shape, i.e. if
/// it's less than the minimum of the number of rows and number of columns.
fn random_with_rank<A, Sh>(shape: Sh, rank: usize) -> Array2<A>
where
 A: Scalar,
 Sh: ShapeBuilder<Dim = Ix2>,
{
 let mut out = Array2::zeros(shape);
 assert!(rank <= usize::min(out.nrows(), out.ncols()));
 for _ in 0..10 {
 let left: Array2<A> = random([out.nrows(), rank]);
 let right: Array2<A> = random([rank, out.ncols()]);
 general_mat_mul(A::one(), &left, &right, A::zero(), &mut out);
 if let Ok(out_rank) = out.rank() {
 if out_rank == rank {
 return out;
 }
 }
 }
 unreachable!("Failed to generate random matrix of desired rank within 10 tries. This is very unlikely.");
}

This implementation uses general_mat_mul to write the result of the matrix multiplication into out, so that we have easy control over its layout.

@emiddleton, I had the same thought, but you beat me to it! I wrote fully-tested implementation here, but had not yet covered the in-place use case.

The main difference between our implementations is that my algorithm attempts QR decomposition first and then falls back to singular value decomposition. My reasoning is that the most common use case for computing the Moore-Penrose pseudoinverse is performing least-squares regression, where we would like to know the pseudoinverse of the design matrix x. Most of the time x will have full column rank. When x is full rank, QR decomposition is about 4 times faster than singular value decomposition with dgesvd() and at least twice as fast as degsdd(). There is a separate pinv_svd() method when the caller believes the matrix is rank-deficient.

What are your thoughts?