Pith. sign in

REVIEW 3 major objections 5 minor 23 references

The Stability of Block Eliminations and Additive Modifications

T0 review · 3 major / 5 minor · reviewed 2026-08-04 · deepseek-v4-flash

Pith's one-line read Block LU factorization, as used in the BEAM solver, has backward error bounded by the square of a block-aligned growth factor rather than the cube of the pointwise growth factor, extending stability analysis to matrices that pointwise LU ca

desk verdict Solid improvement of the Demmel-Higham-Schreiber stability bound for block LU; the BEAM transfer is not as clean as the headline claims, but the core theorem holds. read the letter →

arxiv 2509.07305 v1 pith:LVXDK5NP submitted 2025-09-09 math.NA cs.NA

classification math.NAcs.NA MSC 15A2365F05
keywords blockLUfactorizationgrowthfactorbackwarderroranalysisBEAMadditivemodificationsGaussianeliminationsingularvaluedecompositionconditionnumber
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

This paper proves that block LU factorization is more stable than previous analyses indicated: its backward error grows at most quadratically in a block-aligned growth factor, not cubically in Wilkinson's classic pointwise growth factor. The improvement matters because block LU is the core of the BEAM solver, a communication-avoiding alternative to LU with partial pivoting, and the new measure lets the theory cover matrices for which ordinary pointwise LU does not even exist, such as a matrix whose leading 2x2 block is [0 1; 1 0]. The proof route is a general backward-error theorem for block LU combined with a norm bound for the computed L and R factors when the diagonal blocks of L are unitary, which is exactly the SVD-based situation in BEAM. The paper separately analyzes BEAM's additive modifications, bounding how far they can inflate the condition number and showing that Zielke's matrix still produces exponential growth. It closes with simple thresholds, based on diagonal dominance and related quantities, under which BEAM applies no modifications at all, so the behavior of ordinary block LU is preserved.

What carries the argument

The central object is the block-aligned growth factor P^I_alpha, defined as the maximum over elimination steps of the norm of the current Schur complement divided by the norm of the original matrix; it replaces Wilkinson's pointwise growth factor and is independent of how diagonal blocks are factored. The second main ingredient is the blockwise backward-error Theorem 2.3, which reduces the error of the whole factorization to assumed normwise errors of each block operation. Theorem 2.6 then bounds ||L|| ||R|| in Schatten p-norms when diagonal blocks of L are unitary, using the Schur-complement identity for subdiagonal blocks of L; the product is O(n_t^2 P_p^2 kappa_p(A) ||A||), which is what

What would settle it

Compute BEAM's backward error ||A - Lhat Rhat||_2 on a matrix with moderate blockwise growth but with a diagonal block whose smallest singular value sits just below tau and whose singular vectors interact badly with the trailing update; if the per-unit-roundoff error grows faster than P^2 as the block count increases, the quadratic bound does not transfer to the implemented clamped-SVD step.

Watch

Extended reading notes

Core claim

The paper's central claim is that the backward error of block LU is controlled by a blockwise growth factor P^I_alpha = max_k ||A^(k)||_alpha / ||A||_alpha over the Schur complements produced during elimination, and that when the diagonal blocks of L are unitary (as they are in BEAM), the bound becomes O(u P^I ||A||) + O(u n_t^2 (P^I)^2 kappa(A) ||A||) up to polynomial factors. This is quadratic in the blockwise growth factor, improving on the previous cubic-in-pointwise-growth bound, and it applies to matrices whose pointwise LU factorization is undefined. The argument proceeds by first proving a general backward-error theorem for Algorithm 2.1 under per-block normwise assumptions, then bou

Load-bearing premise

The load-bearing premise is that every block operation in Algorithm 2.1, including BEAM's SVD factorization of a diagonal block with small singular values clamped to tau, has roundoff error bounded by a small constant times machine precision; the paper assumes this for the clamped-SVD step rather than proving it.

Editorial extensions

