Matrix Integral Bootstrap Techniques

• Matrix Integral Bootstrap is a rigorous method that uses loop (Schwinger–Dyson) equations and moment positivity to bound observables in large‑N matrix models.
• It employs a relaxation strategy via semidefinite programming to convert non-convex polynomial constraints into tractable convex problems.
• Algorithmic innovations and symmetry-based block-diagonalization enable high-precision computation of eigenvalue distributions and phase structures.

The matrix integral bootstrap is a rigorous, semi-algorithmic approach to solving large-$N$ matrix models by combining loop (Schwinger–Dyson) equations with positive semidefiniteness constraints on moment matrices. This methodology systematically bounds observables, reconstructs eigenvalue distributions, and can in many cases achieve higher precision than traditional techniques such as Monte Carlo, especially in models where analytic solutions are unavailable or where measure positivity enables exact optimizations.

1. Fundamental Principles and Loop Equations

In matrix models, observables of interest are typically normalized single-trace moments, e.g., $m_k = \langle (1/N)\,\mathrm{Tr}\,M^k \rangle$ for Hermitian matrices. The central objects are large-$N$ (planar) limits of polynomial matrix integrals, which admit a factorization property wherein products of traces factorize up to corrections of order $1/N^2$. The Schwinger–Dyson equations, sometimes called loop equations, arise from invariance of the matrix integral under suitable infinitesimal shifts: $$0 = \int dM\,\partial_M (M^k\,e^{-N\,\mathrm{Tr}\,V(M)}) \implies \langle\, (1/N)\,\mathrm{Tr}\,V'(M)\,M^k\,\rangle = \sum_{\ell=0}^{k-1} m_\ell\,m_{k-1-\ell}$$ These equations recursively relate all higher moments to lower ones and closures are enabled in the large-$N$ limit by factorization, reducing the system to a finite-dimensional, often nonlinear, set of polynomial constraints among the moments. For multi-matrix systems, similar loop equations are derived for words built from several matrices, closing after accounting for all operator combinations up to a prescribed cutoff length (Kazakov et al., 2021, Li et al., 8 Oct 2025).

2. Positivity Constraints and Moment Matrices

The bootstrap imposes, in addition to the loop equations, positivity constraints derived from the measure positivity of the matrix ensemble. For Hermitian matrices this means every real linear combination $f = \sum_j v_j M^j$ yields

$\langle \mathrm{Tr}\,f^\dagger f \rangle \geq 0$

This translates to the requirement that all finite truncations of the Hankel matrix, $C^{(K)}_{ij} = m_{i+j-2}$, are positive semidefinite (PSD): $C^{(K)} \succeq 0$. For multi-matrix systems, the analog is a tracial moment matrix labeled by words, also required to be PSD. For models with $O(D)$ symmetry, block-diagonalization with respect to this symmetry dramatically reduces the computational cost, as only small blocks corresponding to irreducible representations must be checked for positivity (Kazakov et al., 2021, Li et al., 8 Oct 2025, Lin et al., 28 Jul 2025).

A key result is that, in the infinite cutoff ($K \to \infty$), positivity of the moment matrix is mathematically equivalent (by the Hamburger moment problem) to the existence of a positive eigenvalue measure: the resolvent

$G(z) = \sum_{k=0}^\infty m_k z^{-k-1}$

must be the Stieltjes transform of a measure supported on the real line with non-negative weight. Therefore, enforcing positivity at the matrix level ensures the physicality of reconstructed distributions (Kazakov et al., 2021, Kováčik et al., 19 Sep 2025).

3. Relaxation Bootstrap and Semidefinite Programming Formulations

The nonlinearity of the large-$N$ loop equations, stemming from bilinear terms $m_\ell m_m$, renders the system non-convex. The relaxation bootstrap replaces these products by lifting to new matrix variables $X_{ij} \geq m_i m_j$, introducing an auxiliary “relaxation” block: $$\begin{pmatrix} 1 & x^\top \ x & X \end{pmatrix} \succeq 0$$ This converts the search for feasible moments into a convex semidefinite program (SDP) in the variables $(x, X)$ subject to linearized loop equations and positivity constraints on both the truncated moment matrix and the relaxation block. This framework enables global bounds for any linear observable, e.g., $m_2$, and systematically converges to the true solution as the cutoff increases (Kazakov et al., 2021, Lin et al., 28 Jul 2025, Li et al., 8 Oct 2025).

