¹

1 Introduction↩︎

Random-matrix theory (RMT) provides universal descriptions of spectra and eigenvectors in complex quantum systems, ranging from nuclear spectra to quantum chaotic dynamics and disordered many-body systems [1]–[5]. In the standard forward formulation, one samples a matrix from an ensemble and then studies its eigenvalues, eigenvectors, level spacings, or spectral form factors. In many microscopic problems, however, matrix elements are not independent primitive variables. They arise from overlaps of hidden amplitudes, projected wave functions, effective modes, or constrained internal degrees of freedom. This motivates the inverse question addressed here: given a desired set of RMT-like element statistics, can one reconstruct a structured set of vectors whose overlaps generate those statistics?

This question is distinct from the usual diagnosis of chaos through Lyapunov exponents, out-of-time-order correlators, or sensitive dependence on initial conditions. The present construction is static and algebraic. It does not compare two nearby dynamical trajectories, nor does it establish exponential separation of initially close states. Instead, it asks how far low-order RMT constraints determine an underlying vector representation. In the numerical implementation used below, the algebraic parameters are allowed to be complex. The resulting matrix is complex symmetric, \(M=M^T\), but not Hermitian, \(M\ne M^\dagger\). The relevant diagnostics are therefore non-Hermitian branch diagnostics rather than positive-semidefinite Gram checks. The problem is therefore best viewed as an inverse random-matrix reconstruction problem, or equivalently as a moment-constrained latent-vector learning problem.

The present problem is naturally viewed as a finite-dimensional constrained inverse problem. In this setting, the role of the imposed moment equations is analogous to that of measurement constraints in quantum-state estimation: they restrict the admissible manifold but do not, by themselves, guarantee a unique or physically preferred solution. This motivates the use of branch diagnostics, loss functions, and regularity criteria rather than a bare direct inversion alone [6]. A second issue is specific to finite Hilbert spaces. Global constraints can strongly reshape the distribution of weights and may produce either extended vectors or condensation-like concentration in a small number of components [7]. For this reason, the reconstructed eigenvectors should be monitored not only through their eigenvalues but also through inverse participation ratios (IPR), maximal weights, and related localization measures. This viewpoint is consistent with many-body studies in which fluctuations of few-body observables across eigenstates are controlled by the delocalization of eigenvectors in an appropriate basis [8]. Recent work on Floquet-MBL variational simulation emphasizes that Haar randomness is a useful benchmark but not necessarily the desired endpoint: a state ensemble can be useful precisely because it remains structured and non-Haar-like while avoiding excessive localization [9]. The moment constraints are therefore supplemented by branch-level, eigenvector-level, and low-weight-correlation diagnostics.

2 Vector ansatz↩︎

