pith. sign in

arxiv: 2605.21679 · v1 · pith:KQLSQFWDnew · submitted 2026-05-20 · 🧮 math.NA · cs.NA

Component-wise accurate computation of the square root of an M-matrix

Dario A. Bini , Bruno Iannazzo , Beatrice Meini , Jie Meng This is my paper

Pith reviewed 2026-05-22 09:01 UTC · model grok-4.3

classification 🧮 math.NA cs.NA
keywords M-matrixprincipal square rootcomponent-wise stabilitytriplet representationcyclic reductionNewton iterationnumerical stability
0
0 comments X

The pith

M-matrices admit triplet representations that let adapted iterations compute their principal square roots with component-wise stability independent of singularity or condition number.

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

The paper shows how to compute the principal square root of an M-matrix so that the error in each entry stays small relative to the entry itself. It represents the matrix A by a triplet (P, u, v) where P holds the negated off-diagonal entries, u is positive, and A u equals v. New versions of the cyclic reduction and incremental Newton methods are written to work directly with these triplets. If the original matrix has such a representation then its square root also has one, and the iterations preserve the structure while delivering component-wise accuracy that does not degrade when the matrix is singular or ill-conditioned. Numerical tests confirm the stability claims.

Core claim

If an M-matrix A possesses a triplet representation (P, u, v) with u positive and v nonnegative, then its principal square root exists and likewise possesses a triplet representation. Reformulations of the cyclic reduction and incremental Newton iterations that operate on these triplets therefore compute the square root while remaining component-wise numerically stable, with the stability independent of whether A is singular and independent of the condition number of A.

What carries the argument

The triplet representation (P, u, v) of an M-matrix, which encodes off-diagonal structure in P and the relation A u = v, allowing structure-preserving reformulations of cyclic reduction and incremental Newton iterations.

If this is right

  • The computed square root is guaranteed to be an M-matrix that itself admits a triplet representation.
  • Component-wise relative errors are bounded independently of the condition number of the input matrix.
  • The same stability holds for both singular and nonsingular M-matrices.
  • The adapted iterations preserve the nonnegativity and structural properties required by the triplet form throughout the computation.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • The same triplet encoding might be adapted to compute other matrix functions such as the logarithm or exponential while retaining component-wise accuracy for M-matrices.
  • Applications that repeatedly form square roots of M-matrices arising from discretised differential operators could obtain more reliable results without extra conditioning safeguards.
  • Sparse implementations of the triplet iterations could be tested on large-scale problems from Markov chain analysis to measure practical speed and accuracy gains.

Load-bearing premise

The input M-matrix admits a triplet (P, u, v) with u positive and v nonnegative, and the principal square root inherits a compatible triplet representation.

What would settle it

Apply the triplet-based cyclic reduction or Newton iteration to a concrete singular M-matrix whose exact principal square root is known analytically, then check whether every component-wise relative error remains bounded by a small multiple of machine epsilon.

Figures

Figures reproduced from arXiv: 2605.21679 by Beatrice Meini, Bruno Iannazzo, Dario A. Bini, Jie Meng.

Figure 1
Figure 1. Figure 1: Log plot of the absolute value of the matrix square root in Test 1 [PITH_FULL_IMAGE:figures/full_fig_p021_1.png] view at source ↗
Figure 2
Figure 2. Figure 2: Log plot of the relative errors in the matrix square root of Test 1 with [PITH_FULL_IMAGE:figures/full_fig_p021_2.png] view at source ↗
Figure 3
Figure 3. Figure 3: Log plot of the absolute value of the matrix square root in Test 2 with [PITH_FULL_IMAGE:figures/full_fig_p023_3.png] view at source ↗
Figure 4
Figure 4. Figure 4: Log plot of the relative errors in the matrix square root of Test 3 with [PITH_FULL_IMAGE:figures/full_fig_p023_4.png] view at source ↗
Figure 5
Figure 5. Figure 5: Log plot of the absolute value of the matrix square root in Test 3 with [PITH_FULL_IMAGE:figures/full_fig_p024_5.png] view at source ↗
read the original abstract

