# Faster than forward substitution?

I have a matrix of the form:

$$M:=\begin{pmatrix} S_1 & & & \\\ Q_1 & S_2 & & \\\ & ... & ... & \\\ & & Q_n & S_n\end{pmatrix}$$

where the blocks are square. The $$S_i$$ are asymmetric, the $$Q_i$$ are symmetric and positive definite. Would you have a reference/hints for efficient solution algorithms for systems of the type $$Mx=b$$?

What further properties are in case needed on the $$Q_i$$ blocks? We can assume that the resulting matrix $$M$$ is non singular.

I encountered such structure while solving parabolic PDEs on time dependent domains.

The suggested method should be faster than solving $$S_ix_i-Q_{i-1}x_{i-1}=b_i$$ recursively.

Moreover, we can assume that the number of blocks is large (e.g. $$1e4$$) and that the size of each block is also large, e.g. $$1e3$$.

Parallel strategies are very welcome.

• Do you care about all the xi's (timesteps), or just the final one (xn)? Mar 3 at 17:48
• They are all known from the start @rchilton1980 (I'm curious about how your answer would change otherwise: in fact, we are inside a Newton iteration here, so that the matrices do not depend on the unknown: but before linearization, they do, so that this is interesting for me as well) Mar 4 at 7:28
• Are the Qs and Ss banded? If so, a BandedBlockedBanded matrix solver might be appropriate (see github.com/JuliaLinearAlgebra/BlockBandedMatrices.jl). Mar 5 at 4:36
• I don't know of a paper, but it looks like they use a QR for linear solve (and also support Cholesky which won't help you) The source is here github.com/JuliaLinearAlgebra/BlockBandedMatrices.jl/blob/… and there are docs for the data layout (but not the algorithms) here julialinearalgebra.github.io/BlockBandedMatrices.jl/dev/… Mar 6 at 3:02
• Perhaps a simple suggestion: you could factorize the $S_i$ matrices in parallel, these factorizations do not depend on each other so this should parallelize well. Afterwards, the backwards substitiution would essentially be serial, but you might get a significant speedup. Mar 6 at 8:35

Let's assume that the size $$m$$ of the individual blocks is fixed, but that the number of blocks $$n$$ grows. Then the one-step-at-a-time algorithm takes $$O(m^3n)$$ operations if you chose to invert the diagonal blocks via Gauss elimination of LU decomposition. That is, the effort is $$O(n).$$ At least asymptotically, it doesn't get any faster than that. Your speedup will have to come from solving the diagonal blocks faster; for example, if you can use a multigrid method for the $$S_i$$ blocks, then you can get the complexity down to $$O(mn)=O(N)$$ where $$N$$ is the size of the matrix $$M$$ -- this is now truly something you can no longer beat as far as the complexity is concerned.
The considerations above were for counting operations. In practice, one may want to count wallclock time instead, and for that it is worth considering inverting the diagonal blocks in parallel if an explicit inversion is possible. In that case, you can do all of the inverting of the $$S_i$$ on background threads and the complexity of the main thread drops from $$O(m^3n)$$ to $$O(m^2n)$$ where the $$O(m^2)$$ factor results from either (i) applying $$S_i^{-1}$$ or (ii) the forward-backward substitution when applying an LU decomposition of $$S_i$$.
• Thanks for the insightful answer. Parallel in time methods seem also to need many CPUs to yield a speedup, e.g. $50, 100$. If this is not the case, is there a way out here, while keeping the solution algorithm for the $S_i$ block fixed? E.g. to your experience, would solving $Mx=b$ directly, with an iterative algorithm (e.g. GMRES), bring any wall clock time additional speed up? Mar 7 at 6:01
• You can run GMRES on the whole system, of course (assuming that you can actually store all of it in memory at once), but GMRES will only converge quickly if you have good preconditioners. A good preconditioner would like do a triangular forward-substitution sweep through all variables of the form $x_i = \widehat{S_i^{-1}} (b_i - Q_{i-1}x_{i-1})$ with an approximation $\widehat{S_i^{-1}}$ of $S_i^{-1}$. This is the equivalent of a block-SOR algorithm that utilizes the block-triangular structure of the matrix. Whether that is faster than running GMRES on each time step is unclear to me. Mar 7 at 18:47