PaperlayerView exact revision on arXiv

Source-linked AI summary

Efficient computation of the asymptotics of extensive-rank HCIZ integrals

Antoine Maillard, Jean-Christophe Mourrat

arXiv:2609.07771v1math.PRcond-mat.dis-nncs.NAmath-phmath.MPmath.NA

TL;DR

In the extensive-rank regime, HCIZ asymptotics are linked to a hydrodynamical boundary-value problem, while behavior as one density’s rank vanishes remains an open question. The paper develops a particle discretization for computing the asymptotic value and establishes its validity under stated boundary-density conditions.

  • Problem

    The behavior of hydrodynamical trajectories when one density has vanishingly small rank remains an interesting question in the transition between low-rank and extensive-rank HCIZ regimes.

    Source: §1.1.3 · p. 5

  • Method

    The paper designs a particle-based discretization to compute J(µ0, µ1) and establishes its validity for compactly supported absolutely continuous boundary measures with L3 densities.

    Source: §1.2 · p. 7, §1.2.2 · p. 8

  • Results

    The discretization asymptotically approaches the HCIZ value for boundary measures satisfying the theorem’s stated conditions.

    Source: §1.2.2 · p. 8

  • Takeaways & Limitations

    Nearest-neighbor spacings make the discretization local in space, simplifying its implementation and analysis.

    Source: §1.3 · p. 12

    The theorem’s L3 boundary-density assumption may be weakenable, but the tradeoffs involved in regularizing singular densities remain unquantified.

    Source: §1.3 · p. 12

Abstract

from arXiv · show

We study the high-dimensional asymptotics of Harish-Chandra-Itzykson-Zuber (HCIZ) integrals in the extensive-rank regime. The limit of these integrals is governed by a one-dimensional boundary-value hydrodynamical problem originally derived by Matytsin (1994) and rigorously proved by Guionnet and Zeitouni (2002). Despite its wide-ranging applications, explicit solutions to this problem are known only in a few specific cases. In this work, we introduce an efficient numerical scheme based on a particle discretization and prove its convergence to the continuous boundary-value problem for generic boundary densities. We validate our approach against known analytical solutions and apply it to generic densities, uncovering interesting dynamical phenomena. The high-dimensional limit of HCIZ integrals appears in various contexts, from the large deviations of random matrix spectra to the limiting free energy of disordered systems, high-dimensional statistics, and machine learning. As such, our contribution opens the way towards the numerical exploration of a wide range of high-dimensional models that were previously intractable.

1 Introduction and main results

The paper develops a convergent particle-based numerical scheme for extensive-rank HCIZ asymptotics governed by Matytsin’s hydrodynamical boundary-value problem, targeting generic boundary densities. It compares the scheme with known analytical predictions and explores trajectories and methodological extensions, including support merging and singular boundaries.

1.1.1 Notation and definition

  • The introduction defines the spherical Harish-Chandra–Itzykson–Zuber integral for matrix pairs and retains a scaling parameter θ that can be varied.The integral is introduced with references.Source: §1.1.1 · p. 2

1.1.2 Related literature

  • High-dimensional HCIZ asymptotics arise in random-matrix spectral large deviations, disordered-system and spin-glass free energies, and high-dimensional statistics and machine learning.The extensive-rank regime is linked to growing-rank spectral problems and to high-rank matrix denoising, factorization, and neural-network learning dynamics.Source: §1.1.2 · p. 3

1.1.3 Matytsin’s solution

  • Under weak convergence and common compact support assumptions, the extensive-rank HCIZ limit is independent of β and continuous in its limiting boundary measures; with finite logarithmic energies, its variational infimum has a unique path that is smooth in the interior.The extensive-rank regime is governed by the hydrodynamical boundary-value problem derived by Matytsin and proved by Guionnet and Zeitouni.Source: §1.1.2 · p. 3, §1.1.3 · p. 4
  • Because generic boundary densities remain largely intractable analytically, the authors introduce a convergent numerical scheme and aim to enable exploration of previously intractable high-dimensional models.They benchmark the algorithm against known analytical cases and anticipate applications beyond Bayes-optimal estimation or extreme-temperature free energies.Source: §1.1.3 · pp. 6–7

1.2.2 Validity of the discretization scheme

  • The particle discretization is proved to approximate the HCIZ asymptotic functional to arbitrary accuracy, with validity established for compactly supported absolutely continuous boundary measures with L3 densities.The implemented scheme additionally discretizes time; the paper states that this time-step discretization converges to the particle problem as the time steps shrink.Source: §1.2 · p. 7, §1.2.2 · p. 8

