Gaussian belief propagation solvers for nonsymmetric systems of linear equations

TL;DR

This paper extends Gaussian belief propagation to nonsymmetric linear systems, introduces a new matrix inversion algorithm, and demonstrates its advantages over classical methods in terms of robustness and efficiency.

Contribution

It adapts Gaussian belief propagation for nonsymmetric matrices and proposes a new belief propagation-based matrix inversion algorithm.

Findings

Algorithms outperform classical solvers like Gauss-Seidel and BiCGSTAB.

Gaussian belief propagation serves as an effective multigrid smoother.

Proposed methods show increased robustness in solving linear elliptic equations.

Abstract

In this paper, we argue for the utility of deterministic inference in the classical problem of numerical linear algebra, that of solving a linear system. We show how the Gaussian belief propagation solver, known to work for symmetric matrices can be modified to handle nonsymmetric matrices. Furthermore, we introduce a new algorithm for matrix inversion that corresponds to the generalized belief propagation derived from the cluster variation method (or Kikuchi approximation). We relate these algorithms to LU and block LU decompositions and provide certain guarantees based on theorems from the theory of belief propagation. All proposed algorithms are compared with classical solvers (e.g., Gauss-Seidel, BiCGSTAB) with application to linear elliptic equations. We also show how the Gaussian belief propagation can be used as multigrid smoother, resulting in a substantially more robust solver than the one based on the Gauss-Seidel iterative method.

Tables

Table 1. Table 1: Computational complexity of GaBP (error correction scheme) with precomputed Λ Λ \Lambda messages in comparison with the classical Gauss-Seidel relaxation schemes. Line GaBP refers to the partition presented in Figure   3a and x / y 𝑥 𝑦 x/y -GS is a classical line smoother. As one can see from the theory of generalized GaBP, it is not possible to apply line GaBP for 9 9 9 points stencil, because large regions do not cover all edges of the original graph. However, one indeed can construct line GaBP smoothers for this case, too, but we do not consider them here.

Sweeps	Stencil	GaBP	line GaBP	GS	$x/y$-GS
$1$	5 points	$18N$	$38N$	$9N$	$14N$
	9 points	$32N$	N/A	$17N$	$21N$
$2$	5 points	$30N$	$65N$	$18N$	$28N$
	9 points	$56N$	N/A	$34N$	$42N$
$M\ge 3$	5 points	$12N\cdot M+6N$	$28N\cdot M+9N$	$9N\cdot M$	$14N\cdot M$
	9 points	$24N\cdot M+8N$	N/A	$17N\cdot M$	$21N\cdot M$

Equations

$$\mathbf{Ax}=\mathbf{b}$$

$${\mathbf{x}}^{\star }=\arg \underset{\mathbf{x}}{\min }\left(\frac{{\mathbf{x}}^T\mathbf{Ax}}{2}-{\mathbf{x}}^T\mathbf{b}\right).$$

$$p(\mathbf{x})\sim \exp \left(-\frac{{\mathbf{x}}^T\mathbf{Ax}}{2}+{\mathbf{b}}^T\mathbf{x}\right).$$

$${\mathbf{A}}_{22}\leftarrow {\mathbf{A}}_{22}-\frac{{\mathbf{A}}_{21}{\mathbf{A}}_{12}}{A_{11}},{\mathbf{b}}_2\leftarrow {\mathbf{b}}_2-\frac{{\mathbf{A}}_{21}{b}_1}{A_{11}}.$$

$$\mathbf{A}=\left(\begin{array}{ccccc}{A}_{11}& 0& A_{13}& 0& 0\\ 0& A_{22}& A_{23}& 0& 0\\ A_{31}& A_{32}& A_{33}& 0& A_{35}\\ 0& 0& 0& A_{44}& A_{45}\\ 0& 0& A_{53}& A_{54}& A_{55}\end{array}\right),\mathbf{b}\in {\mathbb{R}}^5.$$

$$\underset{={\tilde{A}}_{33}}{\underbrace{\left(A_{33}-\frac{A_{31}{A}_{13}}{A_{11}}-\frac{A_{32}{A}_{23}}{A_{22}}\right)}}{x}_3+A_{35}{x}_5=\underset{={\tilde{b}}_3}{\underbrace{b_3-\frac{A_{31}{b}_1}{A_{11}}-\frac{A_{32}{b}_2}{A_{22}}}};$$

$$\left(A_{55}-\frac{A_{53}{A}_{35}}{{\tilde{A}}_{33}}-\frac{A_{54}{A}_{45}}{A_{44}}\right)x_5=b_5-\frac{A_{53}{\tilde{b}}_3}{{\tilde{A}}_{33}}-\frac{A_{54}{b}_4}{A_{44}}.$$

$$\underset{={\tilde{A}}_{55}}{\underbrace{\left(A_{55}-\frac{A_{54}{A}_{45}}{A_{44}}\right)}}{x}_5+A_{53}{x}_3=\underset{={\tilde{b}}_5}{\underbrace{b_5-\frac{A_{53}{b}_4}{A_{44}}}};$$

$$\left(A_{33}-\frac{A_{35}{A}_{53}}{{\tilde{A}}_{55}}-\frac{A_{31}{A}_{13}}{A_{11}}-\frac{A_{32}{A}_{23}}{A_{22}}\right)x_3=b_3-\frac{A_{31}{b}_1}{A_{11}}-\frac{A_{32}{b}_2}{A_{22}}-\frac{A_{35}{\tilde{b}}_5}{{\tilde{A}}_{55}}.$$

$$\left(A_{33}+{\Lambda }_{31}+{\Lambda }_{32}\right)x_3+A_{35}{x}_5=b_3+{\Lambda }_{32}{\mu }_{32}+{\Lambda }_{31}{\mu }_{31}.$$

$$\left(A_{55}+{\Lambda }_{53}+{\Lambda }_{54}\right)x_5=b_5+{\Lambda }_{53}{\mu }_{53}+{\Lambda }_{54}{\mu }_{54},$$

$${\Lambda }_{53}=-\frac{A_{35}{A}_{53}}{A_{33}+{\Lambda }_{31}+{\Lambda }_{32}},{\mu }_{53}=\frac{b_3+{\Lambda }_{31}{\mu }_{31}+{\Lambda }_{32}{\mu }_{32}}{A_{35}}.$$

$$x_j=\frac{b_j+\sum _{k\in \text{neighbours of }j}{\Lambda }_{jk}{\mu }_{jk}}{A_{jj}+\sum _{k\in \text{neighbours of }j}{\Lambda }_{jk}}.$$

$$p(x)=\frac{1}{Z}\prod _{i\in \mathcal{V}}{\varphi }_i(x_i)\prod _{(ij)\in \mathcal{E}}{\psi }_{ij}(x_i,x_j)\equiv \frac{1}{Z}\exp \left(-E(\mathbf{x})\right).$$

$$p(\mathbf{x})=\frac{1}{Z}\exp \left(-\frac{{\mathbf{x}}^T\mathbf{Ax}}{2}+{\mathbf{b}}^T\mathbf{x}\right)\equiv \mathcal{N}\left(\mathbf{x}|\underset{\text{mean}}{{\mathbf{A}}^{-1}\mathbf{b}},\underset{\begin{array}{c}\text{covariance}\\ \text{matrix}\end{array}}{{\mathbf{A}}^{-1}}\right),$$

$$p_r({\mathbf{x}}_r)=\sum _{\mathbf{x}\backslash {\mathbf{x}}_r}p(\mathbf{x}).$$

$$p_i(x_i)=\mathcal{N}\left(x_i|{\left({\mathbf{A}}^{-1}\mathbf{b}\right)}_i,{\left({\mathbf{A}}^{-1}\right)}_{ii}\right)\equiv \mathcal{N}\left(x_i|{\mu }_i,{\beta }_i\right).$$

$$m_{ij}(x_j)\leftarrow \sum _{x_i}{\varphi }_i(x_i){\psi }_{ij}(x_i,x_j)\prod _{k\in N(i)\backslash j}{m}_{ki}(x_i),$$

$$b_i(x_i)\sim {\varphi }_i(x_i)\prod _{k\in N(i)}{m}_{ki}(x_i).$$

$$m_{ji}(x_i)\sim \exp \left(-\frac{{\Lambda }_{ji}{\left(x_i-{\mu }_{ji}\right)}^2}{2}\right),$$

$$\begin{array}{rl}& {\mu }_{ji}^{(n+1)}=\frac{b_j+\sum _{k\in N(j)\backslash i}{\Lambda }_{kj}^{(n)}{\mu }_{kj}^{(n)}}{A_{ji}},{\Lambda }_{ji}^{(n+1)}=-\frac{A_{ij}{A}_{ji}}{A_{jj}+\sum _{k\in N(j)\backslash i}{\Lambda }_{kj}^{(n)}};\\ & {\mu }_i^{(n)}=\frac{b_i+\sum _{j\in N(i)}{\Lambda }_{ji}^{(n)}{\mu }_{ji}^{(n)}}{A_{ii}+\sum _{j\in N(i)}{\Lambda }_{ji}^{(n)}},{\beta }_i^{(n)}=A_{ii}+\sum _{j\in N(i)}{\Lambda }_{ji}^{(n)}.\end{array}$$

$$m_{ji}^{(n)}\equiv\mu_{ji}^{(n)}\Lambda_{ji}^{(n)},~{}\widetilde{\Lambda}_{ji}^{(n)}\equiv\Lambda_{ji}^{(n)}\big{/}A_{ji}.$$

