4
$\begingroup$

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.

Improve this question
$\endgroup$
8
  • $\begingroup$ Do you care about all the xi's (timesteps), or just the final one (xn)? $\endgroup$
    – rchilton1980
    Mar 3 at 17:48
  • 1
    $\begingroup$ 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) $\endgroup$
    – Leonardo
    Mar 4 at 7:28
  • 1
    $\begingroup$ Are the Qs and Ss banded? If so, a BandedBlockedBanded matrix solver might be appropriate (see github.com/JuliaLinearAlgebra/BlockBandedMatrices.jl). $\endgroup$
    – Oscar Smith
    Mar 5 at 4:36
  • 1
    $\begingroup$ 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/… $\endgroup$
    – Oscar Smith
    Mar 6 at 3:02
  • 1
    $\begingroup$ 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. $\endgroup$
    – Thijs Steel
    Mar 6 at 8:35

2 Answers 2

Reset to default
3
$\begingroup$

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$.

@FedericoPolini points out the idea of parallel-in-time methods which are indeed used for exactly this kind of purpose. They can, in some circumstances, reduce the wallclock time for the solution of equations such as yours, but generally at a very large cost in terms of CPU time: They are quite wasteful with CPU time to achieve modest wallclock-time speedups.

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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? $\endgroup$
    – Leonardo
    Mar 7 at 6:01
  • 2
    $\begingroup$ 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. $\endgroup$
    – Wolfgang Bangerth
    Mar 7 at 18:47
3
$\begingroup$

I am not an expert in the field, but since there are no other answers at least I can suggest you a keyword to start a literature search: "parallel-in-time integration methods". This is a family of methods that aim to solve this exact problem, essentially applying iterative methods to your linear system.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Interesting suggestion, thanks $\endgroup$
    – Leonardo
    Mar 5 at 17:44

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.