We construct a low-parameter ansatz for vector families whose induced Gram-type matrices reproduce selected low-order statistical properties of the target random matrix. Based on the structured vector ansatz, we determine its free parameters by matching selected low-order moment constraints of the target random-matrix ensemble. In contrast to a real-GOE reconstruction, the branch reported here keeps the original complex algebraic solutions. We focuses on \(N=8\), which is the largest size treated exactly in the present algebraic implementation. We define \[c=\cos\theta,\qquad s=\sin\theta,\] and take \(\theta=\pi/4\) in the numerical implementation. The unknown parameters are complex algebraic variables \[{\color{black}\{X_i,\,\widetilde{X}_i\}_{i=1}^{N/2}\subset \mathbb{C}.}\] For each \(i=1,\ldots,N/2\), we introduce four length-\(N\) vectors, \[\begin{align} \boldsymbol{p}^{(1)}_i &= (c,c,\ldots,c) \Big|_{2i-1\rightarrow X_i c,\;2i\rightarrow \widetilde{X}_i s}, \tag{1} \\ \boldsymbol{p}^{(2)}_i &= (-s,-s,\ldots,-s) \Big|_{2i-1\rightarrow -s,\;2i\rightarrow -c}, \tag{2} \\ \boldsymbol{q}^{(1)}_i &= (s,s,\ldots,s) \Big|_{2i-1\rightarrow s,\;2i\rightarrow -c}, \tag{3} \\ \boldsymbol{q}^{(2)}_i &= (c,c,\ldots,c) \Big|_{2i-1\rightarrow X_i c,\;2i\rightarrow -\widetilde{X}_i s}. \tag{4} \end{align}\] The notation means that all entries are first set to the displayed uniform background and then the two pair-resolved components \((2i-1,2i)\) are replaced. The two vector families are then assembled as \[\begin{align} P &= \frac{1}{\sqrt{N^2/2}} \begin{pmatrix} \boldsymbol{p}^{(1)}_1\\ \vdots\\ \boldsymbol{p}^{(1)}_{N/2}\\ \boldsymbol{p}^{(2)}_1\\ \vdots\\ \boldsymbol{p}^{(2)}_{N/2} \end{pmatrix}, \tag{5} \\ Q &= \frac{1}{\sqrt{N^2/2}} \begin{pmatrix} \boldsymbol{q}^{(1)}_1\\ \vdots\\ \boldsymbol{q}^{(1)}_{N/2}\\ \boldsymbol{q}^{(2)}_1\\ \vdots\\ \boldsymbol{q}^{(2)}_{N/2} \end{pmatrix}. \tag{6} \end{align}\] The corresponding bilinear overlap matrices are \[{\color{black}G_P=PP^{T},\qquad G_Q=QQ^{T}.}\] Since the transpose is not a conjugate transpose, these objects are complex symmetric rather than Hermitian once \(X_i\) and \(\widetilde{X}_i\) are complex. The reconstructed matrix is \[{\color{black}M=a_1G_P+a_2G_Q,} \label{eq:Mdef}\tag{7}\] so that \(M=M^T\) but generally \(M\ne M^\dagger\). Thus the branch used here is a complex-symmetric non-Hermitian reconstruction, not a real symmetric GOE reconstruction. In the implementation discussed here, \(a_1=3\) and \(a_2=-3\). A Hermitian complex variant would instead require \(G_P=PP^\dagger\) and \(G_Q=QQ^\dagger\), which is not the case in this work.

2.0.0.1 Moment-matching constraints.

The unknown parameters \(\{X_i,\widetilde{X}_i\}\) are not fitted to a specific target matrix element by element. Instead, they are determined by ensemble-level moment constraints. Denote the full matrix-element average by \[\overline{A_{\alpha\beta}}= \frac{1}{N^2}\sum_{\alpha,\beta=1}^{N}A_{\alpha\beta},\] the diagonal average by \[\overline{A_{ii}}_{\rm d}= \frac{1}{N}\sum_{i=1}^{N}A_{ii},\] and the off-diagonal average by \[\overline{A_{ij}}_{\rm od}= \frac{1}{N(N-1)} \sum_{i\ne j}A_{ij}.\] The algorithm solves the following algebraic constraints: \[\begin{align} \overline{M_{\alpha\beta}^{2}} &= \frac{1}{N^2} \left[ \frac{2(a_1^2+a_2^2)}{N} + \frac{(N^2-N)(a_1^2+a_2^2)}{N^2} \right], \tag{8} \\ \overline{M_{ii}}_{\rm d} &= \frac{a_1+a_2}{N}, \tag{9} \\ \overline{M_{ii}^{2}}_{\rm d} &= \frac{2(a_1^2+a_2^2)}{N^2}, \tag{10} \\ \overline{M_{ij}}_{\rm od} &=0, \tag{11} \\ \overline{M_{ij}^{2}}_{\rm od} &= \frac{a_1^2+a_2^2}{N^2}. \tag{12} \end{align}\] The individual Gram sectors are additionally constrained by \[\begin{align} \overline{(G_P)_{ii}}_{\rm d} &= \overline{(G_Q)_{ii}}_{\rm d} = \frac{1}{N}, \tag{13} \\ \overline{(G_P)_{ij}}_{\rm od} &= \overline{(G_Q)_{ij}}_{\rm od} = 0. \tag{14} \end{align}\] Finally, to avoid a systematic mixed contribution from the two Gram sectors, we impose \[\overline{F_{ij}}_{\rm od}=0,\qquad F_{ij}=2a_1a_2(G_P)_{ij}(G_Q)_{ij}. \label{eq:FFconstraint}\tag{15}\] Equations 8 –15 form a closed nonlinear algebraic system for the \(N\) unknowns \(\{X_i,\widetilde{X}_i\}\) at \(N=8\).

