nonsquare (underdetermined) LQ solves

[stevengj]

Fixes JuliaLang/LinearAlgebra.jl#689.

Seems to be about 2x faster than pivoted QR for a large underdetermined solve, with similar accuracy (assuming full row rank):

julia> A = rand(1000,1500); b = rand(size(A,1));

julia> @btime lq($A) \ $b;
  101.955 ms (14 allocations: 19.37 MiB)

julia> @btime qr($A, Val{true}()) \ $b;
  250.584 ms (8046 allocations: 31.98 MiB)

julia> x = lq(A) \ b; norm(x - qr(A, Val{true}()) \ b) / norm(x)
3.5357944054772334e-15

[dkarrasch]

One of us is getting over- and under-determined wrong, perhaps worth a second thought. The other issue I couldn't see resolved is ldiv!, which promises solve in-place of B, but must necessarily widen/stretch B to fit the longer X in A*X = B. It's not easy to track the changes in the test, but I couldn't find tests for ldiv!.

[stevengj]

For non-square problems with more cols n than rows m, ldiv! only looks at the first m rows of its input (which should have n rows to store the output solution)

(As noted above, ldiv! is exercised by the tests because it is actually called by \.)

You’re right that I was mixing up over/under determined some of the time — it's fixed now.

[stevengj]

win32 CI failure seems to be an unrelated Error in testset SharedArrays

[stevengj]

win32 failure is now another unrelated glitch cp: cannot stat 'dist-extras/7z.*': No such file or directory.

Should be good to merge?

[stevengj]

I wonder if we should change the default A \ b polyalgorithm to use LQ for undetermined problems ("wide" A)? It still seems to be twice as fast on Julia 1.11.

[dkarrasch]

Currently, we only distinguish square from non-square, right? Sounds like a good idea.

[stevengj]

Currently, we don't have row-pivoted lq, so using lq(A) \ b as-is might be problematic for badly conditioned "wide" A (i.e. nearly dependent rows).

However, we could equivalently use pivoted QR on $A^T$ via qr(A', ColumnNorm())' \ b. That seems to be the fastest option of all, actually:

julia> A = rand(100, 1000); b = rand(100);

julia> @btime $A \ $b; @btime qr($A, ColumnNorm()) \ $b; @btime lq($A) \ $b; @btime qr($A', ColumnNorm())' \ $b;
  6.802 ms (42 allocations: 1.89 MiB)
  6.851 ms (42 allocations: 1.89 MiB)
  3.847 ms (21 allocations: 926.41 KiB)
  2.454 ms (29 allocations: 1008.34 KiB)

julia> A \ b ≈ qr(A, ColumnNorm()) \ b ≈ lq(A) \ b ≈ qr(A', ColumnNorm())' \ b  # same answer
true

Technically, the fastest option of all (by a factor of 10!) is to use $A^T (A A^T)^{-1} b$ via Cholesky, but this is potentially inaccurate because it squares the condition number:

julia> @btime (cholesky($A*$A') \ $b);
  287.375 μs (8 allocations: 157.31 KiB)

julia> A \ b ≈ A' * (cholesky(A*A') \ b)
true

(which is also why we don't use the analogous Cholesky-based algorithm in the "tall" overdetermined case).