$$\begin{array}{rl}& m_{ji}^{(n+1)}={\tilde{\Lambda }}_{ji}^{(n+1)}\left(b_j+\sum _{k\in N(j)\backslash i}{m}_{kj}^{(n)}\right),{\tilde{\Lambda }}_{ji}^{(n+1)}=-\frac{A_{ij}}{A_{jj}+\sum _{k\in N(j)\backslash i}{\tilde{\Lambda }}_{kj}^{(n)}{A}_{kj}};\\ & {\mu }_i^{(n)}=\frac{b_i+\sum _{j\in N(i)}{m}_{ji}^{(n)}}{A_{ii}+\sum _{j\in N(i)}{\tilde{\Lambda }}_{ji}^{(n)}{A}_{ji}},{\beta }_i^{(n)}=A_{ii}+\sum _{j\in N(i)}{\tilde{\Lambda }}_{ji}^{(n)}{A}_{ji}.\end{array}$$

$$\text{error}=\max \left({\Vert {\tilde{\mathbf{\mu }}}^{(n+1)}-{\tilde{\mathbf{\mu }}}^{(n)}\Vert }_{\infty },{\Vert {\tilde{\mathbf{\Lambda }}}^{(n+1)}-{\tilde{\mathbf{\Lambda }}}^{(n)}\Vert }_{\infty }\right).$$

$$c_r=1-\sum _{i\in A(r)}{c}_i,$$

$$\sum _{r\in \mathcal{R},v\in r}{c}_r=\sum _{r\in \mathcal{R},e\in r}{c}_r=1.$$

$${\mathbf{m}}_{a\to b}({\mathbf{x}}_b),$$

$${\mathbf{b}}_r({\mathbf{x}}_r)\sim \prod _{i\in {\mathcal{V}}_r}{\varphi }_i(x_i)\prod _{(ij)\in {\mathcal{E}}_r}{\psi }_{ij}(x_i,x_j)\prod _{a\in B(r),b\in S(r)}{\mathbf{m}}_{a\to b}({\mathbf{x}}_b),$$

$$\forall l,r\in \mathcal{R},l\subset r\Rightarrow \sum _{x_r\backslash x_l}{\mathbf{b}}_r({\mathbf{x}}_r)={\mathbf{b}}_l({\mathbf{x}}_l).$$