Component-wise accurate algorithms for computing the principal square root of an M-matrix are designed in terms of triplet representations. A triplet representation of an M-matrix $A$ is the triple $(P, {\bf u},{\bf v})$, where the matrix $P$ is such that $p_{ij}=-a_{ij}$ for $i\ne j$, $p_{ii}=0$, and ${\bf u}>0$, ${\bf v}\ge 0$ are two vectors such that $A{\bf u}={\bf v}$. It is shown that if $A$ is an M-matrix representable by a triplet, then its principal square root exists and is an M-matrix represented by a triplet as well. New versions of the Cyclic Reduction and the Incremental Newton iterations are provided in terms of triplets, to compute the principal matrix square root of $A$. It is shown that these algorithms are component-wise numerically stable independently of the singularity of $A$ and of its condition number. Numerical experiments are shown to confirm the component-wise stability.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

2 major / 2 minor

Summary. The paper defines a triplet representation (P, u, v) for an M-matrix A with u > 0, v ≥ 0 and A u = v. It proves that if A admits such a representation then its principal square root X also admits a triplet representation, and rewrites the Cyclic Reduction and Incremental Newton iterations to operate directly on the triplet. The central claim is that these adapted algorithms are component-wise numerically stable with relative error O(ε) in each entry of the computed square root, independently of whether A is singular and independently of cond(A). Numerical experiments are presented to illustrate the claimed stability.

Significance. If the component-wise stability result holds, the work supplies a practical and theoretically grounded method for accurate square-root computation on M-matrices that arise in Markov chains, discretizations of elliptic PDEs, and other applications where matrices are often ill-conditioned or singular. The structure-preserving triplet approach and the independence from conditioning constitute a genuine technical contribution; the explicit adaptation of two standard iterations is a clear strength.

