{"id":"27ff410a-ded2-413a-80eb-15fe02c328ac","arxiv_id":"2501.04935","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":3,"one_line_summary":"A joint Kronecker-structured Inverse-Wishart variational approximation, optimized with an orthogonalized pullback metric, outperforms mean-field variational Bayes for Gaussian multiway covariances.","lead":"The authors develop a Riemannian geometry for variational Bayesian approximations to the covariance of Gaussian multiway data, using a pullback of the affine-invariant metric on Kronecker-structured covariance matrices. They show a joint Kronecker-structured approximation converges faster and fits the exact Inverse-Wishart posterior better than a mean-field approximation, and they derive efficient gradients and a Monte Carlo sampler.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Convexity claim in §2.1 is unsupported and the empirical speed comparison does not control for step-size tuning; the main risk is that the reported joint-vs-mean-field advantage is optimization artifact, not method advantage.","rationale":"I agree with the reader's weak assumption: §2.1's unsupported leap from Wiesel's convexity to full ELBO convexity is the most load-bearing gap, because both the algorithm's justification and the speed comparison depend on it. The reader also flagged the Supplement/main-text contradiction and the hand-selected step sizes, which I fold into the concrete test. The geometric derivations (Propositions 2.1–2.3) appear internally consistent, and the construction of the orthogonalized pullback metric is a genuine contribution. The fidelity claim for the joint approximation is qualitatively supported, but without a convexity guarantee or controlled step-size comparison, the efficiency claim is not established. These are addressable, so CONDITIONAL is appropriate; no basis for outright rejection.","tokens_in":21541,"tokens_out":1579,"duration_ms":13990,"concrete_test":"Derive and verify the Hessian of the joint ELBO from §3.1 along a geodesic of Proposition 2.3 on a small non-diagonal example, e.g., D=2, d1=d2=2, checking convexity numerically over a grid of endpoints. If any negative second derivative is found, the convexity claim fails. Independently, rerun Figures 1–3 with both methods using their respective largest stable step sizes (for example, ϵ_joint=10^{-4} and ϵ_MF=10^{-5.5}, matching the stability limits already observed), and report the iteration counts and wall-clock times with error bars over multiple random seeds; if mean-field then converges in comparable iterations, the efficiency claim is an artifact of step-size choice.","verdict_should_be":"CONDITIONAL","load_bearing_attack":"The paper's central claim—that the joint Kronecker variational family outperforms mean-field in both optimization efficiency and fidelity—rests on two load-bearing pillars. First, §2.1 asserts without proof that 'all ELBO computations considered later in this paper are geodesically convex,' citing Wiesel [26] for the multivariate normal likelihood. That reference covers maximum likelihood on Kronecker-structured SPD matrices, not the variational ELBO, which adds entropy terms, inverse-Wishart prior terms, and the trace term tr(⊗_i A_i^{-1}(S+Λ)). No theorem establishes geodesic convexity on the orthogonalized product manifold under the pullback metric, nor for the mean-field ELBO under the product metric. If the ELBO is not geodesically convex, Algorithm 2's iterates could stall at local optima or saddle points, and the convergence differences in Figures 1–3 would not be attributable to the geometry or family. Second, the empirical speed comparison is not controlled: the joint model used one global step size ϵ=10^{-4.4} while mean-field used ϵ=10^{-6}, and the text admits mean-field was 'not necessarily run until convergence.' A one-parameter step-size choice can dominate iteration counts in Riemannian gradient descent, so the claimed 'superiority' in optimization efficiency is not demonstrated. The fidelity comparison (Figure 2) is more defensible, but its caption concedes the mean-field penalty may be 'unsurprising.' The Supplement's LBMF derivation also contradicts the main-text mean-field ELBO (main text has νvj/2 S; Supplement omits it), so the mean-field objective actually optimized in the experiments is ambiguous. These issues are addressable but currently leave the headline claim under-supported.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper develops a Riemannian optimization framework for variational Bayes on Kronecker-structured covariance matrices arising from tensor-normal (multiway) data. It constructs the pullback of the affine-invariant metric under the Kronecker product map, proves that the naive pullback metric is degenerate, and proposes an orthogonalization condition (unit determinants on all but one factor) that renders the pullback metric positive definite and block-diagonal. On this basis, the paper derives Euclidean and Riemannian gradients for ELBOs of a joint Inverse-Wishart variational family with Kronecker scale and of a mean-field family, and it reports simulation experiments comparing the two families in terms of optimization iterations and posterior predictive Mahalanobis distance, plus an application to UN commodity trade data.","tokens_in":21830,"tokens_out":13755,"duration_ms":115588,"significance":"If the geometric and variational claims survive scrutiny, the paper would make a useful contribution to Bayesian inference for multiway covariance models: the pullback-metric construction with its explicit degeneracy resolution is a natural and nontrivial geometric result, and the efficient computation of the trace term in Proposition 1.3 is practically valuable. The paper also ships a concrete algorithm with gradient derivations, which strengthens its reproducibility potential. However, the current manuscript contains load-bearing weaknesses in the convexity justification, the mean-field ELBO derivation, and the experimental comparisons, so the significance can only be assessed after these are repaired.","major_comments":[{"comment":"The assertion that 'all ELBO computations considered later in this paper are geodesically convex' is not supported. Wiesel [26] establishes geodesic convexity for maximum likelihood estimation of Kronecker-structured covariance matrices, not for a variational ELBO that includes an entropy term, an Inverse-Wishart prior term, and a trace term, all under a pullback metric on an orthogonalized product manifold. No theorem or proof is given for the joint or mean-field ELBOs. This is load-bearing because Algorithm 2's convergence and the speed comparisons in Section 4 rely on the optimization being well-behaved; without a convexity proof, iterates could be converging to local optima, and the reported gap between joint and mean-field would not measure what the paper claims. Please either prove geodesic convexity for the specific ELBOs or substantially weaken the claim and discuss the risk of local optima.","section":"Section 2.1"},{"comment":"The supplement's derivation of the mean-field ELBO is internally inconsistent with Section 3.2. The final displayed equation in Supplement Section 6.2 is exactly the joint ELBO from Supplement Section 6.1, with a single degree-of-freedom parameter ν_v and the same trace term tr([⊗_i A_i^{-1}](S+Λ)). In contrast, the main-text mean-field bound in Section 3.2 uses separate ν_{vi} and a product ∏_j (ν_{vj}/2) in the trace term. Because the mean-field gradients in Section 3.2 and all mean-field experiments derive from the mean-field ELBO, this duplication means the mean-field results lack a correct derivation. Please fix the supplement to derive the factored mean-field bound consistently with Section 3.2.","section":"Supplement, Section 6.2"},{"comment":"The optimization-efficiency claim is not demonstrated because the comparison is not controlled. The joint approximation uses a global step size ε_joint = 10^{-4.4} while the mean-field uses ε_MF = 10^{-6}, and the text states that the mean-field was 'not necessarily run until convergence.' A single step-size choice can dominate iteration counts in Riemannian gradient descent, so the reported gap may be an optimization artifact rather than a property of the variational family. Please report wall-clock times, run both methods to a common convergence criterion, and perform a step-size sensitivity analysis across a grid for both methods, rather than selecting the 'most successful' step per method post hoc.","section":"Section 4.1, Figures 1-3 and Table 1"},{"comment":"The fidelity comparison is not apples-to-apples. The joint approximation targets the unstructured Inverse-Wishart posterior, while the mean-field baseline approximates a different Bayesian model with independent priors (Section 3.2, based on Hoff [12]). Therefore the abstract's claim that the joint approximation provides 'a superior approximation to an unstructured Inverse-Wishart posterior' is not established: the mean-field family is not attempting to approximate that posterior. The caption itself concedes the mean-field penalty may be 'unsurprising.' A fair comparison would use a mean-field family targeting the same unstructured posterior (for example, a factorized approximation over the modes or entries of the covariance), or the claim should be restated as a comparison of two different model/posterior choices rather than a statement about approximation fidelity to a fixed target.","section":"Section 4.1, Figure 2 and Abstract"}],"minor_comments":[{"comment":"The proof's Schur-complement computation appears incorrect in the diagonal blocks: g_ii - g_{i1} g_{11}^{-1} g_{1i} is not generally zero because the outer product v(Σ_i^{-1})v(Σ_i^{-1})^T differs from Σ_i^{-1}⊗Σ_i^{-1}. The degeneracy nevertheless holds via the scaling-invariance direction (V_i = c_i Σ_i with ∑ c_i = 0), so the proposition's conclusion is correct, but the proof should be revised.","section":"Section 2.2, Proposition 2.2"},{"comment":"The Euclidean gradient formula for ∇_{A_i} LB includes a term '-(n+ν)/2 d_{-1} A_1^{-1}' for a general index i. Under the orthogonalization condition, the prior log-determinant term contributes only to A_1; for i>1 the term should be absent (or, in a non-orthogonalized formulation, it would be d_{-i} A_i^{-1}). As written, the formula is dimensionally inconsistent for i>1.","section":"Section 3.1"},{"comment":"In the sentence 'we would instead have ∑_{j=1}^D d_{-j} log|A_j| = d_{-1} A_1', the right-hand side should be 'd_{-1} log|A_1|'; the log is missing.","section":"Section 3.1"},{"comment":"The gradient for ν_{vi} contains the undefined symbol 'C_p(S+Λ)'; please define this quantity or correct the notation.","section":"Section 3.2"},{"comment":"The paper consistently misspells 'ELBO' as 'EBLO' in the supplement (e.g., 'non-orthogonalized EBLO'); please correct throughout.","section":"Supplement, Sections 6.1 and 6.2"},{"comment":"The definition of T^{(k)}(S, Σ_{-i}) is difficult to parse; the formula for C_k(S, Σ_{-i}) has a missing closing parenthesis and the indexing of S_Γ is not immediately transparent. A small worked example would improve readability.","section":"Section 1.3, Proposition 1.3"},{"comment":"The data-generation details are incomplete: no random seeds, prior hyperparameter values for Λ, or exact model specification are given for the simulations, which hinders reproducibility. Please provide these details and, preferably, the code used for the experiments.","section":"Section 4.1"},{"comment":"The trade-data analysis claims an 'identical eigenvalue structure' and 'similar clustering pattern' compared to Hoff [12], but no quantitative comparison (e.g., Procrustes alignment or eigenvector correlation) is provided; the claim is currently based on visual inspection.","section":"Section 4.2"}],"recommendation":"major_revision","confidential_remarks":"The paper's core geometric contribution — the pullback metric and its degeneracy resolution — appears sound and is likely of interest to the statistical computing audience. However, the manuscript in its current form has a serious internal inconsistency in the supplement's mean-field derivation, an unsupported geodesic-convexity assertion, and uncontrolled empirical comparisons. These issues are fixable within the scope of a major revision, but they need to be addressed carefully before the paper can be considered for publication. I would also encourage the authors to publish the code and data-generation details, as the empirical claims are central to the paper's impact."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Dear colleague,\n\nThe thing to know about this paper: the geometry contribution is real and the empirical headline is not yet supported. The core derivation—that the naive pullback of the affine-invariant metric under the Kronecker product is degenerate, and that fixing |Σ_i|=1 for i>1 makes it block-diagonal and positive definite—is coherent and, as far as I can tell, correct. That is worth having.\n\nWhat else is new: the joint Inverse-Wishart variational ELBO with Kronecker scale, the corresponding Euclidean and Riemannian gradients, and the efficient trace computations via Proposition 1.3. The trade data analysis is a nice demonstration of the intended use.\n\nThe soft spots are real. Section 2.1 asserts without proof that all ELBO computations are geodesically convex, citing Wiesel's result for the multivariate normal likelihood. That reference covers maximum likelihood, not the variational objective with entropy and prior terms. I could not find a theorem in the paper that establishes geodesic convexity on the orthogonalized product manifold. It may be true, but it is unproven, and Algorithm 2's convergence claims depend on it.\n\nThe empirical speed comparison is also not clean. The joint model used one global step size, the mean-field another, and the text admits the mean-field was not necessarily run until convergence. With Riemannian gradient descent, a one-parameter step-size choice can dominate iteration counts. So the \"optimization efficiency\" claim is not demonstrated. The fidelity comparison is more defensible, but it compares a joint approximation to an unstructured prior against a mean-field approximation of a different model with independent priors. Some of the gap is by construction.\n\nThe Supplement has a concrete error: the mean-field ELBO in Section 6.2 is identical to the joint ELBO, contradicting the main text's mean-field ELBO, which includes the νvj/2 S term. That needs fixing. Also no code or data is provided, and there are no error bars in the figures.\n\nWho is this for? Computational Bayesian statisticians working on tensor normal models and anyone doing optimization on Kronecker-structured SPD manifolds. The geometry section is worth reading now; the empirical claims need revision before they can be relied on.\n\nI would send this to peer review: the derivations deserve scrutiny by a knowledgeable referee, and the issues are addressable. A serious referee should demand the convexity proof (or a careful statement of where it fails), a corrected Supplement, a fairer empirical protocol, and code/data release.","headline":"Real geometry result, unsupported empirical headline; deserves review but not acceptance as is.","tokens_in":22409,"tokens_out":2358,"would_cite":true,"duration_ms":22454,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["62F15","62H12","53B21","90C26"],"pacs":[],"model":"deepseek-v4-flash","headline":"A joint Kronecker-structured variational family, optimized on a pullback metric made positive definite by determinant normalization, converges faster and tracks the unstructured Inverse-Wishart posterior better than mean-field variational…","keywords":["affine-invariant metric","pullback metric","Kronecker covariance","tensor normal model","variational Bayes","Inverse Wishart","mean-field approximation","Riemannian gradient"],"falsifier":"Take a small array, say $D=3$ with dimensions $(5,6,4)$, choose two feasible points in the orthogonalized product manifold, evaluate the ELBO along the affine-invariant geodesic connecting them, and test the convexity inequality $f(\\gamma(t)) \\leq (1-t)f(p)+t f(q)$. A single violation would refute the geodesic-convexity premise; equivalently, running Algorithm 2 from many random starts and seeing divergent optima would indicate the reported convergence is not to a global optimum.","tokens_in":21316,"feed_emoji":"📐","tokens_out":7713,"duration_ms":68550,"temperature":0.7,"pith_summary":"This paper tries to prove that variational Bayesian inference on the covariance of Gaussian multiway data is best done jointly on the Kronecker-structured space rather than with the usual mean-field factorization. The authors construct a Riemannian geometry on the product of mode-wise covariance manifolds by pulling back the affine-invariant metric from the full Kronecker space, and show the naive pullback degenerates until one imposes the normalization $|\\Sigma_i|=1$ on all but one mode. With that normalization, the pullback metric becomes a positive-definite direct sum, so each mode can be moved along its own geodesic during optimization. They then argue that a joint Inverse-Wishart variational family with a Kronecker scale matrix, optimized on this geometry, converges in fewer iterations than the mean-field approximation and gives a closer match to the unstructured Inverse-Wishart posterior under average Mahalanobis distance, while preserving the multiway interpretation of the covariance.","feed_headline":"Joint Kronecker VB beats mean-field on multiway covariances","feed_subtitle":"A determinant-normalized pullback metric turns separable covariance VB into fast geodesic optimization.","key_machinery":"The central object is the pullback of the affine-invariant metric under the Kronecker product map, together with a determinant projection. The affine-invariant metric $g_{\\Sigma}(U,V)=\\mathrm{tr}(\\Sigma^{-1}U\\Sigma^{-1}V)$ is the Hessian of the Boltzmann entropy and connects to Fisher information. Its pullback to the product manifold is degenerate; Proposition 2.3 repairs this by requiring $|\\Sigma_i|=1$ for $i>1$, which zeroes out cross terms and yields the direct-sum metric $G(\\Sigma_i|\\Sigma_{-i})=(1/d_{-i})\\Sigma_i^{-1}\\otimes\\Sigma_i^{-1}$. Optimization then uses geodesic updates $\\Sigma_i(t)=\\Sigma_i^{1/2}\\exp(t\\Sigma_i^{-1/2}V_i\\Sigma_i^{-1/2})\\Sigma_i^{1/2}$, with tangent vectors projected by $P_{\\Sigma_i}(V)=V-\\mathrm{tr}(V\\Sigma_i^{-1})/d_i\\,\\Sigma_i$. The trace identity of Proposition 1.3 supplies the efficient gradient terms $T^{(k)}(S,\\Sigma_{-k})$ that avoid sample-size-scaled tensor contractions.","core_discovery":"On its own terms, the paper's discovery is that the degeneracy of the pullback metric has a simple fix, and the fix is what makes joint variational inference practical. The map $\\phi: \\times_{i=1}^D P(d_i) \\to P(\\otimes_{i=1}^D d_i)$ sending $(\\Sigma_1,\\ldots,\\Sigma_D)$ to $\\otimes_{i=1}^D \\Sigma_i$ induces a metric whose off-diagonal blocks wipe out the determinant: under the affine-invariant metric, the pullback $\\Phi^*G$ is not positive definite. Proposition 2.3 gives the resolution: impose $|\\Sigma_i|=1$ for $i>1$, so the tangent space condition $\\mathrm{tr}(d\\Sigma_i \\Sigma_i^{-1})=0$ kills all cross terms and leaves $\\Phi^*G = \\oplus_{i=1}^D (1/d_{-i}) \\Sigma_i^{-1}\\otimes\\Sigma_i^{-1}$, which is positive definite. On this orthogonalized manifold, the Riemannian gradients of the joint ELBO reduce to simple expressions, the geodesic updates factor mode by mode, and the joint variational family $\\mathrm{IW}(\\nu, \\otimes_{i=1}^D A_i)$ has a single variational degree of freedom $\\nu$. The empirical claim, demonstrated in simulations, is that this joint approximation is both faster and more faithful to the unstructured posterior than the mean-field approximation of the independently-prior model, and that it reproduces the mode-wise eigenstructure of a full sampler on commodity-trade data.","pith_inferences":["If the ELBO's geodesic convexity is confirmed, Algorithm 2 would come with a global-optimality guarantee, making the comparison with mean-field a fair test rather than a race to different local optima.","The determinant normalization singles out mode 1 as the carrier of overall scale; an equivariant counterpart that spreads the normalization across modes could make the geometry invariant to how the tensor is indexed.","Because mean-field samples are always perfectly Kronecker-separable while joint Inverse-Wishart samples are not, the joint approximation should also be preferred when the goal is uncertainty about nonseparable structure, not just the separable mean.","A natural extension is to test the two approximations by posterior predictive coverage on held-out tensor entries under misspecification, since average Mahalanobis distance on separable truth is only one scoring rule."],"forward_implications":["In the truly separable case, a single joint Inverse-Wishart family with Kronecker scale replaces $D$ independent families, so variational inference has one degrees-of-freedom parameter instead of $D$ and converges in fewer iterations.","The orthogonalized pullback metric gives positive-definite Riemannian gradients and geodesic updates that factor mode by mode, making the separable covariance variational problem a well-defined geometric optimization.","The trace identity lets the ELBO and its gradients be computed in a way that does not scale with sample size, removing a computational bottleneck of Gibbs samplers for tensor data.","On commodity-trade data, the joint approximation recovers the same exporter, importer, and commodity eigenstructure as a full sampler, so the approximate posterior remains interpretable.","Under misspecification with low-rank noise added to a separable covariance, the joint approximation stays stable while the mean-field approximation needs dramatically more iterations or fails to converge."],"supporting_citations":[{"why":"Defines the tensor normal model, the independent-prior Bayesian model used as the mean-field baseline, and the trade-data application.","marker":"[12]"},{"why":"Supplies the geodesic-convexity result for Kronecker-structured covariance likelihoods that the paper relies on for ELBO convexity.","marker":"[26]"},{"why":"Provides the closed-form affine-invariant geodesic used in the position updates of Algorithm 2.","marker":"[18]"},{"why":"Derives the Riemannian metric as the Hessian of the Boltzmann entropy and the associated gradient formula for SPD matrices.","marker":"[19]"},{"why":"Provides the Bartlett decomposition used for efficient Monte Carlo sampling of the multiway Inverse-Wishart.","marker":"[15]"},{"why":"Gives the lower bound on the error of the best Kronecker approximation, used to argue the mean-field approximation's limitation.","marker":"[22]"}],"fun_headline_variants":["Geodesic VB on Kronecker covariances beats mean-field","Joint Kronecker VB outperforms mean-field via geodesic metric","Affine-invariant metric speeds up variational multiway covariances","Geometric fix enables efficient joint multiway covariance VB","New VB geometry makes joint multiway covariance inference faster"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The load-bearing premise is that the full ELBO, including the variational entropy and trace terms, is geodesically convex on the orthogonalized product manifold; the paper asserts this from a cited convexity result for the Gaussian likelihood but does not prove it for the whole objective.","fun_headline_variants_meta":{"raw":{"variants":["Geodesic VB on Kronecker covariances beats mean-field","Joint Kronecker VB outperforms mean-field via geodesic metric","Affine-invariant metric speeds up variational multiway covariances","Geometric fix enables efficient joint multiway covariance VB","New VB geometry makes joint multiway covariance inference faster"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000773,"raw_usage":{"total_tokens":3437,"prompt_tokens":973,"completion_tokens":2464,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":589,"completion_tokens_details":{"reasoning_tokens":2378}},"tokens_in":589,"tokens_out":2464,"duration_ms":18310,"temperature":1.0,"reasoning_tokens":2378,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-10T21:22:42.023541+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Take a small array, say $D=3$ with dimensions $(5,6,4)$, choose two feasible points in the orthogonalized product manifold, evaluate the ELBO along the affine-invariant geodesic connecting them, and test the convexity inequality $f(\\gamma(t)) \\leq (1-t)f(p)+t f(q)$. A single violation would refute the geodesic-convexity premise; equivalently, running Algorithm 2 from many random starts and seeing divergent optima would indicate the reported convergence is not to a global optimum.","supporting_citations":[{"cited_title":"Separable covariance arrays via the Tucker product, with applications to mul- tivariate relational data","cited_arxiv_id":null,"evidence_quote":"Defines the tensor normal model, the independent-prior Bayesian model used as the mean-field baseline, and the trade-data application."},{"cited_title":"Geodesic convexity and covariance estimation","cited_arxiv_id":null,"evidence_quote":"Supplies the geodesic-convexity result for Kronecker-structured covariance likelihoods that the paper relies on for ELBO convexity."},{"cited_title":"Symmetric positive-definite matrices: From geome- try to applications and visualization","cited_arxiv_id":null,"evidence_quote":"Provides the closed-form affine-invariant geodesic used in the position updates of Algorithm 2."},{"cited_title":"The Riemannian geometry of the space of positive-definite matrices and its application to the regularization of positive-definite matrix-valued data","cited_arxiv_id":null,"evidence_quote":"Derives the Riemannian metric as the Hessian of the Boltzmann entropy and the associated gradient formula for SPD matrices."},{"cited_title":"Bartlett decomposition and Wishart distribution","cited_arxiv_id":null,"evidence_quote":"Provides the Bartlett decomposition used for efficient Monte Carlo sampling of the multiway Inverse-Wishart."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Gives the lower bound on the error of the best Kronecker approximation, used to argue the mean-field approximation's limitation."}],"review_version":1}