QuanVI: Score-based Variational Inference via Quantum Maximally Mixed States

Conference: NeurIPS2026
arXiv: 2609.39164
Area: Optimization & Theory
Keywords: score matching, variational inference, maximally mixed states, quantum tensor networks, local Markov dependence

TL;DR

QuanVI replaces EigenVI's single-eigenvector solution with a maximally mixed state on a low-energy subspace and compresses the density operator through local quantum tensor networks, scaling to 100-dimensional chain-structured synthetic targets while remaining limited by local-window capacity and cost.

Background & Motivation

Variational inference approximates a complex posterior with a computationally tractable distribution. Conventional methods usually minimize KL divergence or maximize the ELBO; score-based variational inference instead compares gradients of log densities, so the target's normalization constant disappears from its score. GSM and BaM use Gaussian variational families, which have limited capacity for multimodal, ring-shaped, or funnel geometry. EigenVI represents a density as the squared amplitude of an orthogonal function expansion and turns its score objective into a minimum-eigenvalue problem. This goes beyond a single Gaussian, but using K basis functions per coordinate requires K to the power D coefficients and K to the power 2D entries in the global operator.

Storage is not the only problem: a single eigenvector is also an unstable representation unit. When the smallest eigenvalue is degenerate, a solver can return different directions within its eigenspace, and squaring their amplitudes can produce different densities. Near degeneracy, sampling noise or numerical perturbations can change the selected direction. Compressing only the coefficient tensor does not automatically resolve this non-uniqueness; introducing only a mixed state still leaves an exponentially large density matrix.

QuanVI therefore changes both what constitutes a solution and how that solution is stored. It replaces one direction with a uniform density operator on a subspace and exploits local scores under chain-structured Markov dependence to avoid constructing a global training matrix. Core Idea: learn a maximally mixed state on a low-energy subspace, compress that subspace with orthogonality-constrained local tensor cores, and use tensor contractions for training and probability queries instead of explicitly solving for a high-dimensional global eigenvector.

Method

Overall Architecture

The inputs are the target score function and orthogonal basis functions for each coordinate; the experiments use Hermite bases. The output is a queryable variational probability density, not an image or a diffusion model. Per-coordinate basis vectors form tensor-product features, and a learned density operator acts on this feature space.

Four interdependent designs define the method: the maximally mixed state specifies the variational family, local score operators specify the decomposition of the training objective, the quantum tensor network specifies compact operator storage, and measurement contractions implement density, marginal, and conditional queries after training. Here, “quantum” refers to complex tensors and orthogonality constraints evaluated on classical computers, not a requirement for quantum hardware; the appendix states that all experiments run on CPUs.

These designs primarily involve operator algebra and tensor contractions rather than serial neural-network modules, so a pipeline diagram would not clarify their matrix mechanism. Training contractions insert local loss operators built from the target score; probability-query contractions insert measurement matrices built from basis functions. These are distinct inputs, not interchangeable forms of the same data flow.

Key Designs

1. Maximally mixed state: replace a non-unique direction with a subspace

EigenVI's pure state can be written as the outer product of a unit vector, with the corresponding squared amplitude defining its density. QuanVI instead learns a complex matrix U with orthonormal columns, each representing one orthogonal direction in a subspace. Equal weights on these directions produce a positive semidefinite operator with rank r and trace 1; sandwiching this operator between tensor-product features gives the density.

\[ \rho=\frac{1}{r}UU^{\dagger},\qquad U^{\dagger}U=I_r,\qquad q_{\rho}(x)=\Phi(x)^{\dagger}\rho\Phi(x). \]

Here, U has \(K^D\) rows and r columns, and \(\Phi(x)\) is the tensor product of the per-coordinate orthogonal basis vectors. Positive semidefiniteness makes the density nonnegative; basis orthonormality and unit trace make its integral equal to 1. “Maximally mixed” means uniform mixing within the selected subspace, not a default identity matrix on the entire feature space and not an arbitrary learned set of mixture weights.

Replacing U with another orthonormal basis of the same subspace leaves the sum of outer products unchanged, removing dependence on a particular degenerate eigenvector. This has an important boundary: without tensor-network restrictions, the fixed-rank linear trace objective selects the average of the lowest r eigen-directions, rather than ensuring that every direction attains the global minimum eigenvalue. It can represent the entire lowest-energy eigenspace stably when r matches that degenerate subspace and it is separated from higher energies; near degeneracy and rank mismatch still affect the result. Restricted QTN expressiveness and nonconvex optimization further limit whether this ideal solution is reached.

A second theoretical boundary must also remain explicit: Appendix A.1 derives the quadratic Fisher objective for a real pure state, after which the main text extends that quadratic form to a linear trace objective over mixed states. The supplied full text does not prove that this equals the Fisher divergence of the resulting mixed density itself, so the two should not be equated unconditionally.