$${\left({\Pi }_k\right)}_{ij}=\{\begin{array}{l}{\delta }_{ij},i,j\in k,\\ 0,\text{otherwise}.\end{array}$$

$${\mathbf{\Pi }}_k^T\mathbf{A}{\mathbf{\Pi }}_k={\left[\mathbf{A}\right]}_k,{\mathbf{\Pi }}_k^T\mathbf{b}={\left[\mathbf{b}\right]}_k;$$

Full text

1 Introduction

A basic problem of numerical linear algebra is to solve a linear equation

[TABLE]

with an invertible matrix $\mathbf{A}$. The textbook technique is the LU decomposition, equivalent to Gaussian elimination [11, ch. 3]. However, when $\mathbf{A}$ is large and sparse, algorithms that exploit sparsity are used instead of the direct elimination [8]. Among iterative methods for sparse systems, one can mention classical relaxation techniques such as Gauss-Seidel (GS), Jacobi, Richardson, and projection methods such as conjugate gradients, generalized minimal residuals, biconjugate gradients, and others [25, 24].

An easy way to understand the projection methods is to reformulate the original equation as an optimization problem [27]. For example, for a symmetric positive-definite matrix $\mathbf{A}$, one has

[TABLE]

Such a reformulation allows one to apply new techniques and leads to methods of steepest descent, conjugate directions, and cheap and efficient conjugate gradients [14].

Another reformulation of the problem is known, but is less explored. It also goes back to Gauss and his version of elimination. To derive an LU solution of (1), one can consider the probability density function $p(\mathbf{x})$ of multivariate normal distribution (see also equation (13) below)

[TABLE]

We can consider the first component $x_{1}$ of $\mathbf{x}$ and integrate it out in $(\ref{multivariate_normal})$ (a process called “marginalization”). The resulting marginal distribution for the remaining components $x_{2},x_{3},\dots$ is again multivariate normal, but with the covariance matrix given by the Schur complement of $A_{11}$111In the article we use boldface for matrices or matrix blocks, and regular font for scalar values and matrix components. In this case $A_{11}$ is an element of the matrix $\mathbf{A}$ in the first row and the first column, and $\mathbf{A_{22}}$ is a square matrix that contains all elements of $\mathbf{A}$ excluding the first row and the first column. and the mean vector modified accordingly, i.e.

[TABLE]

It is well known that the LU decomposition consists of the very same steps [28]. When $x_{1}$ is not a scalar, but a subset of variables, marginalization of multivariate normal distribution results in a block LU decomposition.

Thus, the most popular direct technique for the solution of linear equations with dense matrices is intimately connected with the marginalization problem, which belongs to the class of inference problems. Recently, many other intriguing connections between statistical inference and linear algebra have been pointed out. For instance, in [7], [12], and [4], the authors provide a method to recover the Petrov-Galerkin condition from the Bayesian update and construct a Bayesian version of the conjugate gradients. Paper [19] constructs a state-of-the-art multigrid solver using the game theory and statistical inference. These works demonstrate that ideas from statistical inference allow for new and useful insights into problems of linear algebra. It is thus reasonable to explore how other inference algorithms are translated to the realm of numerical linear algebra. Among them are expectation propagation [18], Markov chain Monte Carlo, mean field, other variational Bayesian approximations [6, chapters 8, 10], [32], and belief propagation with its generalized counterparts. The latter two are the focus of the present work.

The first comparison between classical methods and belief propagation appeared in [33]. Then in [26], authors argued explicitly for the belief propagation as a solver and later, Bickson [5] presented a more systematic treatment of the Gaussian belief propagation (GaBP) in the same context. Among other proposed methods was an algorithm that treats nonsymmetric matrices through diagonal weighting [5, 5.4] and the usual trick from linear algebra, $\mathbf{A}\rightarrow\mathbf{A}^{T}\mathbf{A}$. Both techniques are of limited use because of slow convergence in the first case and fill-in in the second. We improve on these results and propose several new algorithms.

In particular, in this work we:

• •

explain how belief propagation can be applied to nonsymmetric matrices with no computation overhead compared to the original belief propagation (which was limited to symmetric matrices) (Algorithm 1);

• •

design a family of linear solvers based on the generalized belief propagation (Algorithm 2) and relate them to the block LU decomposition (see Section 3.3);

• •

introduce a two-layer region graph and derive generalized belief propagation rules (36) that are substantially less demanding computationally compared to the basic generalized belief propagation algorithm (see (44));

• •

show how proofs of sufficient condition for convergence and consistency for the original GaBP can be modified to hold for the new algorithms (for GaBP see A, B, for generalized GaBP – C, D);

• •

explain how one can speed up GaBP using multigrid methodology which results in a robust solver (see Figure 4);

• •

implement the new algorithms [1] and benchmark them against several classical multigrid solvers.

The rest of the paper is organized as follows. In Section 2.1, we start with the intuitive explanation of GaBP based on the connection between the algorithm and Gaussian elimination. The general description of how to treat the problem (1) as an inference problem together with the basic terminology and main facts about Gaussian belief propagation are introduced in Section 2. In Section 2.3 we introduce GaBP that can be used for nonsymmetric matrices, prove consistency of the proposed algorithm, and establish a sufficient condition for convergence in appendices A, B. In Section 3, extensions of the belief propagation are given: in 3.1, we describe the generalized belief propagation algorithm (GaBP) (parent-to-child in [39]); in 3.2, we derive message update rules for region graphs with two layers, and in 3.3, we discuss the generalized GaBP from the elimination perspective and explain why the algorithm can be applied to the case $\mathbf{A}^{T}\neq\mathbf{A}$; the resulting algorithm is introduced in Section 3.4 and analyzed in C, D. Section 4 explains how to use GaBP within a multigrid scheme, we discuss smoothing properties, complexity and describe a rather unusual behaviour for singularly perturbed elliptic equations. Section 5 contains numerical examples. In Section 6, we summarize the main results and discuss possible future research.

2 Gaussian belief propagation

2.1 GaBP from the elimination perspective

To build intuition about GaBP, we consider connection with the Gaussian elimination first. Ideas of this section are similar to those in [22], but the presentation is more straightforward and after appropriate modifications (see Section 2.3), applies to non-symmetric matrices as well.

To illustrate the main ideas, consider a linear problem $\mathbf{A}\mathbf{x}=\mathbf{b}$ with the matrix and right-hand side defined as

[TABLE]

For simplicity, require that $\mathbf{A}$ be positive definite and that all elements of $\mathbf{A}$, not explicitly indicated as zeros, are nonzero.

To obtain GaBP rules, we introduce a graph of the matrix (5). For this section, it is sufficient to associate the set of vertices with the set of diagonal terms and the collection of edges with nonzero entries $A_{ij},~{}i\neq j$. One can find the resulting graph in figure 1. The correspondence between graphs and matrices is discussed in detail in two subsequent subsections 2.2, 2.3.

Suppose one wants to calculate variable $x_{5}$. To do that, we exclude variables $x_{1}$, $x_{2}$ from the equation for $x_{3}$ and then eliminate variables $x_{3}$, $x_{4}$ from the equation for $x_{5}$

[TABLE]

Figure 1a captures this particular elimination order. In the same vein, to find $x_{3}$ one may follow the order presented in figure 1b. The resulting equations are

[TABLE]

From these calculations, one can make two observations:

In the course of elimination one successively changes the diagonal elements $A_{jj}$ and the right-hand side $b_{j}$. 2. 2.

The exclusion schemes in figures 1a and 1b share the same computations. For example, terms A_{31}A_{13}\big{/}A_{11}, A_{32}A_{23}\big{/}A_{22} and A_{54}A_{45}\big{/}A_{44} appear on the way to equation (7b) as well as to (6b). It would be more advantageous to reuse the same computations, not to redo them each time one needs to eliminate a variable.

The first observation suggests that one can introduce corrections to the diagonal terms and $b_{j}$, that come from the elimination of variable $i$. For the sake of convenience, we denote them $\Lambda_{ji}$ and $\mu_{ji}\Lambda_{ji}$, respectively. For example, equation (6a) becomes

[TABLE]

Since corrections are the same for any order of elimination, to reuse them, we can regard $\Lambda_{ji}$ and $\mu_{ji}$ as a message that node $i$ sends to node $j$ along the edge of the graph. Once computed, these messages are in use in expressions like (8) and (9). To complete rewriting the elimination in terms of messages, one needs to introduce the rules to update messages when a new variable is excluded. To derive the rules, we rewrite equation (7b) using the definition of messages

[TABLE]

and use (8) to get

[TABLE]

It is easy to see that one needs to accumulate all messages from neighbors of $i$ except for $j$ to send the message from node $i$ to node $j$. Since the update rule includes only messages from the previous stages of elimination, one can iterate equations like (10) till convergence. And then, when all messages arrive, the solution can be read off as follows

[TABLE]

Note that (7) and (6) have exactly this form. Equations for the update of messages that we deduced in this section coincide with the GaBP update rules given by (19), which are derived from the probabilistic perspective below.

To summarize, the GaBP rules can be understood as a scheme that propagates messages on the graph, corresponding to the matrix of the linear system under consideration. These messages, namely $\Lambda_{ji}$ and $\mu_{ji}$, represent the corrections to the diagonal terms of matrix $\mathbf{A}$ and to the right-hand side $\mathbf{b}$, resulting from the elimination of variable $x_{i}$ from the $j$-th equation, $A_{jj}x_{j}+A_{ji}x_{i}+\dots=b_{j}$. Consistency, convergence, stopping criteria, and other practical matters are discussed in the following two sections.

2.2 Conventional belief propagation

Here we give a more traditional introduction to GaBP as a technique for statistical inference in graphical models. Following [5, 26], we reformulate (1) as an inference problem. For this purpose, consider a small subset of undirected graph models that are known as Gauss-Markov random fields. First, we define a pairwise Markov random field. The graph $\Gamma$ is the set of edges $\mathcal{E}$ and vertices $\mathcal{V}$. Each vertex $i$ corresponds to the random variable $x_{i}$ (discrete or continuous), and each edge corresponds to interactions between variables. The set of non-negative integrable functions $\left\{\phi_{i},\psi_{ij}\right\}$ together with the graph $\Gamma$ completely specifies the form of the probability density function of a pairwise Markov random field

[TABLE]

Here, $Z$ is a normalization constant (known in statistical physics as a partition function). The second equation in (12) (that is, the Gibbs distribution) should be considered as a definition of energy $E(\mathbf{x})$. The Gauss-Markov random field is a particular instance of a pairwise Markov model with a joint probability density function given by a multivariate normal distribution

[TABLE]

where $\mathbf{A}$ is a symmetric positive-definite matrix. The edges of $\Gamma$ correspond to the nonzero $A_{ij}$, and note that the splitting of the product in (12) into $\phi_{i}$ and $\psi_{ij}$ is not unique. A common task in the inference process is a computation of a partial distribution (or a marginalization)

[TABLE]

Integrals replace sums if $\mathbf{x}\in\mathbb{R}^{N}$. For the Gauss-Markov model, marginal distributions are known explicitly. For individual components of the vector $\mathbf{x}$, which is distributed according to (13), one can obtain distributions in closed form

[TABLE]

As the means of marginal distributions for the model (13) coincide with the elements of the solution vector for (1), methods from the domain of probabilistic inference can be applied directly to obtain the solution.

A popular algorithm that exploits the structure of the underlying graph to find the marginal distribution efficiently is Pearl’s belief propagation [20]. Pearl’s algorithm operates with local messages that spread from node to node along the graph edges, and beliefs (approximate or exact marginals) are computed as a normalized product of all incoming messages after the convergence. More precisely, belief propagation consists of (i) the message update rule

[TABLE]

where $m_{ij}$ is a message from node $i$ to node $j$ and $N(i)$ is the set of neighbors of the node $i$, and (ii) the formula for marginals

[TABLE]

Although, for continuous random variables the problem of marginalization and the algorithm of belief propagation are harder in general, it is not the case for the normal distribution. Namely, for the Gauss-Markov model, one can parameterize messages in the form of the normal distribution

[TABLE]

and explicitly derive update rules, means, and precision

[TABLE]

These update rules correspond to the flood schedule such that at the current iteration step, each node sends messages to all its neighbours based on messages received at the previous step. Equations for the mean and precision should be put to use only after saturation according to some criteria, for example $|\mu^{(n+1)}-\mu^{(n)}|\leq\text{tolerance}$, and the same for $\mathbf{\Lambda}$. Rules (19) are collectively known as Gaussian belief propagation.

Belief propagation was designed to give an exact answer if $\Gamma$ has no loops. In the presence of loops, the result appears to be approximate if delivered at all. In the case of GaBP, the situation is more optimistic. We briefly recall some useful facts about GaBP that we discuss later in more detail. If GaBP converges on the graph of arbitrary topology, the means are exact, but variances can be incorrect [33]. The best sufficient condition for convergence of the Gauss-Markov model with symmetric positive-definite matrix can be found in [17], we discuss it later in greater detail. The fixed point of GaBP is unique [13]. On the tree, GaBP is equivalent to the Gaussian elimination [22].

Many different schemes that extend belief propagation and GaBP have been developed [38], [9], [29], [18], [10]. Here, we are mainly interested in generalized belief propagation proposed in [38] and subsequently developed in [37], [36] ,[39]. This new algorithm is significantly more accurate [39, Fig. 15] than Pearl’s algorithm, but at the same time it can be computationally costly. In what follows, we show how to use the generalized belief propagation in the context of numerical linear algebra.

2.3 GaBP for a nonsymmetric linear system

As explained in Section 2.1, the GaBP rules can be understood as corrections to the right-hand side and the diagonal elements of the matrix under successive elimination of variables. It means that in principle, one can apply the rules to solve at least some nonsymmetric systems. However, there is a problem which is specific to nonsymmetric case. Namely, it is possible to have $A_{ij}=0$ and $A_{ji}\neq 0$. In this case, rules (19) lead to singularity as $A_{ij}$ appears in the denominator. Since parametrization of messages is not unique both from the elimination and probabilistic perspectives, it is possible to define new set of messages $\widetilde{\mathbf{\Lambda}}$ and $\mathbf{m}$ as follows

[TABLE]

Then update rules become

[TABLE]

Note that reparametrization (20) has a problem in that it is not one-to-one if $A_{ji}=0$. However, the quick look at the equations (6), (7) makes clear that indeed it is possible to define messages in that way. That is, if $A_{jk}=0$, one does not need to eliminate $x_{k}$ from the second equation so the message $\Lambda_{kj}$ is indeed zero.

For a given matrix $\mathbf{A}\in\mathbb{R}^{N\times N}$, we construct a graph with $N$ vertices $v\in\mathcal{V}$ corresponding to the variables $x_{1},\dots,x_{N}$ and the set of directed edges $\mathcal{E}$. The edge pointing from the vertex $j$ to the vertex $i$ belongs to the set of edges iff $A_{ij}\neq 0$, i.e. $A_{ij}\neq 0\Leftrightarrow e_{ji}\in\mathcal{E}$. This definition fixes the correspondence between directed graphs and nonsymmetric matrices and allows us to use GaBP (see Algorithm 1) beyond its usual domain of applicability.

Algorithm 1 is sequential, but can run in parallel after some modifications. The stopping criteria can be different, for example, it is possible to use different norms, or $\text{error}=\left\|\mathbf{x}^{(n+1)}-\mathbf{x}^{(n)}\right\|_{\infty}$, or

[TABLE]

Note that the update of $\widetilde{\mathbf{\Lambda}}$ decouples from the one for $\widetilde{\boldsymbol{\mu}}$. So it is possible to construct an algorithm that computes only messages $\widetilde{\boldsymbol{\Lambda}}$ and returns diagonal elements for the inverse matrix. Later, these messages can be used in the course of all successive iterations if one resorts to the error correction scheme. We discuss how the algorithm of this kind can be utilized to decrease the number of floating point operations in the context of a multigrid scheme.

One of the central results of the present work is that two classical theorems from GaBP theory, summarized below, can be readily established for nonsymmetric matrices.

Theorem 2.1**.**

If there is $N\in\mathbb{N}$ such that $\widetilde{\mu}^{(N+k)}_{e}=\widetilde{\mu}^{(N)}_{e}$, $\widetilde{\Lambda}^{(N+k)}_{e}=\widetilde{\Lambda}^{(N)}_{e}$ for all $e\in\mathcal{E}$ and for any $k\in\mathbb{N}$, then $\mu^{(N+k)}_{i}=\mu^{(N)}_{i}=\left(\mathbf{A}^{-1}\mathbf{b}\right)_{i}$.

The analogous result for symmetric matrices first appeared in [33]. In A, we show how to extend the proof for the nonsymmetric case.

Theorem 2.2**.**

If $A_{ii}\neq 0~{}\forall i$, $\widetilde{|R|}_{ij}=\left(1-\delta_{ij}\right)\frac{\left|A_{ij}\right|}{\left|A_{ii}\right|}$, and $\rho(\widetilde{|\mathbf{R}|})<1$, then the Algorithm 1 converges to the solution $\mathbf{x}^{\star}=\mathbf{A}^{-1}\mathbf{b}$ for any $\mathbf{b}$.

Sufficient condition for symmetric positive-definite matrices was established in [17]. Appendix B contains the proof with necessary modifications that holds for nonsymmetric matrices.

To make connections with the classical theory of iterative methods, we give another (less general) sufficient condition.

Corollary 2.0.1**.**

If $\mathbf{A}$ is the $M$-matrix (see [24, Definition 1.30, Theorem 1.31]), Algorithm 1 converges to the solution $\mathbf{x}^{\star}=\mathbf{A}^{-1}\mathbf{b}$ for any $\mathbf{b}$.

Proof.

For $M$-matrix $\rho(\mathbf{I}-\mathbf{D}^{-1}\mathbf{A})<1$, $A_{ij}\leq 0,~{}i\neq j$ and $A_{ii}>0$, where $\mathbf{D}$ is a diagonal of $\mathbf{A}$. It means that $\mathbf{\widetilde{R}}=\mathbf{I}-\mathbf{D}^{-1}\mathbf{A}=\left|\mathbf{\widetilde{R}}\right|$ and $\rho\left(\left|\mathbf{\widetilde{R}}\right|\right)<1$. ∎

3 Generalized GaBP solvers

3.1 Parent-to-child schedule

We present a particular version of the parent-to-child schedule from [39] applied to the pairwise Markov graphical models.

First, we define a region $r$ as a connected subgraph of the original graph $\Gamma$. Each vertex $a_{1}$ in the region graph is a region that may be connected by a directed edge with another vertex $a_{2}$ if $a_{2}\subset a_{1}$. The direction of the edge is from a larger region to smaller. If there is a directed edge from $a$ to $b$, then $a$ is a parent of $b$, and $b$ is a child of $a$. If the vertices $a$ and $b$ are connected by a directed path starting from $a$, then $a$ is an ancestor of $b$, and $b$ is a descendent of $a$. The set of all vertices of the factor graph is $\mathcal{R}$, the set of all edges is $\mathcal{E}_{\mathcal{R}}$, the set of all parents, children, ancestors, descendants of $a$ are $P(a)$, $C(a)$, $A(a)$, and $D(a)$, respectively. We supplement each region $r\in\mathcal{R}$ with a counting number,

[TABLE]

and require that each vertex $v\in\mathcal{V}$ and each edge $e\in\mathcal{E}$ of the original graph be counted exactly once,

[TABLE]

The definitions of a counting number (23) and condition (24) are justified in the framework of the cluster variation method [21] (or the Kikuchi approximation [16]). Equation (24) is a result of the Möbius inversion formula applied to the sum over partially ordered sets [3], and (24) can be proven using definitions of Möbius and Zeta functions [3, equation 16].

A sample region graph is shown in Figure 2. For example, by region $5689$, we mean all nodes and all links between them that are present on the original graph. It is easy to see that the counting number condition (24) is satisfied for all the links and nodes.

The last two definitions that we need are the shadow of the region $S(r)=D(r)\cup r$ and a Markov blanket of the region $B(r)=P\left(S(r)\right)\backslash S(r)$. An example of both the shadow and Markov blanket of region $5689$ can be found in Figure 2b.

A parent-to-child algorithm consists of three elements: (i) messages that propagate along the directed edges of the region graph

[TABLE]

where $x_{b}$ corresponds to variables belong to the cluster $b$, (ii) a formula for the cluster beliefs

[TABLE]

and (iii) message update rules that follow from the consistency conditions

[TABLE]

In [39, equation 114], one can find message update rules for the general situation, but for our simple region graph, conditions (27) suffice.

3.2 Two-layer generalized GaBP

In this section, we consider the simplest possible valid region graph that consists of two layers, Figure 3b. The large regions, the horizontal and vertical stripes in 3a, are parents of small regions presented by the individual nodes. To proceed, we need to establish some further notation. First, we define a projector on the region $k$,

[TABLE]

Here, $|j|=|k|$ and $|i|>|k|$ is chosen to be conformable depending on the context. Ordering is global, i.e., it is fixed for the whole graph and maintained the same way in all manipulations. We also introduce brackets,

[TABLE]

where $\mathbf{C}\in\mathbb{R}^{|k|\times|k|}$ and according to our notation, the sizes of the matrix and the vector (29b) depend of the context whereas in (29a) the size of the matrix is $|k|\times|k|$ and the size of the vector is $|k|$.

Let $\left\{L_{i}\right\}$ and $\left\{l_{i}\right\}$ be sets of large and small regions in the two-layer region graph. For messages, we use the following parameterizations

[TABLE]

According to (26), the belief of any region reads

[TABLE]

For the small region $l$, the shadow is $\mathcal{S}(l)=l$ and the Markov blanket is $\mathcal{B}(l)=P(l)$, whereas for the large region $L$, the shadow is $\mathcal{S}(L)=C(L)$ and the Markov blanket is $\mathcal{B}(L)=P(C(L))\backslash L$. Consistency condition (27) allows one to derive update rules for each message sent from the parent region to the child region,

[TABLE]

Using standard results for marginals of the normal distribution, we can deduce

[TABLE]

If region $L$ has many children and $|l|\ll|L|$, the direct use of equation (32c) is not efficient. Instead, we use the formula following from the Woodbury matrix identity,

[TABLE]

and split for each child region $l$,

[TABLE]

Then, for precision and mean of messages, we obtain

[TABLE]

Equations (36) are especially useful in the situation when the small regions consist of the single node, i.e. the case of graph in Figure 3. In this situation, one needs to invert the matrix corresponding to the large region only once whereas the direct application of (33) leads to $|C(a)|$ inversions.

From formulae (36), one can derive GaBP rules. As the large regions in the Bethe approximation consist of two vertices with the edge connecting them, we can rewrite messages as

[TABLE]

Then, the Gaussian belief propagation rules follow from

[TABLE]

The validity of the presented rules for nonsymmetric linear problems does not follow from the derivation above. Nevertheless, one can apply the generalized GaBP rules with no modifications to solve them, as we explain in the next sections.

3.3 Elimination perspective

First, we define a hypergraph $\mathcal{G}$ based on the region graph. The set of nodes coincides with the set of large regions, and each common child corresponds to the edge in the graph. The example of such a hypergraph is shown in Figure 3c. With this definition, messages from equations (32) can be redefined with no reference to the small region as long as a one-to-one correspondence between edges of $\mathcal{G}$ and child regions of the original region graph are established. Now we study the single message from region $j$ to region $i$ with $k=j\cap i$ and $\overline{j}=j\backslash k$. According to (32) and (33), the precision part of the message is

[TABLE]

and the mean is

[TABLE]

So the subset $k$ of $\mathbf{A}_{i}$ and $\mathbf{b}_{i}$ receive a correction from region $j$,

[TABLE]

The corrections above are the same as in the ordinary block LU decomposition,

[TABLE]

Thus, we have a correspondence between the block LU operations and the generalized GaBP for the two-layer region graph. As in the case of the regular GaBP, LU applies in the fashion of dynamic programming, i.e., one does not just solve a smaller subproblem as in the block iterative scheme, but rather forms a recursive procedure that decouples different blocks from each other. Again, one should reparametrize messages for them to be valid for the arbitrary invertible matrix $\mathbf{A}$. Note that if $\mathbf{A}_{jk}=0$, then $\mathbf{\Lambda}_{ji}$ is not invertible. However, this is not a problem because in this case $\mathbf{A}_{i}$ does not receive corrections $\Delta\left(\mathbf{\mathbf{A}}_{i}\right)_{kk}$, and the reparametrization for the mean messages $\mathbf{m}_{Ll}\equiv\mathbf{\Lambda}_{Ll}\boldsymbol{\mu}_{Ll}$ resolves the issue.

3.4 The algorithm

In the case of generalized GaBP, messages propagate on the region graph. To remind, $\left\{L_{i}\right\}$ is the set of large regions, and $\left\{l_{i}\right\}$ is the set of small regions.

The algorithm, as we describe it in this section, acts on a given region graph. We do not provide an algorithm or recommendations on how to build a region graph. Some observations about the influence of a particular choice of the large regions can be found in Section 5. Regarding complexity, for each region, one needs to solve a linear system of $\left|L\right|$ equations, and also find a submatrix of the inverse matrix of size $\left|L\right|\times\left|L\right|$. We do not specify how to do it. But generally, the former task is hard to accomplish asymptotically faster than the whole inverse, so the computational cost of the entire scheme depends mostly on this operation. For this reason, it what follows we mostly consider region graphs with small regions each of which contains a single node. In this situation, one can avoid fill-in and estimate only the diagonal of the inverse matrix which is usually more straightforward. Also, there is no need for an additional inverse step during the message update. Moreover, when all small regions are single nodes, the whole algorithm is either a method that speeds up GaBP or a particular schedule of GaBP depending on the chosen way of finding the inverse and the solution of a linear system.

To perform a worst-case analysis both of Algorithm 2 and update rules (36), for a given region with $N$ variables, we denote the number of children by $M$, the number of variables of each child region by $n_{i}$, the number of parents for each child by $p_{i}$, and use the LU decomposition to find an inverse. For the Algorithm 2, the number of operations is

[TABLE]

The first term in brackets is due to the inverse during the message update stage, the second term in brackets is due to the message update and message accumulation steps, the last term is from the inverse of a matrix for the large region. For update rules (36), we obtain

[TABLE]

Since

[TABLE]

one obtains a speed-up if $M>3$ for an arbitrary region graph. For certain regular partitions, for example the one in Figure 3a, $M$ scales like $N$, and in this cases, the Algorithm 2 performs $O(N^{3})$ operations whereas rules (36) perform $O(N^{4})$ operations for each large region.

It is easy to see that if the LU method is employed, the number of operations for the single sweep is $O(K)$, where $K$ is an overall number of variables, only if the number of variables, for some regular partition for which the limit makes sense, in each large region scales like $O(1)$. Clearly, LU is not the best option for all cases, for example, matrices can have a particular structure (tridiagonal as Figure 3a), or it may be more advantageous to use probing or other techniques of estimation of certain subblocks of the inverse matrix [30] in combination with some iterative scheme for the solution of linear system.

The Algorithm 2 can be justified theoretically on the basis of the following two theorems.

Theorem 3.1**.**

If there is $N\in\mathbb{N}$ such that $\mathbf{m}^{(N+k)}_{e}=\mathbf{m}^{(N)}_{e}$, $\mathbf{\Lambda}^{(N+k)}_{e}=\mathbf{\Lambda}^{(N)}_{e}$ for all $e\in\mathcal{E}_{\mathcal{R}}$ and for any $k\in\mathbb{N}$, then for each large region $\left[\mathbf{x}\right]_{L}\equiv\widetilde{\mathbf{\Lambda}}_{0}^{-1}\widetilde{\mathbf{b}}_{0}=\left[\mathbf{A}^{-1}\mathbf{b}\right]_{L}$ (see Algorithm 2 for details).

That is, the steady state of the message flow, if it exists, corresponds to the exact solution. The proof is given in C.

For the sufficient condition for convergence we need additional definitions. First, based on a given region graph $\left\{\mathcal{R},\mathcal{E}_{\mathcal{R}}\right\}$, we define a set of variable subsets

[TABLE]

where $l_{i}$ is a small region and $L_{j}$ is a large region. Using $F$, we form a partition of the matrix $\mathbf{A}$ and the right-hand-side vector $\mathbf{b}$

[TABLE]

where each diagonal block corresponds to the element of the set $F$. We also define

[TABLE]

Note that the second line in the preceding equation contains a definition of matrix $\left\|\widetilde{\mathbf{R}}\right\|$, which depends on the operator norm $\left\|\cdot\right\|$ (see [15, ch. 5, Definition 5.6.3]).

The following statement gives sufficient conditions for convergence.

Theorem 3.2**.**

If for matrix (47) which is based on partition (46) $\det\mathbf{A}_{ii}\neq 0~{}\forall i$ and $\rho\left(\left\|\widetilde{\mathbf{R}}\right\|\right)<1$ in some operator norm, then two-layer generalized GaBP (algorithm 2) converges to the exact solution $\mathbf{x}=\mathbf{A}^{-1}\mathbf{b}$.

The proof of this theorems appears in D. From the second part of the argument in D.2, one can deduce the following

Corollary 3.0.1**.**

Generalized GaBP (algorithm 2) converges whenever GaBP converges (Theorem 2.2) and all submatrices corresponding to the large blocks are invertible (see equations (46), (47)).

The opposite does not hold. For example, consider a matrix

[TABLE]

In this case, the spectral radius of matrix $\left|\mathbf{\widetilde{R}}\right|$ defined in Theorem 2.2 equals $\sim 1.03$ and GaBP diverges222Note that the divergence of GaBP does not follow from $\left|\mathbf{\widetilde{R}}\right|>1$ as Theorem 2.2 provides only sufficient conditions. For this particular case, the pathological behavior of GaBP follows from the numerical experiment (see [1] for details).. On the other hand, the spectral radius of $\left|\mathbf{\widetilde{R}}\right|$ defined by (48) and the partition given in (49) is smaller than one in $l_{\infty}$ and spectral norms [15, Examples 5.6.5, 5.6.6] (see [1] for further details).

4 GaBP as a smoother for the multigrid method

The most straightforward view on the geometric multigrid is to describe it as an acceleration scheme for classical iterative methods. For completeness, we briefly recall the main ideas.

The multigrid consists of four essential elements: a projection operator $\mathbf{I}_{V}^{V^{{}^{\prime}}}:V\rightarrow V^{{}^{\prime}}$ ($V$, $V^{{}^{\prime}}$ are linear spaces) that reduces the number of degrees of freedom, an interpolation operator $\mathbf{I}_{V^{{}^{\prime}}}^{V}:V^{{}^{\prime}}\rightarrow V$ that acts in the ”inverse” way, a smoothing operator $\mathbf{S}_{V}:V\rightarrow V$ which is usually a classical relaxation method, and a set of linear operators $\mathbf{A}_{V^{{}^{\prime}}}$ that approximate $\mathbf{A}$ on coarse spaces $V^{{}^{\prime}}$. What we describe next is a two-grid cycle.

• •

For the current approximation $\mathbf{x}^{n}$ of solutions of $\mathbf{A}\mathbf{x}=\mathbf{b}$, one performs several relaxation steps $\overline{\mathbf{x}}=\mathbf{S}^{\nu}\mathbf{x}^{n}$.

• •

Then, based on properties of $\mathbf{S}$, the linear space $V^{{}^{\prime}}$ and the projection operator $\mathbf{I}_{V}^{V^{{}^{\prime}}}:V\rightarrow V^{{}^{\prime}}$ are constructed. The purpose of this space is to represent the residual $\mathbf{r}=\mathbf{b}-\mathbf{A}\overline{\mathbf{x}}$ and an error $\mathbf{e}=\mathbf{x}_{\text{exact}}-\overline{\mathbf{x}}$ accurately using fewer degrees of freedom $\left|V^{{}^{\prime}}\right|<|V|$.

• •

Having the space $V^{{}^{\prime}}$, one constructs an operator $\mathbf{A}^{{}^{\prime}}$ that approximates $\mathbf{A}$ and solves the error equation $\mathbf{A}^{{}^{\prime}}\mathbf{e}^{{}^{\prime}}=\mathbf{I}_{V}^{V^{{}^{\prime}}}\mathbf{r}$.

• •

The error, after projection back to $V$, gives the next approximation to the exact solution, $\mathbf{x}^{n+1}=\mathbf{S}^{\mu}\left(\mathbf{x}^{n}+\mathbf{I}^{V}_{V^{{}^{\prime}}}\mathbf{e}^{{}^{\prime}}\right)$.

The multigrid utilizes a two-grid cycle to solve the error equation $\mathbf{A}^{{}^{\prime}}\mathbf{e}^{{}^{\prime}}=\mathbf{I}_{V}^{V^{{}^{\prime}}}\mathbf{r}$ itself. It produces the chain of spaces (grids in the geometric setup), projection operators that allow moving between them, and a set of approximate linear operators. For more details, we refer the reader to other resources: a simple introduction to geometric multigrid can be found in [25, ch. 13], for the algebraic multigrid a recent review [35], physical considerations about algebraic multigrid can be found in the introduction of [23], and among other books on the subject, [31] provides a comprehensive introduction for practitioners.

Here, we consider only linear systems of equations arising from finite difference discretization of elliptic equations with smooth coefficients in two space dimensions. In this case it is possible to use grids in place of linear spaces. Let the finest grid contain $2^{J}+1$ points $\mathcal{G}_{J}=\left\{(i+j)h|h=2^{-J},i,j=\overline{0,2^{J}}\right\}$, then the coarser grid $\mathcal{G}_{J-1}$ contains each second points along both directions. As we are working in the physical space, the restriction operator $\mathbf{I}_{V}^{V^{{}^{\prime}}}$ computes a weighted average of neighbouring points, operator $\mathbf{I}^{V}_{V^{{}^{\prime}}}$ performs interpolation, and $\mathbf{A}$ on the grid $\mathcal{G}_{J^{{}^{\prime}}}$ is a finite difference approximation of the differential operator. In this article we always use full weighting restriction [31, eq. 2.3.3] and bilinear interpolation [31, eq. 2.3.7].

The smoother should be a mapping $\mathbf{S}:\mathbf{x}^{n}\rightarrow\mathbf{x}^{n+1}$. Although GaBP is not of this form, one can use an error correction scheme as explained in Algorithm 3. In the next two subsections we analyze smoothing properties of Algorithm 3 and estimate its computational complexity.

4.1 Local Fourier Analysis

Local Fourier Analysis allows us to compute the spectral radius of the two-grid cycle, the smoothing factor of the relaxation scheme, and the error contraction in a chosen norm [31, ch. 4]. In this subsection, we apply the analysis to the central difference discretization of the Laplace equation in two spatial dimensions

[TABLE]

If after $\mu$ sweeps of Algorithm 1 the solution has a form $\mathbf{S}\mathbf{r}^{n}$, then the output of Algorithm 3 is $\mathbf{x}^{n+1}=\mathbf{S}\left(\mathbf{b}-\mathbf{A}\mathbf{x}^{n}\right)+\mathbf{x}^{n}$. Thus, for the error we have

[TABLE]

Now we need to find a stencil of the operator $\mathbf{S}$. For GaBP it differs for parallel and sequential versions. For one and two sweeps of parallel version on the infinite lattice we have

[TABLE]

It is easy to compute symbols

[TABLE]

As the smoothing factor is

[TABLE]

we can conclude that the parallel version of GaBP shows no smoothing properties. This conclusion is confirmed by our numerical experiments.

In the sequential case, one can deduce the form of $\mathbf{S}$ based on the elimination perspective. When one starts to move along the lattice, messages correspond to the elimination of variables, which means that

[TABLE]

Then, for the smoothing factor we have

[TABLE]

for $\theta_{2}=\frac{\pi}{2}$ and $\theta_{1}=2\arctan\frac{1}{3}$. This means that the smoothing factor for the sequential GaBP coincides with the one for sequential Gauss-Seidel iteration scheme [31, Example 4.3.4]. It is also clear that for the anisotropic problem

[TABLE]

both Gauss-Seidel scheme and sequential GaBP lose their smoothing properties when $\epsilon\rightarrow 0$. However, numerical experiments (figure 4) show that the convergence rate of GaBP does not depend on $\epsilon$. An explanation for this particular case is straightforward. For sufficiently small $\epsilon$ equations for each vertical line (i.e., in $y$ direction) are independent. GaBP is an exact solver for trees. The single multigrid iteration eliminates variables from $2+2=4$ neighbours in case 4a and from $3+3=6$ neighbours in case 4b. When messages cover the whole line of $2^{6}-1=63$ nodes, the system of linear equations for each vertical line is solved exactly. It gives $63/4\sim 15$ iterations for 4a, and $63/6\sim 10$ iterations for 4b.

The effect displayed in Figure 4 is a manifestation of the dynamic nature of GaBP. Even as part of the multigrid it maintains information about all previous iterations. More convergence histories can be found below, in the section with numerical examples. Overall, we conclude that sequential GaBP as part of the multigrid behaves similarly to Gauss-Seidel in the absence of anisotropy, but is substantially more robust in the presence of anisotropy. The behavior captured in Figure 4 also illustrates that Local Fourier Analysis is not an appropriate tool to analyze GaBP.

4.2 Reducing computational complexity

The number of floating point operations per iteration for algorithms 1 and 2 depends on the graph of the matrix $A$. Here, we consider the operator with the dense $9$ point stencil

[TABLE]

which can come from the second order finite difference approximation of a differential operator containing second and first derivatives. The same analysis for the 5 points stencil is straightforward. For convenience, we split Algorithm 1 (sequential version) into three parts:

• •

Accumulation stage. $\Sigma$ and $m$ are computed.

• •

Update stage. New messages $\widetilde{\Lambda}$ and $\widetilde{\mu}$ are constructed from the previous ones.

• •

Termination stage. The final answer $m/\Sigma$ is obtained.

We also neglect all effects from boundaries. Under these assumptions, the number of floating point operations for the single sweep GaBP is

[TABLE]

Here, we perform only a half of accumulation stage and a half of update stage, because we do not need to receive messages from nodes that we have not visited yet, nor we need to send messages to already visited nodes.

For lexicographical Gauss-Seidel scheme, the number of floating point operations is $\#_{\text{LEX GS}}=17N.$ It means that a single sweep of 1 takes slightly fewer floating point operations than three sweeps of Gauss-Seidel smoother $\#_{\text{GaBP}_{1}}\sim 3\#_{\text{LEX GS}}$. For $M$ sweeps of GaBP one has $\#_{\text{GaBP}_{M}}=N(64M-12)$. In the context of multigrid, it is important to have a cheap smoother, but even $52N$ is too expensive. However, it is possible to reduce computational complexity by precomputing all required messages $\widetilde{\Lambda}$, which depend only on the matrix $\mathbf{A}$ and not on the right-hand-side vector. For $M$ sweeps of GaBP with precomputed $\widetilde{\Lambda}$, we have $\#^{\widetilde{\Lambda}}_{\text{GaBP}_{M}}=N(24M+8)$. We summarize all these results regarding the complexity of GaBP in Table 1.

5 Numerical examples

In this section, we present numerical experiments with matrices that arise from second-order finite difference approximations of two-dimensional elliptic differential equations with $(x,y)\in\left[0,1\right]^{2}$. The grid is assumed to be uniform and consists of $2^{6}-1$ inner points along each direction. We use Dirichlet boundary conditions in all the examples. These conditions are not specified directly and should be extracted from the exact solution. In the same vein, the form of the source term $g(x,y)$ (ride-hand side) can be derived from the exact solution and is not given explicitly. In all the experiments, the stopping criterion is $\left\|r\right\|_{\infty}\leq 2\cdot 10^{-4}$. In the tables below, $N_{\text{it}}$ denotes the number of iterations and $N$ is the number of variables. Before we begin the main discussion, we summarize the main properties of the solvers that are used.

5.1 Note about solvers and smoothers

In addition to the number of floating point operations (FLOP) (see Table 1), an important characteristic of a solver is its degree of parallelism, which is provided for various solvers in the following table (for classical methods see e.g. [31, ch. 6]):

[TABLE]

Note that the line version of GaBP possesses a better degree of parallelism than GS versions. When we consider GaBP as a multigrid component, we always use V-cycle, bilinear restriction, and prolongation operators and LU as a coarse-grid solver. In our notation $V(J_{1},J_{2})$ means that the fine grid consists of $2^{J_{1}}-1$ points, the coarse grid of $2^{J_{1}-(J_{2}-1)}-1$ points; numbers $(n,m)$ before the smoother name refer to the number of pre- and post-smoothing steps.

5.2 GaBP as a stand-alone solver

Classical relaxation methods are rarely used outside the AMG (algebraic multigrid) or GMG (geometric multigrid) to solve linear systems. Nevertheless, we present an example of their performance below. As a linear problem we use the following elliptic boundary value problem:

[TABLE]

with $g(x,y)$ and boundary conditions chosen such that $\phi_{\text{exact}}=\cos(\pi x)\cos(\pi y)$ is the exact solution (this method of manufactured solutions is used in the remaining examples as well). The performance of various methods on this problem is shown in the following table:

[TABLE]

The last line in the table corresponds to the three sweeps of the error correction scheme with precomputed messages $\Lambda$ (see Section 4.2). We see that GaBP does not provide particular advantages over classical relaxation methods as a stand-alone solver, even though there are some reports of its excellent performance (see, e.g., [33], [5]). The main bottleneck here is the computational complexity of the scheme. To some extent, one can mitigate this problem by precomputing $\Lambda$ before the iteration process begins. Still, in this particular situation, both the 4-colors GS and Jacobi provide better alternatives due to their low cost and high degree of parallelism. We observe the analogous behavior for other elliptic problems.

5.3 GaBP as a multigrid smoother

As a rule, relaxation solvers become applicable to real large-scale problems and are competitive with projection methods only in the framework of multigrid schemes. We now present several situations that could potentially challenge state-of-the-art GMG smoothers. Note that in this section, we always use the error correction version of GaBP. Since we apply the same version of multigrid, we compute the FLOP score solely for the smoother and on the fine level only.

5.4 Large mixed derivative

The first equation of interest is of the form

[TABLE]

The main problem here is that for small $\epsilon$, the ellipticity is almost lost. Below one can see the table with the best in terms of FLOP $V$-cycle solver of each kind:

[TABLE]

We can see that the performance of GaBP is comparable with the 4-color GS smoother. So GaBP can be considered to be robust for the almost non-elliptic equations. We also stress that both 4-color GS and GaBP are preferable over the line smoothers for this problem because of their better degree of parallelism.

5.5 Boundary layers

The other practically relevant case that is a challenge for standard geometrical smoothers is an advection-diffusion problem in which advection dominates. We take as an example the following problem:

[TABLE]

As one can see, there are two boundary layers near $x=1$ and $y=1$, each of width $\sim\epsilon$. The solution is not large, but the derivative is \sim 1\big{/}\epsilon. The table below shows the performance results.

[TABLE]

We observe two interesting trends. First, for $\epsilon=0.02$, it is not possible to apply the red-black GS. However, it is still possible to construct a smoother from GaBP, using a sufficiently large number of sweeps. Second, the line GaBP smoother with a reasonable amount of steps performs nearly as well as the alternating-zebra GS. Still, since the GABP can process all lines simultaneously, we conclude that it can outperform classical geometric smoothers for such advection-dominated elliptic problems. Furthermore, if the precomputation of $\Lambda$ is affordable from the perspective of the additional storage required, it is far better to use the red-black GaBP smoother.

5.6 Inner layers

We consider another advection-dominated diffusion problem

[TABLE]

This equation differs from (62) in two respects: 1) the solution has two inner layers, and 2) they are not aligned with the coordinate lines. Now the performance is as follows.