If this is right

  • A backward error that scales like (P^I)^2 instead of (P^I)^3 means BEAM and other block LU implementations can tolerate substantially larger element growth before the solution is polluted by roundoff.
  • Matrices that pointwise LU cannot factor, such as those with a leading block [0 1; 1 0], are now covered by a stability guarantee, so the blockwise analysis broadens the class of matrices for which pivotless block elimination is justified.
  • For block diagonally dominant matrices and inverses of row-block-diagonally dominant matrices, the combined bounds give backward error constants independent of the condition number, matching the best known pointwise results without pivoting.
  • If tau-hat kappa_2(A) is at most 0.5, BEAM's modifications change the inverse by at most a factor of 2, so the Woodbury-based correction step is provably stable for moderately conditioned or mildly modified systems.
  • For diagonally dominant, H-, symmetric positive definite, and block diagonally dominant matrices, BEAM takes the modification-free path for sufficiently small tau, making the extra Woodbury work unnecessary on those well-behaved classes.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • Because the growth factor P^I_alpha is defined on Schur complements that are already computed, it could be monitored cheaply at runtime and used to decide when to switch from BEAM to a pivoting strategy.
  • A direct rounding-error analysis of the clamped-SVD diagonal-block factorization would settle whether the assumed bound (2.4) holds with moderate constants; until then, the quadratic bound applies to BEAM's formulation conditional on that backward stability.
  • The Zielke example suggests that additive modifications alone cannot eliminate worst-case growth, so a practical BEAM variant might combine SVD clamping with block-dependent pivot selection, as hinted by the diagonal-block pivoting strategy noted in the paper.
  • The modification-free criteria for H-matrices could be deployed as a cheap precheck that skips Woodbury correction and iterative refinement entirely for matrices whose diagonal dominance factors are large enough.
Share X Bluesky LinkedIn Reddit HN

Signed reviews

No signed human review yet.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

3 major / 5 minor

Summary. The paper studies the numerical stability of block LU factorization and of the BEAM algorithm, which replaces pivoting by additive modifications of small singular values in the diagonal blocks. It defines a block-aligned growth factor P^I_α, proves a backward error bound for block LU that is quadratic in this growth factor rather than cubic in the pointwise growth factor (Theorems 2.3 and 2.6), and applies the framework to several matrix classes (diagonally dominant, block diagonally dominant, SPD). In the second half it analyzes the effect of BEAM's additive modifications on the condition number and growth of the modified matrix, gives an exponential bound for the general case, studies the Zielke matrix as a worst-case example, and provides sufficient conditions, based on diagonal dominance, H-matrix properties, and SPD structure, under which BEAM makes no modifications.

Significance. The block LU result is a genuine theoretical improvement: it is not limited to the standard formulation with identity diagonal blocks of L, it replaces a cubic dependence on the pointwise growth factor with a quadratic dependence on a block-adapted growth factor, and it applies to matrices for which pointwise LU may not even exist, such as a matrix with leading block [0 1; 1 0]. The paper is careful in stating assumptions (2.4)-(2.7) with explicit constants, and it is transparent about limitations, notably that the general bound on the inverse of the modified matrix is exponential. If the BEAM-specific claims are properly qualified, the analysis is a useful contribution to the numerical stability literature.

major comments (3)
  1. [§2.4, Table 2.1, Thms. 2.6–2.7] For BEAM, the matrix actually factored is eA = A + M_U M_Σ M_V^T, as stated in §1.1. Applying Theorem 2.6 to eA replaces P, κ, and ||A|| by P(eA), κ(eA), and ||eA||. The paper does not prove κ(eA) ≤ min(κ(A), τ̂^-1) nor P(eA) ≤ C P(A); Section 3 leaves the general bound on ψ = ||eA^{-1}||/||A^{-1}|| exponential in n (Thm. 3.2). Therefore the Table 2.1 row “BEAM modified (P^I)^2 min(κ(A), τ̂^-1)” is not supported. Either prove the required bounds for eA or state the BEAM row as conditional on the growth and condition of eA.
  2. [§2.2, eqs. (2.4)–(2.7), Algorithm 1.1 lines 6–13] For the clamped-SVD diagonal factorization, the unmodified singular value σ_i ≤ τ is replaced by τ, so U_k Σ_clamped V_k^T differs from the input diagonal block by (τ−σ_i) u_i v_i^T, an O(τ) perturbation rather than O(u)||A_kk||. The paper's escape is that BEAM factors the modified matrix eA, but this transfer is implicit. Please state and prove a backward-stability lemma: with eA defined using the computed SVD factors and clamped values, (2.4) holds for eA with an error of order u||A_kk||. Also clarify that the quantities A^(k) in Theorem 2.3 are exact Schur complements; applying the theorem to a recursion on computed Schur complements requires an additional accumulation argument.
  3. [§3, second paragraph] The display “P_α ≤ 1 + max_k || eA^{-1}_{1:k,1:k} ||_α ||A||_α” uses the original matrix norm ||A||_α, whereas the analogue of (2.3) for the matrix eA being factored must use ||eA||_α. As written this is not a consequence of Theorem 2.2 and it obscures which matrix the growth factor P refers to. This is not a pure typo because it feeds into the overclaim about the BEAM row in Table 2.1.
