# Methods for solving rectangular, full-rank systems of equations — which is best?

Suppose I have a large, sparse, $m \times n$ matrix $A$, with $m \gt n$ and $\text{rank}(A) = n$. I wish to solve $Ax=b$.

Suppose I know that $A$ has the following characteristics:

• $A$ is somewhat diagonally-dominant (it may have a few side-bands offset quite a bit from the main diagonal, and a few rows do not fit the band structure).
• $A$ is not symmetric.
• $A$ is not positive-definite.
• $rank(A) = n$ (a unique solution exists!).
• $\text{cond}(A^T A)$ is quite high: $10^{\gt 20}$.
• $m$ is only slightly larger than $n$.

Suppose I require an accuracy of one part in $10^{\gt 9}$, or roughly 9 digits of precision.

This unique problem is neither under- nor over-determined, but it involves a rectangular matrix. It is not in the realm of either least-squares or minimum length. This places constraints on which approaches would be suitable to solve it.

I know I could use Gaussian elimination, SVD (projecting into the orthogonal complement to the null space), or a few other methods. I also know several other methods (e.g. LU decomposition) require a square matrix, so they could be used to solve $\left(A^T A\right) x = A^T b$, though I'm concerned that the high condition number of $A^T A$ would be prohibitive for these approaches.

What methods would be appropriate to solve this system of equations, and what are their relative advantages/disadvantages?

In particular:

• Which method(s) hold up best against high condition numbers?
• Which method(s) offer the best accuracy?
• Which method(s) are the quickest / most efficient?
• I'm curious- what reason do you have to believe that $b$ is in the range of $A$? If not, then you will need a least squares solution. – Brian Borchers Mar 23 '16 at 23:29
• The problem that inspired me to ask this general question was the singular Poisson equation with Neumann boundary conditions, with an added constraint that endures uniqueness. In this case, $m=n+1$. Ordinarily, the Poisson equation with Neumann boundary conditions has infinitely many solutions which are related via an additive constant. That's why it's singular. Constraining the value at one point (the $(n+1)^{th}$ row) is sufficient to make it non-singular. – jvriesem Mar 23 '16 at 23:35
• @BrianBorchers I think that's sufficient to put $b$ in the range of $A$, though I'm still learning about what that means. – jvriesem Mar 23 '16 at 23:49
• @BrianBorchers I should add: I know the problem is well-posed and a unique solution exists. It is for this reason I believe some method should be able to solve the matrix without resorting to a least squares approach...(I think!) – jvriesem Mar 23 '16 at 23:54
• From the way you've described diagonally dominant, it sounds like you're describing the sparsity/band-structure of the matrix. Can you confirm if indeed it is 'somewhat' diagonally dominant. – namu Mar 25 '16 at 0:52