3 Estimator viewpoint and constrained branch selection↩︎

The algebraic solution can be interpreted as a direct-inversion estimator for the latent parameters: the observed data are the desired moment values, and the inferred object is a vector realization \(\mathcal{S}=\{P,Q\}\). As in quantum-state estimation, direct inversion can satisfy linear or moment equations while still leaving ambiguity among physical branches. We therefore define a reconstruction loss \[\mathcal{L}(\mathcal{S}) =\sum_{\mu} w_{\mu} \left[\mathcal{E}_{\mu}(\mathcal{S})-\mathcal{E}_{\mu}^{\ast}\right]^2, \label{eq:loss}\tag{16}\] where \(\mathcal{E}_{\mu}\) denotes the left-hand side of one of the constraints and \(\mathcal{E}_{\mu}^{\ast}\) denotes its target value. Exact algebraic solving corresponds to imposing \(\mathcal{L}=0\) within numerical tolerance. When exact solving becomes impractical at larger \(N\), the same expression can be used as an optimization objective, possibly supplemented by priors or regularizers. This is directly analogous in spirit to replacing raw direct inversion by constrained maximum-likelihood or Bayesian estimation in finite-dimensional state-estimation problems [6].

The algebraic system may have multiple complex branches. These branches are not equivalent from the viewpoint of the reconstructed vector geometry. Low-order moment matching fixes only certain averages of \(M\), \(G_P\), and \(G_Q\), while leaving residual freedom in the distribution of overlaps and in the complex spectral structure. A solution branch \(\mathcal{S}\) is first required to satisfy the residual condition \[{\color{black}\epsilon(\mathcal{S})=\max_{\mu}\left|{\cal E}_{\mu}(\mathcal{S})-{\cal E}_{\mu}^{\ast}\right|<\epsilon_{\rm tol}.} \label{eq:residual}\tag{17}\] For the present non-Hermitian complex-symmetric branch, positive-semidefinite Gram checks are not appropriate, because \(PP^T\) and \(QQ^T\) are bilinear rather than Hermitian Gram matrices. Instead we report the complex-symmetry score, the Hermiticity violation, and the imaginary spectral weight, \[\begin{align} \eta_T(\mathcal{S}) &= \|M-M^T\|_F,\\ \eta_H(\mathcal{S}) &= \|M-M^\dagger\|_F,\\ \eta_{\rm Im}(\mathcal{S}) &= \left\|{\rm Im}[\mathrm{spec}(M)]\right\|_2. \label{eq:nonhermdiag} \end{align}\tag{18}\] The exact bilinear construction gives \(\eta_T=0\) up to numerical precision, while nonzero \(\eta_H\) and \(\eta_{\rm Im}\) quantify the genuinely non-Hermitian character of the branch. We also retain the approximate row-orthogonality score \[{\color{black}{\cal R}(\mathcal{S})=\left\|{\rm offdiag}\,G_P(\mathcal{S})\right\|_F+\left\|{\rm offdiag}\,G_Q(\mathcal{S})\right\|_F.} \label{eq:orthscore}\tag{19}\] This criterion does not enforce Hermitian orthogonality; it only measures how small the bilinear off-diagonal overlaps are after the moment constraints have been imposed. In a Bayesian variant, one could assign a posterior weight to each branch, \[{\color{black}\Pi(\mathcal{S}|\mathcal{E}^{\ast})\propto \exp\!\left[-\frac{\mathcal{L}(\mathcal{S})}{2\sigma^2}-\lambda_{\rm R}\mathcal{R}(\mathcal{S})-\lambda_{\rm max}W_{\rm max}(\mathcal{S})-\lambda_H\eta_H(\mathcal{S})\right],} \label{eq:posteriorbranch}\tag{20}\] where \(W_{\rm max}\) penalizes anomalously concentrated vector components. The algorithm performs the branch selection and reports all diagnostic scores, so that a later version can replace manual branch choice by an explicit optimization criterion.

