PaperlayerView exact revision on arXiv

Source-linked AI summary

Halley's Method for Rectangular Matrix Variables and the Matrix Schwarzian Derivative

Shintaro Yoshizawa

arXiv:2609.06592v1math.NAcs.NA

TL;DR

Rectangular matrix generalizations of Halley’s method must address division, non-commutativity, and computation without forming a large Hessian operator. The paper develops a matrix Schwarzian framework and proves local cubic convergence; numerical tests examine convergence behavior and practical trade-offs.

  • Problem

    Extending Halley’s method to rectangular matrix variables raises unresolved questions about division, non-commutativity, and feasible computation.

    Source: §1.2 · p. 3

  • Method

    The paper develops a matrix Schwarzian framework for Halley’s method, interprets it through information geometry, and derives a square-root-free construction from the Newton map’s third derivative.

    Source: §9 · p. 28

  • Results

    High-precision tests confirm cubic convergence for Halley and Laguerre and quadratic convergence for Newton, while Oja-system tests find convergence quality governed by target eigenvalue gaps and cubic gains growing with ill-conditioning.

    Source: §7.4 · p. 22, §8.6 · p. 27

  • Takeaways & Limitations

    Algorithm choice should weigh convergence order against per-step cost and the required accuracy.

    Source: §6.4 · p. 16

    A uniform scalar power parameter cannot generally cancel the matrix-valued second-order error, and its naive matrix extension remains essentially quadratic in the reported experiments.

    Source: §4.2 · p. 12

Abstract

from arXiv · show