For multi-matrix or quantum mechanical systems with extensive symmetry, operator basis reduction and SDP block-diagonalization allow practical computations for cutoffs up to length 12–14. At each cutoff level, the allowed region in the space of low moments shrinks rapidly, often producing compact “islands” that converge to a unique solution at large cutoff (Li et al., 8 Oct 2025, Lin et al., 28 Jul 2025).

4. Analytic–Trajectory and Reconstruction Techniques

Parallel to the moment positivity approach, analytic-trajectory bootstraps posit that families of correlators fit explicit analytic patterns, especially in high-symmetry models. For example, in bosonic Yang–Mills integrals, one-length and multi-length trace patterns are captured by shifted and rescaled Gaussian moments, parameterized by a small set of unknowns (such as $z_{\mathrm{max}}$ for the eigenvalue endpoint). Loop equations and group-theoretic constraints fix these ansatz coefficients, yielding closed-form approximations for all moments, remarkably accurate even at moderate $D$ (Li et al., 8 Oct 2025).

Once moments are consistently estimated, one reconstructs the underlying eigenvalue probability distribution using either Stieltjes inversion of the resolvent or orthogonal polynomial expansions. For single-matrix cases, the orthogonal polynomial method yields the eigenvalue density as a finite series in determined moments, approaching exponential accuracy in the large cutoff limit. Physical observables and the free energy are then computed directly from the reconstructed density (Kováčik et al., 19 Sep 2025).

5. Phase Structure and Free-Energy Selection

While the bootstrap identifies all mathematically consistent sets of moments, physical selection among possible solutions (such as phases with different symmetry breaking patterns or multi-cut behaviors) is performed through free-energy minimization. The reconstructed eigenvalue density is used to compute the large-$N$ free energy functional: $$F[\rho] = \int V(\lambda)\,\rho(\lambda)\,d\lambda - \iint \rho(\lambda) \rho(\mu) \log|\lambda-\mu|\, d\lambda\, d\mu$$ The phase with minimal $F$ is thermodynamically preferred. This approach discriminates among multiple convex islands in bootstrap space, mapping out phase diagrams, first-order transitions, and coexistence regions, as realized concretely in multi-trace quartic models (Kováčik et al., 19 Sep 2025).

6. Algorithmic Innovations and High-Precision Results

The practical success of the matrix integral bootstrap relies on algorithmic advances in SDP solvers, moment matrix reduction, and relaxation schemes. In quantum mechanical and high-symmetry multi-matrix systems, block-diagonalization in symmetry irreps enables high cutoffs ($L = 11$–$14$), producing 6–8 digit precision in observables. For matrix models on the circle (unitary ensembles), the "shoestring bootstrap" bypasses SDP entirely by leveraging the rapid decay of Fourier moments, yielding hundreds of digits of precision for small $N$ (Berenstein et al., 28 Feb 2025).

Comparisons with Monte Carlo simulations indicate that the bootstrap routinely achieves higher precision in the large-$N$ limit, especially as cutoff increases. The critical and singular behavior, e.g., multi-cut transitions, spontaneous symmetry breaking, and the emergence of nontrivial saddle points, is robustly and systematically captured (Kazakov et al., 2021, Lin et al., 28 Jul 2025, Li et al., 8 Oct 2025).

7. Extensions and Significance

The matrix integral bootstrap generalizes naturally to Dirac-ensemble models, noncommutative-geometric settings, and models with complex or non-positive weights provided appropriate analytic input can be supplied. It is broadly applicable to bosonic multi-matrix models, Yang–Mills reductions, and even quantum mechanical systems obtained by dimensional reduction of gauge theory. The method not only serves as a rigorous alternative to Monte Carlo for large-$N$ integrals but also provides a conceptual bridge between modern convex optimization, representation theory, and non-perturbative field theory (Kazakov et al., 2021, Li et al., 8 Oct 2025, Hessam et al., 2021, Lin et al., 28 Jul 2025, Berenstein et al., 28 Feb 2025).

Key insights include the mathematical equivalence of moment positivity and positive measures, the efficiency of block-diagonalization for exploiting symmetry, the role of analytic trajectories in circumventing sign problems, and the fusion of algebraic and convex-analytic constraints for robust, high-precision determination of nontrivial observables in matrix models.