Component-wise accurate fixed point iterations for computing the square root of a singular M-matrix

1 Introduction

• Square roots and inverse square roots of discretized differential operators arise in Dirichlet-to-Neumann and Neumann-to-Dirichlet maps, as well as fractional diffusion problems in complex network analysis.Source: §1 · pp. 1–2
• An M-matrix has the form A = αI − B with B ≥ 0 and ρ(B) ≤ α; connected-graph Laplacians are examples of irreducible singular M-matrices satisfying property C.Source: §1 · p. 2
• For singular M-matrices, the Binomial iteration is sublinear with generic initialization despite its component-wise numerically stable nonnegative operations.This makes the iteration of scarce practical use in that setting.Source: §1 · p. 2
• A triplet-based CR approach [5] uses a shift to restore quadratic convergence after singularity reduces CR convergence from quadratic to linear, but the shift requires strong assumptions on A.Source: §1 · p. 2
• Suitable initial approximations make the Binomial and U-based iterations linear for singular A, with rates expressible through A's eigenvalues and component-wise relative accuracy.The cost per step is much lower than for CR, so the iterations may be a stable and fast alternative for singular A.Source: §1 · pp. 2–3
• Experiments illustrate the iterations' effectiveness on sparse or bipartite graph Laplacians, where the accelerated CR version from [5] cannot be applied.Source: §1 · p. 3

2 Preliminaries

• For matrices satisfying property C, the principal square root is an M-matrix; the paper assumes A = I − N with nonnegative N and writes A^1/2 = I − G.Source: §2 · p. 3

2.1 Component-wise accurate solution of a system with an M-matrix

• A triplet representation (−offdiag(A), u, v) with Au = v supports cancellation-avoiding recovery of diagonal entries and a GTH-like component-wise accurate solve for Ax = b.The solve is based on an extension of the GTH algorithm.Source: §2.1 · pp. 3–4

3 Reduction to a quadratic matrix equation

• For an M-matrix with property C, A^1/2 = I − G, where G is the minimal nonnegative solution of G^2 − 2G + N = 0 and is generated by the binomial iteration from G_0 = 0.Source: §3 · p. 4
• When A is singular, ρ(G) = 1 and any positive null vector u of A satisfies Gu = u; if 1 is a simple eigenvalue of N, it is the only eigenvalue of G with modulus 1.Source: §3 · p. 4
• If 1 is a simple eigenvalue of N, spectral mapping shows that all other eigenvalues of G have modulus less than 1.Source: §3 · p. 5
• For a more general splitting A = M − N, the square root can be expressed as M^1/2 − W, with W the minimal nonnegative solution of an algebraic Riccati equation.Source: §3 · p. 5

4 The Natural iteration

• The update Xk+1=1/2(Xk^2+N) preserves nonnegativity, and starts satisfying X0u≤u preserve Xku≤u.When X0u=u, the equality is preserved at every step.Source: §4 · pp. 5–6
• With X0=0, the binomial iteration is linear for nonsingular A but sublinear when A is singular.For singular A, the Jacobian at G has spectral radius ρ(G)=1, explaining the sublinear convergence for generic starts.Source: §4 · p. 5
• With X0u=u, the iteration converges linearly to G at rate σ; if also v⊤X0=v⊤, then σ=|λ2|.The rate is bounded by 1/2 max_{i≥2, j≥1}|λi+λj|; the additional left-eigenvector condition gives |λ2|.Source: §4 · pp. 6, 8
• The simple-eigenvalue condition can be relaxed to semisimplicity if X0 matches all right and left eigenvectors associated with eigenvalue 1; X0=I satisfies these conditions.The theorem’s stated setting uses a positive right eigenvector u and a simple eigenvalue 1 of N.Source: §4 · pp. 6, 8
• The component-wise stable implementation avoids cancellation in the nonnegative updates and reconstructs diagonal entries using the triplet representation.The diagonal formula uses the off-diagonal entries and the positive null vector u.Source: §4 · p. 8

5 The U-based iteration

• An iteration derived from a U-based iteration computes the minimal nonnegative solution, computing the triplet for 2I−Xk+1 from that for 2I−Xk and exiting when |Xk+1−Xk|≤ϵ|Xk+1|.Source: §5 · pp. 9–10
• The U-based iteration is reported to converge faster than the natural iteration and has an entry-wise stable implementation.Its implementation uses triplet representations for the iterates and the square-root approximation.Source: §5 · pp. 9–10
• Under the assumptions of Theorem 2, the U-based iterates converge to G; if X0u=u, then Xku=u is preserved and the convergence rate is expressed using G’s eigenvalues.With the additional left-eigenvector condition, the theorem gives a corresponding rate expression.Source: §5 · p. 9