4 Delocalization, IPR, and finite-dimensional constraints↩︎

A useful eigenvector diagnostic is the inverse participation ratio \[{\rm IPR}_{n} = \sum_{\alpha=1}^{N} |v_{n\alpha}|^{4}, \label{eq:ipr}\tag{21}\] where \(\boldsymbol{v}_n\) is a normalized eigenvector of \(M\). For a maximally delocalized vector, \({\rm IPR}\sim 1/N\), whereas localized states have larger IPR. The effective participation number is \[D_{\rm eff}^{(n)}=\frac{1}{{\rm IPR}_n}. \label{eq:deff}\tag{22}\] This diagnostic is especially relevant because in interacting many-body systems the suppression of eigenstate-to-eigenstate fluctuations of observables is often controlled by eigenvector spreading in an appropriate basis [8]. For a diagonal observable \(A={\rm diag}(A_1,\ldots,A_N)\) in the reconstruction basis, \[A_n=\langle v_n|A|v_n\rangle =\sum_{\alpha}|v_{n\alpha}|^2A_{\alpha}.\] If the weights \(|v_{n\alpha}|^2\) sample many basis states with weak residual correlations, the variance of \(A_n\) across nearby eigenvectors is reduced roughly by the effective number of participating components. A schematic diagnostic is therefore \[\delta A^2 \propto {\rm Var}_{\alpha}(A_{\alpha}) \left\langle {\rm IPR}_n\right\rangle, \label{eq:fluctipr}\tag{23}\] where the proportionality constant depends on correlations among the weights. Equation 23 is not imposed by the reconstruction algorithm; it is a post-reconstruction test of whether the learned matrix behaves like an effectively delocalized finite-dimensional system.

The finite-dimensional nature of the problem also matters. Energy-constrained Haar-random states illustrate that imposing a global constraint can lead either to an extended high-temperature-like regime or to condensation, where one eigenstate carries macroscopic weight [7]. The present constraints are not energy constraints, but the lesson is directly relevant: moment constraints alone do not guarantee uniform delocalization. The maximum weight reads \[W_{\rm max}^{(n)}=\max_{\alpha}|v_{n\alpha}|^2. \label{eq:wmax}\tag{24}\] Small IPR together with small \(W_{\rm max}\) indicates broad spreading, whereas a small set of large weights signals localization or condensation-like concentration. In the present exact algebraic implementation, systematic scaling of Eqs. 21 and 24 is limited by solvability beyond \(N=8\). This limitation motivates replacing exact symbolic solving by numerical optimization or sampling at larger \(N\).

5 Haar-design, entanglement, and low-weight diagnostics↩︎

The MBL-initialization work suggests an additional caution that is particularly relevant for the present reconstruction: delocalization and Haar randomness are not identical. A branch can have moderately small IPR while still retaining low-complexity structure inherited from the pair-resolved ansatz. Conversely, a branch can become nearly Haar-like in its low-order moments, in which case the ansatz has effectively washed out the structure that the inverse reconstruction was meant to expose. We therefore introduce diagnostics that compare the reconstructed eigenvectors with Haar-design benchmarks.