minor comments (5)
  1. [Abstract] “we analyzed the modifications” should be “we analyze” to match the tense of the rest of the abstract.
  2. [§2.2, Theorem 2.4 proof] The summation limits in the first displayed bound after “Continuing the iteration” are swapped relative to the theorem statement: c_L1 should run to n_t and c_L2 to n_t−1. The theorem statement has the correct limits.
  3. [§2.3] “A(k) independent of how the diagonal blocks are factored” is missing a verb; should be “A(k) is independent”.
  4. [§3, after Thm. 3.2] The sentence “all by the last row” should be “all but the last row”.
  5. [§2.6 (Theorem 2.6 proof)] The proof uses implicitly that (A^(k))^{-1} is a trailing principal submatrix of A^{-1}, so that ||(A^(k))^{-1}||_p ≤ ||A^{-1}||_p. This should be stated, since it is the key step relating the L-factor to κ_p(A).

Circularity Check

0 steps flagged · score 2.0 of 10

No significant circularity: central block-LU error bound is derived from exact Schur-complement quantities under explicit stability assumptions; flagged caveats (self-citations and the unproven clamped-SVD stability premise) do not reduce the result to its inputs.

full rationale

The paper's central claim—the quadratic-in-growth backward-error bound (Theorems 2.3 and 2.6)—is not circular. Theorem 2.3 is a conditional bound parameterizing per-block operation errors (eqs. 2.4-2.7) and bounding the residual A - Lhat Rhat in terms of the exact Schur-complement growth factor P (eq. 2.1); P is defined via exact arithmetic quantities (successive Schur complements), not from the rounding errors being bounded, and no parameter is fitted to the target quantity. Theorem 2.6's bound on ||L||_p ||R||_p follows from the exact identity (2.11) and blockwise inversion, with the external identity cited from Demmel et al. [4]; the new growth factor is explicitly related to Wilkinson's classic factor via Theorem 2.1, so the cubic-to-quadratic improvement is not an artifact of re-defining the growth measure. The authorial self-citations are not load-bearing: the dissertation [13] is provenance (the section 2 analysis appears in full in this paper), and the one cited prior bound [14, Thm. 4.3] used in the proof of Theorem 3.1 is re-derived within the paper via the Woodbury formula (eqs. 3.1-3.3). The legitimate limitation is that the applicability to BEAM requires the clamped-SVD step to satisfy the normwise bounds (2.4)-(2.7) for the modified matrix; the paper states these as assumptions ('For the sake of flexibility, we parameterize the error of the blockwise operations') and does not prove backward stability of the clamping step. This is a correctness gap (as is the paper's own admission that 'bounding psi is a difficult question' in Section 3), not a circular reduction: the bound does not assume the smallness of the residual it purports to prove.

Assumptions & free parameters 0 free parameters · 6 assumptions · 0 invented entities

No parameters are fitted to data. The analysis depends on standard norm properties, block nonsingularity, and cited external bounds; the tolerance tau-hat is an algorithmic input, not a fitted constant.

