Block-Wise Jacobi Iteration

• Block-wise Jacobi iteration is a parallelizable fixed-point method that partitions variables into blocks, enabling independent and concurrent updates.
• Its algorithmic structure employs localized matrix inversions and optimized GPU routines to solve large-scale linear systems and nonconvex optimization problems.
• Rigorous convergence criteria and adaptive blocking strategies ensure high computational efficiency and scalability in distributed and high-throughput computing architectures.

A block-wise Jacobi iteration is a parallelizable fixed-point method applied to matrix or operator equations, in which the primary variable is partitioned into blocks, each corresponding to a subdomain or coordinate subgroup. The global system is recast so that each block update depends only on the current (or previous) iterates of neighboring or coupled blocks, and all block updates can, in most variants, be performed independently and concurrently. This concept manifests across linear solvers, domain decomposition, nonconvex optimization, classical and non-orthogonal matrix diagonalization, and modern application-specific iterative refinement, with block-wise Jacobi serving as a backbone for scalable, high-throughput computation in distributed and GPU-centric architectures.

1. Mathematical Foundations and Block Decomposition

The block-wise Jacobi procedure begins with a partitioning of the variable (vector or matrix) into $s$ contiguous blocks, accompanied by a corresponding block structure in the main operator (matrix $A$). For a linear system $A u = f$ arising, for instance, from a constant-coefficient 3D finite-volume or Laplacian stencil, the operator $A \in \mathbb{R}^{m \times m}$ is block-partitioned:

$$A = \begin{bmatrix} A_{11} & A_{12} & \cdots & A_{1s} \ A_{21} & A_{22} & \cdots & A_{2s} \ \vdots & \vdots & \ddots & \vdots \ A_{s1} & A_{s2} & \cdots & A_{ss} \end{bmatrix},$$

where $A_{BB} \in \mathbb{R}^{n \times n}$ and $n = b_x b_y b_z$ for spatial blocks ($b_x$, $b_y$, $b_z$ are block dimensions) (Birke et al., 2012). The block-Jacobi splitting sets $M = \mathrm{blockdiag}(A_{11}, \ldots, A_{ss})$, $N = M - A$, with per-block update:

$$x_B^{(k+1)} = A_{BB}^{-1} \left( b_B - \sum_{C \neq B} A_{BC} x_C^{(k)} \right)$$

Convergence is guaranteed if $\rho(I - M^{-1}A) < 1$, typically enforced by block strict diagonal dominance or symmetry/positive-definiteness (Birke et al., 2012). For nonconvex optimization with block linear coupling, the variable $x=(x_1,\ldots,x_s)$ is split with a constraint $\sum_{i=1}^s A_i x_i = b$, giving rise to augmented Lagrangian block updates (Subramanyam et al., 2021), or, in the dual-ROF context, blockwise subdomain solvers (Lee et al., 2019).

2. Algorithmic Structure and Parallel Implementation

Block-wise Jacobi variants are inherently designed for parallel execution:

• Linear systems/Stencil problems: Each block update requires inverting local $A_{BB}$, often possible analytically or using LAPACK routines (for SPD blocks). In GPU architectures, blocks are mapped to thread-groups with local shared memory, invoking CUDA kernels that load block data, execute multiple Jacobi iterations within shared memory, and synchronize global updates (Islam et al., 2020, Birke et al., 2012).
• Optimization frameworks: In proximal-ADMM schemes, each block forms and solves its own local optimization subproblem, with blockwise Bregman distances or strong convexity regularizers ensuring parallel decoupling (Melo et al., 2017, Subramanyam et al., 2021).
• Matrix diagonalization: For symmetric (real or complex) matrices, each Jacobi step targets one or more off-diagonal block pairs, diagonalizes them via block rotations or Eberlein transformations, and applies the result globally. Strategies range from cyclic/quasi-cyclic pivoting to maximum block improvement scheduling (Hari et al., 2016, Begovic et al., 1 Apr 2025, Demmel et al., 4 Jun 2025).

Parallel efficiency hinges on cache/block size selection, load balancing (multi-level blocking, e.g., inter-node, intra-node, cache-level partitions (Singer et al., 2010)), and efficient aggregation of local updates leveraging BLAS2/BLAS3 operations.

3. Convergence Theory and Iteration Complexity

Block-wise Jacobi methods are fundamentally stationary iterations, with rigorous convergence criteria:

• Spectral radius: For linear systems, global contraction requires $\rho(I-M^{-1}A)<1$, quantifiable via smoothing factor and block size (Birke et al., 2012).
• Operator contraction: In matrix diagonalization, the off-norm $S(A)=\|A-\diag(A)\|_F$ decreases geometrically per sweep, with contraction factor $c<1$ determined by block structure and minimum singular value of rotation blocks (UBC property) (Hari et al., 2016).
• Nonconvex optimization: Blockwise-proximal ADMM and similar frameworks guarantee global convergence to stationary points, with iteration complexity $O(1/\sqrt{k})$ in general nonconvex cases or $O(1/k)$ for specific domain decomposition (dual-ROF) settings (Melo et al., 2017, Lee et al., 2019, Subramanyam et al., 2021).
• Nonoverlapping domain decomposition (image processing): Pre-relaxed and fast-Jacobi variants attain accelerated convergence rates ($O(1/n^2)$), with theoretical bounds depending on the number of color-blocks and subdomain interface length (Lee et al., 2019).

Empirically, blockwise Jacobi is robust but may converge slower than block Gauss–Seidel per sweep; however, Jacobi is parallelizable, and synchronous execution can compensate for slower per-iteration progress in large-scale environments.

4. Communication and Arithmetic Cost Analysis

A central theme in high-performance applications is optimal reduction of arithmetic and communication complexity:

• Blocked Jacobi for eigenproblems: With block size $b$, arithmetic cost per sweep is $O(n^2 b + n^3/b)$; communication is minimized when $b = O(\sqrt{M})$, $M$ fast memory, yielding $O(n^3/\sqrt{M})$ words moved per sweep, matching the matrix multiplication lower bound (Demmel et al., 4 Jun 2025).
• Recursive and parallel blocked algorithms: Multi-level cache-oblivious blocking and $P \times P$ processor grid layouts achieve near-optimal complexity in both flops and bandwidth, with cost per sweep telescoping to $O(n^{\omega_0})$ for Strassen-like fast multiplication, where $\omega_0 < 3$ (Demmel et al., 4 Jun 2025, Singer et al., 2010).
• GPU implementations: Hierarchical subdomain Jacobi iterations performed in shared memory yield $6$–$8\times$ speedup over global memory-based Jacobi for large Poisson problems (Islam et al., 2020). Properly chosen block sizes optimize both local arithmetic intensity and inter-block communication.
• High-level optimization: Distributed ADMM blockwise updates communicate only local block residues, enabling near-linear scaling with thousands of cores/nodes (Subramanyam et al., 2021).

5. Blockwise Jacobi in Advanced Applications

Block-wise Jacobi iteration extends beyond classical numerical analysis:

• Non-Euclidean geometry and Bregman distances: Proximal Jacobi-ADMM can target block-specific geometries, accelerating subproblems through tailored metrics (Melo et al., 2017).
• Joint approximate diagonalization: Non-orthogonal Jacobi-type diagonalization on the special linear group leverages GLU/GQU plane updates, Riemannian manifold gradients, and Łojasiewicz convergence theory to guarantee global convergence to stationary points (Li et al., 2020).
• Sampling in deep generative models: GS–Jacobi hybrid iterations in TarFlow image generation exploit block sensitivity metrics to segment tough inversion blocks and allocate adaptive iteration budgets, achieving $2.5$–$5.3\times$ speedups while retaining model quality (FID) (Liu et al., 19 May 2025).
• Image processing: Fast nonoverlapping block Jacobi on dual-ROF domain decomposition matches optimal energy functional decay rates, with practical acceleration when coupled with FISTA-style momentum (Lee et al., 2019).

6. Block Size, Partitioning, and Strategy Selection

Block size ($b$), partitioning, and pivot strategy are pivotal for both convergence and computational efficiency:

• Convergence rate vs. block size: Larger blocks typically decrease the spectral radius of the error propagator, enhancing smoothing and convergence (Birke et al., 2012), but increase per-block computation cost and may limit parallel occupancy.
• Pivot strategy: Cyclic, quasi-cyclic, generalized serial, and maximum block improvement (MBI) strategies impact convergence rate and parallel scheduling. Generalized serial and modulus pivoting are proven globally convergent, with full diagonalization in a constant number of sweeps (Hari et al., 2016, Demmel et al., 4 Jun 2025, Singer et al., 2010).
• Multi-level blocking: Hierarchical (three-level) partitioning maximizes cache reuse, load balancing, and parallel scalability (Singer et al., 2010).
• Adaptive partitioning in learning or sampling models: In TarFlow, block segmentation is influenced by empirical convergence metrics (CRM, IGM) to ensure both stability and throughput (Liu et al., 19 May 2025).

7. Representative Performance and Empirical Results

Comprehensive benchmarks establish block-wise Jacobi as highly scalable:

• Stencil solvers: On $64^3$ patches with constant-coefficient stencils, increasing block size from $(2,2,2)$ to $(8,8,8)$ yields monotonic error reduction, smoothing factor improvement, and GPU throughput of $\geq 67$ million cells/sec (Birke et al., 2012).
• Parallel matrix diagonalization: Three-level block Jacobi for Hermitian matrices matches relative accuracy of classical Jacobi and attains $30$–$70\%$ speedup on $16$–$40$ cores for $n\geq 10,000$ (Singer et al., 2010).
• Image generation/sampling (TarFlow): GS–Jacobi-based sampling yields $2.5$–$5.3\times$ faster inference at negligible FID degradation, with “tough” blocks identified and adapted by convergence ranking metrics; all code and checkpoints are publicly accessible (Liu et al., 19 May 2025).

• (Birke et al., 2012) Block-Relaxation Methods for 3D Constant-Coefficient Stencils on GPUs and Multicore CPUs
• (Melo et al., 2017) Iteration-complexity of a Jacobi-type non-Euclidean ADMM for multi-block linearly constrained nonconvex programs
• (Lee et al., 2019) Fast Nonoverlapping Block Jacobi Method for the Dual Rudin--Osher--Fatemi Model
• (Begovic et al., 1 Apr 2025) On the block Eberlein diagonalization method
• (Subramanyam et al., 2021) A Globally Convergent Distributed Jacobi Scheme for Block-Structured Nonconvex Constrained Optimization Problems
• (Hari et al., 2016) Convergence of the Cyclic and Quasi-cyclic Block Jacobi Methods
• (Islam et al., 2020) Hierarchical Jacobi Iteration for Structured Matrices on GPUs using Shared Memory
• (Demmel et al., 4 Jun 2025) Minimizing the Arithmetic and Communication Complexity of Jacobi's Method for Eigenvalues and Singular Values
• (Singer et al., 2010) Three-Level Parallel J-Jacobi Algorithms for Hermitian Matrices
• (Liu et al., 19 May 2025) Accelerate TarFlow Sampling with GS-Jacobi Iteration
• (Li et al., 2020) Convergence of gradient-based block coordinate descent algorithms for non-orthogonal joint approximate diagonalization of matrices