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 theproperties of the rank of matrix products and the fact that matrices generated withndarray_linalg::generate::random are almost always full-rank. So, given two random matrices of shapem × r andr × n, wherer <= min(m, n), their product will almost always have rankr.

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.fnrandom_with_rank<A,Sh>(shape:Sh,rank:usize) ->Array2<A>whereA:Scalar,Sh:ShapeBuilder<Dim =Ix2>,{letmut out =Array2::zeros(shape);assert!(rank <=usize::min(out.nrows(), out.ncols()));for _in0..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);ifletOk(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 usesgeneral_mat_mul to write the result of the matrix multiplication intoout, so that we have easy control over its layout.

Codecov Report

Merging#299 (970232c) intomaster (082f01d) willincrease coverage by0.17%.
The diff coverage is97.43%.

❗ Current head970232c differs from pull request most recent head942b7d7. Consider uploading reports for the commit942b7d7 to get more accurate results

@@            Coverage Diff             @@##           master     #299      +/-   ##==========================================+ Coverage   89.01%   89.19%   +0.17%==========================================  Files          71       75       +4       Lines        3578     3656      +78     ==========================================+ Hits         3185     3261      +76- Misses        393      395       +2

Continue to review full report at Codecov.

Legend -Click here to learn more
Δ = absolute <relative> (impact),ø = not affected,? = missing data
Powered byCodecov. Last update082f01d...942b7d7. Read thecomment docs.

@emiddleton, I had the same thought, but you beat me to it! I wrote fully-tested implementationhere, 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 matrixx. Most of the timex will have full column rank. Whenx is full rank, QR decomposition is about 4 times faster than singular value decomposition withdgesvd() and at least twice as fast asdegsdd(). There is a separatepinv_svd() method when the caller believes the matrix is rank-deficient.

What are your thoughts?