For a complex Haar-random state in a \(d\)-dimensional Hilbert space, the expected \(t\)th-order participation moment is \[\mathbb{E}_{\rm Haar}\,{\rm IPR}_{t} = \mathbb{E}_{\rm Haar}\sum_{\alpha=1}^{d}|v_{\alpha}|^{2t} = \frac{d!\,t!}{(d+t-1)!}. \label{eq:complexhaariprt}\tag{25}\] For \(t=2\), this gives \(2/(d+1)\), the same benchmark used to diagnose the onset of state-design behavior in Floquet-MBL circuits [9]. Because the present branch is complex and non-Hermitian, the complex-Haar number is more relevant as a qualitative reference than the real-orthogonal value \(3/(d+2)\). However, right eigenvectors of a non-normal matrix are not distributed as Haar-random orthonormal vectors, so the benchmark should not be overinterpreted. The algorithm therefore reports the raw IPR and the maximum component weight rather than claiming convergence to a GOE or GUE eigenvector ensemble. A useful normalized comparison is \[{\color{black}\mathcal{D}_{\rm IPR}^{(\mathbb{C})}=\left|\frac{\langle {\rm IPR}_{2}\rangle}{2/(N+1)}-1\right|.} \label{eq:dipr}\tag{26}\] A large value indicates localization, eigenvector non-normality, or structured deviation from random-vector behavior. In the current numerical branches \(\langle {\rm IPR}_2\rangle\) is close to unity, so the reconstructed right eigenvectors are highly concentrated rather than Haar-like.

When \(N=2^n\), as in the present \(N=8\) calculation with \(n=3\), each eigenvector of \(M\) can also be interpreted as an \(n\)-qubit pure state in a computational basis. For a bipartition \(A\cup B\), define \[S_A(\boldsymbol{v}_n) = -{\rm Tr}\,\rho_A\log\rho_A, \qquad \rho_A={\rm Tr}_B |v_n\rangle\langle v_n|. \label{eq:ententropy}\tag{27}\] The branch-averaged entropy \(\langle S_A\rangle\) should be compared with the Page value for a random pure state with subsystem dimensions \(d_A\le d_B\), \[S_{\rm Page}(d_A,d_B) = \sum_{j=d_B+1}^{d_Ad_B}\frac{1}{j} - \frac{d_A-1}{2d_B}. \label{eq:pageentropy}\tag{28}\] This comparison separates three possibilities that IPR alone can blur: localized vectors with low participation, structured delocalized vectors with broad support but low entanglement, and Haar-like vectors with both broad participation and near-Page entanglement.

A still more targeted diagnostic is obtained by borrowing the low-weight stabilizer Rényi entropy used in the Floquet-MBL analysis [9]. Let \(\mathcal{P}_{n,k}\) be the set of \(n\)-qubit Pauli strings of weight no larger than \(k\), including the identity and fixing the overall phase to \(+1\). For an eigenvector \(|v_n\rangle\), define \[M_{t,k}(|v_n\rangle) = \frac{1}{1-t} \log\left[ \frac{1}{|\mathcal{P}_{n,k}|} \sum_{P\in\mathcal{P}_{n,k}} \left| \langle v_n|P|v_n\rangle \right|^{2t} \right]. \label{eq:lowweightsre}\tag{29}\] In a fully Haar-like state, fixed-weight Pauli expectation values are small, whereas a structured or localized branch can retain anomalously large low-weight expectations. For the present problem this quantity should not be interpreted as a literal many-body-localization order parameter, because \(M\) is not generated by a Floquet circuit. Its value is methodological: it detects whether the branch keeps low-weight, ansatz-resolved memory that is invisible to entrywise second moments. In practice, \(M_{2,2}\) is the most useful first choice because it probes fourth moments while remaining cheap for \(N=8\).

These diagnostics motivate a refined branch classification. A branch is called localized if both IPR and \(W_{\rm max}\) are large. It is Haar-like if \(\mathcal{D}_{\rm IPR}\) is small, \(S_A\) is close to \(S_{\rm Page}\), and low-weight Pauli expectations approach their Haar scale. It is structured-delocalized if the IPR is close to the random-vector benchmark but \(S_A\) or \(M_{2,2}\) still shows a clear deviation from Haar behavior. The last regime is the most interesting for inverse reconstruction, because it indicates that RMT-like matrix moments coexist with nontrivial latent-vector structure.

6 Algorithm↩︎

The reconstruction algorithm is summarized as follows.

1. Fix \(N\), \(\theta\), \(a_1\), and \(a_2\).

2. Build the pair-resolved vector ansatz \(P,Q\) from Eqs. 5 and 6 .

3. Construct \(G_P=PP^T\), \(G_Q=QQ^T\), and \(M=a_1G_P+a_2G_Q\).