[TABLE]

For this problem we can see the same pattern as for the previous example. The classical red-black solver is unable to smooth the error, whereas the GaBP-based color iteration scheme works fine. Additionally, due to its excellent degree of parallelism, the red-black GaBP significantly outperforms the alternating-zebra GS smoother.

5.7 Stretched grid

Another situation of practical interest is given by the following problem:

[TABLE]

To understand this problem, consider the finite difference discretization of the Laplace equation on a grid which is highly concentrated near the edges of the domain. If one denotes $a_{\text{in}}$ and $a_{\text{bn}}$ to be characteristic scales of the coefficients, related to inner and boundary points respectively, the ratio a_{\text{bn}}\big{/}a_{\text{in}} will be large. We achieve the same effect in equation (64) on the uniform grid by multiplying second derivatives by positive terms of the form \left(1+\left(\left(x-\frac{1}{2}\right)^{2}+\eta\right)^{p}\big{/}\epsilon\right) which are approximately equal to $1$ inside the domain, but grow rapidly to large values near the boundaries. The table below shows thee results.

[TABLE]

In this table, the first set of parameters p,\eta,\epsilon=20,1\big{/}2,10^{-6} corresponds to the linear stretching of the grid by a factor of $\sim 40$ and the second p,\eta,\epsilon=20,1\big{/}2,8\cdot 10^{-8} to the linear stretching by a factor of $\sim 160$. We conclude that there is a version of the red-black GaBP with excellent convergence rate and good degree of parallelism. The same is true for the line GaBP smoother. Both of them can be used as an alternative to the classical alternating-zebra smoother.