2. Local score operators: avoid an exponentially large global matrix

The paper uses a path graph along an ordering of variables as a representative source of locality: the target density is a product of pairwise factors between adjacent variables. An interior coordinate's score needs only its left neighbor, itself, and its right neighbor; endpoint scores need only two neighboring variables. Expressing the score residual in the orthogonal basis makes each coordinate's operator act on the corresponding local feature space, with orthogonal integration yielding identities elsewhere.

Training thus contracts local operators term by term rather than storing a complete \(K^D\times K^D\) matrix. Embedding the local terms in the full space and summing them recovers the main text's global trace objective. Section 5 suppresses this embedding notation; a small local matrix should not be mistaken for a matrix already matching the global density operator's dimensions.

This mechanism depends on actual locality of the target score, not merely on strong tensor-network compression. A general dense dependency graph does not automatically admit the same local-operator decomposition, and variable ordering also affects compressibility. Two different windows must be distinguished: an interior path-graph score involves three variables, whereas the experimental default L=2 means each learnable QTN core acts on two wires, not that the target score depends on only two variables.

3. Quantum tensor network: confine exponential dependence to the local window

QuanVI represents U and its conjugate transpose with an MPO-style chain of sliding-window tensor cores. D system wires correspond to data coordinates, and E ancilla wires provide environmental degrees of freedom for mixed states; every wire has local dimension K. Local cores have orthogonality or unitary constraints, and fixed, non-trainable boundary vectors close the contraction. Training must therefore preserve these constraints rather than modify a whole unconstrained matrix arbitrarily.

E=0 gives the pure-state variant; the main text states that the effective mixed-state rank is \(r=K^E\). This relation must be interpreted together with column orthonormality: r cannot exceed the system feature-space dimension \(K^D\), so adding ancilla wires is not an unlimited free source of effective rank. Increasing E changes both mixing freedom and the number of network cores, so improvements with E cannot be attributed exclusively to resolving degeneracy.

L is the number of wires acted on by each local core, the core count is \(D+E-L+1\), and entries per core grow as \(K^{2L}\). The paper gives an effective bond dimension growing as \(K^{L-1}\). Parameter growth in D and E is thus polynomial for fixed K and L, but increasing L to accommodate dense or long-range dependence brings exponential costs back.

Complex tensors expand the local parameterization; they are not equivalent to higher floating-point precision. The appendix's float-versus-cfloat comparison primarily changes real versus complex representation. Complex double precision in the posterior experiments is a separate experimental setting and should not be conflated with a general precision improvement.

4. Measurement contractions: one density operator supports multiple probability queries

After training, evaluating a coordinate value inserts the outer product of its basis vector; a full-density query inserts these measurement matrices on all system wires. Integrating a variable out replaces its measurement matrix with the identity. This is not approximate marginalization: the orthogonal-basis integral identity implements marginalization within the learned variational distribution.

A conditional density is obtained by dividing a joint-density contraction by a marginal-density contraction over observed variables, provided the denominator is nonzero. Queryability does not imply that the variational approximation equals the target; it means that integration and conditioning of the learned q share a computational interface.

Sampling applies the chain rule coordinate by coordinate: fix previously sampled coordinates, replace future coordinates with identities, obtain the current one-dimensional conditional density, and sample it with numerical inverse-CDF sampling. Both network contraction and numerical inversion require computation, so queryability does not imply constant-time generation of a high-dimensional sample.

A Worked Example

Take the main text's D=4, E=1, L=2 setting and assume K basis functions per coordinate, so the ancilla wire corresponds to effective rank K. To compute \(q(x_2,x_4)\), insert measurement matrices for the specified coordinates on system wires 2 and 4, insert identities on wires 1 and 3, and contract the ancilla wire and network cores.

To obtain the conditional density of coordinates 1 and 3 given \(x_2=z_2,x_4=z_4\), the numerator fixes all four coordinates, while the denominator fixes only coordinates 2 and 4 and integrates out coordinates 1 and 3. Their contraction ratio gives the conditional density without training an additional conditional model.

This also shows why orthogonal bases are part of the probability semantics rather than merely an optimization convenience: without their integral-to-identity property, replacing measurement matrices with identities would no longer automatically implement marginalization.

Loss & Training

Training minimizes the sum of traces of per-coordinate local score terms. The tildes below explicitly denote embedding each local term in the full space, avoiding dimensional ambiguity from the main text's suppressed embedding notation.

\[ \mathcal{L}(\Theta)=\sum_{d=1}^{D}\operatorname{tr}\!\left(\rho_{\Theta}\widetilde M^{(d)}\right),\qquad \rho_{\Theta}=\frac{1}{r}U_{\Theta}U_{\Theta}^{\dagger}. \]