major comments (2)
  1. [§5] §4 (adapted iterations) and §5 (stability analysis): the update formulas for the new triplet (P', u', v') contain divisions and subtractions whose denominators are entries of u or linear combinations involving u and v. When A is singular, v lies in the range of A and u may approach a kernel direction, so these denominators can be as small as O(ε)‖u‖ while the mathematical quantities remain well-defined. The perturbation argument given does not supply a uniform bound showing that the relative perturbation in each entry of the computed square root remains O(ε) under these conditions.
  2. [Theorem 3.2] Theorem 3.2 (triplet preservation for the square root): the existence proof is algebraic and relies on the positivity properties of M-matrices, but the subsequent numerical-stability claim in §5 is not shown to be independent of the linear dependence that appears between u and v precisely when A is singular. A concrete counter-example or a refined error-propagation estimate that remains valid in the singular limit is needed to support the strongest claim.
minor comments (2)
  1. [Numerical experiments] The numerical experiments section would benefit from explicit reporting of the smallest singular values of the test matrices and of the observed relative errors in the singular cases, so that the independence claim can be directly verified from the tables.
  2. [Notation] Notation: the vectors u and v are sometimes written in bold and sometimes not; consistent boldface throughout would improve readability.

Simulated Author's Rebuttal

2 responses · 0 unresolved

We thank the referee for the detailed and constructive report. The comments highlight important aspects of the stability analysis in the singular case, and we address each point below. We will revise the manuscript to provide a refined error analysis that explicitly treats the singular limit.

read point-by-point responses

  1. Referee: [§5] §4 (adapted iterations) and §5 (stability analysis): the update formulas for the new triplet (P', u', v') contain divisions and subtractions whose denominators are entries of u or linear combinations involving u and v. When A is singular, v lies in the range of A and u may approach a kernel direction, so these denominators can be as small as O(ε)‖u‖ while the mathematical quantities remain well-defined. The perturbation argument given does not supply a uniform bound showing that the relative perturbation in each entry of the computed square root remains O(ε) under these conditions.

    Authors: We agree that the perturbation argument in §5 needs to be strengthened for the singular case. Although the current analysis relies on positivity and component-wise bounds derived from the M-matrix structure, it does not explicitly derive a uniform relative-error bound when denominators in the triplet updates become small. In the revised manuscript we will add a refined perturbation analysis in §5 that tracks the error propagation through the divisions and subtractions by exploiting the fact that the triplet updates preserve the M-matrix property and that the accumulated rounding errors remain relatively small in each entry. This will establish the O(ε) component-wise bound independently of singularity. revision: yes

  2. Referee: [Theorem 3.2] Theorem 3.2 (triplet preservation for the square root): the existence proof is algebraic and relies on the positivity properties of M-matrices, but the subsequent numerical-stability claim in §5 is not shown to be independent of the linear dependence that appears between u and v precisely when A is singular. A concrete counter-example or a refined error-propagation estimate that remains valid in the singular limit is needed to support the strongest claim.

    Authors: Theorem 3.2 is an algebraic existence result that holds for both singular and nonsingular M-matrices; it does not address numerical stability. The stability claim is made for the adapted algorithms and is analyzed in §5. We acknowledge that the present error-propagation argument does not explicitly address the linear dependence between u and v that arises at singularity. In the revision we will supply a refined component-wise error-propagation estimate that remains valid in the singular limit, using the limiting behavior of the triplet updates and the nonnegativity properties of M-matrices. We will also augment the numerical experiments with additional singular and near-singular test cases to illustrate the claimed independence from cond(A). revision: yes

Circularity Check

0 steps flagged

No circularity: structure preservation and stability proved directly from triplet relations and positivity

full rationale

The paper first proves (via direct algebraic verification using the defining relation A u = v and off-diagonal nonnegativity) that an M-matrix given by triplet (P, u, v) has a principal square root that also admits a triplet representation. It then rewrites the standard Cyclic Reduction and Incremental Newton iterations to act componentwise on the triplet entries while preserving the M-matrix structure. Component-wise backward stability independent of singularity and cond(A) is obtained by explicit perturbation bounds that exploit the positivity of u and the exact relation A u = v to control relative errors entrywise; these bounds are derived from the iteration formulas themselves rather than from any fitted parameter, self-referential prediction, or load-bearing self-citation. The derivation chain is therefore self-contained and does not reduce any claimed result to its own inputs by construction.

Axiom & Free-Parameter Ledger

0 free parameters · 1 axioms · 1 invented entities

The approach rests on standard properties of M-matrices and the newly introduced triplet representation; no fitted numerical parameters are used.

axioms (1)
  • domain assumption An M-matrix representable by a triplet has a principal square root that is also an M-matrix representable by a triplet.
    Invoked to guarantee existence and structure preservation for the algorithms.
invented entities (1)
  • Triplet representation (P, u, v) no independent evidence
    purpose: To encode M-matrix structure for component-wise accurate iteration
    New representational device introduced to enable the stable algorithms; no independent falsifiable evidence outside the paper is provided.

pith-pipeline@v0.9.0 · 5719 in / 1230 out tokens · 59858 ms · 2026-05-22T09:01:21.774101+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Lean theorems connected to this paper

Citations machine-checked in the Pith Canon. Every link opens the source theorem in the public Lean library.

What do these tags mean?
matches
The paper's claim is directly supported by a theorem in the formal canon.
supports
The theorem supports part of the paper's argument, but the paper may add assumptions or extra steps.
extends
The paper goes beyond the formal theorem; the theorem is a base layer rather than the whole result.
uses
The paper appears to rely on the theorem as machinery.
contradicts
The paper's claim conflicts with a theorem or certificate in the canon.
unclear
Pith found a possible connection, but the passage is too broad, indirect, or ambiguous to say the theorem truly supports the claim.

Reference graph

Works this paper leans on

26 extracted references · 26 canonical work pages · 1 internal anchor

  1. [1]

    Alefeld and N

    G. Alefeld and N. Schneider. On square roots ofM-matrices.Linear Algebra Appl., 42:119–132, 1982

    work page 1982

  2. [2]

    A. S. Alfa, J. Xue, and Q. Ye. Accurate computation of the smallest eigenvalue of a diagonally dominantM-matrix.Math. Comp., 71(237):217– 236, 2002. 25

    work page 2002

  3. [3]

    Benzi, D

    M. Benzi, D. Bertaccini, F. Durastante, and I. Simunec. Non-local network dynamics via fractional graph Laplacians.J. Complex Netw., 8(3):cnaa017, 29, 2020

    work page 2020

  4. [4]

    Berman and R

    A. Berman and R. J. Plemmons.Nonnegative matrices in the mathematical sciences, volume 9 ofClassics in Applied Mathematics. Society for Indus- trial and Applied Mathematics (SIAM), Philadelphia, PA, 1994. Revised reprint of the 1979 original

    work page 1994

  5. [5]

    Berman, R

    A. Berman, R. S. Varga, and R. C. Ward.ALPS: matrices with nonposi- tive off-diagonal entries.Linear Algebra Appl., 21(3):233–244, 1978

    work page 1978

  6. [6]

    D. A. Bini, G. Latouche, and B. Meini.Numerical methods for structured Markov chains. Numerical Mathematics and Scientific Computation. Ox- ford University Press, New York, 2005

    work page 2005

  7. [7]

    Chen, R.-C

    C. Chen, R.-C. Li, and C. Ma. Highly accurate doubling algorithm for quadratic matrix equation from quasi-birth-and-death process.Linear Al- gebra Appl., 583:1–45, 2019

    work page 2019

  8. [8]

    Deadman, N

    E. Deadman, N. J. Higham, and R. M. Ralha. Blocked schur algorithms for computing the matrix square root.Lecture Notes in Computer Sci- ence, 7782 LNCS:171 – 182, 2013. All Open Access; Green Accepted Open Access; Green Open Access

    work page 2013

  9. [9]

    Druskin and L

    V. Druskin and L. Knizhnerman. Gaussian spectral rules for the three-point second differences: I. a two-point positive definite problem in a semi-infinite domain.SIAM Journal on Numerical Analysis, 37(2):403–422, 1999

    work page 1999

  10. [10]

    Evard and F

    J.-C. Evard and F. Uhlig. On the matrix equationf(X) =A. volume 162/164, pages 447–519. 1992. Directions in matrix theory (Auburn, AL, 1990)

    work page 1992

  11. [11]

    Frommer, G

    A. Frommer, G. Ramirez-Hidalgo, M. Schweitzer, and M. Tsolakis. Poly- nomial preconditioning for the action of the matrix square root and inverse square root.ArXiv preprint arXiv:2401.06684, 2024

    work page arXiv 2024

  12. [12]

    W. K. Grassmann, M. I. Taksar, and D. P. Heyman. Regenerative analysis and steady state distributions for markov chains.Operations Research, 33(5):1107 – 1116, 1985

    work page 1985

  13. [13]

    G. Gu, L. Wang, and R.-C. Li. Highly accurate Latouche-Ramaswami logarithmic reduction algorithm for quasi-birth-and-death process.J. Math. Study, 55(2):180–194, 2022

    work page 2022

  14. [14]

    C.-H. Guo. On algebraic Riccati equations associated withM-matrices. Linear Algebra Appl., 439(10):2800–2814, 2013

    work page 2013

  15. [15]

    Guo and D

    C.-H. Guo and D. Lu. On algebraic Riccati equations associated with regular singularM-matrices.Linear Algebra Appl., 493:108–119, 2016. 26

    work page 2016

  16. [16]

    N. J. Higham.Functions of matrices. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2008. Theory and computation

    work page 2008

  17. [17]

    Iannazzo

    B. Iannazzo. A note on computing the matrix square root.Calcolo, 40(4):273 – 283, 2003

    work page 2003

  18. [18]

    Iannazzo and B

    B. Iannazzo and B. Meini. The palindromic cyclic reduction and related algorithms.Calcolo, 52(1):25–43, 2015

    work page 2015

  19. [19]

    Liu, W.-G

    C. Liu, W.-G. Wang, J. Xue, and R.-C. Li. Accurate numerical solution for structuredM-matrix algebraic Riccati equations.J. Comput. Appl. Math., 396:Paper No. 113614, 16, 2021

    work page 2021

  20. [20]

    C. Liu, J. Xue, and R.-C. Li. Accurate numerical solution for shiftedM- matrix algebraic Riccati equations.J. Sci. Comput., 84(1):Paper No. 15, 27, 2020

    work page 2020

  21. [21]

    B. Meini. The matrix square root from a new functional perspective: the- oretical results and computational issues.SIAM J. Matrix Anal. Appl., 26(2):362–376, 2004/05

    work page 2004

  22. [22]

    G. T. Nguyen and F. Poloni. Componentwise accurate fluid queue compu- tations using doubling algorithms.Numer. Math., 130(4):763–792, 2015

    work page 2015

  23. [23]

    G. T. Nguyen and F. Poloni. Componentwise accurate Brownian motion computations using cyclic reduction.arXiv preprint arXiv:1605.01482, 2016

    work page internal anchor Pith review Pith/arXiv arXiv 2016

  24. [24]

    Schneider

    H. Schneider. An inequality for latent roots applied to determinants with dominant principal diagonal.J. London Math. Soc., 28:8–20, 1953

    work page 1953

  25. [25]

    Xue and R.-C

    J. Xue and R.-C. Li. Highly accurate doubling algorithms forM-matrix algebraic Riccati equations.Numer. Math., 135(3):733–767, 2017

    work page 2017

  26. [26]

    J. Xue, S. Xu, and R.-C. Li. Accurate solutions ofM-matrix Sylvester equations.Numer. Math., 120(4):639–670, 2012. 27

    work page 2012