Implement solve (\) for general rectangular static matrices

[cgarling]

This should close #1292. Results are much faster now that it no longer falls back to LinearAlgebra.jl.

Note that in the calculation for the underdetermined case on line 77, q.R' * ((q.R * q.R') \ y) is an approximation for the pseudoinverse pinv which in general may be more accurate, but it also seemed a lot slower in my testing so I did not use it here. It looks like there might a better-conditioned way to calculate it using the results of the QR decomposition but I am not a linear algebra expect so some help there would be appreciated.

[timholy]

Suggested change

	q.R' * ((q.R * q.R') \ y)
	q.R' * (q.R' \ (q.R \ y))

should be a lot faster, if q.R is annotated as UpperTriangular (you could wrap it if not).

[cgarling]

Thanks for the review, couple questions:

q.R has dimensions (Sa[1], Sa[2]) with Sa[1] < Sa[2]. I don't see how to use UpperTriangular on a non-square matrix as, for example, UpperTriangular(rand(2,3)) gives ERROR: DimensionMismatch: matrix is not square: dimensions are (2, 3). Am I missing something?

Using the expression q.R' * (q.R' \ (q.R \ y)) you provided results in a stack overflow error, I believe because of the (q.R \ y) which will recurse infinitely because q.R has the same dimensions as argument a and y has the same dimensions as argument b. Any suggestions to improve this are welcome

[mateuszbaran]

It's a good idea but why don't just copy the logic from https://github.com/JuliaLang/LinearAlgebra.jl/blob/16f64e78769d788376df0f36447affdb7b1b3df6/src/qr.jl#L700 and adapt to StaticArrays.jl? Also, I'd prefer to have the logic for wide/thin case in (\)(Q::QR, b::StaticVecOrMat), and _solve_general only calling qr and using (\) on that QR decomposition.

I also support the use of the UpperTriangular wrapper. You can also safely use view instead of constructing a new SMatrix, it does the right thing for SMatrix, MMatrix and SizedMatrix.

[cgarling]

Updates implementing suggestions have been pushed

[mateuszbaran]

Thanks for the updates. Could you also address the last remaining point from my post? To be more specific, I think we could use the logic of _wide_qr_ldiv! from LinearAlgebra.jl instead of q.R' * ((q.R * q.R') \ y).

[cgarling]

I spent some time look at the code in _wide_qr_ldiv! and do not see an easy way to adapt it. I am not an expert in linear algebra or the internals of StaticArrays and I struggle adapting algorithms that write to mutable buffers (like τ in this case) and rewriting them to work with immutable StaticArrays. Basically, I think this is somewhat beyond my present ability (or would, at least, take me more time than I am willing to commit at the moment).

Is your concern one of accuracy/correctness or performance? If my solution is correct, I would ask that it be accepted as it is many times faster than the current LinearAlgebra fallback and does not allocate, which solves issues like #1292. I could add a TODO annotation in the code to implement the version with improved performance in the future.

[mateuszbaran]

I've ported that method myself and added some tests which revealed a problem. The good news is that it works fast, without any allocations (Julia can stack-allocate non-escaping mutable structs). The bad news is that LinearAlgebra by default uses pivoted QR which is necessary for rank-deficient problems (see the new test case 4). Which is slower and currently not handled correctly by solve in QR.

So, the two remaining issues:

1. If we decide to use QR by default in \, it needs to be pivoted QR.
2. If we decide to implement \ for QR objects, it needs to handle pivoting.

Could you address point 2? It's actually a matter of correctness as you can see in the failing tests.

[cgarling]

Thanks for the help! One new problem is that with the pivoted QR we need to permute the result from _solve, I just pushed a commit that does that. It uses invperm from utils.jl that presently allocates for me -- I will open another PR to try to fix that. Your new test case 4 still fails, unfortunately.

Edit: It's our invperm definition here, and the solution I thought I had to remove the allocation doesn't seem to work :(

[timholy]

One risk of calling a non-public function is that it could change in future Julia releases. Better to copy the definition here.