Paper deep dive
Representation Learning with Quantum Signal Processing
Junqi Wang, Junyu Liu
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 95%
Last extracted: 9/1/2026, 5:03:02 AM
Summary
This paper establishes Quantum Signal Processing (QSP) as a solvable model for representation learning in quantum neural networks. It derives the exact mean and variance of the Quantum Neural Tangent Kernel (QNTK) at arbitrary depth, revealing an input-dependent angular geometry that remains non-self-averaging even under Haar randomness. The authors prove a sparse-data guarantee where the full nonlinear gradient flow converges to an integrable scalar flow with explicit convergence times, and demonstrate via numerical results that higher data density leads to coupled evolution beyond frozen-kernel or scalar descriptions.
Entities (8)
Relation Signals (6)
Junqi Wang → affiliatedwith → University of Pittsburgh
confidence 99% · Affiliation: Department of Computer Science, University of Pittsburgh
Junyu Liu → affiliatedwith → University of Pittsburgh
confidence 99% · Affiliation: Department of Computer Science, University of Pittsburgh
Quantum Signal Processing → isusedas → Representation Learning
confidence 95% · We establish quantum signal processing (QSP) as a solvable quantum model of the representation-learning regime.
Quantum Neural Tangent Kernel → hasproperty → Input-dependent angular geometry
confidence 93% · revealing an input-dependent angular geometry whose diagonal remains non-self-averaging
gradient flow → convergesto → Integrable scalar flow
confidence 92% · the realized dynamics converges to an integrable scalar flow
Quantum Neural Tangent Kernel → hasspectrum → Chebyshev
confidence 90% · its spectrum is Chebyshev structured
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Representation learning begins when training changes the features that define similarity between data. A frozen-kernel model only reweights a fixed geometry. We establish quantum signal processing (QSP) as a solvable quantum model of the representation-learning regime. At arbitrary depth, we compute the exact mean and variance of its quantum neural tangent kernel, revealing an input-dependent angular geometry whose diagonal remains non-self-averaging even when the underlying unitary approaches Haar randomness. We also prove a sparse-data guarantee for the full nonlinear gradient flow without freezing or ensemble-averaging the kernel: the realized dynamics converges to an integrable scalar flow with a time-dependent kernel closure and explicit convergence times. A finite-depth speed limit holds for every data set and trajectory. At higher data density, numerical results show coupled evolution beyond both the scalar and frozen-kernel descriptions. These results give a controlled theory of learned quantum data geometry with provable training dynamics beyond the frozen limit.
Tags
Links
- Source: https://arxiv.org/abs/2608.28828v1
- Canonical: https://arxiv.org/abs/2608.28828v1
Trouble viewing inline? Open PDF directly →
Full Text
87,743 characters extracted from source content.
Expand or collapse full text
Representation Learning with Quantum Signal Processing Junqi Wang Email: junqi.wang@pitt.edu Affiliation: Department of Computer Science, University of Pittsburgh, Pittsburgh, PA 15260, USA Junyu Liu Email: junyuliu@pitt.edu Affiliation: Department of Computer Science, University of Pittsburgh, Pittsburgh, PA 15260, USA August 28, 2026 Abstract Representation learning begins when training changes the features that define similarity between data. A frozen-kernel model only reweights a fixed geometry. We establish quantum signal processing (QSP) as a solvable quantum model of the representation-learning regime. At arbitrary depth, we compute the exact mean and variance of its quantum neural tangent kernel, revealing an input-dependent angular geometry whose diagonal remains non-self-averaging even when the underlying unitary approaches Haar randomness. We also prove a sparse-data guarantee for the full nonlinear gradient flow without freezing or ensemble-averaging the kernel: the realized dynamics converges to an integrable scalar flow with a time-dependent kernel closure and explicit convergence times. A finite-depth speed limit holds for every data set and trajectory. At higher data density, numerical results show coupled evolution beyond both the scalar and frozen-kernel descriptions. These results give a controlled theory of learned quantum data geometry with provable training dynamics beyond the frozen limit. I Introduction Learning can adjust both an output function and the representation in which data are compared. For a differentiable model, the tangent feature vector at an input is the gradient of the output with respect to its trainable parameters. Its Gram matrix is the neural tangent kernel (NTK). In an infinite-width or lazy limit these features barely move, the NTK is effectively frozen, and training becomes kernel regression [9, 11, 3]. Its fixed features prevent the model from changing its own notion of similarity during training. The quantum neural tangent kernel (QNTK) brings this learning-theory viewpoint to variational quantum circuits [15, 14, 1, 22]. A running QNTK distinguishes representation learning from the frozen regime: when the circuit parameters move, the data-dependent tangent features and their pairwise geometry move with them [15]. The same language has clarified how quantum laziness differs from barren plateaus and noise, how symmetry changes effective QNTK dimension and trainability, and how tangent kernels can diagnose practical QNN training [12, 24, 20]. Most existing guarantees assume wide, large-qubit, perturbative, or lazy limits [14, 6, 1]. Here we vary circuit depth, a different scaling resource from width or Hilbert-space dimension. A deep circuit may randomize its state without concentrating its tangent geometry. This work asks whether a provable training law remains possible when the realized kernel is random at initialization and continues to evolve during learning. Quantum signal processing (QSP) makes this question concrete. QSP is a depth-only, single-qubit polynomial model and the primitive underlying quantum singular-value transformation and many quantum algorithms [17, 18, 7, 19]. Its phases are trainable, while circuit width and Hilbert-space dimension remain fixed. Lin and collaborators made high-degree QSP phase synthesis numerically practical [5], characterized the nonconvex landscape of symmetric phase optimization and proved a locally strongly convex basin [23], and established a controlled d→∞d→∞ limit for phases representing smooth targets through infinite QSP [4]. These works study inverse synthesis around structured, target-dependent solutions. Here, “representation learning” has the learning-theory meaning: the tangent feature map induced by those phases changes under gradient flow. We initialize all phases independently and uniformly, placing the problem far from a prescribed solution and outside a deterministic frozen-kernel description. This setting exposes two linked difficulties. First, the initial kernel is a quenched random function of the data indices, so replacing it by its ensemble mean can erase fluctuations that survive at large depth. Second, the same realization evolves with the residuals, so even a fixed random-kernel model is insufficient. Architecture-specific quantum kernels already show how physical nonlinearities can modify convergence and generalization [16]. Finite-width tangent hierarchies, differential QNTKs, and nonlinear mean-field theories organize corrections beyond frozen kernels, and quantum neural networks can show time-dependent kernels and depth-driven dynamical transitions [8, 21, 14, 25, 26]. QSP complements these approaches with an exactly tractable, realization-resolved limit. At every depth, we solve the initialization geometry: the exact mean kernel is local in the angular data coordinate, its spectrum is Chebyshev structured, and its normalized diagonal converges to a nondegenerate beta law. For sufficiently sparse angular data, we prove that the full nonlinear QSP gradient flow converges on finite rescaled-time intervals to independent integrable scalar equations. A higher-density numerical example then shows how overlapping tangent neighborhoods produce coupled, time-dependent dynamics beyond both the scalar law and the initial-kernel prediction. I QSP as a representation-learning model For x∈[0,1]x∈[0,1], a degree-d QSP sequence and its real response are Ud(x,ϕ) U_d(x, φ) =eiϕ0Z∏j=1d[W(x)eiϕjZ], =e^i _0Z _j=1^d [W(x)e^i _jZ ], (1) gd(x,ϕ) g_d(x, φ) =Re[Ud(x,ϕ)]11, =Re[U_d(x, φ)]_11, where W(x)=xI+i1−x2XW(x)=xI+i 1-x^2X, and ϕ=(ϕ0,…,ϕd) φ=( _0,…, _d) contains M=d+1M=d+1 phases. For fixed phases, gdg_d extends to a real polynomial of degree at most d with parity d. At initialization, all phases are drawn independently and uniformly from [0,2π)[0,2π). Given n training pairs (xa,fa)a=1n\(x_a,f_a)\_a=1^n, with xa∈[0,1]x_a∈[0,1] and fa∈[−1,1]f_a∈[-1,1], we define the residuals and mean-squared loss by ea(ϕ)=gd(xa,ϕ)−fae_a( φ)=g_d(x_a, φ)-f_a and L(ϕ)=(2n)−1∑a=1nea2L( φ)=(2n)^-1 _a=1^ne_a^2. The QNTK is defined as Kab(ϕ)=∇ϕgd(xa,ϕ)⋅∇ϕgd(xb,ϕ).K_ab( φ)= _ φg_d(x_a, φ)· _ φg_d(x_b, φ). (2) Under continuous gradient flow, ϕ˙=−∇ϕL φ=- _ φL, the residuals obey the exact dynamics ˙=−1nK(ϕ(t)). e=- 1nK( φ(t)) e. (3) Here K(ϕ(t))K( φ(t)) is the realized, time-dependent QNTK. To make the learning interpretation explicit, define the tangent representation ϕ(x)=∇ϕgd(x,ϕ) _ φ(x)= _ φg_d(x, φ), with Kab(ϕ)=ϕ(xa)⋅ϕ(xb)K_ab( φ)= _ φ(x_a)· _ φ(x_b). Thus KabK_ab is the inner-product geometry assigned to the data pair (xa,xb)(x_a,x_b) by the current circuit. A frozen-kernel model keeps this geometry fixed at K(ϕ(0))K( φ(0)). QSP representation learning means that ϕ(t) _ φ(t) and hence the full data-indexed matrix Kab(t)K_ab(t) evolve as the phases are trained. I Random quantum data geometry Write x=cosαx= α and y=cosβy= β. Products of QSP responses can be represented as replicated transfer processes. Averaging each independent phase imposes a matching condition between the Z charges of the replicas, reducing the calculation to finite-dimensional transfer sectors. For two replicas, this gives the exact finite-depth mean QNTK Kd(x,y)=M4[cosd(α−β)+cosd(α+β)].EK_d(x,y)= M4 [ ^d(α-β)+ ^d(α+β) ]. (4) The direct term depends on the angular separation α−βα-β and measures local similarity, while the reflected term α+βα+β enforces the parity and boundary structure of QSP polynomials. Thus the data indices of KabK_ab are not interchangeable: their angular placement determines which residuals can communicate during training. Setting y=xy=x gives Kd(x,x)=M[1+(2x2−1)d]/4EK_d(x,x)=M[1+(2x^2-1)^d]/4. The corresponding four-replica construction gives the second moment. The phase average selects balanced four-copy charge sectors, whose transfer matrices yield the variance exactly at every finite depth. The full expression and transfer matrices are given in the Supplemental Material below. For every fixed interior input x∈(0,1)x∈(0,1), Kd(x,x)M EK_d(x,x)M ⟶14,VarKd(x,x)M2⟶1144, 14, VarK_d(x,x)M^2 1144, (5) VarKd(x,x)Kd(x,x) VarK_d(x,x)EK_d(x,x) ⟶13. 13. Thus both the mean and standard deviation grow linearly with depth, and increasing the number of phases does not suppress the relative fluctuations of a diagonal QNTK entry. The QSP unitary still randomizes completely. For every fixed x∈(0,1)x∈(0,1), the random walk generated by W(x)eiϕZW(x)e^iφ Z approaches Haar measure on SU(2)SU(2) [10, 2]. Consequently, Ud(x,ϕ)⇒d→∞dist.H,H∼Haar(SU(2)),U_d(x, φ) [d→∞]dist.H, H (SU(2)), (6) and gd(x,ϕ) g_d(x, φ) ⇒d→∞dist.G≡ReH11, [d→∞]dist.G _11, (7) pG(g) p_G(g) =2π1−g2,−1<g<1. = 2π 1-g^2, -1<g<1. To identify the origin of the surviving fluctuations, write Ud=gdI+id⋅U_d=g_dI+i c_d· σ. In the frame representation, each phase derivative is the projection of d c_d onto an associated Bloch-sphere axis j n_j. It follows exactly that Kd(x,x)=dSddK_d(x,x)= c_d TS_d c_d, where Sd=∑j=0djjS_d= _j=0^d n_j n_j T. At large depth, the normalized frame tensor becomes isotropic, Sd/M→I3/3S_d/M→ I_3/3 in probability. Since ∥d∥2=1−gd2 c_d ^2=1-g_d^2, this gives 3Kd(x,x)M−(1−gd(x)2)→d→∞Prob.0. 3K_d(x,x)M-(1-g_d(x)^2) [d→∞]Prob.0. (8) Combining Eq. (8) with the Haar limit yields 3Kd(x,x)M⇒d→∞dist.Beta(32,12). 3K_d(x,x)M [d→∞]dist.Beta ( 32, 12 ). (9) The frame representation and the proof of its large-depth isotropization are given in the Supplemental Material below. The beta law ties the surviving kernel fluctuations to the realized output amplitude, which is also the quantity moved by learning; the unitary has already mixed to Haar. The finite-depth distribution of the normalized diagonal QNTK and its ratio of standard deviation to mean are shown in Figs. 1(a) and 1(b), respectively. Figure 1: Non-self-averaging initialization and angular locality of the mean QNTK. (a) Distribution of the normalized diagonal QNTK z=3Kd(x,x)/Mz=3K_d(x,x)/M at x=0.6x=0.6. The distributions approach the Beta(3/2,1/2)Beta(3/2,1/2) law (black). The finite-depth weight above z=1z=1 reflects the residual anisotropy of the frame tensor. (b) Exact coefficient of variation VarKd(x,x)/Kd(x,x) VarK_d(x,x)/EK_d(x,x) for three fixed interior inputs. All three curves approach the nonzero asymptotic value 1/31/3 (dashed). (c) Along α+β=π/2α+β=π/2, the normalized mean kernel 4Kd(x,y)/M=cosd(α−β)4EK_d(x,y)/M= ^d(α-β) approaches e−u2/2e^-u^2/2 under the rescaling u=d|α−β|u= d\,|α-β|. The inset shows the corresponding narrowing in the unscaled angular coordinate. The mean kernel also determines the scale over which different inputs are coupled at initialization. On the physical half interval θ∈[0,π/2]θ∈[0,π/2], with normalized measure dμ(θ)=(2/π)dθdμ(θ)=(2/π)dθ, its integral operator is diagonal in the parity-matched cosine basis, equivalently the Chebyshev basis, with eigenvalues λm=M22−d(d(d−m)/2),m=d,d−2,…≥0. _m= M22^-d d(d-m)/2, m=d,d-2,…≥ 0. (10) Its algebraic rank is D=⌊d/2⌋+1D= d/2 +1, consistent with the parity restriction on degree-d QSP responses. Its participation ratio, rpart≡(∑mλm)2/∑mλm2r_part≡( _m _m)^2/ _m _m^2, is asymptotic to πd/2 π d/2. The Chebyshev modes are the global directions seen by the mean tangent representation, and their binomial eigenvalues quantify the initialization bias toward different data variations. This motivates the sample density w=nd,w= n d, (11) which measures how many samples occupy one mean-kernel correlation length. Away from the boundary, cosd(α+β) ^d(α+β) is suppressed at large depth, while expanding cosd(α−β) ^d(α-β) near α=βα=β gives the Gaussian profile shown in Fig. 1(c), with width O(d−1/2)O(d^-1/2). The parameter w is its natural occupancy scale: w≪1w 1 describes isolated tangent neighborhoods, whereas w≳1w 1 allows several data points to share features. The rigorous dynamical limit below requires a stronger condition that also controls the accumulated cross-input terms. IV Provable representation learning for sparse data Before taking a large-depth limit, the QSP dynamics already obeys a universal constraint. Each phase derivative satisfies |∂ϕjgd(xa)|≤1| _ _jg_d(x_a)|≤ 1, and hence Kaa≤MK_a≤ M and λmax(K)≤trK≤nM _ (K) ≤ nM. Using Eq. (3), −L˙=n−2K≤2ML- L=n^-2 e TK e≤ 2ML. Therefore, L(t)≥L(0)e−2Mt,tδ≥12Mlog1δ,L(t)≥ L(0)e^-2Mt, t_δ≥ 12M 1δ, (12) where tδt_δ is the first time at which L(tδ)=δL(0)L(t_δ)=δ L(0), 0<δ<10<δ<1. This bound holds for every depth, sample set, initialization, and gradient-flow trajectory. Physically, increasing depth adds more tangent directions; each direction can contribute to the rate at most one, so their combined effect limits the fastest loss descent to a rate proportional to M. The typical loss dynamics also depends on how the kernel is distributed across inputs. We take equally spaced angular nodes on the training interval, θa=(a−12)π/(2n) _a=(a- 12)π/(2n) and xa=cosθax_a= _a, a=1,…,na=1,…,n. Their angular separation is O(n−1)O(n^-1). When this separation is large compared with the correlation length O(d−1/2)O(d^-1/2), different training inputs are weakly coupled. Let M=d+1→∞M=d+1→∞ and allow n to grow with d. We require n3/2[1+logn]M⟶0. n^3/2[1+ n] M 0. (13) This condition is stronger than w→0w→ 0 because it controls both the accumulation and the dynamical transport of cross-input interactions. We introduce the rescaled time τ=Mt/(3n)τ=Mt/(3n). For compactness, write ga(τ)=gd(xa,ϕ(3nτ/M))g_a(τ)=g_d(x_a, φ(3nτ/M)) and Kaa(τ)=Kaa(ϕ(3nτ/M))K_a(τ)=K_a( φ(3nτ/M)) below. The QSP outputs reduce, on every fixed interval 0≤τ≤T0≤τ≤ T, to the scalar dynamics dg¯adτ=−(1−g¯a2)(g¯a−fa),g¯a(0)=gd(xa,ϕ(0)). d g_adτ=-(1- g_a^2)( g_a-f_a), g_a(0)=g_d(x_a, φ(0)). (14) Sparse representation-learning guarantee.— For every fixed T<∞T<∞, sup0≤τ≤T[1n∑a=1n|ga(τ)−g¯a(τ)|2]1/2→d→∞Prob.0. _0≤τ≤ T [ 1n _a=1^n |g_a(τ)- g_a(τ) |^2 ]^1/2 [d→∞]Prob.0. (15) Here the probability is over the random initial phases. The scalar initial conditions are matched to the same QSP realization, so the reduction preserves its initialization fluctuations. Along the same trajectory, the evolving diagonal kernel obeys sup0≤τ≤T[1n∑a=1n|3Kaa(τ)M−(1−ga(τ)2)|2]1/2→d→∞Prob.0. _0≤τ≤ T [ 1n _a=1^n | 3K_a(τ)M- (1-g_a(τ)^2 ) |^2 ]^1/2 [d→∞]Prob.0. (16) In the sparse regime, coupling between inputs becomes subleading, while the diagonal tangent kernel is directly determined by the output, Kaa≃M(1−ga2)/3K_a M(1-g_a^2)/3. Substitution into the exact output dynamics gives Eq. (14) for each gag_a. The loss decreases faster near ga=0g_a=0 and slows as the response approaches the boundaries ga=±1g_a=± 1. Unless a sample begins at its target, gag_a changes by an order-one amount on the rescaled time scale, so Eq. (16) makes Kaa(t)K_a(t) change at leading order as well. The resulting closure is nonlinear and describes an evolving representation even though different inputs become asymptotically independent. Equation (14) is integrable. The time required for one sample to move from ga,0g_a,0 to gag_a is determined by τ=∫gaga,0u/[(1−u2)(u−fa)]τ= _g_a^g_a,0du/[(1-u^2)(u-f_a)]. For example, when fa=0f_a=0, define the relative-error convergence time by g¯a(τq)2=qg¯a(0)2 g_a( _q)^2=q g_a(0)^2, with 0<q<10<q<1. One obtains τq=12log[1−qg¯a(0)2q[1−g¯a(0)2]],tq=3nMτq. _q= 12 [ 1-q g_a(0)^2q[1- g_a(0)^2] ], t_q= 3nM _q. (17) Thus Eq. (17) gives an explicit convergence time for a single sample. In the sparse regime, each sample evolves independently, so their contributions to the multipoint loss add directly. For fa=0f_a=0 at all training points, the scalar loss is L¯(τ)=12n∑a=1ng¯a(0)2g¯a(0)2+[1−g¯a(0)2]e2τ. L(τ)= 12n _a=1^n g_a(0)^2 g_a(0)^2+[1- g_a(0)^2]e^2τ. (18) The multipoint convergence time is therefore obtained as the unique solution of the monotonic equation L¯(τδ)=δL¯(0) L( _δ)=δ L(0), with tδ=3nτδ/Mt_δ=3n _δ/M. Detailed proofs of Eqs. (15) and (16), as well as the sparse condition Eq. (13), are given in the Supplemental Material below. Figure 2: Learning with isolated and overlapping quantum data neighborhoods. QSP loss for a representative initialization (blue) and the interquartile band over 3030 initializations, compared with the realization-matched scalar flow (red dotted), the realized kernel frozen at K(0)K(0) (yellow dashed), and the continuous-flow speed-limit envelope e−2Mte^-2Mt (gray). (a) For (d,n,w)=(1024,4,0.125)(d,n,w)=(1024,4,0.125), both reduced descriptions track the QSP descent, with median ratios tsc/tQSP=0.97t_sc/t_QSP=0.97 and tK(0)/tQSP=1.04t_K(0)/t_QSP=1.04. (b) For (d,n,w)=(128,64,5.66)(d,n,w)=(128,64,5.66), overlapping tangent neighborhoods drive the QSP trajectory away from both predictions. The numerical protocol is given in the Supplemental Material. V Coupled quantum data geometry at higher density The scalar law relies on angular locality. Since w=n/dw=n/ d measures the number of training points within one mean-kernel correlation length, increasing w increases the overlap of tangent neighborhoods. The exact cross-input contribution to the output dynamics is −n−1∑b≠aKab(t)eb(t)-n^-1 _b≠ aK_ab(t)e_b(t). When this contribution is no longer subleading, changes at one input modify the outputs at other inputs through the shared tangent geometry, and the dynamics cannot be reduced to independent scalar equations. Figure 2 compares the QSP training trajectory with the scalar flow matched to the same random initialization and with the evolution generated by the frozen kernel K(0)K(0). In the sparse case, the scalar flow follows the nonlinear QSP descent. In the representative high-density case, the QSP trajectory departs from both the scalar prediction and the frozen-K(0)K(0) reference, showing that neither independent scalar dynamics nor a fixed realized kernel describes that case. The same QSP model therefore exhibits two distinct forms of representation learning. Under Eq. (13), the realization-dependent kernel closes on the outputs and the multipoint problem reduces to independent nonlinear scalar equations. At higher density, the numerical example exhibits coupled, time-dependent training: the learned representation changes both the self-similarity of each input and the communication between inputs. Here w quantifies geometric overlap, whereas Eq. (13) is the sufficient asymptotic condition proved for scalar closure. VI Discussion and outlook QSP separates three notions that are often conflated: randomization, concentration, and representation learning. Increasing QSP depth randomizes the unitary at every fixed interior input. The normalized diagonal QNTK remains realization dependent. The frame picture explains this persistence: the tangent axes become isotropic, while the output amplitude remains random and controls the size of the tangent feature vector. The state approaches Haar randomness while presenting a random metric to the learning problem. The surviving randomness affects the training trajectory. In the sparse regime, Eqs. (15) and (16) give a trajectory-level guarantee for the full QSP learner. The parameters move, the outputs move, and the diagonal kernel follows them through Kaa(t)≃M[1−ga(t)2]/3K_a(t) M[1-g_a(t)^2]/3. The dynamics remains solvable: each datum follows an integrable nonlinear flow with explicit convergence times. This trajectory-level closure establishes provable representation learning beyond the frozen limit. The universal speed limit supplies a complementary statement that holds without sparsity, randomness, or asymptotics. The exact dependence of KabK_ab on its data indices also suggests a useful theory of quantum data geometry. The QSP encoding turns x=cosθx= θ into an angular space with resolution d−1/2d^-1/2, Chebyshev eigenmodes, and a depth-dependent spectral bias. The density w=n/dw=n/ d then has a direct interpretation: it counts how many observations occupy one tangent neighborhood. Sparse observations primarily change their own local metric; overlapping observations can reshape shared features and transmit error across the data set. The kernel serves as a convergence-rate matrix and as a dynamical description of how the quantum model organizes data. The angular geometry organizes the next questions. Varying w probes the crossover from isolated to overlapping tangent neighborhoods and asks whether the off-diagonal frame sums admit a coupled continuum closure. The Chebyshev spectrum provides a mode-resolved route from kernel evolution to generalization, while the scalar sector provides an exact setting in which regularization can be followed through relaxation rates and solution selection. Extending the transport argument to irregular grids and higher-dimensional encodings would determine which aspects of the closure follow from localized quantum feature geometry. Theorem-level analyses of quantum gradient descent likewise tie controlled learning dynamics to explicit structural hypotheses such as sparsity and dissipation [13]. Together, these examples place geometry and transport structure at the center of quantum-learning guarantees. QSP is a minimal quantum model with a random, time-dependent, and provably tractable learned data geometry. It provides a bridge from frozen quantum kernels to a first-principles theory of nonlinear quantum representation learning. Data and code availability.— The plotting code for Figs. 1 and 2, together with the numerical data underlying Fig. 2, is available at https://github.com/AlexWangPhy233/qsp-representation-learning. Acknowledgements. Acknowledgments.— JW and JL are supported in part by the University of Pittsburgh, School of Computing and Information, Department of Computer Science, Pitt Cyber, Pitt Momentum fund, PQI Community Collaboration Awards, John C. Mascaro Faculty Scholar in Sustainability, Switzerland NSF award 2000-1-243053, NSF award 2535915, 2610010, DOE Genesis Program, Thinking Machines Lab and Cisco Research. This research used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR22725. References [1] E. Abedi, S. Beigi, and L. Taghavi (2023) Quantum lazy training. Quantum 7, p. 989. External Links: Document Cited by: §I. [2] D. Applebaum (2014) Probability on compact lie groups. Probability Theory and Stochastic Modelling, Vol. 70, Springer, Cham. External Links: Document Cited by: §S3.B, §I. [3] L. Chizat, E. Oyallon, and F. Bach (2019) On lazy training in differentiable programming. Advances in neural information processing systems 32. Cited by: §I. [4] Y. Dong, L. Lin, H. Ni, and J. Wang (2024) Infinite quantum signal processing. Quantum 8, p. 1558. External Links: Document Cited by: §I. [5] Y. Dong, X. Meng, K. B. Whaley, and L. Lin (2021) Efficient phase-factor evaluation in quantum signal processing. Physical Review A 103, p. 042419. External Links: Document Cited by: §I. [6] D. García-Martín, M. Larocca, and M. Cerezo (2025) Quantum neural networks form gaussian processes. Nature Physics 21, p. 1153–1159. External Links: Document Cited by: §I. [7] A. Gilyén, Y. Su, G. H. Low, and N. Wiebe (2019) Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics. In Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing, p. 193–204. External Links: Document Cited by: §I. [8] J. Huang and H. Yau (2020) Dynamics of deep neural networks and neural tangent hierarchy. In Proceedings of the 37th International Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 119, p. 4542–4551. Cited by: §I. [9] A. Jacot, F. Gabriel, and C. Hongler (2018) Neural tangent kernel: convergence and generalization in neural networks. Advances in neural information processing systems 31. Cited by: §I. [10] Y. Kawada and K. Itô (1940) On the probability distribution on a compact group. I. Proceedings of the Physico-Mathematical Society of Japan, 3rd Series 22, p. 977–998. Cited by: §S3.B, §I. [11] J. Lee, L. Xiao, S. Schoenholz, Y. Bahri, R. Novak, J. Sohl-Dickstein, and J. Pennington (2019) Wide neural networks of any depth evolve as linear models under gradient descent. Advances in neural information processing systems 32. Cited by: §I. [12] J. Liu, Z. Lin, and L. Jiang (2024) Laziness, barren plateau, and noises in machine learning. Machine Learning: Science and Technology 5 (1), p. 015058. External Links: Document Cited by: §I. [13] J. Liu, M. Liu, J. Liu, Z. Ye, Y. Wang, Y. Alexeev, J. Eisert, and L. Jiang (2024) Towards provably efficient quantum algorithms for large-scale machine-learning models. Nature Communications 15 (1), p. 434. External Links: Document Cited by: §VI. [14] J. Liu, K. Najafi, K. Sharma, F. Tacchino, L. Jiang, and A. Mezzacapo (2023) Analytic theory for the dynamics of wide quantum neural networks. Physical Review Letters 130, p. 150601. External Links: Document Cited by: §I, §I. [15] J. Liu, F. Tacchino, J. R. Glick, L. Jiang, and A. Mezzacapo (2022) Representation learning via quantum neural tangent kernels. PRX Quantum 3, p. 030323. External Links: Document Cited by: §I. [16] J. Liu, C. Zhong, M. Otten, A. Chandra, C. L. Cortes, C. Ti, S. K. Gray, and X. Han (2023) Quantum kerr learning. Machine Learning: Science and Technology 4 (2), p. 025003. External Links: Document Cited by: §I. [17] G. H. Low and I. L. Chuang (2017) Optimal hamiltonian simulation by quantum signal processing. Physical Review Letters 118, p. 010501. External Links: Document Cited by: §I. [18] G. H. Low and I. L. Chuang (2019) Hamiltonian simulation by qubitization. Quantum 3, p. 163. External Links: Document Cited by: §I. [19] J. M. Martyn, Z. M. Rossi, A. K. Tan, and I. L. Chuang (2021) Grand unification of quantum algorithms. PRX Quantum 2, p. 040203. External Links: Document Cited by: §I. [20] F. Scala, C. Zoufal, D. Gerace, and F. Tacchino (2026) Towards practical quantum neural network diagnostics with neural tangent kernels. npj Quantum Information. External Links: Document Cited by: §I. [21] M. Seleznova and G. Kutyniok (2022) Neural tangent kernel beyond the infinite-width limit: effects of depth and initialization. In Proceedings of the 39th International Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 162, p. 19522–19560. Cited by: §I. [22] N. Shirai, K. Kubo, K. Mitarai, and K. Fujii (2024) Quantum tangent kernel. Physical Review Research 6, p. 033179. External Links: Document Cited by: §I. [23] J. Wang, Y. Dong, and L. Lin (2022) On the energy landscape of symmetric quantum signal processing. Quantum 6, p. 850. External Links: Document Cited by: §I. [24] X. Wang, J. Liu, T. Liu, Y. Luo, Y. Du, and D. Tao (2023) Symmetric pruning in quantum neural networks. In The Eleventh International Conference on Learning Representations, External Links: Link Cited by: §I. [25] X. You, S. Chakrabarti, B. Chen, and X. Wu (2023) Analyzing convergence in quantum neural networks: deviations from neural tangent kernels. In Proceedings of the 40th International Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 202, p. 40199–40224. Cited by: §I. [26] B. Zhang, J. Liu, X. Wu, L. Jiang, and Q. Zhuang (2024) Dynamical transition in controllable quantum neural networks with large depth. Nature Communications 15, p. 9354. External Links: Document Cited by: §I. Supplemental Material for “Representation Learning with Quantum Signal Processing” Junqi Wang and Junyu Liu Department of Computer Science, University of Pittsburgh, Pittsburgh, Pennsylvania 15260, USA Conventions For x∈[0,1]x∈[0,1] we use Ud(x,ϕ)=eiϕ0Z∏j=1d[W(x)eiϕjZ],W(x)=xI+i1−x2X,U_d(x, φ)=e^i _0Z _j=1^d [W(x)e^i _jZ ], W(x)=xI+i 1-x^2\,X, (S1) where factors are ordered from left to right with increasing j. The number of trainable phases is M=d+1M=d+1, and at initialization ϕ0,…,ϕd _0,…, _d are independent and uniform on [0,2π)[0,2π). We write gd(x,ϕ)=Re[Ud(x,ϕ)]11,Jj(x,ϕ)=∂ϕjgd(x,ϕ),Kd(x,y,ϕ)=∑j=0dJj(x,ϕ)Jj(y,ϕ).g_d(x, φ)=Re[U_d(x, φ)]_11, J_j(x, φ)= _ _jg_d(x, φ), K_d(x,y; φ)= _j=0^dJ_j(x, φ)J_j(y, φ). (S2) For training data (xa,fa)a=1n\(x_a,f_a)\_a=1^n, ea=gd(xa,ϕ)−fa,L(ϕ)=12n∑a=1nea2,ϕ˙=−∇ϕL.e_a=g_d(x_a, φ)-f_a, L( φ)= 12n _a=1^ne_a^2, φ=- _ φL. (S3) Thus ˙=−K/n e=-K e/n and L˙=−K/n2 L=- e TK e/n^2 exactly. The phase average is governed by the following matching rule. Label matrix indices by a,b∈0,1a,b∈\0,1\ and define the Z charge χ(0)=+1χ(0)=+1, χ(1)=−1χ(1)=-1. Since [W(x)eiϕZ]ab=W(x)abeiϕχ(b),[W(x)e^iφ Z]_ab=W(x)_abe^iφχ(b), (S4) for signs sr∈+1,−1s_r∈\+1,-1\ one has ϕ∏r=1peisrϕχ(br)=∑r=1psrχ(br)=0.E_φ _r=1^pe^is_rφχ(b_r)= 1\! \ _r=1^ps_rχ(b_r)=0 \. (S5) The signs sr=−1s_r=-1 label conjugated copies. Equation (S5) leaves two balanced states in a two-copy average and six in a four-copy average with two ket and two bra factors. Writing x=cosθx= θ, the response extends from the physical interval to a degree-d parity polynomial, gd(−x,ϕ)=(−1)dgd(x,ϕ),gd(cosθ,ϕ)∈spancos(dθ),cos[(d−2)θ],….g_d(-x, φ)=(-1)^dg_d(x, φ), g_d( θ, φ) \ (dθ), [(d-2)θ],…\. (S6) Unless stated otherwise, all large-d limits hold for fixed x∈(0,1)x∈(0,1). The endpoint formulas at x=0x=0 and x=1x=1 are evaluated separately at finite depth. S1 Two-copy transfer and the exact mean kernel S1.A Derivative insertions Let B0=eiϕ0ZB_0=e^i _0Z and Bj=W(x)eiϕjZB_j=W(x)e^i _jZ for j≥1j≥ 1. With U≤j=B0⋯BjU_≤ j=B_0·s B_j and U>j=Bj+1⋯BdU_>j=B_j+1·s B_d, the commutation of Z with the phase gate gives ∂ϕjUd=U≤j(iZ)U>j. _ _jU_d=U_≤ j(iZ)U_>j. (S7) Equivalently, if Rj=e1U≤jR_j=e_1 TU_≤ j and Cj=U>je1C_j=U_>je_1, Jj(x)=−Im[(Rj)1(Cj)1−(Rj)2(Cj)2].J_j(x)=-Im\! [(R_j)_1(C_j)_1-(R_j)_2(C_j)_2 ]. (S8) Let ux=[Ud(x)]11u_x=[U_d(x)]_11, so that gx=(ux+ux¯)/2g_x=(u_x+ u_x)/2. In the product Jj(x)Jj(y)J_j(x)J_j(y), the holomorphic–holomorphic and antiholomorphic–antiholomorphic terms are removed by the ϕ0 _0 average. The remaining two mixed terms are complex conjugates and give an overall factor 1/21/2. S1.B Charge projection For one ket copy at x and one bra copy at y, Eq. (S5) gives ϕ[(WxeiϕZ)aa′(WyeiϕZ)bb′¯] _φ [(W_xe^iφ Z)_a (W_ye^iφ Z)_b ] =(Wx)aa′(Wy)bb′¯χ(a′)=χ(b′). 51.21495pt=(W_x)_a (W_y)_b 1\! \χ(a )=χ(b ) \. (S9) The balanced subspace is therefore ℬ2=span|00⟩,|11⟩B_2=span\|00 ,|11 \. Put sx=1−x2s_x= 1-x^2 and sy=1−y2s_y= 1-y^2. Between two consecutive phase projections the signal gate acts on ℬ2B_2 as Txy=(xysxsysxsyxy).T_xy= pmatrixxy&s_xs_y\\ s_xs_y&xy pmatrix. (S10) For x=cosαx= α and y=cosβy= β, its normalized eigenvectors are (1,±1)/2(1,± 1)/ 2 and spec(Txy)=xy+sxsy,xy−sxsy=cos(α−β),cos(α+β).spec(T_xy)=\xy+s_xs_y,xy-s_xs_y\=\ (α-β), (α+β)\. (S11) At the differentiated phase, the ket derivative contributes iχ(bx)iχ(b_x) and the bra derivative contributes −iχ(by)-iχ(b_y). On the balanced states bx=byb_x=b_y, their product is χ(bx)χ(by)=1χ(b_x)χ(b_y)=1. The marked position therefore only separates two powers of the same transfer matrix. For every j=0,…,dj=0,…,d, [Jj(x)Jj(y)]=12(TxyjTxyd−j)00,00=12(Txyd)00,00.E[J_j(x)J_j(y)]= 12 (T_xy^jT_xy^d-j )_00,00= 12(T_xy^d)_00,00. (S12) Summing the M=d+1M=d+1 possible derivative positions and diagonalizing TxyT_xy gives Kd(x,y)=M4[cosd(α−β)+cosd(α+β)].EK_d(x,y)= M4 [ ^d(α-β)+ ^d(α+β) ]. (S13) The same unmarked chain gives [gd(x)gd(y)]=1MKd(x,y),Kd(x,x)=M4[1+(2x2−1)d].E[g_d(x)g_d(y)]= 1MEK_d(x,y), _d(x,x)= M4 [1+(2x^2-1)^d ]. (S14) For fixed interior inputs, both powers in Eq. (S13) decay exponentially unless their angular arguments approach a multiple of π. In the local regime |α−β|=O(d−1/2)|α-β|=O(d^-1/2), cosd(α−β)=exp[−d2(α−β)2+O(d|α−β|4)], ^d(α-β)= \! [- d2(α-β)^2+O\! (d|α-β|^4 ) ], (S15) so the direct term has angular correlation length O(d−1/2)O(d^-1/2). S2 Four-copy transfer and the exact variance S2.A Six phase-matched sectors The raw second moment is Kd(x,y)2=∑i,j=0d[Ji(x)Ji(y)Jj(x)Jj(y)].EK_d(x,y)^2= _i,j=0^dE[J_i(x)J_i(y)J_j(x)J_j(y)]. (S16) Expanding the four real derivatives into holomorphic and antiholomorphic parts produces 242^4 conjugation patterns. Label a pattern by the subset C⊂1,2,3,4C⊂\1,2,3,4\ of conjugated copies, ordered as (x,y,x,y)(x,y,x,y). The ϕ0 _0 average imposes zero total charge on the all-zero boundary state. Hence only the six choices C∈12,13,14,23,24,34C∈\12,13,14,23,24,34\ (S17) survive, where 1212 abbreviates 1,2\1,2\, and similarly for the other pairs. The sectors C=13,24C=13,24 give the two XXYY contractions, whereas C=12,14,23,34C=12,14,23,34 give the four XYXY contractions. Both classes are generated by the same construction below. Set z1=z3=xz_1=z_3=x, z2=z4=yz_2=z_4=y, and sc=−1s_c=-1 for c∈Cc∈ C and +1+1 otherwise. For four-bit indices ,∈0,14 a, b∈\0,1\^4, define the unprojected signal step [(C)],=∏c=14C,c(W(zc)acbc),qC()=∑c=14scχ(bc),[A^(C)]_ a, b= _c=1^4 c_C,c\! (W(z_c)_a_cb_c ), q_C( b)= _c=1^4s_cχ(b_c), (S18) where C,c c_C,c complex conjugates its argument exactly when c∈Cc∈ C. Equation (S5) projects onto ℬC=span|⟩:qC()=0,dimℬC=6.B_C=span\\,| b :q_C( b)=0\,\, _C=6. (S19) If PCP_C denotes this projector, every averaged segment is generated by the 6×66× 6 compression MC=PC(C)PC.M_C=P_CA^(C)P_C. (S20) The two differentiated pairs become diagonal insertions on ℬCB_C, D12|⟩=χ(b1)χ(b2)|⟩,D34|⟩=χ(b3)χ(b4)|⟩.D_12| b =χ(b_1)χ(b_2)| b , D_34| b =χ(b_3)χ(b_4)| b . (S21) The four factors of ±i± i multiply to +1+1 in every surviving sector. The factor 1/161/16 below comes from the four real-part decompositions. S2.B The explicit six-dimensional transfer As one representative of the six sector transfer matrices, take C=34C=34 and use the ordered balanced basis |00;00⟩,|01;01⟩,|01;10⟩,|10;01⟩,|10;10⟩,|11;11⟩.|00;00 ,|01;01 ,|01;10 ,|10;01 ,|10;10 ,|11;11 . (S22) Introduce A=x2y2,B=(1−x2)(1−y2),ℳ=xysxsy,L=x2(1−y2),R=y2(1−x2).A=x^2y^2, B=(1-x^2)(1-y^2), =xy\,s_xs_y, L=x^2(1-y^2), R=y^2(1-x^2). (S23) Directly applying Eq. (S18) gives M34=(ALℳRBLA−ℳ−ℳBRℳ−ℳAB−ℳ−ℳBA−ℳRB−ℳ−ℳALBRℳLA).M_34= pmatrixA&L&M&M&R&B\\ L&A&-M&-M&B&R\\ M&-M&A&B&-M&M\\ M&-M&B&A&-M&M\\ R&B&-M&-M&A&L\\ B&R&M&M&L&A pmatrix. (S24) Its six eigenvalues are λ1 _1 =A−B, =A-B, λ2 _2 =A−B+L−R=2x2−1, =A-B+L-R=2x^2-1, λ3 _3 =A+B+L+R=1, =A+B+L+R=1, λ4 _4 =A−B−L+R=2y2−1, =A-B-L+R=2y^2-1, λ5,6 _5,6 =12[2A+2B−L−R±32ℳ2+(L+R)2]. = 12 [2A+2B-L-R± 32M^2+(L+R)^2 ]. (S25) On the diagonal x=yx=y, define r1=2x2−1,r2=6x4−6x2+1.r_1=2x^2-1, r_2=6x^4-6x^2+1. (S26) Then spec(M34)|x=y=1,1,r1,r1,r1,r2.spec(M_34) |_x=y=\1,1,r_1,r_1,r_1,r_2\. (S27) The other three XYXY sectors are related to C=34C=34 by permutations of the copies and have the same transfer spectrum. The two XXYY sectors are evaluated by the same projected-transfer construction. S2.C Double sum Let e0=|0000⟩e_0=|0000 , which belongs to every ℬCB_C. For fixed C and insertion sites i,ji,j, define ij(C)=MCiD12MCj−iD34MCd−j,i≤j,MCjD34MCi−jD12MCd−i,j<i.C_ij^(C)= casesM_C^iD_12M_C^j-iD_34M_C^d-j,&i≤ j,\\ M_C^jD_34M_C^i-jD_12M_C^d-i,&j<i. cases (S28) The exact finite-depth second moment is Kd(x,y)2=116∑|C|=2∑i,j=0d⟨e0|ij(C)|e0⟩.EK_d(x,y)^2= 116 _|C|=2 _i,j=0^d e_0|C_ij^(C)|e_0 . (S29) This expression contains all phase matching and all derivative positions. If MC=∑aλaPaM_C= _a _aP_a is its spectral resolution, every ordered site sum reduces to hd(a,b,c)=∑r+s+t=darbsct.h_d(a,b,c)= _r+s+t=da^rb^sc^t. (S30) For distinct a,b,ca,b,c, hd(a,b,c)=ad+2(a−b)(a−c)+bd+2(b−a)(b−c)+cd+2(c−a)(c−b),h_d(a,b,c)= a^d+2(a-b)(a-c)+ b^d+2(b-a)(b-c)+ c^d+2(c-a)(c-b), (S31) with coincident eigenvalues obtained by continuous limits. Equations (S29) and (S31) give the finite-depth evaluation. In particular, the two unit eigenvalues in Eq. (S27) generate a quadratic site sum and hence an O(d2)O(d^2) raw second moment. Put Δ1=r1d _1=r_1^d and Δ2=r2d _2=r_2^d. Evaluating the six sector sums gives Kd(x,x)2=1432[ _d(x,x)^2= 1432 [ 3d2(27Δ1+Δ2+10)+6d(27Δ1+5Δ2+6) 3d^2(27 _1+ _2+10)+6d(27 _1+5 _2+6) +(8d+16)(Δ2−1)x2(x2−1)+81Δ1+75Δ2+6]. + (8d+16)( _2-1)x^2(x^2-1)+81 _1+75 _2+6 ]. (S32) Although written as a rational expression, the right-hand side is a polynomial after the removable limits are taken. The endpoint branches are Kd(0,0)2=316[1+(−1)d]M2,Kd(1,1)2=38M2.EK_d(0,0)^2= 316[1+(-1)^d]M^2, _d(1,1)^2= 38M^2. (S33) For every fixed x∈(0,1)x∈(0,1), |r1|<1|r_1|<1 and |r2|<1|r_2|<1, so Kd(x,x)M⟶14,Kd(x,x)2M2⟶572,VarKd(x,x)M2⟶1144. EK_d(x,x)M 14, EK_d(x,x)^2M^2 572, VarK_d(x,x)M^2 1144. (S34) More precisely, VarKd(x,x)Kd(x,x)=13+4(3x4−3x2+1)9x2(1−x2)1d+O(d−2). VarK_d(x,x)EK_d(x,x)= 13+ 4(3x^4-3x^2+1)9x^2(1-x^2) 1d+O(d^-2). (S35) For fixed x∈(0,1)x∈(0,1), this expansion holds as d→∞d→∞, but not uniformly up to the endpoints. S3 Frame coordinates, Lax evolution, and the beta law S3.A Exact frame dictionary For an input xax_a, write the SU(2)SU(2) matrix as Ua=gaI+ia⋅,ga2+‖a‖22=1.U_a=g_aI+i c_a· σ, g_a^2+\| c_a\|_2^2=1. (S36) Define the axes n and their outer product sums by Aaj=Ua,≤jZUa,≤j†=aj⋅,Sab<m=∑j=0m−1ajbj,Sab=Sab<M.A_aj=U_a,≤ jZU_a,≤ j = n_aj· σ, S_ab^<m= _j=0^m-1 n_aj n_bj T, S_ab=S_ab^<M. (S37) Every aj n_aj is a unit vector, so trSaa=MtrS_a=M. Moving the derivative in Eq. (S7) to the left gives ∂ϕjUa=iAajUa _ _jU_a=iA_ajU_a. Using (⋅)(⋅)=(⋅)I+i(×)⋅(a· σ)(b· σ)=(a·b)I+i(a×b)· σ, we obtain Jj(xa)=−aj⋅a,Kab=aSab.J_j(x_a)=- n_aj· c_a, K_ab= c_a TS_ab c_b. (S38) The corresponding second derivatives are ∂ϕjϕl2ga=−(aj⋅al)ga+(a,min(j,l)×a,max(j,l))⋅a.∂^2_ _j _lg_a=-( n_aj· n_al)g_a+( n_a, (j,l)× n_a, (j,l))· c_a. (S39) Under gradient flow, ϕ˙j=1n∑b=1neb(bj⋅b). φ_j= 1n _b=1^ne_b( n_bj· c_b). (S40) Introduce the strict-prefix and full fields a,<j=∑k<jϕ˙kak=1n∑bebSab<jb,a=a,<M. h_a,<j= _k<j φ_k n_ak= 1n _be_bS_ab^<j c_b, h_a= h_a,<M. (S41) The QSP dynamics can be rewritten in the frame coordinates: g˙a=−a⋅a,˙a=gaa−a×a,˙aj=2aj×a,<j. g_a=- h_a· c_a, c_a=g_a h_a- h_a× c_a, n_aj=2 n_aj× h_a,<j. (S42) In matrix form, the last equation is A˙aj=i[Ha,<j,Aaj],Ha,<j=a,<j⋅. A_aj=i[H_a,<j,A_aj], H_a,<j= h_a,<j· σ. (S43) This is a Lax form. Each AajA_aj evolves by conjugation, making its fixed spectrum +1,−1\+1,-1\ and the constraint ‖aj‖=1\| n_aj\|=1 manifest. It also makes the exact preservation of adjacent bond angles transparent. dt(aj⋅a,j+1)=2(aj×a,<j)⋅a,j+1+2aj⋅(a,j+1×a,<j+1). ddt( n_aj· n_a,j+1)=2( n_aj× h_a,<j)· n_a,j+1+2 n_aj·( n_a,j+1× h_a,<j+1). (S44) Using a,<j+1=a,<j+ϕ˙jaj h_a,<j+1= h_a,<j+ φ_j n_aj, the two terms cancel and hence dt(aj⋅a,j+1)=0. ddt( n_aj· n_a,j+1)=0. (S45) Evaluating the conserved inner product at initialization gives aj⋅a,j+1=12tr(AajAa,j+1)=2xa2−1. n_aj· n_a,j+1= 12tr(A_ajA_a,j+1)=2x_a^2-1. (S46) S3.B Fixed-input Haar limit Fix x∈(0,1)x∈(0,1) and let T=eiϕZT=\e^iφ Z\. The support of one random step is the circle W(x)TW(x)T. Products of support elements and their inverses generate both T and W(x)TW(x)−1W(x)TW(x)^-1. At an interior input the axes of these two tori are neither parallel nor antiparallel; their infinitesimal generators and commutator generate (2) su(2). The support is also not contained in a coset of a proper closed normal subgroup. The step law is therefore adapted and aperiodic, and the Itô–Kawada theorem [10, 2] gives Ud(x,ϕ)⇒d→∞dist.H,H∼Haar(SU(2)).U_d(x, φ) [d→∞]dist.H, H (SU(2)). (S47) The independent left factor eiϕ0Ze^i _0Z does not alter the limit. For Haar H∈SU(2)H∈ SU(2), H11H_11 is uniform in the unit disk. Hence gd(x,ϕ)⇒d→∞dist.G≡ReH11,pG(g)=2π1−g2,−1<g<1.g_d(x, φ) [d→∞]dist.G _11, p_G(g)= 2π 1-g^2, -1<g<1. (S48) S3.C Frame isotropization and contraction At this fixed input, set q=2x2−1q=2x^2-1 and define Xj=jj−13I,Bm=∑j=0m−1Xj,λ2=P2(q)=6x4−6x2+1.X_j= n_j n_j T- 13I, B_m= _j=0^m-1X_j, _2=P_2(q)=6x^4-6x^2+1. (S49) XjX_j measures the anisotropy of the frame at site j and BmB_m measures the accumulated anisotropy over the first m sites. By Eq. (S46), j⋅j+1=q n_j· n_j+1=q. Conditioned on j n_j, the uniform phase makes the azimuth of j+1 n_j+1 uniform, so j+1=qj+1−q2j, n_j+1=q n_j+ 1-q^2\, u_j, (S50) where j u_j is uniform on the unit circle in the plane orthogonal to j n_j. Hence [j∣j]=0,[jj∣j]=12(I−jj).E[ u_j n_j]=0, [ u_j u_j T n_j]= 12 (I- n_j n_j T ). (S51) It follows that [j+1j+1∣j] [ n_j+1 n_j+1 T n_j] =q2jj+1−q22(I−jj) =q^2 n_j n_j T+ 1-q^22 (I- n_j n_j T ) =13I+λ2Xj, = 13I+ _2X_j, (S52) and therefore [Xj+1∣j]=λ2Xj.E[X_j+1 n_j]= _2X_j. (S53) Iterating this Markov relation and using tr(XjXj)=2/3tr(X_j TX_j)=2/3 gives tr(XjXk)=23λ2|j−k|.Etr(X_j TX_k)= 23 _2^|j-k|. (S54) For a real matrix A, let ‖A‖F=[tr(AA)]1/2\|A\|_F=[tr(A TA)]^1/2 denote the Frobenius norm and ‖A‖op=sup‖v‖2=1‖Av‖2\|A\|_op= _\|v\|_2=1\|Av\|_2 the operator norm, where ‖v‖2=(vv)1/2\|v\|_2=(v Tv)^1/2. Summing the covariance gives the exact finite-m identity ‖Bm‖F2=23[m+2∑r=1m−1(m−r)λ2r].E\|B_m\|_F^2= 23 [m+2 _r=1^m-1(m-r) _2^r ]. (S55) For fixed x∈(0,1)x∈(0,1), |λ2|<1| _2|<1, so the right-hand side grows only linearly with m. Therefore, SxxM=I3+BMM→d→∞L2I3. S_xM= I3+ B_MM [d→∞]L^2 I3. (S56) Using Eq. (S38) and =1−gd(x)2≤1 c T c=1-g_d(x)^2≤ 1, |3Kd(x,x)M−(1−gd(x)2)|≤3‖SxxM−I3‖op≤3‖SxxM−I3‖F. | 3K_d(x,x)M-(1-g_d(x)^2) |≤ 3 \| S_xM- I3 \|_op≤ 3 \| S_xM- I3 \|_F. (S57) Equations (S56) and (S57) show that the difference on the left tends to zero in mean square, and hence in probability. Together with gd(x)⇒Gg_d(x) G, this gives 3Kd(x,x)M⇒1−G2. 3K_d(x,x)M 1-G^2. (S58) Let Y=1−G2Y=1-G^2. The two inverse branches are G=±1−YG=± 1-Y, with |dGdY|=121−Y. | dGdY |= 12 1-Y. Using pG(g)=21−g2/πp_G(g)=2 1-g^2/π, we obtain pY(y)=2pG(1−y)|dy1−y|=2πy1/2(1−y)−1/2,0<y<1.p_Y(y)=2p_G ( 1-y ) | ddy 1-y |= 2πy^1/2(1-y)^-1/2, 0<y<1. (S59) This is the density of Beta(3/2,1/2)Beta(3/2,1/2), and therefore 3Kd(x,x)M⇒d→∞dist.Beta(32,12),x∈(0,1)fixed. 3K_d(x,x)M [d→∞]dist.Beta\! ( 32, 12 ), x∈(0,1)\ fixed. (S60) Because 0≤Kd/M≤10≤ K_d/M≤ 1, the first two moments also converge, reproducing Eq. (S34). At x=1x=1, KdM=sin2(∑jϕj)∼Beta(12,12) K_dM= ^2 ( _j _j ) \! ( 12, 12 ) exactly at every depth. At x=0x=0, odd d gives Kd=0K_d=0, whereas even d again gives the beta-(1/2,1/2)(1/2,1/2) law. S4 Spectrum of the ensemble-mean kernel We now regard the ensemble-mean kernel as the kernel of an integral operator and determine its eigenvalues and eigenfunctions. Write x=cosθx= θ and work on the physical half interval θ∈[0,π/2]θ∈[0,π/2] with normalized measure dμ(θ)=2πdθ.dμ(θ)= 2πdθ. (S61) The ensemble-mean kernel defines a self-adjoint positive finite-rank operator ¯d K_d on L2([0,π/2],dμ)L^2([0,π/2],dμ) by (¯dh)(α)=∫0π/2Kd(cosα,cosβ)h(β)μ(β).( K_dh)(α)= _0^π/2EK_d( α, β)h(β)dμ(β). (S62) Thus the spectrum considered below is the operator spectrum of ¯d K_d. We seek nonzero functions h and scalars λ satisfying the eigenvalue equation ¯dh=λh,equivalently∫0π/2Kd(cosα,cosβ)h(β)μ(β)=λh(α). K_dh=λ h, _0^π/2EK_d( α, β)h(β)dμ(β)=λ h(α). (S63) For the parity fixed by d, use the orthonormal basis of L2([0,π/2],dμ)L^2([0,π/2],dμ) ψm(θ)=1,m=0,2cos(mθ),m>0,m≥0,m≡d(mod2). _m(θ)= cases1,&m=0,\\ 2 (mθ),&m>0, cases m≥ 0, m≡ d 2. (S64) Using cosdδ=2−d∑k=0d(dk)cos[(d−2k)δ] ^dδ=2^-d _k=0^d dk [(d-2k)δ] (S65) together with cosm(α−β)+cosm(α+β)=2cos(mα)cos(mβ), m(α-β)+ m(α+β)=2 (mα) (mβ), (S66) the ensemble-mean kernel has the finite spectral expansion Kd(cosα,cosβ)=∑0≤m≤dm≡d(mod 2)λmψm(α)ψm(β),EK_d( α, β)= _ subarrayc0≤ m≤ d\\ m≡ d\ ( mod\ 2) subarray _m _m(α) _m(β), (S67) where λm=M2 2−d(d(d−m)/2),m=d,d−2,…≥0. _m= M2\,2^-d d(d-m)/2, m=d,d-2,…≥ 0. (S68) Orthogonality therefore gives ¯dψm=λmψm,m=d,d−2,…≥0. K_d _m= _m _m, m=d,d-2,…≥ 0. (S69) These are all the nonzero eigenvalues of ¯d K_d. In the parity-matched orthonormal basis, the higher-frequency modes ψd+2,ψd+4,… _d+2, _d+4,… have eigenvalue zero. Equivalently, on the full Hilbert space, ker¯d=spanψd,ψd−2,…⟂. K_d=span\ _d, _d-2,…\ . (S70) The operator consequently has rank D=⌊d2⌋+1.D= d2 +1. (S71) To quantify how many of the D nonzero eigenvalues carry appreciable spectral weight, define the participation ratio rpart=(∑mλm)2∑mλm2,r_part= ( _m _m)^2 _m _m^2, (S72) where both sums run over the nonzero spectrum in Eq. (S68). Let pc=2−d(d/2)for even d,qc=4−d(2d).p_c=2^-d dd/2 even d, q_c=4^-d 2dd. (S73) Using ∑k=0d(dk)2=(2d) _k=0^d dk^2= 2dd gives rpart=4d2(2d),dodd,(1+pc)22(qc+pc2),deven.r_part= cases 4^d2 2dd,&d\ odd,\\[5.69054pt] (1+p_c)^22(q_c+p_c^2),&d\ even. cases (S74) Stirling’s formula k!∼2πk(k/e)k! 2π k(k/e)^k yields rpart∼πd2.r_part π d2. (S75) Thus, although the ensemble-mean integral operator has exact rank D=O(d)D=O(d), its binomial spectral weight is concentrated on only O(d)O( d) eigenfunctions. This effective spectral bandwidth is the operator counterpart of the O(d−1/2)O(d^-1/2) angular correlation length and motivates the sampling density w=n/dw=n/ d. S5 Scalar flow and convergence times In the sparse regime, introduce the rescaled time τ=Mt3n.τ= Mt3n. (S76) In the sparse limit, isolated inputs satisfy the independent scalar equations dg¯adτ=−(1−g¯a2)(g¯a−fa),g¯a(0)=gd(xa,ϕ(0)). d g_adτ=-(1- g_a^2)( g_a-f_a), g_a(0)=g_d(x_a, φ(0)). (S77) For |f|<1|f|<1, define Ff(u)=−log(1−u)2(1−f)−log(1+u)2(1+f)+log|u−f|1−f2.F_f(u)=- (1-u)2(1-f)- (1+u)2(1+f)+ |u-f|1-f^2. (S78) Separation of variables gives the exact implicit solution Ff(g¯(τ))−Ff(g¯(0))=−τ.F_f( g(τ))-F_f( g(0))=-τ. (S79) For f=0f=0 this becomes g¯(τ)=g0g02+(1−g02)e2τ. g(τ)= g_0 g_0^2+(1-g_0^2)e^2τ. (S80) For the case with a single zero-target input, let τq _q be the first time at which its squared residual is reduced by the factor q∈(0,1)q∈(0,1). Then τq=12log1−qg02q(1−g02),tq=3nMτq. _q= 12 1-qg_0^2q(1-g_0^2), t_q= 3nM _q. (S81) For several inputs, Eq. (S79) determines a loss-threshold time implicitly. A uniform interior condition gives explicit bounds. Put L¯=(2n)−1∑a(g¯a−fa)2 L=(2n)^-1 _a( g_a-f_a)^2. If maxamax|g¯a(0)|,|fa|≤m∗<1, _a \| g_a(0)|,|f_a|\≤ m_*<1, (S82) each trajectory stays between its initial value and target, and −2L¯≤L¯′≤−2(1−m∗2)L¯.-2 L≤ L ≤-2(1-m_*^2) L. (S83) Consequently, the scalar convergence time defined by L¯(τq)=qL¯(0) L( _q)=q L(0) obeys 12log1q≤τq≤12(1−m∗2)log1q. 12 1q≤ _q≤ 12(1-m_*^2) 1q. (S84) Multiplying by 3n/M3n/M converts these bounds to physical time. For an endpoint target f=s∈−1,1f=s∈\-1,1\, set z=(1+sg¯)/(1−sg¯)z=(1+s g)/(1-s g). The exact solution is z+logz=z0+logz0+4τ,g¯=sz−1z+1,z+ z=z_0+ z_0+4τ, g=s z-1z+1, (S85) and the endpoint residual decays algebraically rather than exponentially. S6 Rigorous sparse limit and convergence-time control S6.A Universal speed limit The parameter-shift identity, or directly unitarity in Eq. (S7), gives |Jj(x)|≤1|J_j(x)|≤ 1. Hence Kaa≤M,λmax(K)≤trK≤nM.K_a≤ M, _ (K) ≤ nM. (S86) Together with the exact loss identity this implies −L˙=1n2K≤2ML,L(t)≥L(0)e−2Mt.- L= 1n^2 e TK e≤ 2ML, L(t)≥ L(0)e^-2Mt. (S87) Therefore every relative-loss convergence time satisfies tq≥12Mlog1q.t_q≥ 12M 1q. (S88) S6.B Control equation and midpoint grid For the sparse theorem, use the half-Chebyshev midpoint grid θa=(2a+1)π4n,xa=cosθa,a=0,…,n−1. _a= (2a+1)π4n, x_a= _a, a=0,…,n-1. (S89) Use the rescaled time in Eq. (S76), with primes denoting d/dτd/dτ, and attach input and time indices to the frame fluctuations in Eq. (S49): Xa,j(τ)=aj(τ)aj(τ)−13I,Ba,m(τ)=∑j=0m−1Xa,j(τ),s(τ)=max0≤a<n0≤m≤M‖Ba,m(τ)‖FM.X_a,j(τ)= n_aj(τ) n_aj(τ) T- 13I, B_a,m(τ)= _j=0^m-1X_a,j(τ), s(τ)= _ subarrayc0≤ a<n\\ 0≤ m≤ M subarray \|B_a,m(τ)\|_FM. (S90) Here s(τ)s(τ) is the largest normalized anisotropy among all same-input frame prefixes. Substituting Saa=MI/3+Ba,MS_a=MI/3+B_a,M into the exact frame flow separates the scalar contribution from the diagonal-frame fluctuation and the off-diagonal coupling: ga′ g_a =−(1−ga2)ea+ρadiag+ρaoff, =-(1-g_a^2)e_a+ _a^diag+ _a^off, (S91) ρadiag _a^diag =−3eaMaBa,M(τ)a,ρaoff=−3M∑b≠aebaSab(τ)b. =- 3e_aM c_a TB_a,M(τ) c_a, _a^off=- 3M _b≠ ae_b c_a TS_ab(τ) c_b. (S92) The first term is the scalar flow in Eq. (S77). The proof bounds ρdiagρ^diag through the prefix anisotropies Ba,mB_a,m and ρoffρ^off through the mixed-input sums Sab<mS_ab^<m. Each bound has an initialization estimate and a deterministic transport estimate along the exact training trajectory. S6.C Diagonal-frame anisotropy Initialization.— Write s0=s(0)s_0=s(0), λ2,a=P2(cos2θa) _2,a=P_2( 2 _a), and γa=1−λ2,a=32sin2(2θa). _a=1- _2,a= 32 ^2(2 _a). (S93) All frame quantities in this paragraph are evaluated at τ=0τ=0. Equation (S53) gives [Xa,j∣Xa,j−1]=λ2,aXa,j−1E[X_a,j X_a,j-1]= _2,aX_a,j-1. For 1≤j≤M−11≤ j≤ M-1, define Da,j=Xa,j−λ2,aXa,j−1,ℳa,r=∑j=1rDa,j,ℳa,0=0.D_a,j=X_a,j- _2,aX_a,j-1, _a,r= _j=1^rD_a,j, _a,0=0. (S94) Then ℳa,rM_a,r is a matrix-valued martingale in the Frobenius inner-product space, and summing the definition of Da,jD_a,j gives, for 1≤m≤M1≤ m≤ M, γaBa,m=Xa,0−λ2,aXa,m−1+ℳa,m−1. _aB_a,m=X_a,0- _2,aX_a,m-1+M_a,m-1. (S95) Conditional Jensen gives [‖ℳa,r+1‖F∣ℱa,r]≥‖[ℳa,r+1∣ℱa,r]‖F=‖ℳa,r‖F,E[\|M_a,r+1\|_F _a,r]≥\|E[M_a,r+1 _a,r]\|_F=\|M_a,r\|_F, (S96) where ℱa,rF_a,r is the natural filtration. Thus ‖ℳa,r‖F\|M_a,r\|_F is a nonnegative scalar submartingale. The L2L^2 Doob maximal inequality gives [maxr≤M−1‖ℳa,r‖F2]1/2≤2[‖ℳa,M−1‖F2]1/2. [E _r≤ M-1\|M_a,r\|_F^2 ]^1/2≤ 2 [E\|M_a,M-1\|_F^2 ]^1/2. (S97) Using Eq. (S54), ‖Da,j‖F2=23(1−λ2,a2).E\|D_a,j\|_F^2= 23(1- _2,a^2). (S98) The zero conditional means make distinct martingale differences orthogonal in L2L^2, and hence ‖ℳa,M−1‖F2=∑j=1M−1‖Da,j‖F2=2(M−1)3(1−λ2,a2)≤2M3(1−λ2,a2).E\|M_a,M-1\|_F^2= _j=1^M-1E\|D_a,j\|_F^2= 2(M-1)3(1- _2,a^2)≤ 2M3(1- _2,a^2). (S99) Moreover, since λ2,a∈[−1/2,1] _2,a∈[-1/2,1], ‖Xa,0−λ2,aXa,m−1‖F2=23(1+λ2,a2)−2λ2,a[(a0⋅a,m−1)2−13]≤2.\|X_a,0- _2,aX_a,m-1\|_F^2= 23(1+ _2,a^2)-2 _2,a [( n_a0· n_a,m-1)^2- 13 ]≤ 2. (S100) Applying ‖A+B‖F2≤2‖A‖F2+2‖B‖F2\|A+B\|_F^2≤ 2\|A\|_F^2+2\|B\|_F^2 to Eq. (S95) and then using Eqs. (S97)–(S100) gives maxm≤M‖Ba,m‖F2≤163M1+λ2,a1−λ2,a+4(1−λ2,a)2.E _m≤ M\|B_a,m\|_F^2≤ 163M 1+ _2,a1- _2,a+ 4(1- _2,a)^2. (S101) The midpoint grid satisfies the exact identities ∑a=0n−1γa−1=2n23,∑a=0n−1γa−2=4(n4+2n2)27. _a=0^n-1 _a^-1= 2n^23, _a=0^n-1 _a^-2= 4(n^4+2n^2)27. (S102) Set Za=maxm≤M‖Ba,m‖FZ_a= _m≤ M\|B_a,m\|_F. Since s0=M−1maxaZas_0=M^-1 _aZ_a, a union bound over the inputs followed by Markov’s inequality gives Prob.(s0>εs)≤∑aProb.(Za2>M2εs2)≤1M2εs2∑aZa2.Prob.(s_0> _s)≤ _aProb.(Z_a^2>M^2 _s^2)≤ 1M^2 _s^2 _aEZ_a^2. (S103) Using (1+λ2,a)/(1−λ2,a)=2γa−1−1(1+ _2,a)/(1- _2,a)=2 _a^-1-1 and Eq. (S102), we obtain Prob.(s0>εs)≤1εs2[649n2M−163nM+1627n4+2n2M2].Prob.(s_0> _s)≤ 1 _s^2 [ 649 n^2M- 163 nM+ 1627 n^4+2n^2M^2 ]. (S104) Transport.— The phase-energy identity follows directly from gradient flow in rescaled time: dLdτ=−M3n‖ϕ′(τ)‖22,∫0T‖ϕ′(s)‖22s=3nM[L(0)−L(T)]. dLdτ=- M3n\| φ (τ)\|_2^2, _0^T\| φ (s)\|_2^2ds= 3nM[L(0)-L(T)]. (S105) Define the phase-derivative budget ΓT=43nTL(0)≤46nT, _T=4 3nTL(0)≤ 4 6nT, (S106) where the second inequality uses 0≤L(0)≤20≤ L(0)≤ 2. For each a and j, let Raj(τ,s)∈SO(3)R_aj(τ,s)∈ SO(3) be the fundamental solution ∂τRaj(τ,s)=Ωaj(τ)Raj(τ,s) _τR_aj(τ,s)= _aj(τ)R_aj(τ,s), Raj(s,s)=IR_aj(s,s)=I, where Ωaj(τ)v=−2[∑k<jϕk′(τ)ak(τ)]×v. _aj(τ)v=-2 [ _k<j _k (τ) n_ak(τ) ]× v. (S107) Equation (S43) and uniqueness give aj(τ)=Raj(τ,0)aj(0) n_aj(τ)=R_aj(τ,0) n_aj(0). Write Raj(τ)=Raj(τ,0)R_aj(τ)=R_aj(τ,0). For 0≤j≤M−20≤ j≤ M-2, the two generators differ by the single channel ϕj′(τ)aj(τ) _j (τ) n_aj(τ) and satisfy ‖Ωa,j+1−Ωaj‖op=2|ϕj′|\| _a,j+1- _aj\|_op=2| _j |. Variation of constants gives Ra,j+1(τ)−Raj(τ)=∫0τRa,j+1(τ,s)[Ωa,j+1(s)−Ωaj(s)]Raj(s)s.R_a,j+1(τ)-R_aj(τ)= _0^τR_a,j+1(τ,s)[ _a,j+1(s)- _aj(s)]R_aj(s)ds. (S108) Both matrices multiplying the generator difference are orthogonal, so ‖Ra,j+1(τ)−Raj(τ)‖op≤2∫0τ|ϕj′(s)|s.\|R_a,j+1(τ)-R_aj(τ)\|_op≤ 2 _0^τ| _j (s)|ds. (S109) Define ωj(τ)=4∫0τ|ϕj′(s)|s _j(τ)=4 _0^τ| _j (s)|ds for 0≤j≤M−20≤ j≤ M-2. Cauchy–Schwarz first in the phase channels and then in time gives, uniformly for τ≤Tτ≤ T, ‖(τ)‖1 \| ω(τ)\|_1 =4∫0τ∑j|ϕj′(s)|s≤4Mτ[∫0τ‖ϕ′(s)‖22s]1/2≤ΓT, =4 _0^τ _j| _j (s)|ds≤ 4 Mτ [ _0^τ\| φ (s)\|_2^2ds ]^1/2≤ _T, M‖(τ)‖2 M\| ω(τ)\|_2 ≤4Mτ[∫0τ‖ϕ′(s)‖22s]1/2≤ΓT. ≤ 4 Mτ [ _0^τ\| φ (s)\|_2^2ds ]^1/2≤ _T. (S110) At fixed τ, define ℒaj(τ)(X)=Raj(τ)XRaj(τ)L_aj^(τ)(X)=R_aj(τ)XR_aj(τ) T. For a sequence XjX_j with prefixes Am=∑j<mXjA_m= _j<mX_j and a sequence of linear maps ℒjL_j, summation by parts gives ∑j<mℒjXj=ℒm−1Am+∑j<m−1(ℒj−ℒj+1)Aj+1. _j<mL_jX_j=L_m-1A_m+ _j<m-1(L_j-L_j+1)A_j+1. (S111) These maps preserve the Frobenius norm, and Eq. (S109) gives ‖(ℒaj(τ)−ℒa,j+1(τ))(X)‖F≤ωj(τ)‖X‖F.\|(L_aj^(τ)-L_a,j+1^(τ))(X)\|_F≤ _j(τ)\|X\|_F. (S112) Applying Eq. (S111) with Am=Ba,m(0)A_m=B_a,m(0) therefore yields ‖Ba,m(τ)‖F \|B_a,m(τ)\|_F ≤‖Ba,m(0)‖F+∑j<m−1ωj(τ)‖Ba,j+1(0)‖F ≤\|B_a,m(0)\|_F+ _j<m-1 _j(τ)\|B_a,j+1(0)\|_F ≤Ms0[1+‖(τ)‖1]≤Ms0(1+ΓT). ≤ Ms_0[1+\| ω(τ)\|_1]≤ Ms_0(1+ _T). (S113) Consequently, s(τ)≤(1+ΓT)s0,0≤τ≤T.s(τ)≤(1+ _T)s_0, 0≤τ≤ T. (S114) S6.D Off-diagonal frame sums The off-diagonal remainder in Eq. (S92) is ρaoff(τ)=−3M∑b≠aeb(τ)a(τ)Sab(τ)b(τ). _a^off(τ)=- 3M _b≠ ae_b(τ) c_a(τ) TS_ab(τ) c_b(τ). (S115) Since ‖a‖22=1−ga2≤1\| c_a\|_2^2=1-g_a^2≤ 1, ‖f‖∞≡maxa|fa|≤1\|f\|_∞≡ _a|f_a|≤ 1, and therefore |eb|≤1+‖f‖∞|e_b|≤ 1+\|f\|_∞, the definition of the operator norm gives |ρaoff(τ)|≤3(1+‖f‖∞)M∑b≠a‖Sab(τ)‖op.| _a^off(τ)|≤ 3(1+\|f\|_∞)M _b≠ a\|S_ab(τ)\|_op. (S116) Thus ρaoff _a^off is controlled by operator norms of the mixed-input frame sums. Although Eq. (S115) contains only Sab=Sab<MS_ab=S_ab^<M, the transport estimate also requires the strict prefixes. Set Qa,m(τ)=∑b≠a‖Sab<m(τ)‖op,0≤m≤M,Q_a,m(τ)= _b≠ a\|S_ab^<m(τ)\|_op, 0≤ m≤ M, (S117) so that Qa,M(τ)=∑b≠a‖Sab(τ)‖opQ_a,M(τ)= _b≠ a\|S_ab(τ)\|_op, and define ℭpre(τ)2 C_pre(τ)^2 =9M3n∑a=0n−1∑m=0M−1Qa,m(τ)2, = 9M^3n _a=0^n-1 _m=0^M-1Q_a,m(τ)^2, ℭfull(τ)2 C_full(τ)^2 =9M2n∑a=0n−1Qa,M(τ)2,ℭinit2=ℭpre(0)2+ℭfull(0)2. = 9M^2n _a=0^n-1Q_a,M(τ)^2, C_init^2= C_pre(0)^2+ C_full(0)^2. (S118) These definitions give the deterministic source bounds [1nM∑a=0n−1∑m=0M−1‖3M∑b≠aebSab<m(τ)b‖22]1/2 [ 1nM _a=0^n-1 _m=0^M-1 \| 3M _b≠ ae_bS_ab^<m(τ) c_b \|_2^2 ]^1/2 ≤(1+‖f‖∞)ℭpre(τ), ≤(1+\|f\|_∞) C_pre(τ), [1n∑a=0n−1|ρaoff(τ)|2]1/2 [ 1n _a=0^n-1| _a^off(τ)|^2 ]^1/2 ≤(1+‖f‖∞)ℭfull(τ). ≤(1+\|f\|_∞) C_full(τ). (S119) Initialization.— To estimate Eq. (S118), first compute the second moment of one mixed-input prefix. From Sab<m=∑j=0m−1ajbj,S_ab^<m= _j=0^m-1 n_aj n_bj T, (S120) we obtain ‖Sab<m‖F2 \|S_ab^<m\|_F^2 =tr[(Sab<m)Sab<m] =tr [(S_ab^<m) TS_ab^<m ] =∑j,k<m(aj⋅ak)(bj⋅bk). = _j,k<m( n_aj· n_ak)( n_bj· n_bk). (S121) The required second moment is therefore determined by the joint axis correlation Cab(r)=[(a,k+r⋅ak)(b,k+r⋅bk)].C_ab(r)=E [( n_a,k+r· n_ak)( n_b,k+r· n_bk) ]. (S122) The same random phase acts at inputs a and b, so the joint correlation must be evaluated before averaging. The required paired transfer is constructed as follows. Let aW_a denote the Bloch-space rotation induced by the signal gate W(xa)W(x_a), with ax=x,ay=cay−saz,az=say+caz,ca=cos2θa,sa=sin2θa.W_ax=x, _ay=c_ay-s_az, _az=s_ay+c_az, c_a= 2 _a, s_a= 2 _a. (S123) Averaging the common phase projects a paired vector onto the subspace invariant under simultaneous rotations about z: Π=12π∫02πRz(φ)⊗Rz(φ)φ. = 12π _0^2πR_z( ) R_z( )\,d . (S124) With uvuv denoting ^⊗ u v, this invariant subspace is spanned by E0=zz,E1=xx+yy2,E2=xy−yx2.E_0=z, E_1= x+y 2, E_2= xy-yx 2. (S125) One averaged paired step is ab=Π(a⊗b)Π.T_ab= (W_a _b) . (S126) Using Eq. (S123) and Π(zz)=zz,Π(xx)=Π(yy)=xx+yy2,Π(yz)=Π(zy)=0, (z)=z, (x)= (y)= x+y2, (yz)= (zy)=0, (S127) we find abE0 _abE_0 =cacbE0+sasb2E1, =c_ac_bE_0+ s_as_b 2E_1, abE1 _abE_1 =sasb2E0+1+cacb2E1. = s_as_b 2E_0+ 1+c_ac_b2E_1. (S128) The antisymmetric vector E2E_2 is invariant under the phase projection but is decoupled from E0E_0 and E1E_1. Since the correlation in Eq. (S122) begins and ends in the direction E0=zzE_0=z, only the two-dimensional subspace spanE0,E1span\E_0,E_1\ contributes. In this ordered basis, Eq. (S128) gives Tab=(cacbsasb/2sasb/2(1+cacb)/2).T_ab= pmatrixc_ac_b&s_as_b/ 2\\ s_as_b/ 2&(1+c_ac_b)/2 pmatrix. (S129) Let e0=(1,0)e_0=(1,0) T be the coordinate vector of E0E_0. In the co-moving frames at layer k, both reference axes ak n_ak and bk n_bk are represented by z. Rotational covariance and the independence of subsequent phases therefore give Cab(r)=⟨E0,abrE0⟩=e0Tabre0.C_ab(r)= E_0,T_ab^rE_0 =e_0 TT_ab^re_0. (S130) Substituting Eq. (S130) into Eq. (S121), the j=kj=k terms contribute m, while for each separation r≥1r≥ 1 there are 2(m−r)2(m-r) ordered pairs. Hence ‖Sab<m‖F2=m+2∑r=1m−1(m−r)e0Tabre0.E\|S_ab^<m\|_F^2=m+2 _r=1^m-1(m-r)e_0 TT_ab^re_0. (S131) The eigenvalues of TabT_ab are λ±=1+3cacb±(1−cacb)2+8sa2sb24, _±= 1+3c_ac_b± (1-c_ac_b)^2+8s_a^2s_b^24, (S132) and (1−λ+)(1−λ−)=2sin2(θa−θb)sin2(θa+θb).(1- _+)(1- _-)=2 ^2( _a- _b) ^2( _a+ _b). (S133) Direct inversion gives e0(I−Tab)−1e0=14[csc2(θa−θb)+csc2(θa+θb)].e_0 T(I-T_ab)^-1e_0= 14 [ ^2( _a- _b)+ ^2( _a+ _b) ]. (S134) To bound Eq. (S131), decompose e0Tabre0e_0 TT_ab^re_0 into the two eigenmodes. Since TabT_ab is symmetric, both spectral weights of e0e_0 are nonnegative. For λ≥0λ≥ 0, ∑r=1m−1(m−r)λr≤m∑r=1∞λr≤m1−λ. _r=1^m-1(m-r)λ^r≤ m _r=1^∞λ^r≤ m1-λ. (S135) For −1<λ<0-1<λ<0, the weighted sum is alternating, begins with a negative term, and has decreasing term magnitudes, so it is nonpositive. Applying these bounds to the two nonnegative spectral weights and using Eq. (S134) gives ‖Sab<m‖F2 \|S_ab^<m\|_F^2 ≤mAab, ≤ mA_ab, Aab A_ab =1+12[csc2(θa−θb)+csc2(θa+θb)]. =1+ 12 [ ^2( _a- _b)+ ^2( _a+ _b) ]. (S136) It remains to sum this pairwise bound over the midpoint grid. Let Hk=∑r=1kr−1H_k= _r=1^kr^-1 and H0=0H_0=0. Since |θa−θb|=|a−b|π/(2n)| _a- _b|=|a-b|π/(2n) and sinu≥2u/π u≥ 2u/π on [0,π/2][0,π/2], ∑b≠a1|sin(θa−θb)|≤n(Ha+Hn−1−a)≤2nHn. _b≠ a 1| ( _a- _b)|≤ n(H_a+H_n-1-a)≤ 2nH_n. (S137) For the reflected separation, θa+θb=(a+b+1)π/(2n) _a+ _b=(a+b+1)π/(2n) and sinu=sin(π−u) u= (π-u), so ∑b=0n−11sin(θa+θb)≤2nH2n. _b=0^n-1 1 ( _a+ _b)≤ 2nH_2n. (S138) Using Hk≤1+logkH_k≤ 1+ k and 1+(u2+v2)/2≤1+(u+v)/2 1+(u^2+v^2)/2≤ 1+(u+v)/ 2, we obtain, for an absolute constant C, max∑b≠aAab≤Cn(1+logn). _a _b≠ a A_ab≤ Cn(1+ n). (S139) At initialization, Minkowski’s inequality, ‖S‖op≤‖S‖F\|S\|_op≤\|S\|_F, and Eqs. (S136) and (S139) give [Qa,m(0)2]1/2 [EQ_a,m(0)^2 ]^1/2 ≤∑b≠a[‖Sab<m(0)‖op2]1/2 ≤ _b≠ a [E\|S_ab^<m(0)\|_op^2 ]^1/2 ≤∑b≠a[‖Sab<m(0)‖F2]1/2 ≤ _b≠ a [E\|S_ab^<m(0)\|_F^2 ]^1/2 ≤m∑b≠aAab≤Cmn(1+logn). ≤ m _b≠ a A_ab≤ C m\,n(1+ n). (S140) Squaring this inequality, summing over a and m, and using ∑m=0M−1m=M(M−1)/2 _m=0^M-1m=M(M-1)/2, we find ℭpre(0)2≤C1n2M(1+logn)2,ℭfull(0)2≤C2n2M(1+logn)2.E C_pre(0)^2≤ C_1 n^2M(1+ n)^2, C_full(0)^2≤ C_2 n^2M(1+ n)^2. (S141) Consequently, for an absolute constant C0C_0, ℭinit2≤C0n2M(1+logn)2,Prob.(ℭinit>εc)≤C0n2(1+logn)2Mεc2.E C_init^2≤ C_0 n^2M(1+ n)^2, .( C_init> _c)≤ C_0n^2(1+ n)^2M _c^2. (S142) Transport.— For the mixed-input dyads, define ℒabj(τ)(X)=Raj(τ)XRbj(τ).L_abj^(τ)(X)=R_aj(τ)XR_bj(τ) T. (S143) Because the matrices RajR_aj and RbjR_bj are orthogonal, this map preserves the operator norm: ‖ℒabj(τ)(X)‖op=‖X‖op.\|L_abj^(τ)(X)\|_op=\|X\|_op. (S144) Moreover, aj(τ)bj(τ)=ℒabj(τ)(aj(0)bj(0)), n_aj(τ) n_bj(τ) T=L_abj^(τ) ( n_aj(0) n_bj(0) T ), (S145) and hence Sab<m(τ)=∑j=0m−1ℒabj(τ)(aj(0)bj(0)).S_ab^<m(τ)= _j=0^m-1L_abj^(τ) ( n_aj(0) n_bj(0) T ). (S146) For any matrix X, adding and subtracting Ra,j+1XRbjR_a,j+1XR_bj T gives (ℒabj(τ)−ℒab,j+1(τ))(X) (L_abj^(τ)-L_ab,j+1^(τ))(X) =(Raj−Ra,j+1)XRbj =(R_aj-R_a,j+1)XR_bj T +Ra,j+1X(Rbj−Rb,j+1). +R_a,j+1X(R_bj-R_b,j+1) T. (S147) All rotations in Eq. (S147) are evaluated at τ. Using ‖AXB‖op≤‖A‖op‖X‖op‖B‖op\|AXB\|_op≤\|A\|_op\|X\|_op\|B\|_op, the fact that every rotation has operator norm one, and Eq. (S109) for inputs a and b, we obtain ‖(ℒabj(τ)−ℒab,j+1(τ))(X)‖op \|(L_abj^(τ)-L_ab,j+1^(τ))(X)\|_op ≤(‖Raj−Ra,j+1‖op+‖Rbj−Rb,j+1‖op)‖X‖op ≤ (\|R_aj-R_a,j+1\|_op+\|R_bj-R_b,j+1\|_op )\|X\|_op ≤4∫0τ|ϕj′(s)|s‖X‖op=ωj(τ)‖X‖op. ≤ 4 _0^τ| _j (s)|ds\,\|X\|_op= _j(τ)\|X\|_op. (S148) At initialization, aj(0)bj(0)=Sab<j+1(0)−Sab<j(0),Sab<0(0)=0. n_aj(0) n_bj(0) T=S_ab^<j+1(0)-S_ab^<j(0), S_ab^<0(0)=0. (S149) Substituting this identity into Eq. (S146) and summing by parts gives Sab<m(τ) S_ab^<m(τ) =ℒab,m−1(τ)(Sab<m(0)) =L_ab,m-1^(τ) (S_ab^<m(0) ) +∑j=0m−2(ℒabj(τ)−ℒab,j+1(τ))(Sab<j+1(0)). + _j=0^m-2 (L_abj^(τ)-L_ab,j+1^(τ) ) (S_ab^<j+1(0) ). (S150) Equations (S144), (S148), and (S150) imply ‖Sab<m(τ)‖op≤‖Sab<m(0)‖op+∑j=0m−2ωj(τ)‖Sab<j+1(0)‖op.\|S_ab^<m(τ)\|_op≤\|S_ab^<m(0)\|_op+ _j=0^m-2 _j(τ)\|S_ab^<j+1(0)\|_op. (S151) Summing over b≠ab≠ a gives Qa,m(τ)≤Qa,m(0)+∑j=0m−2ωj(τ)Qa,j+1(0).Q_a,m(τ)≤ Q_a,m(0)+ _j=0^m-2 _j(τ)Q_a,j+1(0). (S152) For the strict prefixes, Cauchy–Schwarz in j gives [∑j=0m−2ωjQa,j+1(0)]2≤‖22∑j=0m−2Qa,j+1(0)2. [ _j=0^m-2 _jQ_a,j+1(0) ]^2≤\| ω\|_2^2 _j=0^m-2Q_a,j+1(0)^2. (S153) Summing this inequality over m<Mm<M yields ∑m=0M−1[∑j=0m−2ωjQa,j+1(0)]21/2≤M‖2[∑m=0M−1Qa,m(0)2]1/2. \ _m=0^M-1 [ _j=0^m-2 _jQ_a,j+1(0) ]^2 \^1/2≤ M\| ω\|_2 [ _m=0^M-1Q_a,m(0)^2 ]^1/2. (S154) Applying the triangle inequality in the prefix ℓ2 ^2 norm, then summing over a and using Eq. (S110), gives ℭpre(τ)≤(1+ΓT)ℭpre(0). C_pre(τ)≤(1+ _T) C_pre(0). (S155) For the full sums, Eq. (S152) with m=Mm=M gives Qa,M(τ)≤Qa,M(0)+∑j=0M−2ωjQa,j+1(0).Q_a,M(τ)≤ Q_a,M(0)+ _j=0^M-2 _jQ_a,j+1(0). (S156) Cauchy–Schwarz in j, followed by the sum over a, gives ∑a=0n−1[∑j=0M−2ωjQa,j+1(0)]21/2≤‖2[∑a=0n−1∑m=0M−1Qa,m(0)2]1/2. \ _a=0^n-1 [ _j=0^M-2 _jQ_a,j+1(0) ]^2 \^1/2≤\| ω\|_2 [ _a=0^n-1 _m=0^M-1Q_a,m(0)^2 ]^1/2. (S157) After applying the normalizations in Eq. (S118), the factor multiplying ℭpre(0) C_pre(0) is M‖2≤ΓT M\| ω\|_2≤ _T. Therefore, ℭfull(τ)≤ℭfull(0)+ΓTℭpre(0). C_full(τ)≤ C_full(0)+ _T C_pre(0). (S158) S6.E Forcing and finite-time trajectory control Define the initialization event M,n=s0≤εs,ℭinit≤εc.G_M,n=\s_0≤ _s,\ C_init≤ _c\. (S159) Equations (S104) and (S142) give Prob.(M,nc)≤1εs2[649n2M−163nM+1627n4+2n2M2]+C0n2(1+logn)2Mεc2.Prob.(G_M,n^c)≤ 1 _s^2 [ 649 n^2M- 163 nM+ 1627 n^4+2n^2M^2 ]+ C_0n^2(1+ n)^2M _c^2. (S160) The diagonal-frame remainder in Eq. (S92) satisfies ‖diag(τ)‖2n≤3(1+‖f‖∞)s(τ)≤3(1+‖f‖∞)(1+ΓT)s0. \| ρ^diag(τ)\|_2 n≤ 3(1+\|f\|_∞)s(τ)≤ 3(1+\|f\|_∞)(1+ _T)s_0. (S161) The full-prefix bound in Eq. (S119), together with Eq. (S158), gives ‖off(τ)‖2n≤(1+‖f‖∞)[ℭfull(0)+ΓTℭpre(0)]. \| ρ^off(τ)\|_2 n≤(1+\|f\|_∞) [ C_full(0)+ _T C_pre(0) ]. (S162) Set =diag+off ρ= ρ^diag+ ρ^off and define βT=(1+‖f‖∞)[3(1+ΓT)s0+ℭfull(0)+ΓTℭpre(0)], _T=(1+\|f\|_∞) [3(1+ _T)s_0+ C_full(0)+ _T C_pre(0) ], (S163) Then supτ≤T‖(τ)‖2n≤βT. _τ≤ T \| ρ(τ)\|_2 n≤ _T. (S164) Let g¯a g_a solve Eq. (S77) with the same initial output as the full QSP flow. For Fa(u)=−(1−u2)(u−fa)F_a(u)=-(1-u^2)(u-f_a), |Fa′(u)|=|3u2−2fau−1|≤2+2‖f‖∞≡Λ≤4,−1≤u≤1.|F_a (u)|=|3u^2-2f_au-1|≤ 2+2\|f\|_∞≡ ≤ 4, -1≤ u≤ 1. (S165) Subtracting the two output equations and writing the Grönwall estimate in integral form gives ‖full(τ)−¯(τ)‖2n \| g^full(τ)- g(τ)\|_2 n ≤∫0τeΛ(τ−s)‖(s)‖2ns ≤ _0^τe (τ-s) \| ρ(s)\|_2 nds ≤eΛτ−1ΛβT≤τeΛτβT. ≤ e τ-1 _T≤τ e τ _T. (S166) Consequently, sup0≤τ≤T‖full(τ)−¯(τ)‖2n≤TeΛTβT. _0≤τ≤ T \| g^full(τ)- g(τ)\|_2 n≤ Te T _T. (S167) S6.F Sparse-limit theorem For fixed T, Eq. (S106) gives ΓT=O(n) _T=O( n). We impose n3/2(1+logn)M⟶0. n^3/2(1+ n) M 0. (S168) This condition implies (n2M+n4M2)1/2=o(n−1/2),n(1+logn)M=o(n−1/2). ( n^2M+ n^4M^2 )^1/2=o(n^-1/2), n(1+ n) M=o(n^-1/2). (S169) For the n4/M2n^4/M^2 term, one uses n5/2M=1n(n3/2M)2⟶0. n^5/2M= 1 n ( n^3/2 M )^2 0. (S170) We may therefore choose εs _s and εc _c so that the two initialization scales in Eq. (S169) are respectively o(εs)o( _s) and o(εc)o( _c), while nεs→0 n\, _s→ 0 and nεc→0 n\, _c→ 0. Equation (S160) then gives Prob.(M,n)→1Prob.(G_M,n)→ 1, while Eq. (S163) gives βT→0 _T→ 0 on M,nG_M,n. We obtain the following statement. Theorem.— Let M=d+1→∞M=d+1→∞. Use the grid in Eq. (S89), independent uniform initial phases, deterministic targets with ‖f‖∞≤1\|f\|_∞≤ 1, and exact continuous gradient flow. If Eq. (S168) holds, then for every fixed T<∞T<∞, sup0≤τ≤T‖full(τ)−¯(τ)‖2n→M→∞Prob.0. _0≤τ≤ T \| g^full(τ)- g(τ)\|_2 n [M→∞]Prob.0. (S171) The same transport estimate controls the time-dependent diagonal kernel. Equation (S38) gives 3Kaa(τ)M−(1−ga(τ)2)=3Ma(τ)Ba,M(τ)a(τ), 3K_a(τ)M-(1-g_a(τ)^2)= 3M c_a(τ) TB_a,M(τ) c_a(τ), (S172) whose right-hand side has absolute value at most 3s(τ)3s(τ). Equations (S104), (S114), and (S168) therefore imply sup0≤τ≤T[1n∑a=0n−1|3Kaa(τ)M−(1−ga(τ)2)|2]1/2→M→∞Prob.0. _0≤τ≤ T [ 1n _a=0^n-1 | 3K_a(τ)M-(1-g_a(τ)^2) |^2 ]^1/2 [M→∞]Prob.0. (S173) The leading diagonal of the tangent kernel therefore changes with the learned output throughout the controlled time interval. The factor n3/2(1+logn)n^3/2(1+ n) in Eq. (S168) combines the midpoint-grid sum of mixed-input prefixes with their transport over a fixed rescaled-time interval. S6.G Transfer of convergence times Let τqsc<T _q^sc<T be the first time at which the scalar loss reaches qL(0)qL(0), where 0<q<10<q<1 and the common initial loss satisfies L(0)>0L(0)>0. Define ϵT=supτ≤T‖full(τ)−¯(τ)‖2n. _T= _τ≤ T \| g^full(τ)- g(τ)\|_2 n. (S174) Since all outputs and targets lie in [−1,1][-1,1], |Lfull(τ)−L¯(τ)|≤2ϵT.|L^full(τ)- L(τ)|≤ 2 _T. (S175) For the relative losses rfull(τ)=Lfull(τ)/L(0)r^full(τ)=L^full(τ)/L(0) and r(τ)=L¯(τ)/L(0)r(τ)= L(τ)/L(0), Eq. (S175) gives |rfull(τ)−r(τ)|≤2ϵTL(0).|r^full(τ)-r(τ)|≤ 2 _TL(0). (S176) Suppose τqsc _q^sc has a neighborhood contained in (0,T)(0,T) on which −r′(τ)≥νq>0-r (τ)≥ _q>0. For sufficiently small ϵT _T, continuity, monotonicity, and Eq. (S176) imply that the full QSP loss also crosses qL(0)qL(0) within this neighborhood. Denoting its first crossing time by τqfull _q^full gives |τqfull−τqsc|≤2ϵTνqL(0).| _q^full- _q^sc|≤ 2 _T _qL(0). (S177) Under the interior condition in Eq. (S82), Eq. (S84) places τqsc _q^sc in a fixed rescaled-time interval. Choosing T strictly above its upper endpoint gives the required neighborhood, while −r′(τqsc)≥2(1−m∗2)q.-r ( _q^sc)≥ 2(1-m_*^2)q. (S178) By continuity, the slope condition holds locally for any νq<2(1−m∗2)q _q<2(1-m_*^2)q. Hence, if the initial loss stays bounded away from zero in probability, Eqs. (S171) and (S177) imply |τqfull−τqsc|→M→∞Prob.0.| _q^full- _q^sc| [M→∞]Prob.0. (S179) Multiplication by 3n/M3n/M converts Eq. (S177) to physical time. S7 Numerical protocol Figure 2 uses (d,n,w)=(1024,4,0.125)(d,n,w)=(1024,4,0.125) and (128,64,5.66)(128,64,5.66), where w=n/dw=n/ d. For each parameter choice, we draw 3030 independent uniform-phase initializations on the grid in Eq. (S89). Each initialization is paired with an independently drawn degree-88 QSP teacher whose response on the same grid defines the targets. With Jacobian matrix J(ϕ)aj=Jj(xa,ϕ)J( φ)_aj=J_j(x_a, φ), training uses the nonlinear update ϕ(r+1)=ϕ(r)−ηnJ(ϕ(r))(ϕ(r)),η=0.1M, φ^(r+1)= φ^(r)- ηnJ( φ^(r)) T e( φ^(r)), η= 0.1M, (S180) with the forward model and Jacobian reevaluated at every step, and physical time t=rηt=rη. The displayed trajectory has the finite QSP convergence time closest to the median among the converged runs; the shaded band is the pointwise interquartile range over all 3030 trajectories. The scalar reference solves Eq. (S77) from the outputs of the displayed QSP initialization. For the frozen-kernel reference, let K0=Vdiag(λr)VK_0=Vdiag( _r)V T and wr=(Vr0)2w_r=(V_r T e_0)^2. Its continuous-flow loss is LK0(t)LK0(0)=∑rwre−2λrt/n∑rwr. L_K_0(t)L_K_0(0)= _rw_re^-2 _rt/n _rw_r. (S181) The remaining reference is the speed-limit envelope e−2Mte^-2Mt. Convergence is defined by L/L(0)=10−4L/L(0)=10^-4, with logarithmic interpolation between adjacent iterations and a cap of 5×1055× 10^5 steps. Two higher-density runs remain above the threshold at the cap; they remain in the trajectory statistics, and no higher-density convergence-time ratios are reported.