Local terms are estimated from mini-batches, and the loss is evaluated by contracting U, a local operator, and conjugate U. The appendix mentions Stiefel projection and retraction to preserve orthogonality, but the full text does not provide complete optimizer pseudocode or all default hyperparameters, so no specific update formula is invented here.

Main synthetic experiments use K=5, E of 0 or 2, default L=2, 2000 training steps, and training sample count B=500D. Posterior experiments use K=5, E of 0 or 2, B=5000, 1000 training steps, learning rate 0.1, momentum 0.9, and complex double precision. They first apply Stan's unconstraining transformation with Jacobian-adjusted density evaluation, then precondition using a full covariance fitted by GSM or BaM. Reference posterior draws are not used to construct the preconditioner.

Complexity should follow the detailed accounting in Appendix B.1.1. The main text summarizes per-step cost as \(O(DBK^{2L})\), but the appendix includes direct contraction through the whole network for each local term and gives total training cost \(O(TD(B+D+E)K^{O(L)})\). A single-point density query costs \(O((D+E)K^{O(L)})\), and N complete sequential samples cost \(O(ND(D+E)K^{O(L)})\). These expressions still hide contraction-order and local-layout constants and do not separately account for inverse-CDF grid accuracy.

Key Experimental Results

Main Results

Synthetic experiments report forward KL, \(\mathrm{KL}(p\|q)\), over five runs, estimated with \(10^4\) target samples; lower is better. Non-Gaussian targets start with two-dimensional geometry and add a nonlinear Markov chain, so 100 dimensions does not mean arbitrary dense 100-dimensional posteriors were tested. The following selection from original Table 1 retains scalability evidence and failure cases.

Dimension / target	ADVI	EigenVI	MoG	QuanVI E0	QuanVI E2
5 / X-shape	1.6816 ± 0.0716	0.4007 ± 0.0743	1.3039 ± 0.0269	0.2613 ± 0.0035	0.2661 ± 0.0054
5 / Funnel	2.8012 ± 0.0653	1.3435 ± 0.1065	0.8183 ± 0.0742	4.3944 ± 1.8299	2.6412 ± 0.4579
20 / Ring	8.9786 ± 0.1421	OOM	10.6938 ± 0.3243	9.1621 ± 1.5400	9.0459 ± 1.0407
100 / Gaussian	0.1308 ± 0.0020	OOM	2.2290 ± 0.0119	≤ 0.01	≤ 0.01
100 / X-shape	47.2649 ± 0.2927	OOM	72.7630 ± 6.0097	42.2971 ± 1.2528	40.5473 ± 4.0462
100 / GMM3	47.1609 ± 0.2688	OOM	68.3951 ± 7.4535	43.0097 ± 2.3348	43.0798 ± 2.9883
100 / Funnel	48.8944 ± 0.2271	OOM	57.0316 ± 0.4351	46.4010 ± 3.1020	47.0005 ± 3.8384

Table 1 labels high-dimensional EigenVI as OOM, but Appendix C.2.2 says D>5 configurations are skipped and marked infeasible. This does not establish that every OOM cell was actually executed until memory exhaustion. Low-dimensional comparisons also use different basis sizes: synthetic EigenVI uses K=4, while QuanVI uses K=5.

Posterior experiments measure forward Fisher divergence with reference samples: target and variational scores are compared through squared error under p, with lower values indicating better approximation. It is not on the same numerical scale as either the extended training trace objective or the forward KL above. Original Table 2 reports four runs; the following preserves its means and standard deviations, with all metrics evaluated in common unconstrained coordinates.

Posterior / dimension	ADVI	EigenVI	QuanVI E0	QuanVI E2
gpregr / 3	1.1509 ± 0.0177	0.0828 ± 0.0681	0.2180 ± 0.0463	0.2330 ± 0.0369
hmm / 4	168.9795 ± 32.6976	4.0295 ± 1.5069	5.0783 ± 0.6214	4.8047 ± 0.4491
hmm_bball_0 / 6	746.7151 ± 2.0978	18.5057 ± 3.0322	13.7126 ± 0.5198	13.4460 ± 0.3597

Ablation Study

Original Table 5 evaluates local windows over five runs using forward KL. Three representative rows retain the complete window sequence below. These ablation runs differ from the main comparison, so their L=2 values are not replaced with main-table values.

Dimension / target	L=1	L=2	L=3	L=4	L=5
5 / X-shape	6.7725 ± 0.0399	0.2645 ± 0.0065	0.2478 ± 0.0050	0.2552 ± 0.0083	0.2460 ± 0.0120
10 / Funnel	18.9833 ± 3.1152	7.9293 ± 1.1948	6.3671 ± 0.4352	4.3833 ± 0.5269	2.8119 ± 0.0503
20 / Funnel	38.5510 ± 5.3004	13.7661 ± 0.5832	11.7725 ± 0.1846	7.0066 ± 0.2342	—

