---
title: Parallel Banding Algorithm Plus (PBA+)
url: https://www.emergentmind.com/topics/parallel-banding-algorithm-plus-pba
type: topic
---

# Parallel Banding Algorithm Plus (PBA+)

The Parallel Banding Algorithm Plus (PBA+) is an $O(n)$ direct solver for large, sparse banded linear systems $A x = b$, where $A$ is an $n \times n$ matrix with nonzero entries confined within a band about the diagonal, and possibly with a small number $T_{\max}$ of extraband nonzero entries. PBA+ extends the core Parallel Banding Algorithm (PBA) framework by addressing these additional extraband elements without compromising linear scaling, and is designed for high accuracy and numerical efficiency, as required in computational physics and chemistry, such as finite-difference discretizations and holonomic constraint problems [1009.0938].

## 1. Problem Definition and Matrix Structure

PBA+ targets systems where the matrix $A$ exhibits upper ($m_u$) and lower ($m_\ell$) half-bandwidths, such that
$$
A_{i,j} = 0 \qquad \text{whenever}\ j - i > m_u \ \text{or}\ i - j > m_\ell.
$$
Extra-band nonzero entries are incorporated as
$$
A' = A + \sum_{t=1}^{T_{\max}} \delta_{i,R_t} \delta_{j,S_t} \Delta_t,
$$
where each pair $(R_t, S_t)$ lies outside the original band, extending applicability to cases dominated by band structure with sparse, but non-negligible, extraband couplings [1009.0938].

Key parameter notation includes:
- $b = m_u$, $\ell = m_\ell$,
- $\mu' = \min(b, \ell)$,
- $\mu_{i,j} = \min\{ b + (i-j),\, \ell - (i-j) \}$.

## 2. Core Algorithmic Principle: QP-Decomposition and Recurrences