6 On the component-wise numerical stability

• Each step avoids cancellation, so its component-wise relative error is bounded by θu, with θ depending on matrix size and u the machine precision.Errors may accumulate across steps, but the iterations are self-correcting because later steps reduce previously accumulated errors.Source: §6 · p. 10
• The error bound depends on a per-step reduction factor ρ<1.The bound increases as ρ approaches 1.Source: §6 · p. 11

7 On the existence of u > 0 such that Au = 0

• A positive vector u with Au=0 exists exactly when the Frobenius form has a singular diagonal block, each singular block row has no off-diagonal blocks, and each nonsingular block row has at least one.The singular block condition is stated for Aij with j≠i; for nonsingular blocks, at least one Aij with j≠i must be nonzero.Source: §7 · pp. 11–12
• Positive left and right null vectors exist together exactly when the Frobenius form is block diagonal with irreducible singular diagonal blocks.Source: §7 · p. 12
• The component-wise accurate square root can then be computed blockwise from the diagonal blocks.For the stated left-and-right-null-vector condition, the square root is block diagonal with the same block partition.Source: §7 · p. 12

7.1 The case of a singular matrix with triplet

• A right triplet, requiring Au=u′ for u>0 and u′≥0, permits a component-wise stable algorithm under a milder assumption than Au=0.The condition for a right triplet is always fulfilled for M-matrices arising in Markov chains and network models.Source: §7.1 · p. 13
• Using Frobenius form, the procedure computes diagonal block roots and solves for upper off-diagonal blocks; blocks above a singular diagonal block are zero.The resulting problem reduces to square roots of nonsingular or singular irreducible M-matrices and linear systems with M-matrix coefficients.Source: §7.1 · p. 14

8 Numerical tests

• The graph Laplacians are scaled as L = α(I−N), with α equal to the largest node degree, and the experiments compute (I−N)^1/2 = I−G.The four test matrices include banded and bipartite graph constructions.Source: §8 · p. 15
• The tests use graph adjacency matrices with at least one zero in every column, a setting where the accelerated CR technique cannot be applied and standard component-wise stable CR has rate 1/2.The tested minimum modulus entries of G range from 10^-14 to 10^-2; smaller values correspond to greater inaccuracy from standard square-root algorithms.Source: §8 · p. 15
• For singular-M-matrix CR, the reported convergence rate is 1/2; Table 1 reports rates for Algorithms 2 and 3 across four tests with X0 = I and n = 100.The table also reports minimum and maximum entry moduli of the square root.Source: §8 · p. 16, Table 1 · p. 16

8.1 Comparing natural and U-based iteration to CR

• S-CR needs the fewest iterations through quadratic convergence, but its component-wise errors remain large because it is not component-wise stable.Its errors are still smaller than those from Matlab sqrtm.Source: §8.1 · pp. 16–17, Table 2 · p. 17
• Algorithms 2 and 3 take substantially less CPU time than triplet-based CR, while Algorithms 2, 3, and CR reach component-wise relative errors close to machine precision.The comparison includes maximum component-wise relative errors and CPU time across four tests.Source: Table 2 · p. 17, §8.1 · p. 17, Figure 2 · p. 19
• In Test 3, natural and U-based iterations converge much faster than CR; in the other tests, natural is slowest and U-based is close to CR.Figure 1 shows the error dynamics for Tests 1–4.Source: §8.1 · p. 17, Figure 1 · p. 18
• For Tests 1, 2, and 4, Algorithms 2 and 3 require more iterations than CR, but their lower per-step cost yields substantially lower CPU time.Algorithm 2 uses one matrix multiplication and one addition per step; Algorithm 3 alternates inversion methods to lower overall CPU time while preserving convergence and component-wise accuracy.Source: §8.1 · pp. 17, 20

• The two fixed-point iterations compute the square root of a singular M-matrix in a component-wise numerically stable way.Source: §9 · p. 20
• A suitable initial approximation preserves linear convergence, with a convergence factor explicitly expressible using the matrix eigenvalues.Source: §9 · p. 20
• Numerical experiments demonstrate the effectiveness of the iterations for certain classes of problems.Source: §9 · p. 20