1.2.1 The discretization scheme

  • The proposed approximation is particle-based: it uses N ordered particle trajectories with continuous time to compute J(µ0, µ1), under finite-log-energy endpoint measures.The particles have fixed endpoint positions and satisfy ordering and trajectory constraints.Source: §1.2.1 · p. 7

1.2.3 First applications

  • For the free-convolution benchmark, the paper compares particle samples with analytical density and velocity predictions across times for a Marchenko–Pastur input with ratio κ = 2.The plotted computation uses N = 212 particles and T = 210 time steps.Source: §1.2.3 · p. 8, Figure 1 · p. 8
  • The solver agrees with free-convolution predictions and semicircle-law results across large and small scaling parameters, even at relatively small N and T.The semicircle comparison covers a wide range of θ values and different (N, T).Source: §1.2.3 · p. 9, Figure 2 · p. 9
  • The solver also performs well for a boundary density outside the stated L3 assumption, while regularization can extend it to singular measures, including one with an atom of mass 1/2.The paper notes that the L3 condition may be weakenable and calls for more quantitative study of regularization tradeoffs.Source: §1.2.3 · pp. 9, 11, Figure 5 · p. 11, §1.3 · p. 12
  • For generic boundary measures, the solver computes trajectories showing initially disconnected support merging near t ≃ 0.08 and singular velocity behavior, where no analytical solution is known.The example uses symmetrized Marchenko–Pastur initial data and a semicircle terminal measure.Source: §1.2.3 · pp. 9, 11, Figure 4 · p. 10

1.3 Discussion and consequences

  • The discretization makes the scheme local in space by approximating density with nearest-neighbor spacings while preserving the coefficient π^2/3.The authors explain this choice heuristically through the dominance of singular near-neighbor interactions in the continuum limit [18].Source: §1.3 · pp. 11–12
  • The numerical results show support components merging but do not determine whether velocity remains continuous then, leaving solution regularity and convergence rates open.The authors identify regularity near critical times as a likely ingredient in sharp quantitative convergence estimates.Source: §1.3 · p. 12
  • A free-boundary extension would optimize the initial particle positions while fixing the terminal configuration, with discretized boundary terms added to the objective.The proposed class of problems includes spectral large deviations, disordered-system free energies, and high-dimensional statistics and machine learning.Source: §1.3 · p. 12
  • Direct shooting is expected to be unstable, characteristics do not cover general boundary densities, and Eulerian grid methods may struggle with nonsmooth densities and velocity fields.The authors present Eulerian optimization as a promising alternative and leave its implementation and quantitative comparison for future work.Source: §1.3 · p. 13

2 Efficient minimization algorithm

The section develops a particle-based convex minimization algorithm with spacing variables, tailored time discretization, and preconditioning. It reports convergence rates for the discretized objective and describes the Newton–conjugate-gradient implementation.

2.1 Parametrization by the spacings

  • Rescaled adjacent-particle gaps simplify the linear constraints to nonnegativity, while optimizing first-particle positions becomes an analytically solvable quadratic problem.The objective is extended to unrestricted gaps using a convex continuation; checking nonnegative gaps at the output ensures the minimum is the original global minimum.Source: §2.1 · pp. 15–16

2.2 Time discretization

  • A nonlinear time grid, symmetric about the midpoint, is motivated by the free-Brownian-bridge description and its short-time additive free-convolution approximation.The paper notes a trapezoidal discretization with a global time-discretization error of O(1/T^2).Source: §2 · p. 14, §2.2 · p. 16

2.3.1 An ill-conditioned minimization problem

  • The spacing-variable Hessian is expected to have condition number Θ(T^2); velocity variables still yield Θ(N^2) when spacings stay order one and N≳T.Velocity parametrization also requires handling non-crossing constraints, which motivates retaining spacing variables.Source: §2.3.1 · p. 17

2.3.2 Preconditioning procedure

  • To address ill-conditioning, the preconditioner approximates the inverse Hessian and solves its associated linear system in O(NT log N) operations.The solve uses independent positive-definite tridiagonal systems for spatial modes and a discrete sine transform.Source: §2.3.2 · p. 18

2.4 Global scaling and the θ parameter

  • Keeping the global scaling parameter θ explicit can avoid numerical issues at θ≪1 or θ≫1; the rescaled time interval is (0, θ^-1).The variational principle and minimization algorithm are adapted to this rescaling.Source: §2.4 · p. 18

