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.

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)