4. Solve the nonlinear moment equations 8 –15 , or minimize the loss in Eq. 16 for larger systems.

5. Reject branches with large residuals or anomalously large bilinear-overlap, non-Hermiticity, conditioning, or concentration diagnostics.

6. Select the branch minimizing the Gram-regularity score in Eq. 19 , or average over branches using Eq. 20 .

7. Analyze the reconstructed vectors, matrix-element statistics, eigenvector IPR, maximum component weight, spectral diagnostics, and, when \(N=2^n\), entanglement and low-weight Pauli-moment diagnostics.

This procedure turns a random-matrix matching problem into a constrained latent-vector reconstruction problem.

7 Numerical diagnostics↩︎

The reconstruction algorithm used for the present numerical run returns six complex algebraic branches from a direct machine-precision algebraic solve. The least-squares fallback is therefore not used in this run. The report contains the constraint residuals, the bilinear orthogonality score \(\mathcal{R}\), the imaginary part of the reconstructed matrix, the complex-symmetry score \(\|M-M^T\|_F\), the Hermiticity violation \(\|M-M^\dagger\|_F\), the imaginary spectral weight \(\|{\rm Im}[\mathrm{spec}(M)]\|_2\), the average right-eigenvector IPR, and the average maximum component weight. The diagnostics used in the table are \[\begin{align} \langle {\rm IPR}\rangle &= \frac{1}{N}\sum_{n=1}^{N}{\rm IPR}_{n},\\ \langle W_{\rm max}\rangle&= \frac{1}{N}\sum_{n=1}^{N}W_{\rm max}^{(n)},\\ \eta_H &= \|M-M^\dagger\|_F,\\ \eta_{\rm Im} &= \|{\rm Im}[\mathrm{spec}(M)]\|_2. \end{align}\] Entanglement entropy and low-weight stabilizer Rényi entropy remain useful optional diagnostics for the \(N=2^n\) interpretation. Since the selected branch has \(N=8=2^3\), a Schmidt-spectrum diagnostic across a qubit bipartition is reported separately in Table 3.

The fifth branch is not selected because it is mathematically unique or globally optimal. All branch-specific quantities reported below use branch 5. If the goal is to define a fully automated workflow, then the selection rule should be changed explicitly: branch 6 gives the smallest displayed residual, branch 4 gives the smallest orthogonality score, and a combined selection rule could minimize a weighted sum of residual, \(\mathcal{R}\), non-Hermiticity, and eigenvector conditioning.

The large values of \({\rm IPR}_j\) and \(W_{\rm max}^{(j)}\) in Table 2 show that the right eigenvectors of the selected branch are strongly component-concentrated in the chosen basis. Therefore the branch should not be described as Haar-like or fully delocalized. It is better described as a structured complex-symmetric non-Hermitian branch that satisfies the imposed low-order moment constraints while retaining strong eigenvector concentration. Rather than listing the individual matrix elements, we report branch-level and eigenpair-resolved diagnostics that are invariant under trivial reorderings of the basis and are directly relevant to the reconstruction problem.

8 Schmidt-spectrum diagnostics inspired by many-body swapping↩︎