5.8 Comparison with a projection method

For the sake of completeness, we also give an example of the performance of BiCGSTAB and GaBP-based multigrid for problem (60).

[TABLE]

As one may have anticipated, the projection method without a suitable preconditioner cannot outperform the geometric multigrid.

Overall, based on the presented result, we conclude that different versions of GaBP perform either comparably well or better than the state-of-the-art smoothers for GMG. The main disadvantage of GaBP is its computational complexity. Even though with a precomputed $\Lambda$, one can substantially decrease the number of FLOPs, the cost of a single iteration is still higher than for the classical smoothers. However, its clear advantages are the robustness and the degree of parallelism. The former allows one to construct new point-based relaxation schemes for situations where classical point-based relaxation methods fail. And the latter enables the line GaBP to outperform the alternating-zebra GS smoother.

6 Conclusions

In this paper, we have introduced a new class of solvers for linear systems that are based on the generalized belief propagation algorithm. The solvers work for both symmetric and nonsymmetric matrices. We show how to reduce the complexity of the resulting algorithm in comparison to that of the straightforward application of the generalized belief propagation. A clear connection between the block LU decomposition and the new algorithm is established. Existing proofs for symmetric systems are generalized to nonsymmetric systems, and two new proofs for a block version of the GaBP are given. Furthermore, we show how to use the geometric multigrid to accelerate the GaBP, which with a precomputed $\widetilde{\Lambda}$ results in a robust solver with the same computational complexity as the one based on the Gauss-Seidel smoother.

