4
$\begingroup$

Let's suppose I'm numerically solving the Poisson equation for a delta function source:

$$ \nabla^2 f(x) = \delta(x-x') $$

I can represent the Laplacian $\nabla^2$ using the finite difference method as a tri-diagonal matrix $A$. I can represent the delta function as a vector $b$ that is $0$ everywhere except at one point where it is $\frac{1}{dx}$ (here $dx$ denotes the length of the mesh).

Thus, my resulting $Ax=b$ system of equations is very sparse -- the matrix has only $3N$ nonzero elements and the $b$ vector has only $1$ nonzero element.

Is there a numerical method that is optimized for solving $Ax=b$ for this rather specific case? I have seen some solvers that utilize the sparseness of the matrix, but I haven't found one that also utilizes the sparseness of the $b$ vector. I might be naive but I feel like a sparse $b$ vector should simplify the solving a lot.

Improve this question
$\endgroup$
1
  • 6
    $\begingroup$ If you solve $A\backslash{}b$ with LU/GE (which keeps the $LU$ factors tridiagonal), then the zero pattern of $b$ implies that during the first forward substitution step $L\backslash{}b$, some of the work can be omitted, for the leading zeros at the start of $b$. This is the same thing that happens when computing the matrix inverse explicitly from its LU factorization ($U\backslash{}(L\backslash{}I)$). In general, though, sparseness of the r.h.s. is just not all that important, because most of the work is done on the matrix. $\endgroup$
    – Kirill
    Commented Apr 7, 2017 at 0:31

2 Answers 2

Reset to default
9
$\begingroup$

When looking at the solution of your system, you will find that almost all entries of $x$ are nonzero although the right-hand side is "sparse". Hence, whatever algorithm you use, it'll have to visit each and every entry at least once, so one wouldn't expect that you can save a lot of time using the sparsity of $b$.

Right-hand sides where you can save a lot of computation are, for example, multiples of eigenvectors of $A$ since the solution of $$ Ax = \alpha v_i $$ is clearly $\frac{\alpha}{\lambda_i} v_i$.

That said, every full-rank tridiagonal equation system can be solved in $O(n)$. Using multigrid, the same goes for the Poisson equation in any dimension.

Improve this answer
$\endgroup$
5
  • $\begingroup$ Thank you. Follow up question: my specific matrix actually has two off-tridiagonal elements -- what is the name of the algorithm that solves tridiagonal matrices in $\mathcal{O}(N)$ time, and does it also work in linear time for a matrix with a few elements more than the tridiagonal matrix? $\endgroup$
    – alexvas
    Commented Apr 7, 2017 at 17:59
  • 2
    $\begingroup$ In general, you want to limit the bandwidth of your matrix. If it's small, $LU$ does the trick without having to worry about too much fill-in. Ergo: It depends where your two entries are. $\endgroup$
    – Nico Schlömer
    Commented Apr 7, 2017 at 18:24
  • $\begingroup$ @alex, if your matrix is indeed tridiagonal with a few off-diagonal elements, you might want to think about using Sherman-Morrison-Woodbury in tandem with the tridiagonal solver you already have. $\endgroup$
    – J. M.
    Commented Apr 8, 2017 at 14:47
  • $\begingroup$ @J.M. Thanks for the suggestion. Just to be clear: you are suggesting that I write my matrix M as $M = A + \Delta$, where $A$ is tridiagonal and $\Delta$ contains the two off-tridiagonal entries. I would then find $u$ and $v$ such that $\Delta = uv^T$ and use the Sherman-Morrison formula to compute $M^{-1}$. Do you think this is faster than performing an LU decomposition on the original $M$? $\endgroup$
    – alexvas
    Commented Apr 8, 2017 at 23:18
  • $\begingroup$ @alex, "Do you think this is faster than performing an LU decomposition on the original $\mathbf M$" - well, that's an experiment you need to try yourself, as I do not have your matrices or your software at hand. It doesn't necessarily have to be $\mathbf u\mathbf v^\top$; if indeed you have a "cyclic tridiagonal" system, I believe you can express $\mathbf \Delta$ as a product of two rank-2 matrices. Also, you don't need to perform an inversion in any case; SMW is still applicable if you can maintain an LU decomposition of the "original" matrix. $\endgroup$
    – J. M.
    Commented Apr 9, 2017 at 2:19
4
$\begingroup$

It is unclear to me from your question whether the answer is of theoretical or practical interest. I'll address both.

For this to be of practical importance, either your system would need to have a very large number of equations or you need to solve it many times. I suspect that simply calling the Lapack function dptsv, which is specific to solving a symmetric tridiagonal system will be significantly faster than invoking MATLAB backslash on a sparse matrix. And it will be close enough to the optimal-time approach that spending further effort would not be worthwhile. It does not exploit a sparse right-hand-side, however. This algorithm is described in section 4.3.6 of Golub and Van Loan. A total of $8n$ floating point operations (flops) are required; $3n$ flops for the factorization of the matrix and $5n$ flops to compute the solution given the right-hand-side.

Chapter 3 of Davis' book Direct Methods for Sparse Linear Systems deals with solving a linear system with a triangular matrix on the left-hand side. He spends considerable time on this topic because most of the methods for general sparse matrices in this book rely on solving these triangular systems $n$ times. He states that these methods would be impractical if each triangular solve assumed that the right-hand-side vector was dense. So he describes a linear solve algorithm that allows for a general, sparse right-hand-side vector and implements this in function cs_spsolve. To solve a symmetric positive definite system with the software described in this book would require a call to cs_chol to factor the sparse matrix followed by two calls to cs_spsolve. This is interesting theoretically but I'm doubtful it would result in less computational time than the Lapack approach unless you are solving for many right-hand-sides (e.g. many cases with a delta function at a different point in the mesh).

Improve this answer
$\endgroup$

Your Answer

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

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