Ways to solve $Ax=b$ for a sparse (banded) $A$ with updates

Question

I want to solve the time-dependent Schrodinger Equation using the Crank-Nicolson scheme. I end up with the following matrix equation

A_fixed = ...  # the matrix that does the finite difference method
B_fixed = ...  # the matrix that does the finite difference method
A = A_fixed - diags(V)  # Add the potential term V to the diagonal of A_fixed
B = B_fixed + diags(V)  # Add the potential term V to the diagonal of A_fixed  

# solve for psi_new at each time step:
A.dot(psi_new) = B.dot(psi_old)   # psi_old is the solution at the previous time-step

And $A$ and $B$ are sparse matrices with all the terms along the main/off-diagonals (is this called 'banded'?). What are the best ways to solve for the above matrix equation?

I am using the scipy library sp.sparse.linalg, where it provides a solve method which uses LU decomposition. If the potential term $V$ is time-independent, then both $A$ and $B$ are actually fixed throughout the whole computation, so I can use sp.sparse.linalg.factorized to pre-factorize $A$. This gives me like a 2x speed up.

But what if $V$ is time-varying, so $A$ and $B$ changes at each time step? Are there are algorithms which work for this specific problem (where both $A$ and $B$ consist of fixed terms and some time-varying terms which are only added to the main diagonal)?

Answer 1

If the only non-zero entries of $A_{ij}$ have $j$ in $\{i - 1, i, i + 1\}$, then $A$ is a banded matrix with bandwidth 1. More generally, you can talk about matrices of bandwidth $k$ where $k$ is any integer. For example, if you were using a higher-order finite difference discretization that used more points to calculate a derivative, you'd get higher-bandwidth matrices.

It's possible to solve a system with a banded matrix in $O(\text{bandwidth}^2\cdot n)$ time, where $n$ is the size of the matrix; this is essentially the Thomas algorithm. Since you're using scipy already, it has a built-in banded matrix solver, although you'll probably have to change how you store the matrices to the banded format. A banded solver should be substantially faster than a general sparse factorization because it avoids having to store the factored form of the matrix.

As far as the time-varying terms are concerned, your best bet is to just call the banded solver on the updated matrix on every timestep. If you were using iterative methods and the time-varying terms were small in magnitude compared to the time-independent terms, there might be some advantage to be had if you were using an iterative linear solver.

While your matrix is banded, many people have asked about updating the LDU factorization of a general sparse or dense matrix after adding, say, a multiple of the identity matrix or a diagonal matrix. For general matrices and full-rank updates, you are out of luck and need to use iterative methods if you want to leverage information from past solves.