assumptions (6)
  • domain assumption The dimension-invariant matrix norm family satisfies the partition inequality (1.1), including submultiplicativity for the norms used in the main theorems.
    Introduced in Section 1.2 to handle submatrices of different sizes; needed for bounds such as ||AB||_alpha <= ||A||_alpha ||B||_alpha in Theorems 2.2-2.9.
  • domain assumption A is block strongly nonsingular with respect to the blocking I, so Algorithm 2.1 completes.
    Section 2, before subsection 2.1: 'this definition is critical as Algorithm 2.1 completes if and only if A is block strongly nonsingular.'
  • domain assumption The diagonal blocks of L are unitary, or the norm is invariant to them, for the ||L|| ||R|| bounds.
    Theorems 2.6-2.9 assume unitary diagonal blocks of L, matching BEAM's SVD-based diagonal factorization.
  • standard math Schatten p-norm submultiplicativity and the principal submatrix norm inequality ||P^* M P||_p <= ||M||_p.
    Used in Theorem 2.6's proof to bound ||((A^{(k)})^{-1})_{k+1:n_t,k}||_p by ||A^{-1}||_p via the blockwise inversion identity.
  • standard math The Woodbury formula and invertibility of the capacitance matrix I + F for the BEAM solve.
    Section 3 uses the Woodbury formula in equation (3.1) to express A~^{-1}; the associated bounds require I + F to be nonsingular.
  • domain assumption Known external results on diagonal dominance, H-matrices, Varah's bound, and the Zielke matrix spectrum.
    Section 4 and Theorem 3.3 rely on Plemmons, Varah, Varga, Turing, and Elkies for existing bounds; these are cited but not proved in the paper.

how reviews work

0 comments
Cite this review

Pith. "Pith review of The Stability of Block Eliminations and Additive Modifications." pith.science (2026). https://pith.science/paper/LVXDK5NP

@misc{pith2026250907305,
  author       = {Pith},
  title        = {Pith review of: The Stability of Block Eliminations and Additive Modifications},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/LVXDK5NP}},
  note         = {Machine review of arXiv:2509.07305}
}
read the original abstract

The block elimination with additive modifications (BEAM) method was recently proposed as a alternative to LU with partial pivoting requiring less communication. Because of the novelty of BEAM, the existing theoretical analysis is lacking. To that end, we analyze both the numerical stability of the underlying block LU factorization and the effects of additive modifications. For the block LU factorization, we are able to improve the previous results of Demmel et al. from being cubic in the element growth to merely quadratic. Furthermore, we propose an alternative measure of element growth that is better aligned with block LU; this new measure of growth allows our analysis to apply to matrices that cannot be factored with pointwise LU. In the second part, we analyzed the modifications produced by BEAM and the effect they have on the condition number and growth factor. Finally, we show that BEAM will not apply any modifications in some cases that regular block LU can safely factor.

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

23 extracted references · 16 canonical work pages

  1. [1]

    Amodio and F

    P. Amodio and F. Mazzia , A new approach to backward error analysis of LU factorization , BIT Numerical Mathematics, 39 (1999), pp. 385--402, https://doi.org/10.1023/A:1022358300517

  2. [2]

    J. L. Barlow and H. Zha , Growth in Gaussian elimination, orthogonal matrices, and the 2-norm , SIAM Journal on Matrix Analysis and Applications, 19 (1998), pp. 807--815, https://doi.org/10.1137/S0895479896309912

  3. [3]

    J. W. Demmel and N. J. Higham , Stability of block algorithms with fast level-3 BLAS , ACM Trans. Math. Softw., 18 (1992), pp. 274--291, https://doi.org/10.1145/131766.131769

  4. [4]

    J. W. Demmel, N. J. Higham, and R. S. Schreiber , Stability of block LU factorization , Numerical Linear Algebra with Applications, 2 (1995), pp. 173--190, https://doi.org/10.1002/nla.1680020208

  5. [5]

    F. M. Dopico and J. M. Molera , Perturbation Theory for Factorizations of LU Type through Series Expansions , SIAM Journal on Matrix Analysis and Applications, 27 (2005), pp. 561--581, https://doi.org/10.1137/040612142

  6. [6]

    I. S. Duff, A. M. Erisman, and J. K. Reid , Direct Methods for Sparse Matrices , Oxford University Press, Oxford, UK, second ed., Mar. 2017

  7. [7]

    N. D. Elkies , Answer to ``2-norm of the upper triangular ``all-ones'' matrix'' , Aug. 2011, https://mathoverflow.net/a/72383

  8. [8]

    George and K

    A. George and K. D. Ikramov , Gaussian elimination is stable for the inverse of a diagonally dominant matrix , Mathematics of Computation, 73 (2004), pp. 653--657, http://www.jstor.org/stable/4099794, https://arxiv.org/abs/4099794