2.5 Initialization, and final algorithm

  • The algorithm initializes spacings by linearly interpolating boundary quantiles, then uses Newton updates with PCG and stops when the Newton decrement falls below tolerance.Reported implementations ran on a single Nvidia GPU, with instance times from a few seconds to a few hours; the authors found Newton–conjugate gradient outperformed Adam and LBFGS empirically, while quantitative comparisons remain future work.Source: §2.5 · p. 19, §2.6 · p. 20

2 Efficient minimization algorithm

  • Newton with preconditioned conjugate gradient uses Hessian-vector products costing O(NT) and preconditioner solves costing O(NT log N) per iteration.The iteration cost is governed primarily by the number of PCG iterations and the preconditioner solve.Source: §2.3.2 · p. 18, §2.5 · p. 19

2.6 Practical implementation and final remarks

  • At fixed N, objective values generally converge as O(T^-2); at fixed T, error is O(N^-1), improvable by Richardson extrapolation.The O(1/N) accuracy is commensurate with the discretization of boundary measures and quantile-based linear-statistic evaluation.Source: §2.6 · p. 20

3 Proof of Theorem 1.2

This section proves Theorem 1.2 by establishing convergence through regularization, preliminary estimates, and upper and lower bounds. It also outlines a conjectural alternative formulation of the limiting functional using inverse-CDF interpolants.

3 Proof of Theorem 1.2

  • For finite-Σ endpoint measures, G has a unique minimizer whose intermediate-time densities are supported on a common compact set [49].Uniqueness is understood up to equality of (ρ, ρv) almost everywhere.Source: §3 · p. 20

3.1 Small-regularization limit

  • The small-regularization limit establishes Jε(µ0, µ1) → J(µ0, µ1) as ε → 0 by passing from approximate minimizers to a continuum candidate.The limiting argument uses convexity of G in (ρ, ρv) and obtains a candidate satisfying the continuity equation with no greater limiting cost.Source: §3.1 · pp. 21, 23

3.2 Upper bound

  • For finite-Σ endpoint measures, the discrete-to-continuous construction gives an upper bound on J by the liminf of the discrete costs as N → ∞.The proof takes N → ∞ and ε → 0 together, since the constructed paths do not have the exact initial measure at finite N.Source: §3.2 · pp. 23, 25

3.3 Lower bound

  • For absolutely continuous endpoints with L3 densities, Proposition 3.4 establishes the converse bound for every ε > 0.Its proof smooths near-optimal continuum paths and bridges the endpoint discrepancy caused by convolution.Source: §3.3 · p. 25, §3.3.1 · p. 28

3.3.1 Proof strategy and preliminary results

  • The lower-bound strategy uses moving quantile particles for regular velocity fields, while smoothing is needed because the continuum minimizer’s velocity need not be regular near the density support boundary.Proposition 3.5 supplies the quantile-particle construction; the proof smooths a near-optimal pair and bridges the endpoints after convolution.Source: §3.3.1 · pp. 26–28
  • The preliminary proof estimates culminate in the announced result, using bounded-variation integrability and bounds involving the Wasserstein distance and Jensen’s inequality.The inverse cumulative distribution functions are bounded and of bounded variation, making their product Riemann-integrable; the subsequent estimates control the remaining terms.Source: §3.3.1 · pp. 28–29

3.3.2 Proof of Proposition 3.4

  • A regularized trajectory with a uniform density lower bound yields particle paths whose approximation error vanishes as N grows, followed by limits in the regularization and endpoint parameters.The construction smooths the trajectory, adjusts its endpoints by bridging short time intervals, and takes limits in the order N, δ, then η.Source: §3.3.2 · pp. 30–31

3.4 Remark on an alternative formulation

  • The authors speculate that J(µ0, µ1) may admit an alternative formulation over continuous inverse-CDF interpolants that are non-decreasing in their second argument and differentiable in the interior.The proposed admissible functions match the endpoint inverse CDFs; the speculation is stated for compactly supported endpoint measures with finite Σ.Source: §3.4 · pp. 31–32

A.1 Scaling and invariance properties of the HCIZ integral

The HCIZ hydrodynamical system admits a scaling reformulation for θ ≥ 0. Rescaling can account for θ through either interpolation time or the space variable.

A.1 Scaling and invariance properties of the HCIZ integral

  • The section presents scaling and invariance properties of the HCIZ integral and its hydrodynamical system.Source: §A.1 · p. 37

