2
$\begingroup$

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?
Improve this question
$\endgroup$
6
  • 1
    $\begingroup$ 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. $\endgroup$
    – Brian Borchers
    Mar 23, 2016 at 23:29
  • $\begingroup$ 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. $\endgroup$
    – jvriesem
    Mar 23, 2016 at 23:35
  • $\begingroup$ @BrianBorchers I think that's sufficient to put $b$ in the range of $A$, though I'm still learning about what that means. $\endgroup$
    – jvriesem
    Mar 23, 2016 at 23:49
  • $\begingroup$ @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!) $\endgroup$
    – jvriesem
    Mar 23, 2016 at 23:54
  • $\begingroup$ 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. $\endgroup$
    – namu
    Mar 25, 2016 at 0:52

1 Answer 1

Reset to default
1
$\begingroup$

I suggest Tim Davis' SPQR which is part of his SuiteSparse package: http://faculty.cse.tamu.edu/davis/suitesparse.html I believe it will do a reasonable job of satisfying your three criteria.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks! I'd still like to know what the pragmatic differences between different approaches are. (That's more what I'm looking for.) $\endgroup$
    – jvriesem
    Mar 23, 2016 at 23:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.