We have demonstrated the performance of the new algorithm with several examples of boundary-value problems of varying complexity. The numerical experiments show that the GaBP is a competitive alternative to classical relaxation schemes. The reason for the good performance of GaBP is that it retains some information about all the previous stages, whereas the Gauss-Seidel, Jacobi, and Richardson solvers do not. Even though large computational overhead is a disadvantage of GaBP, the problem can be alleviated within the framework of the multigrid scheme at the expense of additional storage and precomputation of some messages. Moreover, as part of the multigrid, the GaBP not only smooths high-frequency components of the error, but also effectively decreases the low-frequencies. Our numerical experiments indicate that this feature promotes additional robustness. For example, the convergence rate of GaBP for anisotropic model elliptic problem given by $\left(\partial_{x}^{2}+\epsilon\partial_{y}^{2}\right)u(x,y)=0$ does not depend on $\epsilon$, which is a somewhat unexpected result. Moreover, in the case of sharp inner and boundary layers, stretched grids, and large mixed derivatives, GaBP retains smoothing properties and performs better than the state-of-the-art geometrical smoothers.

The generalized GaBP introduced in the present work is in some sense a block version of the regular GaBP. Our numerical experiments show that the convergence rate increases with the size of the blocks such that after certain scale the generalized GaBP can compete with Krylov subspace methods (see a numerical example at [1]). However, the computational cost increases as well. In practical applications, one should balance these two tendencies to construct an optimal solver. Some considerations about the choice of the blocks can be found in [34], [37], but the issue is not yet fully resolved.