A.1.1 Scaling parameter θ ≥0

  • After rescaling, the endpoint measures are µ at u = 0 and ν at u = θ^-1.Source: §A.1.1 · p. 37
  • The scaling parameter θ ≥ 0 can be incorporated by rescaling either interpolation time or the space variable.The continuity equation remains satisfied under the stated rescaling, with changed boundary conditions.Source: §A.1.1 · p. 37

A.1.2 Galilean transformation

The hydrodynamical equations are invariant under Galilean transformations. Applying the transformation to the terminal measure recovers a property of the finite-dimensional HCIZ integral.

A.1.2 Galilean transformation

  • The hydrodynamical equations remain satisfied under Galilean transformations.Source: §A.1.2 · p. 37
  • Applying the transformation to the terminal measure ν yields the shifted measure να(x) := ν(x − α).Source: §A.1.2 · p. 37
  • The resulting property follows directly from the finite-dimensional expression of the HCIZ integral.Source: §A.1.2 · p. 37

A.2 Remarks on time discretization

The time discretization error is analyzed using smooth trajectories and the trapezoidal rule. Under the stated regularity conditions, the discrete kinetic energy converges quadratically in the time-step size.

A.2 Remarks on time discretization

  • The time discretization error is O(1/T^2), as validated by numerical experiments in Appendix B.Source: §A.2 · p. 38
  • For a C^2 trajectory and a twice continuously differentiable function near its path, the trapezoidal rule gives the stated convergence order for the discretized kinetic-energy term.The bound's constant depends only on the indicated sup norms of the trajectory and function derivatives.Source: §A.2 · p. 38
  • As the maximum time-step size tends to zero, the discrete kinetic energy converges to its continuous counterpart with error O(||Δt||^2).Source: §A.2 · p. 38
  • The proof derives the error bound by expanding the square, applying Fourier expansion and Parseval's identity, and using n^2 ≥ 1.Source: §A.2 · p. 38

A.3 Derivation of eq. (39)

The derivation reduces a step to a convex quadratic minimization and uses its stationarity conditions to characterize the minimizer. It then completes the derivation of eq. (94) using the relation between M and L.

A.3 Derivation of eq. (39)

  • The step is a simple convex quadratic minimization problem.Source: §A.3 · p. 39
  • Stationarity with respect to x1,k is imposed for each interior time index k.Source: §A.3 · p. 39
  • The resulting expression also holds for the fixed boundary variables at k = 0 and k = T.Source: §A.3 · p. 39
  • The derivation of eq. (94) is completed using M = L^-1 N^-1.Source: §A.3 · p. 40

A.4 Minimization algorithm with a rescaling parameter θ

The rescaled discretization uses a time mesh designed to resolve dynamics on two timescales, while retaining the same constrained minimization objective. The appendix also states the assumptions behind the preconditioned conjugate-gradient procedure.

A.4 Minimization algorithm with a rescaling parameter θ

  • The rescaled objective is minimized over δi,k subject to δi,k ≥ 0 and the original boundary conditions.Source: §A.4 · p. 40
  • The time schedule uses the free Brownian bridge law to motivate its discretization over t ∈ [0, 1/2].Source: §A.4 · p. 40
  • For θ ≪ 1, the evolution has timescales t ∼ θ and t = Θ(1), motivating extra time-grid resolution near t = 0.The discretization maps mesh points through gθ, which is chosen to capture both scales.Source: §A.4 · p. 40
  • The rescaled algorithm generalizes the Hessian and its preconditioner to the nonuniform time steps.Source: §A.4 · pp. 40–41

A.5 Preconditioned conjugate gradient

  • PCG seeks a preconditioner that reduces the system’s condition number while keeping preconditioner solves inexpensive.Source: §A.5 · pp. 41–42

B Runtime and convergence analysis experiments

The section reports experiments on total runtime and convergence properties of the algorithm.

B Runtime and convergence analysis experiments

B.1 Runtime analysis

Runtime experiments find favorable scaling with T and different scaling with N for the two tested values of θ. Joint scaling is summarized by a single parameter for each θ, while PCG iteration counts clarify how conditioning varies with N and T.