Show all 23 references
  1. [9]

    George, K

    A. George, K. D. Ikramov, and A. B. Kucherov , On the growth factor in Gaussian elimination for generalized Higham matrices , Numerical Linear Algebra with Applications, 9 (2002), pp. 107--114, https://doi.org/10.1002/nla.258

  2. [10]

    George and Kh

    A. George and Kh . D. Ikramov , Block LU factorization is stable for block matrices whose inverses are block diagonally dominant , Journal of Mathematical Sciences, 127 (2005), pp. 1962--1968, https://doi.org/10.1007/s10958-005-0154-7

  3. [11]

    N. J. Higham , Accuracy and Stability of Numerical Algorithms , Society for Industrial and Applied Mathematics , Philadelphia, PA, USA, second ed., 2002, https://doi.org/10.1137/1.9780898718027

  4. [12]

    X. S. Li and J. Demmel , Making sparse Gaussian elimination scalable by static pivoting , in SC '98: Proceedings of the 1998 ACM / IEEE Conference on Supercomputing, San Jose, CA, USA, Nov. 1998, IEEE Computer Society, pp. 34--34, https://doi.org/10.1109/SC.1998.10030

  5. [13]

    Lindquist , Reducing Communication in the Solution of Linear Systems , PhD thesis, The University of Tennessee, Knoxville, TN, USA, Aug

    N. Lindquist , Reducing Communication in the Solution of Linear Systems , PhD thesis, The University of Tennessee, Knoxville, TN, USA, Aug. 2023

  6. [14]

    Lindquist, P

    N. Lindquist, P. Luszczek, and J. Dongarra , Using additive modifications in LU factorization instead of pivoting , in Proceedings of the 37th ACM International Conference on Supercomputing , Orlando, FL, USA, June 2023, ACM, pp. 14--24, https://doi.org/10.1145/3577193.3593731

  7. [15]

    R. J. Plemmons , M-matrix characterizations. I ---nonsingular M-matrices , Linear Algebra and its Applications, 18 (1977), pp. 175--188, https://doi.org/10.1016/0024-3795(77)90073-8

  8. [16]

    A. M. Turing , Rounding-off errors in matrix processes , The Quarterly Journal of Mechanics and Applied Mathematics, 1 (1948), pp. 287--308, https://doi.org/10.1093/qjmam/1.1.287

  9. [17]

    J. M. Varah , A lower bound for the smallest singular value of a matrix , Linear Algebra and its Applications, 11 (1975), pp. 3--5, https://doi.org/10.1016/0024-3795(75)90112-3

  10. [18]

    R. S. Varga , On diagonal dominance arguments for bounding A -1 , Linear Algebra and its Applications, 14 (1976), pp. 211--217, https://doi.org/10.1016/0024-3795(76)90067-7

  11. [19]

    J. H. Wilkinson , Error analysis of direct methods of matrix inversion , Journal of the ACM, 8 (1961), pp. 281--330, https://doi.org/10.1145/321075.321076

  12. [20]

    J. H. Wilkinson , Rounding Errors in Algebraic Processes , Prentice-Hall, Princeton, NJ, USA, 1963, https://archive.org/details/roundingerrorsin0000wilk

  13. [21]

    M. A. Woodbury , Inverting Modified Matrices , vol. 42 of Memorandum Report, Statistical Research Group, Princeton, NJ, 1950

  14. [22]

    Zhang, T

    F. Zhang, T. Ando, C. Brezinski, R. A. Horn, C. R. Johnson, J. Liu, S. Putanen, R. L. Smith, and G. P. H. Styan , The Schur Complement and Its Applications , Numerical Methods and Algorithms , Springer Science & Business Media, New York, NY, 2005, https://doi.org/10.1007/b105056

  15. [23]

    Zielke , Testmatrizen mit maximaler Konditionszahl , Computing, 13 (1974), pp

    G. Zielke , Testmatrizen mit maximaler Konditionszahl , Computing, 13 (1974), pp. 33--54, https://doi.org/10.1007/BF02268390

Pith tools

Reviewed August 4, 2026 · model on record in the stance chip above.