The generalized GaBP and GaBP as its particular case come from a domain of variational inference. To the best of our knowledge, there is currently no systematic analysis of the general relationship between deterministic/probabilistic inference and linear solvers. In our opinion, such a link may be useful for new interpretation of known techniques and provide insights that may lead to more efficient algorithms for numerical linear algebra.

Acknowledgement

The author is indebted to Dr. Aslan Kasimov for valuable suggestions.

Appendix A Consistency of GaBP

Here, following [33] we prove See 2.1 That is, if there is a steady state under mapping (21), the solution given by the GaBP rules is exact.

We note that the proof in [33] also holds for the nonsymmetric case. We present a slightly different version of their reasoning, without referencing graphical models for normal distribution.

The first concept that we need is a computation tree, which captures the order of operations under the GaBP iteration scheme. The computation tree contains copies of vertices and edges of the graph corresponding to $\mathbf{A}$. The matrix is supposed to be fixed so the computation tree depends on the root node $i$ and the number of steps $n$. We denote it by $T_{n}(x_{i})$. To obtain $T_{n}(x_{i})$ from $T_{n-1}(x_{i})$, we consider each vertex $m\in\mathcal{V}_{T_{n-1}(x_{i})}$ that has no incidence edges, find the corresponding variable on the graph of $A$, add to $T_{n-1}(x_{i})$ copies of each neighbour $k$ of $m$ such that $e_{km}\in\mathcal{E}_{A}$ except for $l$ for which $e_{ml}\in\mathcal{E}_{T_{n-1}(x_{i})}$. The example of the tree $T_{3}(x_{1})$ is in figure 5c, the $T_{2}(x_{1})$ in the dashed box exemplifies the recursion process.

By the construction of the computation tree, the following proposition is true.

Proposition A.1**.**

If $x_{i}^{(n)}$ is the solution on the $n$-th step of the algorithm 1, then it coincides with the one obtained after the elimination of all variables but $x_{i}$ (the root) from the computation tree $T_{n}(x_{i})$.

To relate the matrix $\mathbf{B}$ of the computation tree $T_{n}(x_{i})$ to the matrix $\mathbf{A}$, we define the matrix $\mathbf{O}$ [33, eq. 15] that connects original variables with copies

[TABLE]

or, more precisely, $y_{j}$ is a copy of $x_{i}\Rightarrow O_{ji}=1$ and $\sum\limits_{i}O_{ji}=1$. For example, matrix $\mathbf{O}$ for the tree in figure 5c is

[TABLE]

Now it is not hard to establish the connection between $\mathbf{B}$ and $\mathbf{A}$ [33, eq. 17]

[TABLE]

where $\mathbf{E}$ is nonzero only for the subset of variables that $n$ steps away from the root node on the computation tree $T_{n}(x_{i})$. The final part of the proof depends on the following statement [33, Periodic beliefs lemma].

Proposition A.2**.**

If there is $N\in\mathbb{N}$ such that $\widetilde{\mu}^{(N+k)}_{e}=\widetilde{\mu}^{(N)}_{e}$, $\widetilde{\Lambda}^{(N+k)}_{e}=\widetilde{\Lambda}^{(N)}_{e}$ for all $e\in\mathcal{E}$ and for any $k\in\mathbb{N}$, then it is possible to construct an arbitrary large computation tree $T_{M}(x_{i})$ for any root node $x_{i}$ such that $\mathbf{O}\boldsymbol{\mu}^{(N)}=\widetilde{\mathbf{B}}^{-1}\widetilde{\mathbf{d}}$. Where $\widetilde{B}_{ij}\neq B_{ij}$ and $\widetilde{d}_{i}\neq d_{i}$ only for $i=j$ that are $M$ steps away from the root node.

The crucial part here is that not only the solution for the root node coincides with the steady state solution of GaBP, but also the same is true for all the variables on the modified computation tree.

The proof is as follows. First, following the recursion procedure, we construct a computation tree of desired depth $M$. Then we continue to grow the tree till the subtrees of nodes $M$ steps away from the root reach the depth $N$ which corresponds to the steady state of GaBP. Now, elimination of subtrees results in the desired modified tree with the matrix $\widetilde{\mathbf{B}}$ and the right-hand side $\widetilde{\mathbf{d}}$.

Since we can construct an arbitrary modified computation tree, we can always get for arbitrary large $M$

[TABLE]

And we know that by construction of the modified computation tree

[TABLE]

So we conclude that

[TABLE]

Note, that $\mathbf{O}^{T}\mathbf{O}$ is a diagonal matrix that counts the number of copies of each variable, therefore we can always choose $M$ large enough to make $\det\left(\mathbf{O}^{T}\mathbf{O}\right)\neq 0$ and $\mathbf{A}\boldsymbol{\mu}^{(N)}=\mathbf{b}$ which means that the iterative scheme defined by the algorithm 1 is consistent.

Appendix B Convergence of GaBP

Here we present the version of the proof from [17] that extends to nonsymmetric matrices. Our modifications are relatively minor, but for the sake of logical coherence, we reproduce here the minimal set of arguments from [17] tuning definitions and proposition when needed. The main result of this section is

See 2.2

The whole idea of the proof [17] is to relate GaBP operations with the recursive update of the weights of walks on the graph, corresponding to the matrix $\mathbf{A}$. For the start, we define a walk $w$ as a an ordered set of vertices $w=\left(i_{1},i_{2},\dots,i_{l(w)}\right)$ where $l(w)$ is a length of the walk $w$ and $\forall k<l(w)\Rightarrow e_{i_{k}i_{k+1}}\in\mathcal{E}$. Each walk possesses a weight

[TABLE]

Note that the order is backward, which is a consequence of our definition of the directed graph. For the symmetric matrix, when the order is not essential, the equation (71) coincides with the weight defined in [17, 3.1]. Now, if we have a system $\mathbf{A}\mathbf{x}=\mathbf{b}$, we can rescale it using $\widetilde{b}_{j}=b_{j}/A_{jj}$. This procedure is valid for any $\mathbf{A}$ with nonzero diagonal and results in the equivalent system

[TABLE]

It is possible to represent the solution of (72) in the form of Neumann series (see [15, ch. 5]) because $\rho\left(\widetilde{\mathbf{R}}\right)<1$ and therefore

[TABLE]

However, for being able to rearrange terms in the sum as necessary, which is sufficient to rewrite the inverse matrix using walks, one needs to require absolute convergence which is $\rho(\widetilde{|\mathbf{R}|})<1$. Having this condition it is not hard to prove [17, Proposition 1, Proposition 5]

Proposition B.1**.**

If $\rho(\widetilde{|\mathbf{R}|})<1$, then $\widetilde{A}^{-1}_{ij}=\sum\limits_{w:j\rightarrow i}\phi(w)$ and $x^{\star}_{i}\equiv\left(\widetilde{\mathbf{A}}^{-1}\widetilde{\mathbf{b}}\right)_{i}=\sum\limits_{k\in\mathcal{V}}\sum\limits_{w:k\rightarrow i}\phi(w)\widetilde{b}_{k}.$

Here, by $w:j\rightarrow i$ we mean the set of walks which start from the vertex $j$ and end at the vertex $i$. If one defines [17, 3.2] sets of single-visit $k\overset{\backslash i}{\rightarrow}i$ and single-revisit $i\overset{\backslash i}{\rightarrow}i$ walks by all walks which are not visiting the node $i$ in between given start and end points, the sum over walks can be decomposed [17, eq. 12, 13; Proposition 9]