Table 6 shows the steep window cost on Funnel: at D=10, L=2 takes about 15 minutes and L=3 exceeds 12 hours; D=100 with L=2 takes about 5 hours. Dashes in the window table indicate only that no final KL was available, not zero error, OOM, or numerical divergence.

Figure 3's ancilla-wire ablation uses a noisier D=10 X-shape setting: 500 steps, a freshly sampled B=100 batch at each step, and five seeds per E. The text says E=3 has the best mean result, but the cache does not contain the figure's full numeric values, so no E-sweep table is invented. The real/complex ablation on D=5, E=0 Funnel gives mean KL of about 5.08 / 3.58 and runtimes of about 1 minute 10 seconds / 1 minute 30 seconds; these are not universal improvement magnitudes across targets.

Key Findings

• The main gain is avoiding the storage bottleneck of a global eigenvalue problem, not near-zero high-dimensional error: KL on 100-dimensional non-Gaussian targets remains around 40–47.
• Mixed states do not always outperform pure states; E0 has lower means on 100-dimensional GMM3 and Funnel. The experiments support occasional fitting improvements from moderate ancilla freedom, but do not measure degeneracy directly or isolate its stability benefit.
• Funnel benefits more from larger windows, yet costs can jump from minutes to more than a dozen hours. Polynomial dependence on D does not make every capacity setting inexpensive.
• NF should not be dismissed using only failure cases: Appendix Table 3 reports median KL of 0.0319 on 5-dimensional X-shape, lower than QuanVI's main-table means, but with a different summary statistic. Numerous divergent runs in other settings preclude an unconditional SOTA claim.

Highlights & Insights

• Replacing vector selection with subspace representation encodes invariance to basis changes under degeneracy directly in the solution. This is more fundamental than fixing an eigenvector's sign after optimization, but requires explicit rank and spectral-gap conditions.
• Locality of the objective and locality of the representation work together: one reduces operator support, the other reduces learnable parameters. They are complementary, and neither substitutes for the other's assumptions.
• Basis orthogonality connects normalization with marginalization: the same identity-replacement rule supports multiple probability queries. A transferable principle is to design representations with useful integral identities before building tractable inference interfaces.

Limitations & Future Work

• The authors explicitly acknowledge dependence on local structure and contraction costs for larger windows, proposing richer dependency graphs and adaptive tensor-network architectures.
• High-dimensional evidence primarily comes from chain-constructed synthetic distributions; real posteriors have only 3, 4, and 6 dimensions, leaving performance on dense high-dimensional Bayesian posteriors unestablished.
• Ancilla ablations increase rank and core count together, so they cannot identify the causal benefit of maximally mixed states for degeneracy alone. More targeted tests should control parameter count, known spectral degeneracy, and perturbation strength.
• Uniform weights within a fixed subspace restrict the mixture spectrum. Adaptive rank or nonuniform weights are promising directions but would change the current maximally mixed-state semantics and optimization problem.
• The main text's simplified training complexity differs from the appendix's detailed contraction accounting; NF uses medians while the main table uses means, and sample budgets differ across methods. Efficiency comparisons need matched wall-clock or score-evaluation budgets and stability statistics.

Related Work & Insights

• vs EigenVI: Both use orthogonal function expansions and score-induced operators; EigenVI explicitly solves for one lowest eigenvector, whereas QuanVI learns a mixed subspace constrained by a local QTN. Avoiding the global matrix enables scalability at the cost of nonconvex training and finite window capacity.
• vs GSM / BaM: Gaussian families are easier to optimize but cannot express all non-Gaussian geometry; QuanVI uses richer basis features and tensor correlations. Gaussian preconditioning in the posterior experiments demonstrates complementarity rather than complete replacement.
• vs Born-type tensor networks: Pure-state Born representations turn one squared amplitude into probability; QuanVI uniformly mixes several orthogonal amplitudes. The broader lesson is to model symmetries of the represented object, not merely compress its parameters.

Rating

• Novelty: 4/5 — Mixed subspaces and local QTNs have a clear joint motivation, but equivalence between mixed-density Fisher divergence and the training objective needs clarification.
• Experimental Thoroughness: 3/5 — Synthetic targets, posteriors, and several ablations are covered, but real high-dimensional posteriors and direct degeneracy-robustness tests are missing.
• Writing Quality: 3/5 — The main mechanism is clear, while complexity accounting, local-operator embedding, and infeasible-baseline labels need clarification.
• Value: 4/5 — Useful for structured score-based inference and queryable probability representations, provided local dependence and window budgets are checked first.