B.1 Runtime analysis

  • For fixed N, runtime grows more slowly than O(T), although each iteration requires O(T) operations.Source: §B.1 · p. 42, Figure 7 · p. 43
  • For fixed T, runtime scales as O(N) at θ = 1 and approximately O(N^3/2) at θ = 0.01.The passage attributes the latter pattern as suggesting more PCG iterations with N and a more ill-conditioned problem.Source: §B.1 · p. 42, Figure 8 · p. 43
  • Joint runtime is represented relatively well by NT for θ = 1 and N^3/2T for θ = 0.01.Source: §B.1 · p. 42, Figure 9 · p. 43
  • PCG iteration counts appear independent of T; with increasing N they are roughly constant for θ = 1 and scale approximately as N^1/2 for θ = 0.01.Each PCG update requires O(NT log N) operations, so iteration counts provide a rough runtime estimate.Source: §B.1 · p. 44, Figure 10 · p. 44

B.2 Convergence analysis

The convergence study tests two boundary-measure settings and finds evidence that the leading time and particle discretization errors have different orders and are essentially decoupled. Their contributions have opposite signs.

B.2 Convergence analysis

  • For (C1), µ = ν = σs.c. with θ = 0.01 and J is known analytically; for (C2), J is unknown and rates are estimated from successive differences.(C2) uses the symmetrized Marchenko-Pastur law with κ = 2 and ν = σs.c.,0.5, with θ = 1.Source: §B.2 · p. 44
  • The leading time and particle discretization errors appear to have different orders and be essentially decoupled.The experiments are consistent with the stated error model.Source: §B.2 · p. 44
  • Time discretization overestimates J, while particle discretization underestimates it.The positive coefficients a and b depend on (θ, µ, ν), not on (N, T).Source: §B.2 · p. 44

B.2.1 Convergence as T →∞and ∥∆t∥∞→0

At fixed particle count, the study examines time-discretization convergence in two boundary-measure settings. The time error is compatible with a 1/T^2 decay and is only weakly dependent on N.

B.2.1 Convergence as T →∞and ∥∆t∥∞→0

  • Figure 11 presents the convergence analysis as the number of time-discretization steps T grows for (C1).The supplied caption does not state a specific outcome for this figure.Source: Figure 11 · p. 45
  • For (C1), the time-discretization error is essentially N-independent, with a least-squares fit consistent with 1/T^2 decay.The signed total error changes sign at an N-dependent T as the time error gives way to the negative finite-N bias.Source: §B.2.1 · p. 45
  • For (C2), Figure 12 tracks objective values and successive differences as T grows; the objective is within ∼10^-4 of its large-T limit by T ≃ 30.The analysis therefore restricts T to at most 512.Source: Figure 12 · p. 45
  • The time-discretization error is compatible with a 1/T^2 decay and is only weakly dependent on N.The objective is already very close to its large-T limit for moderate T ∼ 10^2, so accuracy is often limited by N rather than T.Source: §B.2.1 · p. 45

B.2.2 Convergence as N →∞

For both boundary-measure settings, the finite-particle error is compatible with 1/N decay and is nearly independent of the time-grid size T.

B.2.2 Convergence as N →∞

  • Figure 13 compares objective values with the exact asymptotic value and plots successive differences and signed total error for (C1).Its fitted slope is compatible with a 1/N leading-order decay.Source: Figure 13 · p. 46
  • Figure 14 plots objective values and successive differences as N grows for (C2); the fitted slope is compatible with a 1/N decay.Successive differences for different T values are superimposed.Source: Figure 14 · p. 46
  • The finite-N error is compatible with a 1/N decay in both cases and is nearly independent of T.This is consistent with the predicted error model in eq. (103).Source: §B.2.2 · p. 46, Figure 13 · p. 46, Figure 14 · p. 46

B.2.3 Richardson extrapolation in N

Richardson extrapolation is applied to reduce the finite-N error in both test settings. The extrapolated estimates reach 10^-4 error with substantially fewer particles in (C1).

B.2.3 Richardson extrapolation in N

  • Richardson extrapolation eliminates the leading 1/N term at fixed T.The procedure is motivated by the observed 1/N decay and is applied to both (C1) and (C2).Source: §B.2.3 · pp. 46–47
  • Figures 15 and 16 compare raw and extrapolated objective estimates and their successive differences for (C1) and (C2), respectively.Both figures show estimates at successive particle counts.Source: Figure 15 · p. 47, Figure 16 · p. 47
  • Reaching an error of 10^-4 in (C1) requires N ≃ 10^4 for the raw sequence versus N ≃ 10^3 for the extrapolated sequence.The fitted exponent is α ∈ [1.6, 1.7].Source: §B.2.3 · p. 47

Related papers

Loading 2609.07771v1…