New user here, trying to find an efficient way to iteratively solve lots of wide matrices in a loop.
I came across LinearSolve.jl , but haven’t seen anything in the documentation about whether it can solve non-square matrices.
Simply doing something like this throws an error:
A=rand(2,4);B=rand(2);prob=LinearProblem(A,B);solve(prob)
ERROR DimensionMismatch “matrix is not square”
Any help would be appreciated. You’re all awesome.
If the matrix is non square your system of equations is either underdetermined (more unknowns than independent equations) or it’s overdetermined.
Neither case fits to what Linearsolve.jl does.
If you have the first case, think again 
If you have the second case, you should look at running a regression or optimization instead. There’s many options for that, for example: GitHub - JuliaNLSolvers/LsqFit.jl: Simple curve fitting in Julia
Edit: welcome to the community!
You can solve consistent underdetermined systems with \, the LQ factorization (if your matrix is dense), or the QR of the transpose (if sparse).
If you are looking for iterative methods, we have several in Krylov.jl that fit the bill: GitHub - JuliaSmoothOptimizers/Krylov.jl: A Julia Basket of Hand-Picked Krylov Methods
We can add this case. Open an issue and it should have a quick solution. Looks to just be a bug in the automatic algorithm choice as a form of QR would work.
To complete the answer of @dpo, you can use the methods lsmr when m > n and craigmr when m < n of Krylov.jl.
m is the number of rows of A and n its number of columns.
If you want to solve multiple least-norm / least-squares in a loop you can use in-place methods if the linear systems have the same size.
A = rand(3, 5)
b = rand(3)
solver = CraigmrSolver(A, b) # you can also use `solver = CraigmrSolver(m, n, S)` where (m,n) = size(A) and S = typeof(b)
craigmr!(solver, A, b)
solver.x # the solution of the least-norm problem
solver.y # Lagrange multipliers
norm(b - A * solver.x)
norm(b - A * A' * solver.y)
norm(solver.x - A' * solver.y)
A = rand(5, 3)
b = rand(5)
solver = LsmrSolver(A, b) # you can also use `solver = LsmrSolver(m, n, S)` where (m,n) = size(A) and S = typeof(b)
lsmr!(solver, A, b)
solver.x # the solution of the least-squares problem
norm(A' * b - A' * A * solver.x)