The many-body swapping protocol expresses the overlap and postselection cost of a shared state in terms of the Schmidt spectrum of a target many-body state [10]. The present work does not implement an entanglement-swapping communication protocol. Nevertheless, since the selected exact branch has \(N=8=2^3\), each right eigenvector can be viewed, in a basis-dependent manner, as a three-qubit state. This makes the Schmidt spectrum a useful additional diagnostic of whether the reconstructed right eigenvectors are bipartite-entangled or nearly product-like. For a normalized right eigenvector \(|v_n\rangle\), we reshape its components into a coefficient matrix \(C^{(n)}_{ab}\) associated with a chosen bipartition \(A|B\) and define the Schmidt weights \(\{\lambda^{(n)}_\ell\}\) from the squared singular values of \(C^{(n)}\). We then compute \[{\color{black} S_q^{(n)}=\frac{1}{1-q}\log\sum_{\ell}\left(\lambda^{(n)}_\ell\right)^q. }\] Following the functional form that appears in the many-body swapping fidelity and success probability, we introduce the diagnostic quantities \[\begin{align} & F_{\rm sw}^{(n)}=\exp\!\left(S_3^{(n)}-S_2^{(n)}\right),\\ &p_{\rm sw}^{(n)}=\exp\!\left[-2S_3^{(n)}\right] =\sum_{\ell}\left(\lambda^{(n)}_\ell\right)^3 . \end{align}\] Here \(F_{\rm sw}^{(n)}\) and \(p_{\rm sw}^{(n)}\) are not used as operational communication fidelities. Instead, they are Schmidt-spectrum shape diagnostics. They must be read together with the maximum Schmidt weight: a value \(F_{\rm sw}\simeq 1\) can occur both for a nearly flat spectrum and for a nearly rank-one spectrum. In Table 3, the maximum Schmidt weights are close to one, showing that the selected right eigenvectors are weakly entangled across the \(1|23\) partition despite the large values of \(F_{\rm sw}\).

9 Reconstructed Ensemble and the Latent Vectors↩︎

The reconstruction does not learn a unique original vector set. In the present complex-symmetric version, bilinear representations are nonunique up to complex rotations, sign/phase choices, branch choices, and ansatz restrictions. What the method learns is more specific: it identifies a structured vector realization whose overlaps reproduce the imposed random-matrix moments. Consequently, the algorithm is most suitable for four types of problems.

First, it is useful for inverse RMT: instead of asking whether a sampled matrix has RMT statistics, one asks what hidden vector geometry can generate those statistics. Second, it probes how much of a random-matrix ensemble is already encoded in low-order moment constraints. Third, it supplies a controlled way to distinguish matrix-level randomness from vector-level structure. Fourth, it provides a starting point for studying eigenvector delocalization by connecting element-level statistics to Gram-overlap geometry. Fifth, after the MBL-inspired extension, it can distinguish RMT-like delocalization from full Haar-like featurelessness by comparing IPR, maximum weight, entanglement, and low-weight moment diagnostics. The current branches in Table 1 are strongly concentrated in their right-eigenvector components, so they should be described as structured non-Hermitian branches rather than Haar-like random-vector branches.

10 Relation to chaos↩︎

Our method can generate matrices with selected RMT-like element moments and can be extended to test spectral diagnostics such as level repulsion, eigenvector delocalization, and spectral form factors. However, it does not by itself demonstrate sensitive dependence on initial conditions. Thus it is a static reconstruction of RMT-compatible vector structure, not as a direct evidence of deterministic chaos which would require comparing the evolution of nearby initial states, computing Lyapunov-type growth rates, or evaluating dynamical correlators such as out-of-time-order commutators.

11 Discussion↩︎

The main conceptual point is that RMT-like statistics need not be imposed at the matrix-entry level. They may emerge from structured vector overlaps after a constrained reconstruction. The pair-resolved ansatz used here is intentionally low-dimensional: at \(N=8\), only \(N\) complex unknowns are used to satisfy a comparable number of complex moment constraints. This makes the inverse problem solvable and interpretable, but it also restricts the accessible matrix manifold. The reconstruction should therefore be understood as a constructive representative within a chosen ansatz class, not as a unique inversion of an arbitrary random matrix. Because the ansatz uses \(PP^T\) rather than \(PP^\dagger\), the resulting branch belongs naturally to complex-symmetric non-Hermitian matrix theory.

The main limitation of the present exact construction is that the moment equations determine a finite set of complex branches but do not by themselves specify which branch is physically most informative. A natural extension is therefore to replace exact inversion by a loss-based reconstruction in which the moment residuals, regularity penalties, and possible physicality constraints are optimized simultaneously. Such a formulation would also make it possible to attach uncertainty estimates or branch weights to different reconstructed solutions. On the diagnostic side, the IPR should be connected not only to the visual delocalization of eigenvectors but also to fluctuations of reconstructed observables across eigenstates. The most informative regime is expected to be intermediate: the vectors should not be strongly localized or condensation-dominated, but they also need not become fully Haar-like. Instead, the relevant branches are those in which RMT-like low moments coexist with persistent structure inherited from the latent-vector ansatz.