Alefeld (1981) recast Halley's cubic convergence as Newton's method on $g=f/\sqrt{f'}$, linked to the Schwarzian derivative. We generalize this to matrix gradient fields $F=\nabla\phi$ ($X\in\R^{m\times n}$) via a matrix Schwarzian derivative interpreted through information geometry ($\alpha$-connections). Simplifying Palmore (1994), we construct this operator square-root-free derivative from the third Fr\'echet derivative of the Newton map. Our main theorem proves local cubic convergence with an explicit error constant, without self-adjointness, commutativity, or gradient-field assumptions. A matrix-free algorithm (Hessian-vector products only) validates the theory. We contrast this with a power-Newton family (whose naive matrix extension fails) and the matrix Laguerre family (which requires an operator square root). A case study on the Oja-type system $\dot X=AXB-XBX^TAX$ reveals that convergence depends on target eigenvalue gaps, cubic gains grow with ill-conditioning, and 70-digit tests confirm exact theoretical orders alongside a working-precision accuracy budget. Finally, we note coupled formulations excel primarily when sequential deflation is inapplicable.

1 Introduction

The introduction traces Halley’s cubic convergence to the scalar Schwarzian formulation, then frames the obstacles to extending it to rectangular matrix variables. It presents a matrix generalization, an information-geometric interpretation, and a matrix-free implementation as the paper’s goals.

1.1 Background: Alefeld’s scalar theory

  • Alefeld recasts Halley’s method as Newton’s method on g=f/sqrt(f'), whose second derivative is controlled by f’s Schwarzian derivative.Alefeld uses a uniform bound on g'' to obtain a Newton–Kantorovich-type convergence result.Source: §1.1 · p. 3

1.2 Three challenges in the matrix generalization

  • Rectangular matrices raise challenges in defining division, handling non-commutativity, and computing without forming an infeasibly large Hessian operator.For large dimensions, the paper identifies Hessian-vector products as essential to a practical matrix-free implementation.Source: §1.2 · p. 3

1.3 Contributions and outline

  • The paper generalizes Halley’s theory to rectangular matrix variables, connects the matrix Schwarzian to information geometry, and develops a matrix-free implementation.Its stated goals include verifying cubic convergence numerically and quantifying trade-offs against Newton’s method.Source: §1.3 · p. 4

2 Resolving the First Two Challenges

The paper resolves the first two matrix-generalization challenges by restricting the construction to gradient fields on a shared matrix space and by replacing scalar-style equalities with operator-norm inequalities when operators do not commute.

2.1 Challenge I: the gradient-field assumption

  • A strongly convex potential makes F a gradient field on the same matrix space, so its Hessian operator is self-adjoint positive and has a well-defined inverse and square root.This shared domain and codomain is what makes the operator construction possible for rectangular X.Source: §2.1 · p. 4

2.2 Challenge II: making explicit where the equality breaks down

  • For non-commuting operators, the Daleckii–Krein formula replaces naive scalar differentiation, and the desired equality becomes an operator-norm inequality.When the relevant operators commute, the scalar quotient-rule identity is recovered; the inequality framework absorbs non-commutativity without an added hypothesis.Source: §2.2 · p. 5

3 The Matrix Halley Method: Construction and Geometric Meaning

The section constructs matrix Halley iteration and a matrix Schwarzian derivative, then proves local cubic convergence with a directional error constant. It also explains the method’s matrix-free advantage over Laguerre and gives an information-geometric interpretation.

3.1 The operator Halley method (a Hilbert-space instance of Döring’s construction)

  • The operator Halley step is a Hilbert-space instance of Döring’s construction, substituting the Newton direction into the Hessian derivative and requiring I + 1/2 L to be invertible.The paper warns that sign conventions differ across references and recommends using its convention consistently in implementation.Source: §3.1 · pp. 5–6

3.2 Auxiliary map and the matrix Schwarzian derivative

  • The first matrix Schwarzian definition uses the Hessian and third derivative; in the non-commuting case its relation to D2G is an inequality with a commutator correction.When the operators commute, the relation is the exact matrix analogue of the scalar identity.Source: §3.2 · p. 6

3.3 A lighter alternative: the third derivative of the iteration map

  • The matrix Schwarzian is extracted directly from the Newton map’s third Fréchet derivative, requiring neither an operator square root, self-adjointness, nor a gradient-field assumption.On the Section 6 test problem, numerical differentiation via a least-squares quintic polynomial fit agreed with predictions of (8) and (9) within maximum absolute errors of 4.7 × 10^-5 and 2.9 × 10^-4, respectively.Source: §3.3 · pp. 6–7

3.4 Main theorem: local cubic convergence of the matrix Halley method

  • Theorem 8 establishes genuine local cubic convergence near a simple root, without self-adjointness, commutativity, or gradient-field assumptions.The third derivative of the Halley map is governed by the matrix Schwarzian; a nonzero directional value gives exact cubic order, while a vanishing full form with F in C5 yields order at least four.Source: §3.4 · pp. 8–10

3.5 The matrix Laguerre family: why Halley admits a natural matrix-free realization

  • Laguerre solves the same quadratic Padé model directly and converges cubically, but requires a square root of a possibly non-self-adjoint operator; Halley linearizes the model instead.Halley’s linear formulation supports a matrix-free implementation using Hessian-vector products, whereas a fully matrix-free Laguerre implementation is not straightforward.Source: §3.5 · p. 10

3.6 An information-geometric interpretation

  • The information-geometric interpretation identifies the Hessian as a metric, the third derivative as the Amari–Chentsov skewness tensor, and Halley’s correction with the e-connection at α=1.In this account, Newton follows a flat 0-connection approximation and Halley corrects it using skewness.Source: §3.6 · p. 11

4 The Power-Newton Family and a Unifying Perspective

The section repairs a scalar power-Newton family to preserve roots, tests its cubic-order jump, and finds that a naive matrix lift fails. It concludes that a consistent multivariate correction of the power trick reduces to Halley’s method.

4.1 A root-preserving power-deformation family (corrected version)

  • The corrected family preserves the root (gp(x)=0 iff f(x)=0), and at p*(α)=0, Newton’s method applied to gp jumps to at least cubic convergence without evaluating f′′′.The family also satisfies g1=f and approaches ln u as p approaches zero.Source: §4.1 · pp. 11–12
  • For f(x)=e^x−2, p*=0.5 gives empirical order ≈3.2, while the other tested fixed exponents remain at ≈2.0.An adaptive exponent using current-iterate data reaches empirical order ≈2.9 as p_n approaches p*.Source: §4.1 · p. 12

4.2 A naive matrix extension, and why it fails

  • A uniform scalar exponent leaves the naive matrix extension at empirical order 2 because it cannot cancel an error term with mn degrees of freedom.The entrywise-adaptive exponent matrix is also nearly indistinguishable from Newton; its design assumes evaluation at the root, unlike use at the current iterate.Source: §4.2 · p. 12
  • The naive matrix lift uses an entrywise Box–Cox deformation, but its scalar parameter cannot generally eliminate the matrix-valued second-order error.The reported obstruction is a mismatch between mn error degrees of freedom and one scalar degree of freedom.Source: §4.2 · p. 12
  • The adaptive matrix-exponent design is structurally mismatched to implementation because its formula assumes the root while the nonlinear function is evaluated at a non-root iterate.The passage characterizes this as a fundamental limitation, not an implementation error.Source: §4.2 · p. 12

4.3 Unifying conclusion: Halley’s method as the multivariate limit of the power trick

  • A consistent multivariate version of the power-deformation correction requires the local, pointwise construction of Halley’s method itself.The section connects this correction to the skewness-tensor interpretation of the matrix Schwarzian derivative and Amari–Nagaoka’s α-connections.Source: §4.3 · pp. 12–13

5 A Matrix-Free Algorithm

The algorithm avoids explicitly forming the mn × mn operator by using directional products. It computes these products from Hessian-vector products and solves the resulting linear systems iteratively.

5 A Matrix-Free Algorithm

  • The matrix-free implementation needs Hessian-vector and third-order directional products rather than explicitly forming the full operator.The third-order product is obtained by differentiating the HVP computation directionally, without forming a matrix or tensor.Source: §5 · p. 13
  • Conjugate gradients solve the positive-definite Newton system, while MINRES handles the symmetric but possibly indefinite system.A Kronecker-factored approximation can precondition the iterative solves, and warm starts can reduce inner iterations.Source: §5 · p. 13
  • Damping is a practical refinement when the iteration operator risks becoming ill-conditioned or singular.The passage identifies regularization of the operator as the remedy.Source: §5 · p. 13

6 Numerical Experiments

The experiments test Newton, Halley, and related matrix methods on a problem with an analytically known root, then verify theoretical predictions and compare convergence behavior and computational trade-offs. Halley shows near-cubic convergence, while its advantage over Newton depends on the accuracy required.

6.1 Test problem

  • The test problem has a closed-form nontrivial root, allowing exact iteration-error evaluation from a randomly perturbed initial point.The root is simple, with A(X*) ≻ 0.Source: §6.1 · pp. 13–14

6.2 Implementation

  • The implementation supplies derivative operators without forming an mn × mn matrix, using CG for Newton directions and MINRES for Halley corrections.The inner tolerance is 10^-10.Source: §6.2 · p. 14

6.3 Results

  • Halley’s empirical order is ≈2.85 versus Newton’s ≈1.98; Halley reaches 8.70 × 10^-13 in its final iteration, while the plotted errors show cubic versus quadratic convergence.At ∼10^-9, Halley reaches in four iterations what Newton reaches in six; round-off prevents measuring the last asymptotic iterations.Source: §6.3 · p. 15, Figure 1 · p. 15

6.4 Discussion: which method is “better”?

  • At a 10^-10 error tolerance, Halley uses 4 outer iterations versus Newton’s 7, but total inner iterations are comparable at 16 versus 14.Halley is advantageous below 10^-10, whereas Newton is cheaper when a 10^-6 tolerance suffices; outer iteration counts alone do not determine cost.Source: §6.4 · pp. 15–16

6.5 Numerical verification of the matrix Schwarzian derivative (second definition)

  • The numerical values for all three verified quantities agree with their analytic predictions within finite-difference truncation error, supporting the square-root-free construction.The Newton map’s derivatives were extracted by fitting a degree-six least-squares polynomial.Source: §6.5 · p. 16, Table 2 · p. 16

6.6 The power-Newton family: success in the scalar case, failure for matrices

  • For the scalar case f(x) = e^x − 2, Table 3 shows that only the predicted critical exponent p* = 0.5 pushes empirical order toward 3.Source: §6.6 · p. 16
  • In the matrix test, varying a uniform scalar p leaves empirical order near 2, and the naive adaptive exponent matrix nearly coincides with Newton rather than showing Halley’s cubic convergence.The reported obstruction is that a scalar parameter cannot generally cancel a second-order error term with mn degrees of freedom.Source: §6.6 · p. 17, Table 4 · p. 17, Figure 2 · p. 17

6.7 Numerical verification of the matrix Laguerre family

  • On the matrix test, Halley and Laguerre converge cubically, but Laguerre requires a genuine operator square root at each iteration.The required square root is of a 30 × 30 matrix and is nontrivial, with the reported ratio approximately 0.89.Source: §6.7 · pp. 17–18, Figure 3 · p. 18

7 Case Study: A Generalized Oja-Type Matrix Dynamical System

The Oja-type case study shows that convergence depends on the eigenvalue gaps used by the target root, and that cubic methods gain more over Newton as ill-conditioning increases.

7.1 Problem setup

  • The target root is the stable equilibrium assigning the largest eigenvalues in decreasing order, while permuting assignments or flipping eigenvector signs yields many other roots.The iterations are algebraic constructions and apply to this flow without passing through its twisted metric.Source: §7.1 · p. 19

7.2 Experiment I: B = diag(3, 2, 1), X ∈R5×3 (only the top three eigenvalues are used)

  • For the top-three target, Newton converged at order ≈2.0 and Halley/Laguerre at ≈2.7–3.0, with almost no difference between A1 and A2.The used eigenvalue gaps were comparable, although A2’s global condition number was 3.0 × 10^8.Source: §7.2 · p. 20, Table 6 · p. 20, Figure 4 · p. 20
  • Across 15 distant random starts, all three methods usually reached some root, but rarely the intended root; A2’s condition number barely affected this global robustness.The system’s large number of roots helps explain why reaching the particular target is uncommon.Source: §7.1 · p. 19, §7.2 · pp. 20–21, Table 7 · p. 21

7.3 Experiment II: B = diag(5, 4, 3, 2, 1), X ∈R5×5 (all five eigenvalues are used)

  • Using all five eigenvalues makes the problematic gaps part of the target, worsening its Jacobian conditioning and slowing convergence in the corresponding subspace.The well-separated top-three block reached machine precision within a few iterations, while Newton needed 18–19 iterations and Halley/Laguerre 8–9 for the remaining block.Source: §7.3 · pp. 21–22

7.4 Discussion

  • When a near-zero eigenvalue makes an eigenvector’s sign numerically indeterminate, convergence should be diagnosed with ∥F(Xk)∥ rather than sign-dependent distance to X∗.The apparent stagnation corresponded to a sign-flipped, equally valid root, not an algorithmic defect.Source: §7.3 · p. 21, §7.4 · p. 22

7.3 Experiment II: B = diag(5, 4, 3, 2, 1), X ∈R5×5 (all five eigenvalues are used)

  • With all five eigenvalues used, Newton required 19 iterations while Halley/Laguerre converged in 6–9, making the cubic advantage most pronounced under extreme ill-conditioning.The slow convergence was localized to the subspace associated with the near-degenerate or near-zero eigenvalues.Source: §7.3 · pp. 21–22

8 Verification at 70-Digit Precision

At 70-digit precision, the experiments verify the predicted convergence orders and distinguish genuine sign bifurcation from rounding artifacts. They also show that attainable accuracy is limited by working precision and eigenvalue conditioning.

8.1 Tool: mpmath

  • The 70-digit runs used mpmath, and the 25 × 25 operator square-root iteration was the most computationally expensive part of the experiments.The general matrix square root required for Laguerre was implemented with Denman–Beavers iteration.Source: §8.1 · p. 23

8.2 Precise verification of the convergence order for A1, A2

  • At 70-digit precision, errors decreased past the double-precision floor of 10^-16 to around 10^-64.The experiments used the same all-eigenvalues setting as Experiment II and a 10^-6 initial perturbation of A2.Source: §8.2 · p. 23, Figure 6 · p. 24
  • Newton’s empirical order is exactly 2 and Halley/Laguerre’s is exactly 3 in the asymptotic regime, matching theory to the displayed precision.The higher precision reveals behavior that double precision could confirm only approximately.Source: §8 · p. 23, §8.2 · p. 24, Table 9 · p. 24

8.3 Confirming that the sign-bifurcation phenomenon is genuine

  • Repeating the same perturbation at 70-digit precision reproduced convergence to a sign-flipped root, establishing that the bifurcation is precision-independent.The flipped column corresponds to the eigenvalue ≈10^-8, whose sign is intrinsically indeterminate in this setting.Source: §8.3 · p. 25, Figure 8 · p. 25

8.4 Matrices A3, A4 with an exactly controlled eigenvalue gap

  • Because A1 and A2 specify their small gaps only to 10–11 digits, the authors constructed A3 and A4 with exactly controlled eigenvalues using 70-digit Givens rotations.A3 has a gap of 10^-20; A4 has a smallest eigenvalue of 10^-25.Source: §8.4 · pp. 25–26

8.5 The precision-budget plateau phenomenon

  • For A3 and A4, Newton and Halley plateaued at accuracy predicted by working precision minus condition-number digits: about 10^-50 and 10^-47, respectively.Both reached the same precision limit; Halley took one fewer iteration for A3 and one versus two iterations for Newton on A4.Source: §8.5 · pp. 26–27, Figure 9 · p. 27

9 Conclusion and Future Work

The conclusion consolidates the matrix Halley theory, its implementation and numerical validation, and its application to an Oja-type system. It also bounds the practical niche of coupled formulations and identifies directions for further work.

9 Conclusion and Future Work

  • The matrix Schwarzian is interpreted through information geometry, and the Newton-map third derivative supplies a square-root-free construction used in the local cubic-convergence theorem.The theorem gives an explicit asymptotic error constant and requires neither self-adjointness, commutativity, nor a gradient-field assumption.Source: §9 · p. 28
  • A matrix-free algorithm avoids forming an mn × mn matrix; the paper also contrasts Halley with a power-Newton family and a Laguerre family requiring an operator square root.The naive matrix extension of power-Newton fails, while Halley is presented as the multivariate limit of the power-deformation trick.Source: §9 · p. 28
  • The Oja-type case study links local convergence quality to eigenvalue gaps and finds larger practical gains from cubic convergence under ill-conditioning.The 70-digit experiments confirm the theoretical orders and identify a precision-budget plateau.Source: §9 · p. 28
  • Coupled Halley formulations are most practically valuable for genuinely coupled nonlinear problems; in spectral problems, sequential deflation can avoid catastrophic divergence and reduce cost.The reported coupled formulation diverged on a 100 × 100 Hilbert-matrix example with 40 columns, whereas sequential deflation approached the precision limit gracefully.Source: §9 · pp. 28–29
  • Future work includes quantifying non-commutativity corrections, connecting to K-FAC/Shampoo, developing higher-order methods from α-connections, and approximating Laguerre’s operator square root matrix-free.The conclusion also proposes investigating value-distribution constraints for the matrix dynamical systems.Source: §9 · p. 29
Loading 2609.06592v1…