[TABLE]

The decomposition follows from ”topological” considerations alone which depend only on the structure of walks and not on the particular definition of the weight. The last result that we need is [17, Lemma 18]

Proposition B.2**.**

For each finite length walk $k\rightarrow j$ on directed graph of the matrix $A$ there is $n$ and unique walk on the computation tree $T_{n}(x_{i})$.

Now, if we can relate update rules (21) with the recursive structure of walks on a tree, the proof of the proposition 2.2 is done.

On the tree, for each vertex $i$, the sum over single-revisit walks splits into sums over subtrees $T_{k\cup i}$, which are maximal connected parts that contain $i$ and among $N(i)$, only $k$. Then

[TABLE]

but the sums over subtrees $T_{k\cup i}$ can be written as a sum over $T_{k\backslash i}\equiv T_{k\cup i}\backslash\left\{i\right\}$,

[TABLE]

where we used [17, eq. 12, Proposition 9]

[TABLE]

Using the definition of $\widetilde{\mathbf{R}}$, it is easy to see that if one denotes

[TABLE]

then the update rule (76) coincides with the one for $\widetilde{\mathbf{\Lambda}}$ in (21). Note that (78) is well defined because if $A_{ki}=0$, there is no contribution from this particular subtree, and we do not need to use the walk from there. In the same vein, the sum in the numerator of (74) can be decomposed

[TABLE]

Again, using the sum over subtrees $T_{k\backslash i}$

[TABLE]

decomposition on single-visit walks [17, eq. 13] and equations (78), (76), we obtain

[TABLE]

The parameterization introduced in (81) leads to the same update rule for $\widetilde{\boldsymbol{\mu}}$ as in (21). With that, the sufficient condition, given in proposition 2.2, is established.

Appendix C Consistency of generalized GaBP

Here we prove that the two-layer generalized GaBP is consistent. See 3.1 The idea of the proof is the same as for the regular GaBP. One needs to relate, considering the operations of generalized GaBP, equations that the generalized GaBP solves during the $N$-th step, with the original system of linear equations, and then to show that those systems coincide for a sufficiently large $N$ if steady state exists.

To do so, we introduce a flat version of the region graph (an example is shown in figure 6b) that provides less detailed information about parent-child structure. The flat region graph is an undirected graph $\left\{\mathcal{V},\mathcal{E}\right\}$, where $\mathcal{V}$ is the set of large regions and $\left(L,L^{{}^{\prime}}\right)\in\mathcal{E}$ if $L$ and $L^{{}^{\prime}}$ has at least one common child (the example is in figure 6b).

Now one can introduce the computation tree exactly in the same way as for GaBP. The only difference is that, because of an overlap between large regions, when we add a leaf node, we include overlapping variables to the root node. An example of the computation tree $T_{3}(B)$ as well as the $T_{2}(B)$ is in figure 6c. By construction of the computation tree, we know that the following is true.

Proposition C.1**.**

Elimination of all the variables on the computation tree $T_{N}(B)$ leads to the solution $\mathbf{x}_{B}$ that coincides with the one on the $N$-th step of generalized GaBP.

The relation between the matrix $\mathbf{B}$, corresponding to the computation tree, and the original matrix $\mathbf{A}$ is the same as in the equation (67) if one introduces the matrix $\mathbf{O}$

[TABLE]

Here, $\mathbf{x}$ are variables on the graph of matrix $\mathbf{A}$, and $\mathbf{y}$ are the ones on the computation tree.

Having the same relation between $\mathbf{A}$ and $\mathbf{B}$, we can repeat the rest of the proof, using the same arguments as in Section A. So it follows that generalized GaBP is consistent and proposition 3.1 is true.

Appendix D Convergence of generalized GaBP

In this section, we present a sufficient condition for the convergence of the two-layer generalized GaBP.

See 3.2

The proof consists of two parts. In the first one, we show that single-visit and single-revisit walks on a tree possess the same update rules as generalized GaBP messages. In the second part, we show that it is always possible to reorganize walks on the graph coming from the partition $F$ (see equations (46) and (47)) to restore each walk on a computation tree.

D.1 Walk structure on a tree

To complete the first part, we define for a given partition $F$ (equation (46)) of a matrix $\mathbf{A}$ the weight of a walk $w=\left(i_{1}i_{2}\dots i_{L}\right)$ by the product of matrices

[TABLE]

In the view of the standard result [2, ch. 8, Theorem 8.9] on absolute convergence in complete finite metric spaces it is possible to rearrange terms of the sum, such that we can formulate the following statement.

Proposition D.1**.**

If $\rho\left(\left\|\widetilde{\mathbf{R}}\right\|\right)<1$, then $\left(\widetilde{\mathbf{A}}^{-1}\right)_{ii}=\sum\limits_{n=0}^{\infty}\left(\widetilde{\mathbf{R}}^{n}\right)_{ii}=\sum\limits_{w:i\rightarrow i}\phi(w)$, $\mathbf{x}_{i}\equiv\sum\limits_{j\in\mathcal{V}}\left(\widetilde{\mathbf{A}}^{-1}\right)_{ij}\widetilde{\mathbf{b}}_{j}=\sum\limits_{j\in\mathcal{V}}\sum\limits_{w:j\rightarrow i}\phi(w)\widetilde{\mathbf{b}}_{j}$.

Here we used the same definition for the set of walks as in the Appendix B. Again, [17, eq. 12, 13] allows us to rewrite the diagonal blocks of the inverse matrix and the solution vector using single-visit and single-revisit walks

[TABLE]

On the tree we can split the sums over contributions from subtrees $T_{k\cup i}$ for each $k\in N(i)$. Therefore, from comparison with algorithm 2, we can deduce that

[TABLE]

Messages in algorithm 2 propagate along edges of the region graph, whereas messages that we have just defined flow along edges of a graph of the matrix $\widetilde{\mathbf{R}}$. To have a more straightforward connection between them, we consider $\widetilde{\mathbf{R}}$ as a matrix originates from the computation tree itself. Under this set of circumstances, there is a one-to-one correspondence between messages (85) and the ones in algorithm 2.

For single-revisit walks, one has

[TABLE]

If we consider the following matrix

[TABLE]

one can immediately see that update rules (86) indeed coincide with (33).

For single-visit walks, we have

[TABLE]

or using (85), we get

[TABLE]

Since $\mathbf{m}_{ki}=\mathbf{\Lambda}_{ki}\boldsymbol{\mu}_{ki}$ and

[TABLE]

we recover update rules (33). So we conclude that on the computation tree update rules for the two-layer generalized GaBP coincide with recursive relations for the single-visit and single-revisit walks.

D.2 Walk-sums and the graph refinement

The second part of the proof establishes the connection between sets of walks on the graph of the matrix $\widetilde{\mathbf{R}}$ and walks on the computation tree. First, for the matrix (47) we split a single region $i$ into two parts $i_{1}$ and $i_{2}$

[TABLE]

The transformation of the graph is in figure 7. We refer to this procedure as to the elementary refinement of the region $i$. From the construction of the refined matrix $\mathbf{A}$, the following proposition holds.

Proposition D.2**.**

There is a one-to-one correspondence between walks on the graph of $\widetilde{\mathbf{R}}$ and the one obtained by the elementary refinement of the region $i$ excluding three situations: 1) walk crosses $i$, 2) walk ends at $i$, 3) walk starts at $i$.

We discuss each of these situations separately. First, we need to introduce a new notation. Let $k\overset{M}{\longrightarrow}l$ be the set of walks, where each walk starts from $k$, ends at $l$ and newer leaves the subset $M$. It is easy to see that on the refined graph

[TABLE]

• •

Walk on $\widetilde{\mathbf{R}}$ that crosses $i$ has a form $w_{\text{cross}}=\left(\dots jik\dots\right)$ (see figure 7a). The weight of this walk is

[TABLE]

On the refined graph we can consider the set of all walks that coincides with $w$ outside $i$. The sum of weight of all these walks is

[TABLE]

We see that due to equation (92), weights are the same.

• •

Walk on $\widetilde{\mathbf{R}}$ that ends at $i$ has a form $w_{\text{end}}=\left(\dots ji\right)$ and a weight

[TABLE]

On the refined graph we have two set of walks

[TABLE]

that can be combined to have the same weight. Namely, using (92) we find that

[TABLE]

• •

Walk on $\widetilde{\mathbf{R}}$ that starts at $i$ has a form $w_{\text{start}}=\left(ij\dots\right)$ and a weight

[TABLE]

It is possible to relate this walk to two sets of walks $w_{1}=(i_{1}j\dots)$, $w_{2}=(i_{2}j\dots)$ on the refined graph multiplying by the corresponding inverse matrices

[TABLE]

where $l=\left\{i_{1},i_{2}\right\}$. The re-weight is needed because the original linear system and the refined one are multiplied by different block diagonal matrices and have different inverses.

We know the following two propositions to be true.

Proposition D.3**.**

Any computation tree can be, by the set of elementary refinements, turned to a computation tree of GaBP under a proper schedule (see discision before [17, Lemma 18]) operating on the graph of the matrix (47) partitioned according to $F$.

Proposition D.4**.**

For each walk on the graph of the matrix (47), there is a unique walk on a sufficiently large computation tree formed by a proper schedule.

Hence for each walk on the computation tree, it is always possible to find a unique set of walks on the graph of the matrix (47) that has the same weight after the multiplication by an appropriate inverse matrix (see 99). It allows us to conclude that if it is possible to define a walk-sum for matrix (47) (see proposition D.1), walk-sum on the computation tree converges too, so the proposition 3.2 is proven.