There are two natural technical extensions. The first is to enlarge the ansatz by allowing additional variables \(Y_i,\widetilde{Y}_i\) in Eqs. 2 and 3 . In practice this creates an underdetermined complex system and direct solving may become slow; it is better handled by a loss-based optimizer or by adding additional constraints. The second is to replace exact algebraic solving with numerical optimization, minimizing a loss function built from the same moment constraints. A separate Hermitian extension would replace \(PP^T\) by \(PP^\dagger\), leading to real eigenvalues but complex eigenvectors; that is a different model from the one used here.

12 Size limitation and larger-\(N\) strategy↩︎

The exact algebraic implementation should presently be regarded as an \(N=8\) construction. This is not because the ansatz is mathematically meaningless for \(N>8\), but because the combination of pair-resolved complex variables, nonlinear moment constraints, and multiple branch choices makes direct exact solving rapidly inefficient for \(N=10,12,\ldots\). Enlarging the ansatz by including \(Y_i,\widetilde{Y}_i\) further increases the number of variables and may turn the algebraic problem into a high-dimensional or effectively underdetermined branch search. Therefore the present manuscript does not claim an exact larger-\(N\) branch family.

A more realistic route to larger dimensions is to replace exact branch solving by complex least-squares moment fitting. For \(N=8,10,12,16,\ldots\), one can minimize \[{\color{black}\min_{\{X_i,\widetilde{X}_i\}\subset\mathbb{C}} \sum_{\mu} \left| \mathcal{E}_{\mu}(X,\widetilde{X}) - \mathcal{E}_{\mu}^{\ast} \right|^2,} \label{eq:largerNloss}\tag{30}\] where \(\mathcal{E}_{\mu}\) denotes one of the moment functionals and \(\mathcal{E}_{\mu}^{\ast}\) is its target value. This larger-\(N\) procedure would produce approximate reconstructed matrices rather than exact algebraic branches. Such approximate matrices would be useful for future studies of spectral clouds, IPR scaling, Loschmidt-echo scaling, and branch stability, but they should not be mixed with the exact \(N=8\) branch data in Tables 1 and 2. Thus, in the present work, all numerical branch data are reported only for the exact \(N=8\) complex branch calculation.

13 Conclusion↩︎

We have formulated a moment-constrained vector reconstruction scheme for random-matrix statistics. Starting from two structured vector families, the present complex branch constructs \(M=a_1PP^T+a_2QQ^T\) and fixes the unknown complex vector parameters by matching diagonal, off-diagonal, and mixed-overlap moment constraints. Because the construction uses a transpose rather than a conjugate transpose, \(M\) is complex symmetric but generally non-Hermitian. The complex algebraic reconstruction at \(N=8\) returns six branches; the selected branch 5 has residual of order \(3.2\times 10^{-5}\) and strongly complex eigenvalues. The multiple algebraic solutions should therefore be interpreted as reconstruction branches, not as unique physical matrices. Some branches may be closer to bilinear orthogonality, while others may have smaller residual or different non-Hermitian spectral structure.

By incorporating risk/loss language from finite-dimensional estimation, constrained-state lessons from energy-conditioned Haar ensembles, IPR-based delocalization diagnostics from ETH studies, and MBL-inspired Haar-design diagnostics, the reconstruction becomes a more systematic protocol for learning hidden vector-space representations behind apparent RMT behavior. As an inverse structural reconstruction of complex RMT-like statistics, our algorithm constructs \(M\) statically and does not directly evaluate sensitive dependence on initial conditions, despite the RMT statistics are often tied to quantum chaos.

The present construction should not be confused with exceptional-point physics. Complex-symmetric branches may have complex spectra and non-orthogonal right eigenvectors, but exceptional points require eigenvector coalescence and defective Jordan structure. Degenerate or nearly degenerate reconstructed branches should therefore be diagnosed through eigenvector condition number, left-right biorthogonality, normality, and Jordan-defect tests rather than through eigenvalue coincidence alone.

References↩︎