A Geometric Algorithm for Solving Linear Systems
Based on the geometric {\it Triangle Algorithm} for testing membership of a point in a convex set, we present a novel iterative algorithm for testing the solvability of a real linear system $Ax=b$, where $A$ is an $m \times n$ matrix of arbitrary rank. Let $C_{A,r}$ be the ellipsoid determined as the image of the Euclidean ball of radius $r$ under the linear map $A$. The basic procedure in our algorithm computes a point in $C_{A,r}$ that is either within $\varepsilon$ distance to $b$, or acts as a certificate proving $b \not \in C_{A,r}$. Each iteration takes $O(mn)$ operations and when $b$ is well-situated in $C_{A,r}$, the number of iterations is proportional to $\log{(1/\varepsilon)}$. If $Ax=b$ is solvable the algorithm computes an approximate solution or the minimum-norm solution. Otherwise, it computes a certificate to unsolvability, or the minimum-norm least-squares solution. It is also applicable to complex input. In a computational comparison with the state-of-the-art algorithm BiCGSTAB ({\it Bi-conjugate gradient method stabilized}), the Triangle Algorithm is very competitive. In fact, when the iterates of BiCGSTAB do not converge, our algorithm can verify $Ax=b$ is unsolvable and approximate the minimum-norm least-squares solution. The Triangle Algorithm is robust, simple to implement, and requires no preconditioner, making it attractive to practitioners, as well as researchers and educators.