Non-square matrices return errors - no default method implemented

[l8l]

Describe the bug 🐞

A \ b works when A is of dimension $n\times p$ and $b$ is of dimension $n$. However, for LinearSolve, the code

using LinearSolve
lin_prob = LinearProblem(A, b)
solution = LinearSolve.solve(lin_prob)

returns the error

ERROR: DimensionMismatch: matrix is not square: ...

In the documentation of LinearSolve.jl, I also did not find examples of non-square matrices, yet it is advertised that LinearSolve.jl is supposed to replace the \ - Operator, which works fine in my example.

This discussion on julia discourse revealed that " A \ b , for a “tall” matrix, automatically switches to finding the least-square solution (by pivoted QR), while for a “wide” matrix, it switches to finding the minimum-norm solution."

I found out, that least square and least norm solutions are at least available with the following Krylov methods:

using LinearSolve
lin_prob = LinearProblem(A, b)
x_least_squares = LinearSolve.solve(lin_prob,LinearSolve.KrylovJL_LSMR())
x_least_norm = LinearSolve.solve(lin_prob,LinearSolve.KrylovJL_CRAIGMR())

but there is no method that does, what \ seems to do in this case (first check whether A is tall or wide, then act accordingly) and maybe there is an even more optimal behavior, or better heuristics?

Expected behavior

The solver should check automatically if A is a square matrix, or not, and if not, should apply appropriate methods, that follow appropriate heuristics. For example, one could follow the strategy, that A \ b implements for non-square matrices. In the end,

using LinearSolve
lin_prob = LinearProblem(A, b)
solution = LinearSolve.solve(lin_prob)

should just work for non-square matrices (possibly with a note on the chosen solver-strategy, or a remark that A is not square, if desired).

Minimal Reproducible Example 👇

See above.

Error & Stacktrace ⚠️

ERROR: DimensionMismatch: matrix is not square: ...

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

  LinearSolve v2.35.0
  LinearAlgebra v1.11.0

• Output of versioninfo()

Julia Version 1.11.0
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: Intel(R) Xeon(R) CPU @ 3.50GHz
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, skylake-avx512)

[ChrisRackauckas]

Please share a runnable MWE. Non-square works just fine in this:

julia> using LinearSolve

julia> A = rand(4,5); b = rand(5)
5-element Vector{Float64}:
 0.6457484096650807
 0.6939469552745942
 0.41480950583430376
 0.5372284506183679
 0.754965753202426

julia> lin_prob = LinearProblem(A, b)
LinearProblem. In-place: true
b: 5-element Vector{Float64}:
 0.6457484096650807
 0.6939469552745942
 0.41480950583430376
 0.5372284506183679
 0.754965753202426

julia> solution = LinearSolve.solve(lin_prob)
retcode: Default
u: 5-element Vector{Float64}:
  0.0868623000795956
 -0.0009048226591591323
  0.2049022635171876
  0.25749407914124417
  0.360911144674564

[l8l]

Okay, thanks for this lightning-fast reply Chris! Interesting. You are right! It seems that the error occurred in my case because the matrix came from an adjoint:

julia> A = rand(2,3)'; b = rand(3);

julia> LinearSolve.solve(LinearProblem(A, b))
ERROR: DimensionMismatch: matrix is not square: dimensions are (3, 2)
...

julia> LinearSolve.solve(LinearProblem(A, b),LinearSolve.KrylovJL_LSMR())
retcode: Success
u: 2-element Vector{Float64}:
 -0.3774481167815746
  1.1609616348557

julia> LinearSolve.solve(LinearProblem(Matrix(A), b))
retcode: Default
u: 2-element Vector{Float64}:
 -0.3774481167815767
  1.1609616348557

julia> A
3×2 adjoint(::Matrix{Float64}) with eltype Float64:
 0.258806  0.274645
 0.446531  0.922608
 0.520927  0.446123

But shouldn't the default method also just work for adjoints?

[ChrisRackauckas]

It should, this seems like a bug.