PBA+ exploits a unique decomposition, writing
$$
Q A P = I,
$$
where $Q$ is lower-triangular with unit diagonal and $P$ is upper-triangular. The inverse is $A^{-1} = P Q$. This structure leads to two substitution phases:
$$
c = Q b, \qquad x = P c.
$$
Constructive recursion builds $P = P_1 \cdots P_n$, $Q = Q_n \cdots Q_1$; each $P_k$ operates on row $k$, $Q_k$ on column $k$, under explicit recursive updates for the coefficients $\xi_{i,j}$ and $\chi_{i,j}$:
- For $j \geq i$ (diagonal and super-diagonal):
  $$
  \xi_{i,i} = \left(A_{i,i} - \sum_{m = \max(1, i - \mu')}^{i-1} \chi_{i,m} \xi_{m,i} \right)^{-1}
  $$
  $$
  \xi_{i,j}\ (j>i) = \xi_{i,i}\left(-A_{i,j} + \sum_{m = \max(1, i-\mu_{i,j})}^{i-1} \chi_{i,m} \xi_{m,j}\right)
  $$
- For $i > j$ (sub-diagonal):
  $$
  \chi_{i,j} = -A_{i,j} + \sum_{m = \max(1, j-\mu_{i,j})}^{j-1} \chi_{i,m} \xi_{m,j}
  $$

In symmetric cases ($A = A^T$), the sub-diagonal recursion can be avoided by setting $\chi_{i,j} = \xi_{j,i} / \xi_{j,j}$.

After all $\xi_{i,j}, \chi_{i,j}$ are assembled, the solution proceeds via:
- Forward substitution for $c$:
  $$
  c_i = b_i + \sum_{m = \max(1, i - \ell)}^{i-1} \xi_{i,m} c_m
  $$
- Backward substitution for $x$:
  $$
  x_i = \xi_{i,i} c_i + \sum_{k = i+1}^{\min(n, i+b)} \xi_{i,k} x_k
  $$
Each phase costs $O(n b \ell)$, linear in $n$ for fixed bandwidths [1009.0938].

## 3. Extension to Extra-Band Nonzeros: The Plus Variant

PBA+ generalizes the core banded approach by extending the recurrences for $\xi, \chi$ to all columns/rows occupied by extraband entries $(R_t, S_t)$. At each such location, the same analytic recurrences are applied—skipping zeros as appropriate—thus retaining linear scaling with a supplementary complexity term $O(k m)$, where $k$ is the number of extra-band entries and $m$ is a characteristic band width. This modular design facilitates accurate direct inversion for matrices that are banded plus a sparse pattern of extraband couplings without resorting to iterative refinement or full matrix storage [1009.0938].

## 4. Parallelization Strategies

The core recursions couple only $b$ or $\ell$ preceding rows or columns, enabling efficient parallel execution via two main approaches:

**a) Shared-memory (wavefront scheduling, Meurant '85 style):**
- At each diagonal step $d$ ($0\leq d \leq \max(b,\ell)$), all elements $\xi_{i,i+d}$ (or $\chi_{i,i-d}$) for $i=1..n-d$ can be processed in parallel.
- Both forward and backward substitutions advance along anti-/diagonals, leading to a critical-path of $O(n)$, with parallel work $O(n)$ per step.
- On $p$ cores: $O(n b \ell/p + \mathrm{sync}\cdot\max(b,\ell)) \approx O(n/p + b)$ [1009.0938].

**b) Distributed-memory (block-strip decomposition, Polizzi '06 style):**
- The matrix is partitioned into $P$ contiguous blocks, each handling local $\xi, \chi$ evaluation up to block edges, and exchanging $O(b+\ell)$ boundary values.
- Each block performs a local factorization, then all blocks collectively solve a small coupling system of size $O(P(b+\ell))$.
- Final forward and backward solves require interface data broadcasts.
- Total time: $O((n/P) b \ell + P(b+\ell)^3 + (n/P)(b+\ell))$; optimal for $P \ll n/b$ [1009.0938].

Near-linear speedup is achieved so long as $b, \ell \ll n$ and the number of processors satisfies $P \ll n/b$.

## 5. Computation of the Full Inverse

The factors $P, Q$ facilitate retrieval of the entire inverse $A^{-1}$ in $O(n^2)$ operations. Exploiting $A^{-1}=P Q$ and the banded nature, the entries are computed as:
$$
(A^{-1})_{i,j} =
\begin{cases}
  \sum_{\ell = i+1}^{\min(i+b, n)} \xi_{i,\ell} (A^{-1})_{\ell,j} & i < j \\
  \xi_{i,i} + \sum_{\ell = i+1}^{\min(i+b, n)} \xi_{i,\ell} (A^{-1})_{\ell, i} & i=j \\
  \sum_{\ell = j+1}^{\min(j+\ell, n)} (A^{-1})_{i,\ell} \xi_{\ell,j} & i>j
\end{cases}
$$
Each loop over the matrix involves $O(n b)$ operations, resulting in $O(n^2)$ total complexity—competitive with dense inversion but leveraging full sparsity for memory efficiency [1009.0938].

## 6. Complexity, Numerical Performance, and Comparative Table

PBA+ achieves direct $O(n b \ell)$ scaling for banded systems, maintaining leading constants that are empirically observed to be smaller than those for standard banded Gaussian elimination. For banded $+$ sparse cases, the complexity is $O(n b \ell + k m)$, with $k$ as the number of extraband nonzeros. Storage requirements remain at $O(n(b+\ell) + k)$. Forming the full inverse is $O(n^2)$ in both operations and storage.

| Algorithm                    | FLOPs            | Storage           | Remarks                              |
|------------------------------|------------------|-------------------|--------------------------------------|
| Gaussian elimination (banded)| $O(n b \ell)$    | $O(n (b+\ell))$   | pivoting $O(n b)$ overhead           |
| PBA                          | $O(n b \ell)$    | $O(n (b+\ell))$   | no LU storage; smaller constant      |
| PBA+                         | $O(n b \ell+k m)$| $O(n (b+\ell)+k)$ | $k$ = extra-band nonzeros            |
| PBA full inverse             | $O(n^2)$         | $O(n^2)$          |                                      |

In extensive tests (up to $n=10^5$, $b=\ell=m\leq 300$) on random banded matrices:
- With partial pivoting, PBA achieves $\mathrm{Error} \sim 10^{-12}$–$10^{-10}$, nearly $n$-independent, with the error growing as $m^{1.4}$; Gaussian elimination yields $10^{-10}$–$10^{-8}$.
- PBA runtime is $0.6\times$ that of Gaussian elimination for $m\approx 10$, and up to $1.8\times$ faster for $m\approx 30$. Without pivoting, both algorithms deteriorate to $\sim10^{-9}$ error, but PBA remains $\sim2\times$ faster.

For 1D Poisson problems ($n=10^6$, $m=1$), PBA completes in $\approx0.02$ s with error $10^{-13}$, compared to $0.10$ s for banded Gaussian elimination (Numerical Recipes), for a speedup of $5\times$.

In Lagrange multiplier enforcement for protein chains ($m=6$, $n\approx10R+2$, $R\leq1000$), PBA is $1.6\times$ faster than Gaussian elimination, and a symmetry-exploiting PBA+ variant further improves this to $2.7\times$, achieving errors down to $10^{-12}$ [1009.0938].

## 7. Implementation Considerations and Scalability

PBA+ relies on contiguous diagonal-oriented storage formats for bands, supporting cache-efficient sweeps of all $O(n (b+\ell))$ elements needed. Extra-band entries are stored separately and invoked only in the extended recurrences. Partial pivoting, when adopted, incurs $O(n b)$ overhead and guards robustness against small $A_{i,i}$ elements.

Wavefront parallelism is effective for shared memory platforms as long as $b+\ell \ll n$; distributed memory strategies only require $O(b+\ell)$ words exchanged at block interfaces, facilitating scalability for large $n$. In typical physical simulation scenarios, $b, \ell \leq 100$ and $n\gg b\ell$, so memory and computational costs remain highly favorable [1009.0938].

Practical performance confirms PBA+'s suitability for large systems with predominantly banded structure and sparse extraband perturbations, offering robust accuracy, high efficiency, and broad parallel scalability.

Source: https://www.emergentmind.com/topics/parallel-banding-algorithm-plus-pba