HoT-VI: Reparameterizable Variational Inference for Capturing Instance-Level High-Order Correlations
Authors: Junxi Xiao
Qinliang Su
Zexin Yuan
Venue: NeurIPS 2025 poster ()
Keywords: Probabilistic Models
Variational Inference
Instance-level Correlation
Reparameterization
High-order Correlation
Primary area: probabilistic_methods

Abstract: Mean-field variational inference (VI), despite its scalability, is limited by the independence assumption, making it unsuitable for scenarios with correlated data instances. Existing structured VI methods either focus on correlations among latent dimensions which lack scalability for modeling instance-level correlations, or are restricted to simple first-order dependencies, limiting their expressiveness. In this paper, we propose High-order Tree-structured Variational Inference (HoT-VI), that explicitly models $k$-order instance-level correlations among latent variables. By expressing the global posterior through overlapping $k$-dimensional local marginals, our method enables efficient parameterized sampling via a sequential procedure. To ensure the validity of these marginals, we introduce a conditional correlation parameterization method that guarantees positive definiteness of their correlation matrices. We further extend our method with a tree-structured backbone to capture more flexible dependency patterns. Extensive experiments on time-series and graph-structured datasets demonstrate that modeling higher-order correlations leads to significantly improved posterior approximations and better performance across various downstream tasks.

HoT-VI: Reparameterizable Variational Inference for Capturing Instance-Level High-Order Correlations
Junxi Xiao
Qinliang Su Zexin Yuan
School of Computer Science and Engineering, Sun Yat-sen University, Guangzhou, China Guangdong Key Laboratory of Big Data Analysis and Processing, Guangzhou, China {xiaojx7, yuanzz)))maillll sysu. edu. sugliang@mail. sysu. edu.
Abstract
Mean-field variational inference (VI), despite its scalability, is limited by the independence assumption, making it unsuitable for scenarios with correlated data instances. Existing structured VI methods either focus on correlations among latent dimensions which lack scalability for modeling instance-level correlations, restricted to simple first-order dependencies, limiting their expressiveness. In this paper, we propose High-order Tree-structured Variational Inference (HoT-VI that explicitly models k-order instance-level correlations among latent variables. By expressing the global posterior through overlapping k-dimensional local marginals, our method enables efficient parameterized sampling via a sequential procedure. To ensure the validity of these marginals, we introduce conditional correlation parameterization method that guarantees positive definiteness of their correlation matrices. We further extend our method with rree-structured backbone capture more flexible dependency patterns. Extensive experiments on time-series and graph- structured datasets demonstrate that modeling higher-order correlations leads significantly improved posterior approximations and better performance across various downstream tasks.
Introduction
Variational inference (VT) is widely used framework for approximating the posterior distribu- tion in latent-variable models (X, Z) pe pz))))))(Z) where X = [x1,X2,, X2, .xv, and |21, .2...... zN] are observed data and the corresponding latent variables, respectively. VI seeks to approximate the intractable posterior p( (Z with a tractable surrogate distribu- tion qs (Z from parametrized distribution family 2Φ. by maximizing the evidence lower (X, Z) log 90 x)] In typical settings where data instances are assumed to be independent, the joint distribution naturally factorizes across instances as (X, Z) II. (x, z,) where and z, denote the i-th data instance and latent variable. Under these scenarios, we can reasonably use the mean-field posterior q (Z[X) = ПI qo (z, |x;) for model inference and training. However, many real-world datasets exhibit complex relationships among instances, making the independence assumption across data instances untenable. For instance, in mul- tivariate time series [59 211. the latent states z; timestamp may depend on those preceding and succeeding timestamps. Similarly in graph-structured data such as social networks [19] and citation graphs [29], the latent representations and zj of connected nodes are typically correlated due to underlying relational structure. In these scenarios, the mean-field posterior clearly inadequate, as ignores crucial dependencies among data points.
bound C(O, Φ) Ezzag [log Pe
"Corresponding author.
"oode available at: http://ghhhlt.oom/mglhsssphhohhh
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Many existing variational inference methods have attempted to incorporate structured dependencies into posterior approximation [68] 3 3T. but most of these approaches focus on modeling correlations among dimensions within a latent variable. Since the number of latent dimensions is typically small (e.g., dozens to hundreds) compared to the size of datasets (e.g., thousands to millions), methods tar- geting dimension-level correlations cannot be applied to capture instance-level correlations, especially for large datasets. Some recent efforts have been devoted to explicitly model instance-level correla- tions. For example, Correlated Variational Autoencoder (CVAE) (56) introduces tree-structured variational distribution that captures pairwise dependencies between neighboring latent variables. Similarly, Tree- stresstruc Variational Inference VI) [63] builds correlation matrix derived from a tree structure over instances, enabling efficient sampling and scalable inference. However, the reliance on capturing pairwise interactions limits these methods to first-order correlations, and precludes the representation of eyclic or higher-order dependency structure among latent variables. Yet in many real-world domains such financial time series [51], Sensor networks [12,, climate data 40] and evolving graphs (24), correlations among data instances are not merely pairwise. These data frequently exhibit high-order dependencies, where the relationship between two latent variables is mediated by the joint influence of multiple others. In time series, for instance, the latent state at a given time point may not only depend on nearby time steps, but also on patterns that occurred further in the past [77. In such cases, methods constrained to pairwise or tree-structured dependencies are fundamentally limited in expressiveness, necessitating a more expressive variational inference framework capable of capturing higher-ordet instance-leve dependencies.
In this work, introduce novel variational inference framework that overcomes the limitations of existing methods by incorporating higher-order dependency structures among latent variables. Rather than restricting attention first-order correlations, our approach is able to capture more expressive k- order dependencies. We theoretically show that, by imposing a k-order dependency structure into the global variational posterior, the high-dimensional global posterior can be expressed in terms of of k-dimensional local marginal distributions. By leveraging these local marginals, a sequential sampling method is developed draw parameterized samples from the high-dimensional global variational posterior, which are then substituted into the evidence lower bound for training. To ensure the validity of this approach, conditional correlations are introduced to re- parameterize the correlation matrices of local marginals. We prove that by using the conditional correlations to represent the correlation matrices, the positive definiteness of correlation matrices are guaranteed, and so does the validity of the developed VI approach. Later, we further show that the approach can be extended to support more structured dependency patterns by generalizing to tree-structured backbones, enabling even richer representations of latent correlations. Extensive experiments on time-series and graph-structured datasets demonstrate that our proposed method outperforms competitive baselines by effectively capturing higher-order correlations among latent variables, leading to improved performance in downstream tasks.
2
Variational Inference with High-order Correlation
2.1 Variational Posterior with Instance-Level Correlation Structure
To have the paper focus on its primary objective of capturing instance-level cor- relations, we assume dimension-leve level in- dependence in the variational posterior by restricting it to the factorized form 90 (Z) I2 II & ([Z)..1.N XX). where [II,1:N'' denotes the d-th row of the latent-variable matrix Z E DxN, Figure 1: N instances with k-order dependency structure.
Ziz
Zikk
Zikty
Zx
although the method can be easily ex- tended to take the dimension-leve correlation into account. To model instance-level correlations, we set the d-th variational posterior g (ZZ),, X) be a correlated Gaussian distribution as 9 [Z),, XX= = . (V &, Maa ). where μa E N and Pa E R.N denote the mean vec- tor and precision matrix, respectively. By noting that dimensions are handled separately and similarly, in the following, for the conciseness of presentation, we omit the subscript d and observed data X, and
simply denote the dl-h dimensional variational posterior 4o 1 (X) as g(z) N2( μ, P-1).
where
[21, 22, , NN]T a [M1, ,PN]" and P E R
2

---

Since the instance number N is often very large, which could be as large as tens of thousands or even millions in many scenarios, if we simply set the precision matrix P as a general matrix, would be computationally intractable. To balance the computational cost and the capability of modeling high-order correlation, we propose to impose a k-order connection structure into the precision matrix P. as shown in Fig. [ which is equivalent to set the j)th) element of P to zero for any ji-j] >> k, that is. Púj = 0 for any i - 5| > k. Note that we here assume the connection structure of latent variables is built upon chain, which is reasonable for the modeling of sequential data by itself. Later, we will show that the chain backbone is not necessary and can be extended to the more general tree topology to accommodate more diverse data.
Despite the o-order connection structure is imposed into if we simply substitute the variational posterior (z) = (z; p., p=) into the lower bound, we will see that the N latent variables are still coupled together, and we still have to handle them simultaneously. To overcome this issue, we notice that the connection structure is comprised of N k overlapping local sub-structures, with each involving only k + consecutive latent variables Z4z++k [2,, 2,+1,"., = for i i= 1,2, .... .N- - k. The k + consecutive variables Z follow the marginal distribution of g(z), which can be expressed as
(Z(s++)) = (Zzstti Pos+k di axtk e diag(o
(1)
where [R] at 1i+t-1,i4t- with [sst denoting the (s, ()-th element of matrix; N++ 1,its is
the correlation coefficient between latent variables 2i+s-1 and Ei+t-1 1 for any 8, E (1.2, k k+11;
and k [o,.... : -.+k), with o; being the standard deviation of zj. Each of these local marginals ((Z,+4)) represents a localized view of the global variational posterior q(z). Below, we show that the global posterior q(z) can be expressed in terms of these local marginals q ((m.i-)).
Theorem 2.1. For any joint distribution q(z) (Z; μ, P-1) with a precision matrix P that has k-order connection structure, it can be equivalently expressed as
N-K+1
N-
((zi44))
q(z)
II
g(Za++k-))
II
g (1 (2((+)12))))
where q(Zei+k-1) and q( ((zzit+k) are the marginals of g(z) over wistk-1 and 2is k-- Moreover, if the (k + 1)-variate marginals ((<si+*) are valid distribution for all i = 1,2, N k, then q(z) will also be valid distribution.
The theorem reveals that instead of parameterizing the posterior (zz) with Pij in P. we can also use the correlation coefficients r & 1 = 1,... .N- k. Is t| < k s * t to parameterize the local marginals ((zn)4)) and then use the local marginals to construct the global posterior. From the distribution (2). we can see that the number of parameters in I is the same as the number of non-zero Púj in the precision matrix P. Thus, without introducing more parameters, the distribution expressed with local marginals enables us to work on k-variate local marginals rather than the N -ariiitte global posterior, significantly reducing the computational demand.
With the availability of local marginals g(Zaa++) for i = 1, N - k. we can use them to draw samples from the high- figh-dimmnsional global posterior g(z). Specifically, according to the properties of multivariate normal distribution, the conditional distribution of 2 given Zzi+k-1 -=s+k-1 can be expressed as q(2++1/ (=s++=1) (2 As+k. r)=x), with the mean and variance equal to
Ainh Mitk [R. (1) [1 1. [R
(5)
Jisk, (Zai4h-1 1 Ø Oiitk-1,
[ (1) 1441I
(
1.k.1.kk R (1) (LA,4))
Then o44
where denotes element-wise division; [A] means the submatrix of A with rows from to j and columns from s to f. Thus, given the samples from the i-th to the (i + k - 1)-th variable
2- = [Z,,... ... S the sample drawn from g(=e4x can be represented as
Z Auta (T() + (T+)( Eitt, 6 N),,
(4)
where we deliberately write dit and M as (T()) and (T() to emphasize the sample Eet function of correlation parameters p() (74+--1441-) |s, 1,... k + 1, * 1). With the newly obtained sample Z.- as well as the previous samples 2 we can further
3

---

draw the next sample 2 from the conditional distribution ( /2 which can be easily derived from the local marginal g(z++15+++)> By repeating this process sequentially for i 1,2, .N k, we can obtain the sample [ ... EN]] g(z).
should be noted that the sample z can be explicitly expressed terms of coefficients Т. Thus, we use the sample 2 to estimate the expectation of the evidence lower bound (ELBO) L(O, Φ X) Ex-- (log po X, z)] + H(g [z) and give rise to
[(0, Φ, X) log pe [X ZTT) + H[qo(z)),
(5)
where the entropy term H-(-) can be expressed in terms of local marginals thanks to the decomposition as depicted in Eq. (2). The exact expression for the ELBO is provided in Appendix C
To boost inference efficiency, rather than training the coefficients T, it common to parameterize a neural network (* to output the coefficient values as "i+s=1, = f 1,,
where the output value of neural network confined within the interval (-1,1) to be consistent with the range of correlation coefficients. By substituting "i+ts-1.+tt11 fø (xi+s-1. X;+t-1) into
the lower bound (5) the neural network parameters & can be optimized adequately. However, we directly parameterize 7i++-1,i4t01 as fø (x X(+-11) the resulting correlation matrix R
could be non positivii definite. which violates the basic requirement of a Gaussian distribution, causing the whole sampling process and ELBO estimation above invalid.
2.2 Re-parameterizing the Correlation Matrix R ( with Positive Definite Guarantee
To ensure the positive definiteness of R(),, instead of using neural networks to directly parameterize its elements 7,+ we propose new way to parameterize them. Specifically, we notice that for any valid multivariate Gaussian distribution g (Z,,i+)), which is equivalent to have ( > 0, we
can always decompose
qZz++k) = q(2,, 2441 Z4+1444 (1
(6)
where the conditional distribution
(12 214 =N=
2ittk]
2
;
1
of
0
ofth
R
of 0 Ri
of 1): (7)
14 = 1| 1] and p5+ = El 1i 1] are the conditional means; of = B[( 1//224 1444 11/2 are the conditional standard
14) (2411144- 1/ and o E[(zi+k
u
24,44 IZZ is the conditional correlation matrix. Here, 7 represents the conditional correlation parameter between 2; and øith given 2 with index set I,,++ A ii 1,11.......... ................................ .it k 1}. which can be specifically expressed as
deviations; and R M1,++ R
(1)
(8)
7
(8)
ri
[r
141
where (1) [ [ [2:k. 1
R 22 and R
R 12: For conciseness, we use without introducing ambiguity.
the notation + the following context to represent 7i
From we can see that there exists a one-to-one mapping M
7f.+t + Ro.+ that maps the & to the correlation parameter Kis+ in R
conditional correlation parameter ff.+..
For valid distribution g(2ea+*). its conditional distribution g(2;, 1 must be valid, too. This suggests that the correlation matrix R positive definite, which is equivalent to the condition |nf+++| < Therefore, for any valid distribution g(Mas+)), its conditional correlation is ensured to satisfy lnf, < 1. Below, we prove that the converse is also true, that is, we confine h < and set & M().,+))) the correlation matrix ( constructed with it is guaranteed to be positive definite under some condition.

---

Theorem 2.2. By writing the correlation matrix (()) as the following partitioned form
7i
Bitk-1
TiitR 741.44
76,+1
(9)
Witt
74tk-111 1,14k
78.4
2 1.2k-122kk 122411 in the dotted frames are both positive definite, 175,4tal < and we set Watk MM then R is positive definite.
if the upper-left and lower-right sub-matrices R 1. and [R (*)
According Theorem 2.2 if & x le sub-matrices [R and [R 12441244 of the (k 1) x (le + 1) correlation matrix R are both positive definite, and we let Y4,44 MM with h < 1, then the correlation matrix (+) constructed in the form of ( is guaranteed to be positive definite. This gives rise to an iterative construction approach, starting from small sub-matrices and expanding step by step. To illustrate this process, let us take the construction of x correlation matrix as an example, whose eventual form is
112
M(1532)
M(
R
712
1
723
M(
(10)
M(Siap)
723
1
Y34
M(7
M
334
If 7/2 < 1 and 7/2 1, matrinesss (R.",11211 and R (2)33, are known to be positive definite.
Then, if we confine Tis/2 according to Theorem 2.2 the sub-matrix (R(11:33,1 is ensured to be
11:11:3
positive definite. Similarly, if 734: 2/4/3 lie within (-1, 1). we can also ensure [R( > 0. Then, combining with the condition 2/1/23 < we can see from Theorem 2.2 that the correlation
matrix (2) is guaranteed to be positive definite. Continued in this way recursively, positive definite correlation matrices of arbitrary size can be constructed, as depicted in the following corollary.
Corollary 2.3. If all correlation paramelers in Ií {nr,4 1 and r, (nf,++ N- for
t = 2,3. .....,,... ,k lie in the interval (-1,1), then the (I: + 1) (k + correlation matrix R( constructed above is guaranteed to be positive definite.
Therefore, to construct correlation matrix ( with positive definite guarantee, only need to pa- rameterize first-order correlations Fj and higher-order conditional correlations F, for t = 2,3, 2, ... k. and ensure them to lie in the interval (-1, 1). For different orders of correlation coefficients, we can use a specific neural network fØp (.,:) ) to parameterize them as
Ra,+1 fΦ, (x,, X141),
i= 1,2, N
(11)
M1.++ fo. X(+1), i= 1,2 -tt
which represent the first-order correlations and l-order conditional correlations, respectively. Once the positive definite correlation matrix R (8) has been constructed, we can then use the method described safely.
Section 2.1 to optimize the ELBO
The exact cost of inference with our proposed k-order precision matrix P involves three parts: (i) the cost of neural network evaluations for re-parameterizing correlation coefficients, (ii) the cost of sampling from the variational posterior, and (iii) entropy calculation. In our method, to define posterior with k-order correlation
N latent variables. we need specify exactly (N 1) first-order, (N - second-order, and so on, up to (N k) k-order correlations, yielding a total of 1 correlation coefficients. our method, each coefficient is parameterized by the output of re- parameeeriaaaioonn network (- Therefore, to re- parameterize these coefficients, we need to run the network fa(.,-) for (kkV times. The sampling and entropy calculations involve operations like inversion in (3) and determinant computation on lo sub-matrices, incurring a cost of O(k3) FLOPs. Considering that is typically much smaller than N and the complexity of evaluating neural networks, the cost of these operations is negligible compared to the cost of neural network evaluations. Therefore, the total cost approximately amounts to the cost of evaluating O(kN) times of the neural network fo per epoch, which is approximately times of the cost of mean-field amortized VI methods, with the order & controlling the trade-off between expressiveness and computational cost.
(N (N +. (N k) k(2N
5

---

2.3 Extensions to Tree-structured Backbones
Although Theorem 2.1 is built on a chain- structured backbone, our proposed high-order correlations could be extended to more general tree-structured backbones. If the first-order de- pendency structure among latent variables is characterized by tree-structured backhone shown by the solid lines in Figure then we can also impose k-order dependencies over the tree-structured backbone. In this case, every k + 1 consecutive variables the tree forms (ks + 1)-vertex clique C E Cx+1 with indices C = {tv. f.. ... .ik] } S 21.2..... NN. tree-structured backbone, can extend Theorem [21 equivalently express the joint distribution of latent variables using its local marginals over k-vertex cliques C₁ and (A: + 1)-vertex cliques C++1 as follows
sampling
Figure 2: N instances with k-order dependency structure based on a tree-structured backbone.
Based on the
g 2,.." 21p-1 Zive) ( Riz (9(2))."' )2xe)
((2)
II C=(hi:" C=(birn-AeeC,
((=
2
II
(12)
C=flo-h;" JeCnc4 l0 <.n
which is fully determined by local marginals ((io: 2i.) over (k + veretxx cliques C Ck+1- The validity of the local marginals can be similarly guaranteed by parame- terizing the correlation matrices with first-order correlations and higher-order conditional correlations, and further confining them within (-1, 1), as the following corollary shows.
{ia, i.. &&
Corollary 2.4. If the first-order correlations I1 and higher-order conditional correlations I for t= 2, .... ...,k are built upon a tree-structured backbone, and all correlation parameters lie in the interval (-1,1), then can use them construct (k + x (I 1) correlation matrix R with
k-order dependency structure.
Given the marginal distribution and samples Zig:" - draw sample 24, from the condi- tional distribution (2t, (/...... ... _ ). By recursively sampling from the conditional distribution starting from the root node, samples from the joint distribution can be obtained.
3 Related Work
Bayesian inference provides a principled framework for uncertainty estimattinn, but exact inference is often intractable. Variational inference addresses this by approximating the true posterior with more tractable distribution. This requires a trade-off between expressiveness and computational efficiency 8]. A widely-used approach is mean-field variational inference (MFVI) [[], which assumes a fully factorized posterior, treating all latent variables as independent. Despite its broad applicability across domains such image analysis (57) and biology (2). MFVI struggles capture posterior correlations. particularly in settings where latent variables are strongly dependent. To address these limitations, structured variational inference (SVI) enriches variational distributions capture dependencies among latent variables while retaining tractability. Common SVI approaches achieve this through deterministic stochastic transformations, such as normalizing flows [O
61 and implicit models
54 41 Other techniques include modeling local-global dependencies
22
60
using mixture
distributions 42 34 copula- based augmentations (26. [33,, non-conjugacy approximations 28 49]
and hierarchical extensions [9]. While these methods enhance expressiveness, they primarily focus on intra-instance correlations, limiting their scalability to capturing correlations across instances. Another related thread of work is neural relational inference [30 [4]. which models the latent interactions among entities or objects across data points using graph-based representations. While effective in discovering relational structures, these methods focus structure learning and do not explicitly leverage iner-instance dependencies to enhance the variational approximation itself.
Higher-order dependencies have emerged as crucial modeling component in complex systems where first-order representations fall short [46) These dependencies, which account for interactions
involving three or more entities, are prevalent in real-world sequential data such as multivariate
time series, clickstreams (65) citation flows 23), and transportation systems [65]. To capture these
higher-order interactions, several modeling paradigms have been developed, including hypergraphs

---

32
9 simplicial complexes [13 51, motif- networks 77,6, and higher-order Markov models
However, these methods often suffer from scalability issues due to the exponential growth of dependencies. Recent efforts have aimed to incorporate instance-level dependencies into variational inference 43
52
18
63 For instance, DC-GMM [X introduces a prior information matrix to promote similar posteriors across instances for weakly supervised clustering. However, it still relies on a mean-field approximation, which limits its ability to fully capture correlated posteriors. Other methods, such as CVAE [56], attempt to address this limitation by constructing tree-structured variational posteriors, effectively modeling pairwise dependencies among instances. But they remain limited to first-order correlations and struggle represent higher-order dependencies. The work of TreeVI [63] is most similar to ours, but is inherently limited to modeling only first-order correlations. This limitation of VI arises from its reliance on an acyclic tree structure to construct its correlation matrix. While this is sufficient for simple pairwise relationships, attempting to model higher-order correlations inherently introduces loops into the underlying correlation structure. The construction of VT depends on the acyclic property of its backbone and is no longer valid when these loops exist. Moreover, simply modeling higher-order correlation coefficients within the framework of VI does not guarantee the correlation matrix to be positive definite. So even though Tree can capture instance-level correlations, cannot be easily extended to model higher-order correlations. In addition, SIDEC [58] takes a different approach by leveraging variational inference to learn latent dynamics and employing high-order correlations for structural reconstruction. However, its focus is on recovering interaction graphs rather than explicitly modeling high-order dependencies in the latent posteriors themselves.
56
4 Experiments
Tasks & Datasets. We evaluate our method on three tasks: time series anomaly detection, time series forecasting and constrained clustering, using a diverse set of benchmark datasets. For time series anomaly detection, we experiment on three datasets: SMD, SMAP, and MSL. For time series forecasting, we use five widely-used datasets: ETTh1, ETTm1, Electricity, Exchange, and Weather. For constrained clustering, we conduct experiments on four standard datasets: MNIST, Fashion MNIST, Reuters, and 101. Detailed descriptions for each dataset are provided in Appendix [D.1
Baselines & Implementation Details. For time series anomaly detection, we compare our method with four state-of-ther arrte--rtt unsupervised approccc for time series anomaly detection based on VAE:
DAGMM 77). SSST-VAE [45), Omommny [55), and SISVAE (37. For time series forecasting, we compare our method with the state-of-the-art end-to-end methods on time series modeling and forecasting tasks, including VRAE [20]. Informer [69]. GRU-NVP (47). and DecpAR [50). For constrained clustering,
compare our approach against traditional algorithms such PCKMeans 4. SDEC 48), C-IDEC [67. and the state-of-the-art DC-GMM 38]; and also benchmark against generative models such as VaDE (27. DGG [66), and TI (63] To evaluate the effectiveness of our method, we conduct experiments with varying levels of o-order dependency structures, specifically using & E {1. 1, 5, 10). Further implementation details are provided in Appendix D.3
4.1 Time Series Modeling
Generative Time Series Modeling aims to learn the underlying prob- ability distribution of time series data and generate new, synthetic time series samples that exhibit similar characteristics to the ob- served data. However, the major- ity of existing approaches often ignore instance-level correlations during posterior inference, thus failing to comprehensively cap- ture the temporal dynamics of time series. To address this limitation, our method incorporates two key adaptations compared to the vanilla VAE, First, temporal dependencies are introduced by integrating Gated Recurrent Units (GRUs) [15 in both the VAE encoder and decoder. Second, the
Table 2: F1-Score and Evidence Lower Bound comparisons.
Dataset
SMAP
MSL
SMD
Metric
Fl
ELBO
F1
ELBO
F1
ELBO
DAGMM
0.7105
-115.2820
0.7007
-277.3380
0.7094
-155.9460
LSTM VAE
0.7298
1116.9500
0.6780
281.3220
0.7842
-14600500
OmmiAnonaly
0.8434
-98.9217
0,8849
-161.0002
0.8857
-72,0419
SISVAE
0.8299
1101.1878
0.8766
182.6060
0.8775
-7258332
I-order
0.8411
-97.6057
0.8883
-165,5004
0.8901
-69 5278
HoT.V
3-order
0.8552
-95,2314
0.8940
-157.2134
0.9153
.67.4001
10-order
0,8636
-92,2948
0.9145
-134.0815
0.9284
-65.0345
7

---

Table I: Multivariate time series forecasting results with horizon H E {24, 48, 168,336, 720). Best performance is highlighted in bold font and the second best results are underlined.
HoT-VI (Ours)
Method
Informer
GRU-NVP
DecpAR
VRAE
1-order
3-order
10-order
MSE
MAE
MSE
MAE
MSE
MAE
MSE
MAE
MSE
MAE
MSE
MAE
MSE
MAE
24
0.577
01,549
3,540
0.733
1.166
0.836
0.743
0.762
0.664
0.570
0.543
0.505
0.363
0.376
ETTml ETThl
48
0.685
0.625
2.549
0.622
1.154
0.827
0.826
0.801
0.705
0,597
0.578
0.528
0.392
0.392
168
0.931
0.752
3,831
0.774
1.083
0.778
1.070
0.938
0,848
0,681
0.721
0.615
0.510
0.464
336
1.128
0.873
6.877
1.008
1.043
0.766
1.199
1.016
0,990
0.755
0.883
0.702
0.616
0.525
720
1,215
0.896
5.377
1.060
1.075
0.795
1.426
1.164
1.129
0.821
1.021
0.781
0.763
0.630
24
0.453
0,444
0.605
0.437
1,360
0.871
0.687
0,646
0,488
0,455
0.409
0,417
0.253
0.298
48
0.494
0,503
2.787
0,701
1.334
0.866
0.817
0.724
0.648
0,544
0.535
0.488
0.330
0.345
168
0.678
0.614
4.212
0.824
1.170
0.838
0.853
0.794
0.686
0.573
0.578
0.521
0.368
0.373
336
1.056
0.786
5.062
1.019
1.249
0.846
1.091
11.975
0.771
0.628
0.641
0.567
0.434
0.415
720
1.192
0,926
5.799
1.075
1.075
0.770
1.165
0.996
0,886
0.692
0.737
0.626
0.528
0.474
Exchange Electriity
24
0,312
0,387
3,514
1. 1,844
0.211
0.330
0.279
0.396
0.326
0,400
0.256
0.346
0.134
0.238
48
0.392
0,431
3,318
1,786
0.332
0,398
0.317
0.410
0.347
0,415
0.277
0,363
0.152
0.255
168
0.515
0.509
3,482
1.833
1.065
0.811
0.366
0.475
0.373
0.433
0,303
0,382
0.174
0.273
336
0.759
0.625
3,921
1.941
1,040
0.795
0.402
0.515
0.388
0.445
0.319
0.395
0.194
0.293
720
0.969
0.788
4.232
2,020
1.048
0.804
0,450
0.556
0.415
0.463
0.348
0.416
0,230
0.323
24
0,611
0.626
1,557
0.877
1.328
0.692
0.140
0,310
0.098
0.227
0.093
0.227
0.033
0.126
48
0.680
0.644
1.589
0,883
1.345
0.701
0.238
0.435
0.155
0.267
0.171
0,306
0.058
0.164
168
1.097
0.825
1.663
0.903
1.434
0.745
0.642
0.703
0.379
0.466
0.368
0.458
0.196
0.326
336
1.672
1.036
1.682
0.905
1,489
0.778
1.050
0.9:53
0.992
0.835
1.165
0.821
0.496
0.515
720
2.478
1.310
L.748
0.928
1,526
0.793
3.003
1.593
1.988
1.063
2.029
1.090
1.508
0.857
24
0.162
0,235
1.222
0.909
0.205
0.250
0.227
0.315
0,206
0.294
0.186
0,281
0.129
0.179
Weather
48
0.348
0,400
2.319
1.287
0.229
0.267
0.449
0.495
0.325
0.385
0.291
0.361
0.186
0.230
168
0,444
0,463
2.174
1.165
0.344
0.343
0.563
0.648
0.466
0.506
0.429
0.486
0.294
0.313
336
0.5788
0.523
2.119
1.221
0.568
0.527
0.781
01.44
0.767
0.645
0.625
0.575
0.550
0.430
720
1.059
0.741
2.621
1,303
0.571
0.533
1.125
1.058
0.998
0.727
0.808
0.653
0.772
0.510
Average
0.819
0.660
3,112
1.122
0.978
0.678
0.796
0.741
0,641
0,555
0.573
0,516
0.386
0.373
posterior is approximated using k-order dependency. The formal representations are as follows. Let X = =x { denote time series comprising N observations, where each x; E RC represents observations across C channels at time step t. Given a window of T observations X the encoder is represented as: μl, σ f ([e,-T, ..., eq]) ,where the hidden state of GRU encoder e, is updated by x; and e;-1. Similarly, the decoder is represented as: X;-T: = fe fe([d-TT TT+1..... di]), where d, is the hidden state of the GRU decoder and updated by z, and d,-1. We evaluate the modeling capacity of our method through its performance on two downstream tasks.
Time Series Anomaly Detection The objective of this task is to determine whether an observation x, is anomalous based on the preceding T observations. Our model can be directly applied to the anomaly detection task by reconstructing data. Trained solely on normal data, the model is expected exhibit low reconstruction loss for normal data while high for anomalies. Consequently, anomalies are identified by comparing the reconstruction loss against a threshold. The ELBO serves as a metric to evaluate the modeling capacity for normal data, while the F1 score assesses anomaly detection performance. As shown in Table our method demonstrates superior ELBO and Fl scores compared to other generative approaches that neglect instance-level correlations during posterior inference. Notably, even with only first-order dependencies, our method achieves comparable performance to Onomal, a complex model integrating VAE, flow, and State Space Models (SSM). Furthermore, increasing the order of dependencies in our model leads to consistently higher ELBO and Fl scores than all baselines. This indicates that modeling higher-order temporal relationships in time series improves data modeling and anomaly detection performance. By modeling k-order dependency, our model captures fine-grained local dynamics and roarser-grained long-term dependency, leading to a more robust and comprehensive understanding of the complex temporal structure.
Time Series Forecasting This task aims to predict the subsequent H observations given L past observations. Formally, this is a mapping f X,- E L 1 X Y we omit the subscripts hereafter. Our approach is decomposed into two steps: first, learn an expressive and predictable representation of the historical observations via generative modeling; second, perform forecasting based on the representation. For generative modeling, we capture instance-level correla- tions using k-order dependency that existing approaches often overlook. For forecasting, we integrate a feed forward network fo Z + Y into the original model. Formally, we aim to optimize the
IR Cx and
8

---

Table 3: Clustering performances (%) of our proposed method compared with baselines. Means and standard deviations are computed across 10 runs with different random initializations.
HoT-VI (Ours)
Datasel
Metric
VaDE
SDEC
C-IDEC
DGG
DC-GMM
VI
3-order
S-urder
10-order
MNIST
ACC
89.0+5.0
86.2+0.1
96,340.2
95.8+0.1
96.5+0.2
97. 37.4
18.14
98.340.44
98.5+0.3
NMI
82.8+3.0
84,2=0,1
91. 91.8+1.0
91.2+0.2
91,4+0.3
13.1
93.8+0.4
94.2+0.3
94.6+0.3
ARI
80,9+5.0
80.1=0.1
92.1+0.4
91 420.3
92.5+0.5
93.7=0.7
94.920.6
95.2+0.5
95.620.5
IMNIST
ACC
55.1+2.2
54.0=0,2
68.
79.9+0.4
80.5+0.8
81,4+0.6
82,9+0.5
83.2+0.5
83,4=0.4
NMI
57.912.7
57.3.0.1
66.742.0
70.
72.010.4
73.9
74.720.6
74.8
75.1
ARI
41.6+3.1
40.2=0.1
52.3+3.0
64, 1.33
66.4+0.5
67.9+0.9
68.920.5
69.1
69.2+0.4
Reuters
ACC
76.010.7
82.1
94.740.6
93.5
95.410.2
95.9
96.820.6
97.2+0.6
97.620,5
NMI
50.1+1.3
62.3.0.1
81.440.7
81.2+0..
82.740.7
83 33.450.5
84.810.6
15.1
05.4
ARI
58,0+1.4
66.7=0.1
87.7=0.9
87.8=0.5
89.0+0.6
90.2+0.4
91.3+0.5
91.6+0.55
92.0=0.4
STL-10
ACC
77.3+0.5
79.2=0.1
81.6+3.8
89. 9+0.3
89.5+0.5
90 00.4
91.8+0.77
92.2+006
92.4=0.4
NMI
70.6+0.4
78.6=0.1
77.3+1.7
80,9=0.5
80.2+0.7
81.3+0.8
82.4+0.7
82.8+0.6
83.1=0.5
ARI
62.710.4
71.0=0.1
71.823.4
79.020.4
78.420.9
79.5
80.920.7
81.310.5
81.5=0.5
( joint model p(X, Y) = p(X(( - 22,
by maximizing the evidence lower bound: logp(X, Y) 2 log 52 Pv ( X + g x) [log (X|Z))
KL (9 (Z X) [z()), where the first term measures prediction accuracy, typically estimated using L1 or L2 loss, while the subsequent terms serve as regularization for representation learning. The results presented in Table show that our approach outperforms all baselines across all datasets and metrics with only two second-best exceptions. Furthermore, the consistent prediction improvement in our method with increasing dependency order underscores that time series forecasting can benefit from better time series modeling.
4.2 Constrained Clustering
Constrained clustering is a task that incorporates instance-level constraints into the clustering process, allowing users to enforce specific relationships between data points based on prior knowledge. These constraints are expressed by correlation graph G = (V.& A), where denotes the set of instances, and the edge set & = Ess U &c consists of must-link constraints &Mr requiring two instances to be in the same cluster, and cannot-link constraints &c. which require them to be in different clusters. The adjacency matrix A E R encodes both the type and strength of each i > 0 if (i,j) E EM. [A]i < if (i,j) E Èc, and [A]is
if no constraint exists. [0, (0) reflects the confidence in the constraint. Following the generative modeling framework of previous work (38], constrained clustering can be formulated as a probabilistic clustering problem with joint probability pe [X, Z, A) pe X Z p (Z c) A), where the data x; is generated from a normal distribution conditioned on z,; the latent embedding z; is drawn from a cluster-dependent normal distribution p(2,(ca) = ./ (Z.; Me, diea, 2 )); and the cluster assignments c {c N follow distribution conditioned on A, defined as p(c A) अ(()) where (,, A) = I exp([A]oj is weighting function with ō representing the indicator function, π are the cluster weights, and ?((T)) = Ee II, Te, hy A) is normalization constant.
constraint: [A
The magnitude [Ally
1 II, h;,,, A),
To perform inference, where g(c/Z) = II. g(c; ) is computed using Bayes' rule. In standard approaches like DC-GMM 38), the posterior q [ (Z is modeled as fully factorized, which ignores dependencies between instances. We address this limitation by introducing a higher-order dependency structure over the latent space. Specifically, we approximate qu (Z using k-order correlations, where first-order dependencies are guided by a tree structure learned from the correlation graph G. We follow the work of (63] to learn the tree structure from data by optimizing a symmetric adjacency matrix. In our experiments, we set & E {3, 5, 10) and compare our model with baselines over 10 independent runs, reporting average Accuracy (ACC), Normalized Mutual Information (NMI), and Adjusted Rand Index (ARI) in Table 3 The results show that our approach outperforms existing methods across all datasets and metrics, demonstrating the effectiveness of incorporating higher-order correlations in constrained clustering. The averaged improvements of our method incorporating third-order dependency structure are 1.93, 2.35 and 2.43 ACC, NMI and ART against DC-GMM and are 1. 13, 1.00 and 1.18 against Tree underscoring the significance of considering dependencies among latent
use variational posterior of the form (Z, c|X) q qe (Z| )gfc Z),

---

Table 4: Additional experiments of HoT-VI with orders & exceeding 10, including time series anomaly detection on SMAP dataset, time series forecasting on ETThl dataset with horizon 24, and constrained clustering on MNIST dataset.
Methods
Mean-field
k =
k =
3
k
10
k = 50
k = 100
Time Series Anomaly Detection
Runtime (s)
1.00
2.44
8.51
27.60
142.58
294.53
F1
0.7774
0.8411
0.8552
0,8636
0.8711
0.8755
ELBO
-109.2182
99.6057
-95.2314
-92.2948
-90.4291
-89.9577
Time Series Forecasting
Runtime (s)
4.80
11.45
36.94
126.44
675.75
1340.22
MSE
0.739
0.664
0,543
0.363
0.348
0.333
MAE
0.716
0.570
0,505
0.376
0.362
0,352
Constrained Clustering
Runtime (s)
0.25
0.59
1.84
6.22
29.53
60.59
ACC (%)
96.50
97.55
98.12
98.52
98.62
98.69
NMI (%)
91.37
93.44
93.80
94.55
94.63
94.85
ARI (%)
92.54
93.89
94.89
95.65
95.85
96.09
posteriors, particularly higher-order dependencies. Furthermore, the performance of our method consistently improves with increasing dependency order, benefiting from the ability of higher-order correlations to jointly link larger set data instances. This facilitates more effective propagation of cluster assignment constraints compared to methods limited to pairwise dependencies, further underscoring the importance of capturing high-order interactions in constrained clustering.
Performances At Higher Orders The choice of the order k is trade-off between model expres- siveness and computational cost. Generally, as increases, the model's performance consistently improves, as demonstrated in our experimental results. However, as seen from Table the perfor- mance gains diminish as the order (e-g. 50, 100) goes higher. However, the computational cost always scales linearly with k. To balance the gains and cost, we set & to moderate values (up to 10) in our main experiments. By setting to moderate value (e.g. 10), we can only model correlation up to 10-th order, losing the ability to model higher-order correlations. But observed from Table the gains become increasingly weak as the order goes higher.
Conclusion
In this work, we introduced novel variational inference framework for modeling higher-order correlations among latent variables, going beyond the limitations of mean-field and first-order methods. By equivalently formulating the posterior as composition of local marginals, our approach enables expressive k-order dependency modeling. To ensure tractability, we proposed an iterative procedure that guarantees positive definiteness the resulting correlation matrix via conditional correlation parameterizations. This formulation enables reparameterized sampling and allows efficient optimization. We further generalized the model support tree-structured backbone dependencies, enabling flexible incorporation of more structured latent correlations. Empirical results across diverse tasks, including time series modeling and constrained clustering, demonstrate the effectiveness of our method in capturing complex dependency structures and improving downstream performance.
Limitations & Future Work The proposed method requires specifying backbone structure to construct higher-order correlations. This limitation is mitigated by generalizing to learnable tree-structured backbones. For future work, we will investigate the combination of instance-level and dimension-level correlation structure, to further enhance the expressivity of posterior approximation.
Acknowledgement This work is supported by the National Natural Science Foundation China (No. 62276280), Guangzhou Science and Technology Planning Project (No. 2024A04J9967).
10

---

References
[1] Agrawal, A. and Domke, (2021). Amortized variational inference for simple hicrarchical models. In Neural Information Processing Systems.
[2] Akram, M., Adnan, M., Ali, S. F., Ahmad, J.. Yousef, A., Alshalali, T. A. N., and Shaikh, Z. A. (2025). Uncertainty-awar diabetic retinopathy detection using deep learning enhanced by Bayesian approaches. Scientific Reports, 15(1):1342.
[3] Ambrogioni, L., Lin, K., Fertig, E., Vikram, S., Hinne, M., Moore, D., and van Gerven, M. (2021). Automatic structured variational inference. In International Conference on Artificial Intelligence and Statistics, pages 676 186464..P7 PMLR.
[4] Basu, S., Banerjee, A., and Mooney, R. J. (2004). Active semi-supervision for pairwise con- strained clustering. In Proceedings of the 2004 SIAM international conference on data mining.
[5] Battiston, F., Cencetti, G., lacopini, I., Latora, V., Lucas, M., Patania, A., Young. J.-G., and Petri, G. (2020). Networks beyond pairwise interactions: Structure and dynamics. Physics Reports, 874:1-92.
[6] Benson, A. R., Abebe, R., Schaub, M. T., Jadhahaie, A., and Kleinberg, J. M. (2018). Simplicial closure and higher-order link prediction. Proceedings of the National Academy of Sciences, 1 115:E11221 E2300.
[7] Benson, A. R., Gleich, D. F., and Leskovec, J. (2016). Higher-order organization of complex networks. Science, 353:163 166.
[8] Blei, D. M., Kucukelbir, A., and McAuliffe, J. D. (2016). Variational inference: A review for statisticians. Journal of the American Statistical Association, 112:859 - 877.
[9] Carletti, T., Fanelli, D., and Nicoletti, S. (2020). Dynamical systems on hypergraphs. Journal of Physics: Complexity, 1(3):035006.
[10] Caterini, A., Cornish, R., Sejdinovic, D., and Doucet, A. (2021). Variational inference with contimuously-indexed normalizing flows. In Uncertainty in Artificial Intelligence.
[11] Chaikin, P. M. and Lubensky, T. C. (1995). Mean-field theory, page 144-212. Cambridge University Press.
[12] Chen, Q., Shi, W., Sui, D., and Leng, S. (2023a). Distributed consensus algorithms in sensor networks with higher-order topology. Entropy, 25.
[13] Chen, Q., Zhao, Y., Li, C., and Li, X. (2023b). Robustness of higher-order networks with synergistic protection. New Journal of Physics, 25(11):113045.
[14] Chen, S., Wang, J., and Li, G. (2021). Neural relational inference with efficient message passing mechanisms. In AAAI, pages 7055-7063. AAAI Press.
[15] Cho, K., van Merrienboer, B., Gülçehre, Ç., Bahdanau, D., Bougares, F., Schwenk, H., and Bengio, Y. (2014). Learning phrase representations using RNN encoder-decoder for statistical machine translation. In Moschitti, A., Pang, B., and Daelemans, W., editors, Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing, EMNLP 2014, October 25-29, 2014, Doha, Qatar, A meeting of SIGDAT, a Special Interest Group of the ACL, pages 1724-1734. ACL.
[16] Coates, A., Ng, A., and Lee, H. (2011). An analysis of single-layer networks in unsupervised feature learning. In International Conference on Artificial Intelligence and Statistics.
[17) Deng, I.,, Chen, X., Zhao, Y., and Zheng, K. (2021). Hifi: Anomaly detection for multivariate time series with high-order feature interactions. Database Systems for Advanced Applications: 26th International Conference, DASFAA 2021, Taipei, Taiwan, April / 1-14, 2021, Proceedings, Part I 26, pages 641-649. Springer.
[18) Edler, D., Bohlin, I.,, and Rosvall, M. (2017). Mapping higher-order network flows in memory and multilayer networks with infomap. Algorithms, 10(4).

---

[19] Ektefaie, Y., Dasoulas, G., Noori, A., Farhat, M., and Zitmik, M. (2022). Multimodal learning with graphs. Nature Machine Intelligence, 5:340-350.
[20] Fabius, O., van Amersfoort, J. J.R., and Kingma, D. P. (2015). Variational recurrent auto-encoders. In ICLR ( (Workshop).
[21] He, H., Zhang, Q. Y., Shi, K., Niu, Z., and Cao, L. (2022). Distributional drift adaptation with temporal conditional variational autoencoder for multivariate time series forecasting. IEEE Transactions on Neural Networks and Learning Systems, 36:7287-7301.
[22] Hoffman, M. D. and Blei, D. M. (2014). Stochastic structured variational inference. In International Conference on Artificial Intelligence and Statistics.
[23] Holmgren, A., Edler, D., and Rosvall, M. (2023). Mapping change in higher- networks with multilevel and overlapping communities. Applied Network Science, 8(1):42.
(24] HUMNABADKAR, A. S., SIKDAR, A., Cave, B., ZHANG, H., BAKAKI, P., and BEHERA, (2024). High-order evolving graphs for enhanced representation of traffic dynamies. In ECCV 2024 2nd Workshop on Assion-Centite Autonomous Driving (VCAD): Poster sessions.
[25] Hundman, K., Constantinou, V., Laporte, C., Colwell, I., and Söderström, T. (2018). Detecting spacecraft anomalies using LSTMs and nonparametric dynamic thresholding. Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining.
[26] Janke, T., Ghanmi, M., and Steinke, F. (2021). Implicit generative copulas. In Neural Informa- tion Processing Systems.
[27] Jiang. Z., Zheng, Y., Tan, H., Tang, B., and Zhou, (2016). Variational deep embedding: An unsupervised and generative approach to clustering. In International Joint Conference on Artificial Intelligence.
[28] Khan, M. E. and Lin, W. (2017). Conjugate-compptation variational inference: Converting variational inference in non-conjugate models to inferences in conjugate models. In AISTATS, volume 54 of Proceedings of Machine Learning Research, pages 878-887. PMLR.
[29] Kim, K., Kim, J., Zaheer, M., Kim, J. S., Chazal, F., and Wasserman, L., A. (2020). Pllay: Efficient topological layer based on persistent landscapes. In Neural Information Processing Systems.
[30] Kipf, N., Fetaya, E., Wang, K., Welling, M., and Zemel, R. (2018). Neural relational inference for interacting systems. In ICMI,, volume 80 of Proceedings of Machine Learning Research, pages 2693-2702. PMLR.
[31] Ko, J., Kim, K., Kim, W., and Gardner, J. R. (2024). Provably scalable black-box variational inference with structured variational families. In ICML. OpenReview.net.
[32] Kovalenko, K., Romance, M., Vasilyeva, E., Aleja, D., Criado, R., Musatov, D., Raigorodskii, A., Flores, J., Samoylenko, I., Alfaro-Bittner, K., Perc, M., and Boccaletti, S. (2022). Vector centrality in hypergraphs. Chaos, Solitons & Fracials, 162:112397.
[33] Lai, G., Chang, W.C.., Yang, Y., and Liu, H. (2017). Modeling long- and short-term temporal patterns with deep neural networks. The 41st International ACM SIGIR Conference on Research & Development in Information Retrieval.
[34] Lambert, M., Chewi, S., Bach, F., Bonnabel, S., and Rigollet, P. (2022). Variational inference via Wasserstein gradient flows. Advances in Neural Information Processing Systems.
[35] LeCun, Y. and Cortes, C. (2010). MNIST handwritten digit database.
[36] Lewis, D. D., Yang, Y., Russell-Rose, T., and Li, F. (2004). Rev1: A new benchmark collection for text categorization research. Journal of machine learning research, pages 361-397.
[37] Li, L., Yan, J., Wang, H., and Jin, Y. (2021). Anomaly detection of time series with smoothness- inducing sequential variational auto-encoder. IEEE Transactions on Neural Networks and Learning Systems, 32(3):1177-1191
12

---

[38] Manduchi, L., Chin- Chin-Cheong, K., Michel, H., Wellmann, S., and Vogt. J. (2021). Deep con- ditional Gaussian mixture model for constrained clustering. Advances in Neural Information Processing Systems.
[39] Manduchi, L.,, Vandenhirtz, M., Ryser, A., and Vogt, J. (2023). Tree variational autoencoders. Advances in Neural Information Processing Systems.
[40] Modak, A. and Mauritsen, T. (2023), Better-constrained climate sensitivity when accounting for dataset dependency on pattern effect estimates. Atmospheric Chemistry and Physics.
[41] Moens, V. Ren, H., Maraval, A., Tutunov, R., Wang. J., and Ammar, H. (2021). Efficient semi-implicit variational inference. CoRR, abs/2101.06070.
[42] Morningstar, W. R., Vikram, S., Ham, C., Gallagher, A., and Dillon, J. V. (2020). Automatic differentiation variational inference with mixtures. In International Conference on Artificial Intelligence and Statistics.
[43] Ou, Z., Su, Q.. Yu, J., Liu, B., Wang. J., Zhao, R., Chen, C., and Zheng, Y. (2021). Integrating semantics and neighborhood information with graph-driven generative models for document retrieval. In ACL/ICCNLP (1), pages 9238-22499
[44] Paninski, L. (2003). Estimation of entropy and mutual information. Neural Computation, 15:1191-1253.
[45] Park, D., Hoshi, Y., and Kemp, C. C. (2017). A multimodal anomaly detector for robot-assisted feeding using an LSTM-based variational autoencoder. IEEE Robotics and Automation Letters, 3:1544-1551.
[46] Peixoto, T. and Rosvall, M. (2017). Modelling sequences and temporal networks with dynamic community structures. Nature Communications, 8(1):582
[47] Rasul, K., Sheikh, A.S Schuster, I., Bergmann, M., and Vollgraf, R. (2021). Multivariate probabilistic time series forecasting via conditioned normalizing flows. In International Conference on Learning Representations.
[48] Ren, Y., Hu, K., Dai, X., Pan, L., Hoi, S. C. H., and Xu, (2019). Semi-supervised deep embedded clustering. Neurocomputingg pages 121-130.
[49] Salimbeni, H., Eleftheriadis, S., and Hensman, J. (2018). Natural gradients in practice: Non- conjugate variational inference in Gaussian process models. In AISTATS, volume 84 of Proceedings of Machine Learning Research, pages 689-697. PMLR.
[50] Salinas, D.,, Flunkert, V., Gasthaus, J., and Januschowski, T. (2020). Deepar: Probabilis- tic forecasting with autoregressive recurrent networks. International journal of forecasting, 36(3):1181-1191
[51] Santoro, A., Battiston, F., Petri, G., and Amico, E. (2022). Higher-order organization of multivariate time series. Nature Physics, 19:221-229.
[52] Scholtes, I., Wider, N., and Garas, A. (2016). Higher-order aggregate networks in the analysis of temporal networks: path structures and centralities. The European Physical Journal B, 89)3):11.
[53] Smith, M. S. and Loaiza-Mayy Maya, R. (2021). Implicit copula variational inference. Journal of Computational and Graphical Statistics.
[54] Song. J., Meng, C., and Ermon, S. (2021). Denoising diffusion implicit models. 9th International Conference on Learning Representations.
[55] Su, Y., Zhao, Y.. Niu, C., Liu, R. Sun, W. and Pei, (2019). Robust anomaly detection for multivariate time series through stochastic recurrent neural network. Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery Data Mining.
[56] Tang, D., Liang, D., Jebara, T., and Ruozzi, N. (2019). Correlated variational auto-encoders. In International Conference on Machine Learning.
13

---

[57] Tölle, M. and Schlaefer, A. (2021). A mean-field variational inference approach to deep image prior for inverse problems in medical imaging. In International Conference on Medical Imaging with Deep Learning.
[58] Wang, A. and Pang, J. (2024). Structural inference with dynamics encoding and partial correlation coefficients. In International Conference on Learning Representations.
[59] Wang, H., Bhattacharya, A., Pati, D., and Yang. Y. (2022a). Structured variational inference in Bayesian state-space models. In International Conference Artificial Intelligence and Statistics.
[60] Wang, Q and Van Hoof, H. (2020). Doubly stochastic variational inference for neural processes with hierarchical latent variables. In International Conference on Machine Learning.
[61] Wang, Y., Liu, F., and Schiavazzi, D. E. (2022b). Variational inference with nofas: Normalizing flow with adaptive surrogate for computationally expensive models. Journal of Computational Physics.
[62] Xiao, H., Rasul, K., and Vollgraf, R. (2017). Fashion-mnist novel image dataset for benchmarking machine learning algorithms. ArXiv, abs/1708.007747
[63] Xiao, J. and Su, Q. (2024). VIe Reparameterizable tree-structured variational inference for instance-level correlation capturing. In The Thirty-eighih Annual Conference Neural Information Processing Systems.
[64] Xie, J., Girshick, R. B., and Farhadi, A. (2016). Unsupervised deep embedding for clustering analysis. In JCML.
[65] Xu, J., Wickramarathne, T., and Chawla, N. (2016). Representing higher-order dependencies in networks. Science Advances, 2:e 600028.
[66] Yang. L., Cheung, N.-M., Li, J., and Fang, J. (2019). Deep clustering by Gaussian mixture variational autoencoders with graph embedding. In Proceedings of the IEEÉ/CVF international conference computer vision.
[67] Zhang, H., Basu, S., and Davidson, (2019). A framework for deep constrained clustering algorithms and advances. In ECML/PKDD.
[68] Zhao, Y. and Linderman, S. W. (2023). Revisiting structured variational autoencoders. In International Conference on Machine Learning.
[69] Zhou, H., Zhang, S., Peng, J., Zhang, S., Li, J., Xiong, H., and Zhang, W. (2021). Informer: Beyond efficient transformer for long sequence time-series forecasting, In Proceedings of the AAAI conference on artificial intelligence, volume 35, pages 11106-11115
[70] Zong, B., Song, Q., Min, M. R., Cheng, W., Lumezanu, C., ki Cho, D., and Chen, H. (2018).
Deep autoencoding Gaussian mixture model for unsupervised anomaly detection. In International Conference on Learning Representations.
14

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The main claims in the abstract and introduction have clearly reflected the paper's main contributions in the following context.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification:
Limitations our work has been discussed in Section the paper.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
the paper. The authors should use their best
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
15

---

Justification: Assumptions for the theoretical results are given in the descriptions of theorem or corollary, and the proofs are included Appendix Theorem is proved in Appendix
A.1. Theorem 2.2 is proved in Appendix A.2. Corollary 2.3 and Comollar 2.4 are proved in
Appendix A.3
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
The datasets, baseline methods and implementation details needed to reproduce the experimental results are included in Section and we refer to Appendix D.3 for more experimental details.
Justification:
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, other means that appropriate to the research performed.
also be provided via detailed
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (c.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (c.g., to registered users), but it should be possible for other researchers to have some path to reproducing verifying the results.
16

---

5. Open access to data and code
Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, described in supplemental material?
Question:
Answer: [Yes]
Justification: The datasets and experimental settings are provided in Section 4 and Appendix D.3 and our code is submitted to the GitHub repository.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://niis.cc.
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (c.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https
/nips. wwtpcc/pcccccpcccceooooooop httee/ggguudee/oeeeuuubbiiiiiiii ieaibbmsssionoo swiniisssiiiiiiyyy Policyy for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: Training details of all our experiments have been specified in Appendix D.3 Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core the paper to a level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, appendix, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: The statistical significances of our experiments that support the main claims of the paper are shown in Table Table a Table respectively. And the related metrics have been depicted in the context.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emor bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
17

---

The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables figures symmetric error bars that would yield results that are out range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
Question:
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Answer: [Yes]
Justification:
The computer resources needed to reproduce our experiments have been
specified in Appendix D.4
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g.. preliminary or failed experiments that didn't make into the paper).
9. Code of ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics http::////eeeeep..pcppppptttteat
Answer: [Yes]
Justification: We have reviewed the NeurIPS Code of Ethics, and make sure that our research follows the Code of Ethics.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make to preserve anonymity (e.g.. if there is special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive socictal impacts and negative socictal impacts of the work performed?
Answer: [NA]
Justification: Our work mainly focuses on basic theory about variational inference with instance-level and higher-order correlation structure, meaning that there is no societal impact to be addressed.
18

---

Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA No, they should explain why their work has societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (c.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., denloymemt of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology being used intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms monitor how system learns from feedback over time, improving the efficiency and accessibility ML).
11. Safeguards
Question:
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Answer: [NA]
Justification: Our experiments are conducted on standard and public datasets available for everyone, and our proposed methods focus on basic theory with regards variational inference without safety risks.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (c.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: Datasets used in experiments are all public, and their related papers have been cited in our paper.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
19

---

The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g.. CC-BY 4.0) should be included for each asset. For scraped data from a particular source (e.g.. website). the copyright and terms of service of that souree should be provided.
assets are released, the license, copyright information, and terms use in package should be provided. For popular datasets, paperssi thoode com com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors encouraged to reach out to the asset's ' creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification:
No new assets are introduced or released in our paper.
Guidelines:
The answer NA means that the paper does release new assets.
Researchers should communicate the details the datasct/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: Our work does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurlPS Code of Ethics, workers involved data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: Our work does not involve crowdsoureing research with human subjects
Guidelines:
20

---

The answer NA means that the paper does not involve crowdsoureing nor research with
human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and expect authors adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Does the paper describe the usage of LLMs if it is important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification: The core methodology, scientific rigorousness, and originality of our paper are unrelated to LLM usage.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer to LLM for what should or should not be described.
21

---

A Proofs
Followings are the details of proof's of the claim from the main text Theorem Theorem 2.2
Corollary 2.3 and 2.4
A.1 Proof of Theorem 2.1
Lemma A.l. Suppose that A and 2g are conditionally independent given ze where A, B,C.C {i,i +1... , k} are mutually exclusive, then in terms of probabilities.
28, 2c) q(ZA, zc ( zc)
g(ZA, 2n,
(13)
g(zc)
Proof. According to the definition of conditional independence.
(Z 2n ac) = g(Z. bc (alcc(((((zz)))) = g(Za, 2c) 2c) zc)
(14)
((cc)2
Multiplying both sides by q(c)), we can obtain
q(ZA. 2c ( 2c)
( ZB, 2c)
(15)
g(zc)
Theorem A.2. For any joint distribution q(z) = N2; μs, P-1) with a precision matrix P that has k-order connection struclure, it can be equivalently expressed
N-k-1
(Ziitt))
((2)
II i=1
((za+t-1)) 1) II
(16)
((Zsst-1)9(
where g(zea4k-1) 1) and ((zzii+)))) are the marginals of q(z) over 2istk-1 and Zish ko Moreover, if the (k + 1)-variate marginals ((zsstk) are valid distribution for all i 1,2, ... N k, then q(z) will also be valid distribution.
Proof. To prove the theorem, turn prove that for any and f1.... N -t. the (marginal) distribution of Ze:t is given by
1
t-k
g(zii+k)
((Za:))
II s=
g g(Zai4R-))
(17)
q (ZZist ()((2++114)))
I=s
The result is trivial for - = k. where the right-hand side becomes
g(Zs:*+k))
g(Zs:st--1)ed (oosst4-14qq22 1)0 alZ ( (2441:444)) ((Zxs+t) = g(Zo:e).
(18)
To start the induction proof, we first prove for
which is
g(zastk+z) ((Z.:+t-1))9( (Zz
(Zs
((2.41946+1)
g 1) Z4+1 q( (( (Z64222444)
(19)
2 (1
(Z,++1stk)
By letting A = (s), B (s k + 1) and C = {s 1,... 8 k). then the equation above is direct conclusion of Lemma A. where we use the condition that 2, and 2x*8+1 are conditionally independent given zg Zetle+* implied by the k-order dependency structure. We proceed by
22

---

induction and from t - 8 /tott ++1.
The induction hypothesis gives
1-k+1
q ((z.++k))
g(zz:))
II
g ((za44-1)
II
(20)
@ (Zzzii (2
for any à: 5-5- S 1, and we want to show that for - ----1 1.
s+I-k+2
s+1-k+1
q ((iatt))
q( ((2s:0(411)
II
((zaatk-1)
(21)
Letting A {s}, B -s (+ 1} and C {s 1,+1,- 1, + and then applying Lemma
A.1 gives
us
g(Zs:4441))
(22)
(2,+11:44)
where we can use the induction hypothesis to obtain
s+1-k+1
841-
((ziitk)
((2x:s+1)=
II i=
(Z
2 k (g(2
i=
9+11-k+2
8+1-k+1
q ((zzittk)
g g(2++1:ot++1>
II f=8-1
q (Z.
(23)
g) (z)4))iit(+
i=++
s+1-k++
s+1-k
((Zni+*))
g(z++1:01)
II
II
ims+1
Leveraging them to simplify the right-hand side of Eq.
(22)
we can obtain
((bss+k)
RHS = ((2s:stk-1) x
((
199
s+1-k+2
8+1 k+1
(Zistt))
II
(Zut+bb1))
(24)
1=
i=+1
*+1 k-+2
st1-k+1
((za44k)
(Zaatk-1)
II
which completes the induction. Notably, the conclusion restricted chain-structured backbones. but also applies to tree- ere-structured backbones. The corresponding k- variaat and (k + 1) variatee local marginals are defined over k-vertex and (k+ + ve-vertex cliques, respectively, and the conclusion is built up by the Hammersley-Clifee -llifforrd Theorem.
A.2 Proof of Theorem 2.2
Lemma A.3. Let a, b. distinct integers {1,2,...,ht and let be a subset of (1, ,....,KKVV,, (a, b, c}. For a correlation matrix R E [[++1)x(++1), we denote D(L) as the deter- minant of the sub-matrix R[L]- A R then
D((a.b.e)
2 Tablet = D(laccuL)D((b.euu))
(25)
Proof. If a, b, are indices not in L, then define
1
Tobi v.
YaclF.
R R|a, 6, c|L &
Taoll Naclt
Tbe
(26)
Ybelk
23

---

and define R(aa R(a, L] R as principal 2 x sub-matrices of R(a, b, c| Since that
TabloL = Table. Yhol Pocl..
(27)
V Vodl Tbolr.
then
(1 - Tae (1 - -bc) ( - (Table - Yacl l
22 Tablet =
(1 Naoi 1. 1.)(1 Voclr.)
22
Loo - TbelL 710i1 + Vico (oell TablL
(28)
(1 Tacl
-acc (1 (oe)e)
det R[abe| L]
det R|ac|L] det R[bc|L]
If I. = 20, then the above becomes
det R [[abc]
D(fa, b()a.)c (((a,c))D((O.c)
(29)
det R.ac] det R. [(be]
since by definition D({c))
Otherwise for L + Ø, let (Zis 2j, zR, 2r.) a mean zero normal random vector with correlation matrix R|{a, b, U L] and unit variances. Let Vabe
diag(Varlza [Var[z kr), Var[zs kz., Var[ze|zg.) so that 2///R[[bbe LVV is the covariance matrix of (Za: Zb. IZz. Since the determinant of positive definite matrix can be decomposed the mul- tiplication of determinant of its principal sub-matrix and determinant of the corresponding Schur complement, then
Me/RRabbl R[abb] L)b1 det Ri b, D(fa,b,ctuz) D(L)
dVt(V1/2
(30)
det R
so that
D(la,bee
det R(abe] L]] = D(L) Varlz. (2) Var [ar(z)/zz) Var 2 Zy
(31)
Similarly,
u U UL)
det R|ac| det L[be| (2 Var 2o Øy Var[2s Var[ze|zz]
(32)
Hence, det L[bb[L]
Varizelzu. D(fa, UL) U))
(33)
det R[ac| L] det R L]
By another application the determinant decomposition, D(L) Var Varlz [2.]] = D((c) L), which completes the proof.
O
Lemma A.4. For correlation matrix R (+++))((++1) with conditional correlations Tig as defined by Eq. 8 with li 5l S k. its determinant is given by
& k+1-1
det [
II F=1
(1 2 (1-
(34)
Proof. The result is known for l: 1. To start the induction proof, we first prove it for k = 2. As special case of conditional correlation as defined by Eq. (
713 712723
713/2
(35)
V1 - 2 /1 1---- 7/3
so that
2
7/3 2/3 + 22/127237/13
det(R)
7/3/2
(36)
(1 - 72 [1 228
(1 (1 723)
24

---

Hence de) (1 - -1- -1732) We proceed by induction and go from &s to k + 1.
711/2)
The induction hypothesis gives
det R((1 |a)l D(11, ,k), II
(1 7
(37)
and we want to show that det R for size k x k is
b k+1-1
II
II (1 (.Ht) D(LL, kb)(1 - -22,642)
(38)
=11
(1 (1- (1-12-16411) ... 71
By Lemma A.3 this is:
D(11, (1) 1)) D 1,, (})
D(fk ((k)) 13) D(k 2, - 1, 11)D((k-
x
D(i 2,
1,,,,, det (2, .... kk]
x...X D(11. , D(k,, ls + 1))
(39)
D 1-
- 1 .. ;A]) 13)
=D(L, ( 13) II 1=
Dik+1-1,... 1,,1--..-. (()(,,,,))(((k
D((H)) D(L1, [K[1,...,,II)
D(LL, (k))D((k,++1)) 1)) D(LL, ,kk] D((kk+1])
det R X D({+)) det R.
Theorem A.5. By writing the correlation matrix R as the following partitioned form
Batk-l
Tithh... 741.44
76,+1
R
(40)
the upper-left sub-matrix [R. and the riger-right sub-matria [R (+) frames are both positive definite, also if 16 < and we set Tii44 M(Tf,e4x) then R( will be positive definite.
2: 1,2 in the dotted
Proof. Without loss of generality, we use the indices 21,2. ,k replace the original indices {i, i 1, + k} of the correlation R E ((++1)x(++1) giving the correlation matrix
7ik
712
12k
72
(41)
72k
Nk.4+1
71.++1
5.72.++1
To prove that is positive definite, we only need to show that det R > O, given the sub-matrix
[R] is positive definite. Since that both the sub-matrices [R] 1:k and [R] 1,2.44 positive definite, then the corresponding marginal distributions g(Z1.)) and g(22x-1) must be valid,
define conditional correlations Yiliz, for any - -l < k. satisfying that hajj, < 1. Combing the provided (I: 1)-th order conditional 75, can leverage Lemma A.4 to
25

---

compute the determinant of the correlation matrix R as follows
k+1-f
det R
II
II (1- 2
II
3 Tysts)s
(42)
1=1
which is guaranteed to be positive, given that 2/2 also lies in (-1, 1). Therefore, det R > and then R. is positive definite.
A.3 Proof of Corollary
2.3 and 2.4
Corollary A.6. If all correlation coefficients in Tı = 1os4t\er and T, = (igit NN-t for
t = 2.3.... ... . lie in the interval (-1, 1), then we can use them to construct a (k + 1) x (h (k+1) +1) correlation matrix R with k-order dependency.
Proof. The result is known for = 1, since the determinant det R( = 1 - 71,41 > 0. To start the induction proof, we first prove for 2, which show that
-
7,,42
R (1)
71,+1
1
741,12
(43)
m1,+4
744,442
is positive definite with 7,,42 M(1((+2) and l7,442l
1. In this case, the 2 upper-left
submatrix 112,122 ( positive definite, since its determinant
det
71i+1 = - 7,1+1 O,
(44)
71,+11
for |1,,41l 1. And similarly, the lower-right submatrix [R is also positive definite. Leveraging the conclusion of Theorem 2.2oorrmm222anddd the condition that 7,442 = M(n(,442 with |7 <
guarantee R (6) to be positive definite. We proceed by induction and go from 1 to k
where we want to show that the correlation matrix
Bitk- 74144
71++1
M411H
R
(45)
W 7ii+h
is positive definite. The induction hypothesis gives us the upper-left submatrix [R() 1isk 1sk and lower-
right submatrix [R| (x) 12:4+1,2:4++ are both positive definite, by ensuring that correlation coefficients r' and rL for & E (2, ... k - lie in the interval (-1, 1). By further ensuring that h MM
with | I < 1, then we can leverage Theorem 2.2 to guarantee the positive definiteness of R which completes the proof.
Corollary A.7. If the first-order correlations P2 and higher-order conditional correlations I for 2, 2,...... are built upon tre-sstructured backbone, and all correlation parameters lie in the interval (-1,1), then we can use them to construct (k + 1) x (k + 11) correlation matrix R!()) with k-order dependency structure.
Proof. This corollary of the tree-structured backbone can be similarly proved above by induction. Notice that every (k + 1)-vertex clique C E C++1 can be decomposed into two k-vertex cliques C and C2 such that C C, U C2 and [C₂ C2l = k 1. So the correlation matrix with respect to C can be partitioned as two submatrices corresponding to C₁ and C2, respectively. By ensuring their positive definiteness and letting the k-order conditional correlation similarly guarantee the correlation matrix of C to be positive definite. Therefore, we can perform induction by starting from & 2, and sequentially expand the correlation matrix by introducing higher-order conditional correlations and ensuring them to lie in the interval (-1, 1). As the induction
in the interval (-1, 1), we can
26

---

completes, we can use these correlation coefficients to construct a (k: + 1) x (k + 1) correlation matrix R(O).
B Procedure of Constructing the Correlation Matrix
Algorithm 1 Algorithm of constructing the correlation matrix R
Input: Conditional parameters T (rr,41 J 1 I2 2 1.. .xx &i 1
NN-K
Output: Full correlation matrix R of size N x N
1:
function CORRELATION BATOOOAAAARIIIIINOTNNCCCOI
2: R I
Identity matrix
3:
k+-1
> Starting from the first-order
4:
for i to N do
5:
R + 7,.+1
6:
R + 1, f 7.+11
7:
end for
8:
for < 2 to do
b LLopp through higher orders
9:
for i < to N 3: do
10:
7 &- inverse_ pnnerse_coodiiionna(((... :i i k))
Inverting Eq. (8)
11:
R i + kl + 7ii+
R + i] + 7.1+
13:
end for end for
15: return R
end function
C Evidence Lower Bound
The evidence lower bound of our proposed method is given by
L( Φ, x) = Logro(( (x, 2) + (z),
(46)
where denotes the re-parameterized latent variables. The first term above can be directly computed by
log 2) = Σ log (x; + logp(2;),
(47)
where 2; is the reparameterization for latent variable 2i, i = 1, .N. And the entropy of the posterior
N-k-1
go (Zii44))
(z)
II
yo (o(2iit)-))
II
(48)
lΦ (Ziz ())
with k-order dependency structure can be factorized entropy terms with respect to k-variate and (k + 1)-variate local marginals
1++1
N Σ
H. [z)]
H[go (Z
(H(a))-((aaaaii)) Hl H[g[
1=1
N-
Σ H(gg(Zz(++))
Σ
(49)
where the entropy of each normally distributed local marginal can be directly computed by its mean and covariance.
27

---

D Experimental Details
D.1 Datasets
The datasets used in the time series anomaly detection task are the followings:
SMAP (Soil Moisture Active Passive): NASA's Soil Moisture Active Passive mission [25 aims to measure global soil moisture and freeze/thaw states to enhance understanding of Earth's water, energy, and carbon eycles. The SMAP dataset comprises multivariate time series telemetry data collected from the SMAP satellite, including a training and a testing subsets. It includes expert-labeled anomalies in testing subsets, making it suitable for benchmarking time series anomaly detection algorithms.
MSL (Mars Science Laboratory): Originates from NASA's Mars Science Laboratory mission (25]. featuring the Curiosity rover, explores Mars' surface to assess its habitability. The MSL dataset contains multivariate time series telemetry data from the Curiosity rover, with expert annotations identifying anomalous events in the testing subsets.
Collected by researchers from large Internet company [55]. SMD comprises a 5-week-loog collection of multivariate time series data from 28 server machines, each monitored by 38 sensors capturing metrics like CPU usage. memory, and network throughput. The dataset includes labeled anomalies, facilitating supervised learning approaches. Due to the high degree of similarity in temporal characteristics across servers, we conducted experiments solely on machine 1-1 for simplicity.
SMD (Server Machine Dataset):
The datasets used in the time series forecasting task are the followings:
ETT (Electricity Transformer Temperature): This dataset includes the target variable "oil temperature" along with six power load features [69]. It is recorded at two different frequencies: hourly (i.e., ETThl and ETTh2) and every 15 minutes (i.e., ETTml and ETTm2), spanning a period of two years.
Sourced from the UCI Machine Learning Repository and preprocessed follow- ing [33), this dataset contains hourly electricity consumption (in kWh) for 321 clients from 2012 to 2014.
Electricity:
Exchange: This dataset comprises daily exchange rates for eight countries, collected from 1990 to 2016 [44].
Weather Includes 21 metcorological indicators (e.g., temperature, humidity), recorded every 10 minutes throughout the year 2020.
The datasets utilized in the constrained clustering task are as follows:
MNIST: A widely used benchmark dataset containing 70,000 grayscale images of handwrit-
ten digits. Each image is represented as 784-dimensional vector by flattening the original 28x28 pixel grid [35].
Fashion MNIST: A collection of Zalando's ' fashion article images (62), this dataset includes a training set of 60,000 images and a test set of 10,000 images.
Reuters: Contains 810,000 English news articles (36]. Following the preprocessing method of DEC [64], we select four root categoriss-corporateee corporate/innuustiall government/social, mar- kets, and economics and exclude documents with multiple labels. The resulting dataset contains 685,071 articles, each represented using tf-idf features over the top 2,000 words. A random subset of 10,000 documents used for experiments.
STL-10: Composed of 96x96 color images across 10 object classes, with 13,000 labeled samples [16]. For feature extraction, we apply a ResNet-50 model as done in VaDE (27]
https://www www.:wbgc-jee... jena.jena... meona.mpg.... wttp..pp.d//wwtterrrrrrr...
28

---

Table 5: Detailed information of datasets used in time series anomaly detection and forecasting tasks.
Tasks
Dataset
Dim
Size (Train, Validation, Test)
Domain
Forecasting
ETTml
7
(34465, 11521, 11521)
Electricity
ETTh1
7
(8545, 2881, 2881)
Electricity
Electricity
321
(18317, 2633, 5261)
Electricity
Weather
21
(36792, 5271, 10540)
Weather
Exchange
8
(5120, 665, 1422)
Exchange rate
Anomaly Detection
SMD
38
(566724, 141681, 708420)
Server Machine
MSL
55
(44653, 11664, 73729)
Spacecraft
SMAP
25
(108146, 27037, 427617)
Spacecraft
D.2 Further Experiments
We also run our model under univariate forecasting settings, where only single feature is considered in each dataset. The experimental results in Table shows that our method outperforms other fundamental time series modeling techniques. The superior capability of our method in capturing temporal dependencies is more pronounced in this setting, as all models are restricted to fully exploiting temporal correlations without leveraging inter-channel information.
Table 6: Univariate time series forecasting comparisons. Best performance highlighted in bold font and the second best results are underlined.
VRAE
Informer
Autoformer
TCN
Ours
Method
1-order
3-order
10-arder
MSE
MAE
MSE
MAE
MSE
MAE
MSE
MAE
MSE
MAE
MSE
MAE
MSE
MAE
24
0.059
0.215
0.098
0.247
0.057
0.189
0.104
0.254
0.054
0.178
0.038
0.149
0.032
0.127
48
0.097
0.279
0.158
0.319
0.070
0.207
0.206
0.366
0.087
0.229
0.061
0.187
0.052
0.163
168
0.402
0.346
0.260
0.462
0.586
0.161
0.316
0.131
0.278
0.088
0.212
336
0.187
0.400
0.222
0.387
0.119
0.281
0.422
0.564
0.170
0.333
0.149
0.303
0.105
0.240
720
0.244
0.471
0,269
0.435
0.109
0.264
0.438
0.578
0.221
0,392
0.172
0,336
0.139
0.285
24
0.021
0.122
0.030
0.137
0,022
0.115
0.027
0.127
0.018
0.101
0.015
0.091
0.012
0.075
ETTml
0.039
0.172
0.203
0.032
0.138
0.154
0.041
0.154
0.027
0.123
0.022
0.105
168
0,060
0.217
0,194
0.372
0.045
0.168
0.097
0,246
0.052
0.173
0.043
0,158
0.034
0.129
336
0.143
0.344
0.401
0.554
0.071
0.207
0.305
0.455
0.131
0.276
0.091
0.229
0.073
0.192
720
0.428
0.644
0.254
0.445
0.576
0.134
0.287
0.135
0,282
0.099
0.227
24
0.370
0.459
0.251
0.275
0.290
0.411
0.243
0.367
0.252
0.278
0.247
0.285
0.166
Electricity
48
0.459
0.519
0.346
0.339
0.310
0.408
0,283
0.397
0.301
0.309
0.298
0,318
0.202
0.277
168
0.547
0,575
0.544
0.424
0.435
0.490
0.357
0.449
0.413
0.384
0.408
0.386
0.270
0.323
336
0.682
0.660
0.713
0.512
0.646
0.606
0.355
0.446
0.551
0,468
0.537
0.468
0.339
0.369
720
0.889
0.790
1.182
0.806
0.609
0.587
0.387
0.477
0.862
0.650
0.812
0.628
0.454
0.448
Average
0,280
0.404
0.345
0.400
0.202
0.306
0.278
0.403
0.230
0.302
0.211
0.281
0.139
0.228
D.3
Implementation Details
Time Series Anomaly Detection We set the input sequence length to 100 and use GRU and dense layers with 500 hidden units each. The latent dimension is fixed at 3. Models are trained with batch size of 50 for up to 20 epochs using early stopping. Optimization is performed using the Adam optimizer with an initial learning rate of 10 3. L2 regularization with a coefficient of 10 is applied to all layers. During training. 30% of the data is reserved for validation.
Time Series Forecasting We adopt a single-layer fully connected network as the feedforward predictor. The latent representation dimension is set to 128. The model is trained using the Adam optimizer with an initial learning rate of 10-3 decayed by a factor of 0.95 after each epoch. Early stopping is applied within 10 epochs to prevent overfitting.
Constrained Clustering. To ensure a fair comparison with baseline methods, we adopt the same encoder-decoder fced-forward architecture: four fully connected layers with sizes 500, 500, 2000,
29

---

Table 7: Hyperparameters setting of constrained clustering task.
MNIST
IMNIST
Reuters
STL-10
Batch size
256
256
256
256
Epochs
1000
500
500
500
Learning rate
0.001
0.001
0.001
0.001
Decay
0.9
0.9
0.9
0.9
Epochs decay
20
20
20
20
and D units, respectively, where D 10 unless otherwise specified. For all VAE-based baselines and our proposed methods built VAE backbones, apply 10 epochs of pretraining. For DEC-based baselines, we follow their standard training procedure, including 50 epochs of layer-wise pretraining and 100 epochs of fine-tuning. Each dataset split into training and test sets; model training is conducted on the training split, while all reported results are evaluated on the test split. Pairwise constraints randomly generated within the training set: must-link is assigned if two sampled instances share the same label, and cannot-link otherwise. To ensure consistent training conditions across methods, we uniformly set the absolute constraint strength [[A] 10 and sample 6000 such constraints for all datasets. Following DC-GMM, we use the same set of hyperparameters across all four datasets, detailed in Table 7 All models are trained with an initial learning rate of 0.001, which decays by factor of 0.9 every 20 epochs.
D.4 Resource Usage
Experiments were conducted an internal computing cluster. Each experiment configuration used one NVIDIA GPU (either a 2080TI or 3090TI), 16 CPUs and a total of 24GB of memory.
30


Reviews:
[{"review_id": "iQJKsnrtfb", "reviewer": "Reviewer_yDoz", "summary": "The authors propose a more accurate variational inference approximation that takes instance level correlation into account. They do this by approximating the N by N covariance matrix with a K by K matrix that takes into account the correlation between $i$th instance upto the $i+k$th instance. The correlation coefficients of this K by K matrix are parametrised through a neural network. The authors then introduce a method to ensure that the resultant covariance matrix is PD as well as an extension to a tree based correlation structure. The authors show improved performance in time series and constrained clustering tasks.", "questions": "- What is the intuition behind the K level correlation modelling? Is using a value of 1 effectively a Markov assumption in the posterior?\n- At what level of K should we recover the same correlation structure as the baselines? In Table 1 I would have expected the 1 order method to perform as well as an LSTM-VAE that models the same correlation (apart form the fact that the authors use a GRU). \n- Why is it necessary to directly parametrise the correlation coefficients $R$? If terms of the covariance matrix are parametrised won't this allow for parametrising the cholesky of the covariance matrix which will ensure that the covariance is PD by construction?\n- What kind of neural network is $f_{\\phi}$? Are the different $x_i$ concatenated across the feature dimension or is attention applied across the $K$ dimensions?\n- Is there any relation to the inducing points framework used in sparse GPs? The problem being solved seems similar, here an N by N matrix is too complex to invert, so it is approximated by a K by K matrix with K < N.\n- What is the cost of the sampling procedure laid out in L121?", "limitations": "Yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 2, "strengths_and_weaknesses": "**Strengths:**\n- The method is interesting and sound.\n- The experiment are fairly wide ranging and over multiple baselines.\n\n**Weaknesses:**\n- The paper could be a lot clearer:\n     - What is the exact cost of inference with P?\n     - The intuition behind modelling only K correlations (apart from better computational cost) is not clear. What are you losing here? There should also be a discussion on the intuition for choosing this parameter.\n    - The baselines are not sufficiently described. Which of these methods actually model instance level correlation? If they do, what is the difference between the proposed method?", "was_revised": true, "final_justification": "The authors have addressed my concerns and I don't have any other major concerns.", "rebuttal": "We thank the reviewer for the comments and suggestions. Below we respond to each point raised by the reviewer.\n\n**Q1: What is the exact cost of inference with $P$?**\n\nA1: The exact cost of inference with our proposed k-order precision matrix $P$ involves three parts: (i) the cost of neural network evaluations for re-parameterizing correlation coefficients, (ii) the cost of sampling from the variational posterior, and (iii) entropy calculation.\n\nIn our method, to define a posterior with k-order correlation over $N$ latent variables, we need to specify exactly $(N–1)$ first-order, $(N–2)$ second-order, and so on, up to $(N–k)$ k-order correlations, yielding a total of $(N-1)+(N-2)+\\cdots+(N-k)=k(2N-k-1)/2$ correlation coefficients. In our method, each coefficient is parameterized by the output of a re-parameterization network $f_{\\phi}(\\cdot)$. Therefore, to re-parameterize these coefficients, we need to run the network $f_{\\phi}(\\cdot)$ for $\\mathcal{O}(kN)$ times.\n\nThe sampling and entropy calculations involve operations like inversion and determinant computation on $k \\times k$ sub-matrices, incurring a cost of $\\mathcal{O}(k^3)$ FLOPs. Considering that k is typically much smaller than $N$ and the complexity of evaluating neural networks, the cost of these operations is negligible compared to the cost of neural network evaluations. \n\nTherefore, the total cost approximately amounts to the cost of evaluating $\\mathcal{O}(kN)$ times of the neural network $f_{\\phi}$ per epoch, which is approximately k times of the cost of mean-field amortized VI methods.\n\n**Q2: What is the intuition behind k-level correlation modeling? How should $k$ be chosen? What are you losing here? Does k=1 imply a Markov assumption in the posterior?**\n\nA2: The choice of k is a trade-off between model expressiveness and computational cost. Generally, as k increases, the model's performance consistently improves, as demonstrated in our experimental results. However, as seen from Table 1 below, the performance gains diminish as the order k (e.g. 50, 100) goes higher. Meanwhile, the computational cost always scales linearly with k according to our analysis in answer A1. To balance the gains and cost, we set k to moderate values (up to 10) in our main experiments. \n\nBy setting k to a moderate value (e.g. 10), we can only model correlation up to 10-th order, losing the ability to model higher-order correlations. But as observed from Table 1, the gains become increasely weak as the order goes higher. Therefore, setting k to a moderate value is a reasonable choice.\n\nYes, when k=1, our variational posterior is reduced to the first-order Markov random field distribution.\n\n---\nTable 1: Time series anomaly detection on SMAP dataset.\n|Methods|Mean-field|k=1|k=3|k=10|k=50|k=100|\n|-|-|-|-|-|-|-|\n|Runtime(s/epoch)|1.00|2.44|8.51|27.60|142.58|294.53|\n|F1|0.7774|0.8411|0.8552|0.8636|0.8711|0.8755|\n|ELBO|-109.2182|-97.6057|-95.2314|-92.2948|-90.4291|-89.9577|\n\n**Q3: Which baselines actually model instance-level correlation, how do they differ from the proposed method?**\n\nA3: Due to the difficulties of modeling instance-level correlation, existing works mostly only consider the modeling of correlation within one instance. In our baselines, only OmniAnomaly and TreeVI explicitly model instance-level correlations in their posterior distributions.\n\nOmniAnomaly imposes a first-order Markov assumption among latent variables and employs an MLP to parameterize the one-step conditional Gaussian, where $(\\mu_{z_t},\\sigma_{z_t})=\\text{MLP}([z_{t-1},e_t])$. This design makes it essentially a recurrent model, just like an RNN learning a transition function for its hidden state. As a result, it is inherently limited to modeling only first-order correlation. More importantly, OmniAnomaly can only be applied to model chain-structured correlation due to its inherent RNN-like structure, making it unable to model correlation with more complicated structure like trees. In contrast, our method not only can model higher-order correlation, but also can be applied to model correlations with tree backbones.\n\nThe key difference between TreeVI and our method is that TreeVI is inherently limited to modeling only first-order correlations, while our method can capture higher-order correlations. This limitation of TreeVI arises from its reliance on an acyclic tree structure to construct its correlation matrix. While this is sufficient for simple pairwise relationships, attempting to model higher-order correlations inherently introduces loops into the underlying correlation structure. The construction of TreeVI depends on the acyclic property of its backbone and is no longer valid when these loops exist. Moreover, simply modeling higher-order correlation coefficients within the framework of TreeVI does not guarantee the correlation matrix to be positive definite. So even though TreeVI can capture instance-level correlations, it cannot be easily extended to model higher-order correlations.\n\n**Q4: At what level of k should the proposed method recover the same correlation structure and performance as baselines? In Table 1, I would have expected the 1 order method to perform as well as LSTM-VAE.**\n\nA4: At k=1, our method is reduced to TreeVI, with both methods modeling first-order instance-level correlation in the same way. But for other baselines, although some of them also attempt to model first-order instance-level correlation, they are fundamentally different from us in terms of correlation modelling methods and the ability of generalizing to more complicated correlation structures.\n\nAs for the performance gap between LSTM-VAE and our method, the reason is explained below. Although LSTM-VAE also attempts to capture instance-level correlation, it is achieved via the deterministic dependency among the hidden states in RNN. Given the hidden states' values, the latent variables $z_i$ of instances at different time steps are independent. By contrast, our method not only captures the deterministic dependency, but also models the correlation among the stochastic variations of latent variables $z_i$ corresponding to different instances. That is,  given the hidden states' values, in our method, the latent variables $z_i$ are still correlated. Thus, our method is able to capture more accurate and complicated correlation among data. That explains why our method  performs better than LSTM-VAE even under the case of k=1.\n\n**Q5: Why is it necessary to directly parameterize the correlation coefficients $R$? Why not parameterize the Cholesky decomposition of the covariance matrix to ensure positive definiteness by construction?**\n\nA5: We first want to emphasize that this work aims to capture instance-level correlation, under which the dimension of correlation matrix is $N\\times N$, with $N$ denoting the number of instances. Parameterizing the Cholesky factor $\\mathbf{L}$ of a correlation matrix $\\mathbf{R} = \\mathbf{L} \\mathbf{L}^\\top$ does ensure its positive definiteness. However, we have to re-parameterize as many as $N(N+1)/2$ elements, i.e., all elements from the lower-triangular positions of $\\mathbf{L}$ need to be re-parameterized. That means we need to evaluate the neural network $f_\\phi(\\cdot,\\cdot)$ by $\\mathcal{O}(N^2)$ times for every epoch, which is computationally unacceptable, especially considering that $N$ could be very large in practice.\n\nTo reduce the cost, a common practice is to restrict $\\mathbf{L}$ to be low-rank, which reduces the number of parameters to be re-parameterized from $\\mathcal{O}(N^2)$ to $\\mathcal{O}(kN)$ for a low rank $k$. However, each element in $\\mathbf{L}$ does not explicitly correspond to the correlation between any two specific data instances. Therefore, to parameterize each element in $\\mathbf{L}$, we have to feed all $N$ data instances into the neural network $f_\\phi$. Obviously, this will result in the network to be huge, and thereby incurring unacceptable computation and storage burden.\n\nIn contrast, although our method also need to re-parameterize $\\mathcal{O}(kN)$ parameters, in our method, the parameter is the correlation coefficient. For each correlation coefficient, such as $\\gamma_{i,i+t}^c$, it directly represents the correlation between a specific pair of instances $(x_i,x_{i+t})$. Therefore, we only need to feed the pair of instances into the neural network $f_\\phi$, making the network to be small.\n\n**Q6: What kind of neural network is $f_{\\phi}$?**\n\nA6: The function $f_{\\phi}$ is implemented as a standard MLP. Inputs $x_i$ and $x_{i+t}$ are concatenated across the feature dimension to form a single input, which is then processed by the MLP with a tanh activation function to constrain the output to the range of (-1,1).\n\n**Q7: Any relation to inducing points used in sparse GPs, given that both aim to approximate large matrices?**\n\nA7: While both our method and the inducing points method aim for scalability, the inducing points framework cannot handle instance-level correlation modeling effectively. The inducing points method works by a low-rank approximation based on a small set of $M$ representative points. Since these points must summarize the properties of the entire dataset, the network used to learn their optimal locations needs to see all $N$ data instances to derive this global summary. This results in a network with an $\\mathcal{O}(N)$ input size, which is computationally intractable for large datasets. Therefore, low-rank methods like inducing points are not usable in our setting. In contrast, in our method, the re-parameterization neural network used to compute each correlation coefficient only requires a pair of instances as input.\n\n**Q8: Computational cost of the sampling procedure laid out in line 121?**\n\nA8: The cost is composed of two parts: $\\mathcal{O}(kN)$ times of neural network evaluations for re-parameterizing correlation coefficients, and $\\mathcal{O}(k^3)$ FLOPs for matrix inversion. For more details, please refer to A1."}, {"review_id": "ghDZmrk5AZ", "reviewer": "Reviewer_oEn6", "summary": "The paper focusses on approximate Bayesian inference in correlated latent variable models using variational inference. The paper proposes a Gaussian variational family with (chain or tree) structured  precision matrices. \n\nThe paper provide a product decomposition theorem (Theorem 2.1) for the density of a multivariate Gaussian, where the precision matrix exhibits k-order connectivity structure. The theorem also works \"the other way\", giving conditions for when the set of marginal form a valid joint density. The theorem enables efficient sampling from the joint distribution in a sequential manner via the set of marginals.  \n\nFor efficiency, an amortized solution is proposed, where the required correlation coefficents of the variatonal distribution are reparametrized via a neural network. Specifically, the correlation matrices are re-parametrized through a combination of (constrained) first order and higher order correlation coefficient to ensure positive definiteness and proves that the resulting correlation matrix is valid (Theorem 2.2 and Corollary 2.3).\n\nThe paper concludes with an empirical evaluation of the proposed method in the context of time series anomaly detection, forecasting and constrained clustering. The results for anomaly detection shows that the proposed method outperforms the included baseline in terms F1 (better outlier detection) and ELBO (better representation of the exact posterior as measured by the KL divergence) and similar strong performance for the forecasting and clustering experiments", "questions": "I believe the paper could be stronger if  the authors could elaborate on the following points:\n\n- The results generally show improved performance as a function of the order of the structured precision matrix. It is therefore natural to wonder whether the performance would be even better if the author went beyond 10? Did the authors experiment with approximation order beyond 10?\n  \n- Moreover, it would be interesting and informative to understand the effect of increasing the order of the approximation, all the way from mean-field to a fully flexible Gaussian (possible on a smaller dataset) and possibly also comparing to other variational families, e.g. low rank Gaussians and flows.\n  \n- It would also be interesting to investigate the trade in computational complexity when the order is increased from mean-field to a full Gaussian. What happens to the run time? When will it become more efficient to simply use a full Gaussian rather than the structured approximation. \n  \n-  There are no error bars/uncertainties associated with the metrics. How robust are the results?\n  \n- What is the experimental set-up for the baselines? Did the authors do an effort to make them appear as strong as possible or did they use them blindly?\n  \n- What is the likelihood of the VAE? A Gaussian? How many posterior samples were used during training and evaluation?", "limitations": "yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "The paper is well-structured, well-written and easy to follow. The proposed method for bridging the gap between a mean-field posterior and a general Gaussian is interesting. There is a lot of prior work on structured precision matrices in general, but I believe that the proposed idea is novel in the context of variational inference (at least to the best of my knowledge) and will be of interest to the community. \n\nGenerally, the experimental design seems sound and shows a clear benefit relative to the included baselines. However, I am not that familiar with the set of included baselines and therefore, it is no clear to me whether they are strong baseline or easy to beat. Moreover, it is unclear how the baseline methods are configured, to what degree they are tuned etc.\n\nMoreover, some of the implementation details are unclear. For example, what is the likelihood/loss in the VAE? It is a standard Gaussian? How many posterior samples is used for 1) training and 2) evaluation? Finally, the authors do not report uncertainty for any of the evaluation metrics, i.e. no error bars etc.", "was_revised": true, "final_justification": "The detailed response from the authors addressed most of my questions/concerns, and therefore, I believe the paper should be accepted, and hence, my final score is 5.", "rebuttal": "We thank the reviewer for the comments and suggestions. Below we respond to each point raised by the reviewer.\n\n**Q1: Are the chosen baselines strong or merely easy to beat? How were they configured and tuned during experimentation?**\n\nA1: Our chosen baselines represent strong and state-of-the-art methods in their respective domains, selected to ensure a fair evaluation of our model. For the time series modeling tasks, we compared our model against both leading generative methods and prominent forecasting models. For constrained clustering, we benchmarked against state-of-the-art algorithms. \n\nTo ensure a fair comparison, we configured and tuned all methods to achieve their best possible performance. Hyperparameters of baseline methods are set consistent with our method as far as possible.\n\n**Q2: What is the specific likelihood or loss function used in the Variational Autoencoder (VAE) component? And how many posterior samples were utilized during both the training and evaluation phases of the model?**\n\nA2: The specific likelihood function $p( \\mathbf{X} |\\mathbf{Z})$ used in our VAE is task-dependent. For the time series tasks, the likelihood is a Gaussian distribution with parameters produced by a GRU-based decoder. For constrained clustering, the likelihood is a normal distribution conditioned on a cluster-dependent latent variable. Across all applications, the loss function is the negative of the ELBO. For all experiments, we utilized a single sample from the posterior distribution during both training and evaluation.\n\n**Q3: Why are there no error bars or uncertainties reported for the evaluation metrics, and how robust are the presented results without them?**\n\nA3: We would like to clarify that we did report error bars in Table 3 across 10 independent runs. For Tables 1 and 2, we omitted the error bars in the main paper due to space constraints, and we will include them in the appendix in our revised manuscript. To address your concern directly, we present the F1-Scores and ELBOs for the time series anomaly detection task with their corresponding standard deviations from 10 runs below:\n\n|      |      |        k=1       |        k=3       |       k=10       |\n| - | - | - | - | - |\n| SMAP |  F1  |   0.8411±0.0078  |   0.8552±0.0070  |   0.8636±0.0065  |\n|      | ELBO |  -97.6057±4.4278 |  -95.2314±3.4668 |  -92.2948±3.835  |\n|  MSL |  F1  |   0.8883±0.0070  |   0.8940±0.0067  |   0.9145±0.0064  |\n|      | ELBO | -165.5004±3.9339 | -157.2134±3.7234 | -134.0815±3.7679 |\n|  SMD |  F1  |   0.8901±0.0084  |   0.9153±0.0081  |   0.9284±0.0087  |\n|      | ELBO |  -69.5278±1.3356 |  -67.4001±1.5207 |  -65.0345±1.2238 |\n\n**Q4: Were experiments conducted with approximation orders exceeding 10? If so, did performance continue to improve at higher orders?**\n\nA4: Yes, we have conducted additional experiments with orders k exceeding 10, specifically for k=50 and k=100, to explore the performance limits of our model, as shown in Tables 1-3 below.\n\n---\n\nTable 1: Time series anomaly detection performance on SMAP dataset.\n| Methods  | Mean-field | k=1 | k=3 | k=10 | k=50 | k=100 |\n| - | - | - | - | - | - | - |\n| Runtime (s/epoch) | 1.00 | 2.44 | 8.51 | 27.60 | 142.58 | 294.53 |\n| F1     | 0.7774 | 0.8411 | 0.8552 | 0.8636 | 0.8711 | 0.8755 |\n| ELBO   | -109.2182 | -97.6057 | -95.2314 | -92.2948 | -90.4291 | -89.9577 |\n\n---\n\nTable 2: Time series forecasting performance on ETTh1 dataset with horizon 24.\n| Methods  | Mean-field | k=1 | k=3 | k=10 | k=50 | k=100 |\n| - | - | - | - | - | - | - |\n| Runtime (s/epoch) | 4.80 | 11.45 | 36.94 | 126.44 | 675.75 | 1340.22 |\n| MSE | 0.739 | 0.664 | 0.543 | 0.363 | 0.348 | 0.333 |\n| MAE | 0.716 | 0.570 | 0.505 | 0.376 | 0.362 | 0.352 |\n\n---\n\nTable 3: Constrained clustering performance on MNIST dataset.\n| Methods  | Mean-field | TreeVI | k=1 | k=3 | k=5 | k=10 | k=50 | k=100 |\n| - | - | - | - | - | - | - | - | - |\n| Runtime (s/epoch) | 0.25 | 0.52 | 0.59 | 1.84 | 3.09 | 6.22 | 29.53 | 60.59 |\n| ACC (\\%)     | 96.50 | 97.43 | 97.55 | 98.12 | 98.39 | 98.52 | 98.62 | 98.69 |\n| NMI (\\%)     | 91.37 | 93.14 | 93.44 | 93.80 | 94.17 | 94.55 | 94.63 | 94.85 |\n| ARI (\\%)     | 92.54 | 93.73 | 93.89 | 94.89 | 95.24 | 95.65 | 95.85 | 96.09 |\n\n---\n\nThe results show that performance does continue to improve as the order k increases beyond 10, which further validates our central claim that capturing higher-order dependencies is beneficial.\n\nHowever, we observed that the marginal performance gains diminish as k becomes excessively large, and concurrently the computational runtime increases linearly with the order k. We focused on k up to 10 in the main paper because this range demonstrates the advantages of our approach while maintaining a practical and efficient balance between performance and computational tractability.\n\n**Q5: What is the effect of systematically increasing the approximation order from a mean-field approach all the way to a fully flexible Gaussian (potentially on a smaller dataset)? How does the proposed method compare to other variational families, such as low-rank Gaussians and normalizing flows?**\n\nA5: Increasing the approximation order k consistently improves model performance by capturing higher-order correlations, as seen from Tables 1-3. To capture all correlations, k would need to equal $N-1$, making our model equivalent to a fully flexible Gaussian. However, the computational cost scales linearly with k, with the runtime being approximately k times that of a mean-field method. This makes a fully flexible Gaussian computationally intractable for instance-level correlation modeling, even on smaller datasets.\n\nSince other variational families like low-rank Gaussians and normalizing flows are not designed for instance-level modeling, we apply our method to model dimension-level correlations within the latent space of a VAE to provide a fair comparison. We configured a VAE with 30 latent dimensions on MNIST and report the ELBO on the test set below.\n\n| Method              | Test ELBO |\n|---------------------|:---------:|\n| Mean-Field Gaussian |   -84.62  |\n| Low-Rank Gaussian (rank=10)  |   -83.15  |\n| Normalizing Flow (10 steps of planar flow)   |   -82.51  |\n| Ours (k=3)          |   -83.54  |\n| Ours (k=5)          |   -82.91  |\n| Ours (k=10)         |   -82.38  |\n| Ours (k=30, Full-Rank Gaussian) |   -81.73  |\n\nThe results show that our method bridges the gap between mean-field approximations and fully expressive models. Even with a small order (k=3), our method outperforms the mean-field approach and becomes competitive with the low-rank Gaussian. As k increases to 10, our method surpasses both the low-rank Gaussian and the normalizing flow in this setup, and approaches the performance of the full-rank Gaussian, demonstrating that it offers a competitive and highly flexible alternative to other variational families.\n\n**Q6: What is the trade-off in computational complexity and runtime when the approximation order is increased from mean-field to a full Gaussian? At what point does it become more computationally efficient to simply use a full Gaussian rather than the structured approximation?**\n\nA6: The choice of the approximation order k is a direct trade-off between model expressiveness and computational cost. As k increases, our model captures higher-order correlations, leading to consistently better performance, as shown in Tables 1-3 above. However, the computational cost also scales linearly with k, with the runtime being approximately k times that of a mean-field model. Because of this trade-off between diminishing performance gains and linearly increasing cost, we selected a moderate k (up to 10) in our main experiments to strike a balance between performance and efficiency.\n\nWhen k is set to $N-1$, where $N$ is the number of instances, our method becomes equivalent to a full Gaussian. However, a full Gaussian would require evaluating $\\mathcal{O}(N^2)$ times of neural networks $f_\\phi(\\cdot)$ per epoch, which is computationally intractable for instance-level correlation modeling, especially considering that the number of instances can be very large in practice."}, {"review_id": "S1DLpTYvWf", "reviewer": "Reviewer_XNGw", "summary": "Within the context of variational inference (VI) for latent-variable models, this paper proposes a variational structure that builds in dependence between k neighbouring latent indicators. In principle, doing so will allow improved dependence (e.g. temporal) modelling within the latents, and better modelling and predictive abilities on the observed data. The paper uses a Gaussian distribution for the variational approximation, but with some conditioning and reparameterising tricks, the ability to handle this computationally is managed, with reasonable scalability. The simulations show improved performance in a number of modelling tasks compared to existing methods with weaker dependence among latents.", "questions": "I was looking for some information about the computational overheads needed for this method (given that this is effectively the problem the paper is solving). Both for the presented simulations in Section 4 (comparative times for analysis of each method), but also an understanding (empirical and general) of how these overheads scale as k increases in a given computational task. Is k=3,5,10 the limit that this method can reasonably tackle, or is k=50,100,250 possible too?\n\n\nSome trivial typos (no response required):\n\nFig 1 - should the last node be Z_N?\nl.97 \"connection structure [is] imposed\" - delete \"is\"\nl.102 \"follows\" -> \"follow\"\nl.138 \"for THE ELBO\" insert \"the\"\nl.146 \"above TO BE valid.\" insert \"to be\"\nl.298 and l.300 \"state of THE GRU\" insert \"the\"\nl.319 \"learn AN expressive\" change \"a\" to \"an\"\nl.327 Write L1, L2 as L_1, L_2, delete the extra space before the comma from \"L_2 loss ,\"\nReferences lots of missing capitalisations. E.g. \"lstm\", \"gaussian\", \"wasserstein\", \"bayesian\", \"Journal of machine learning research\" etc.\nReference [57] and [58] look suspiciously similar. Is there a publication racket going on here?", "limitations": "Yes, but as discussed above, it would be good to have some further discussion on this, as it was only mentioned very briefly.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths: The paper has seemingly found a way to build structure between latent variables, via a normal variational density, while managing to control the computational overheads for increasing dimension of dependence. This structure not only includes a k-Markov correlation structure, but this can be reworked to encompass a \"tree-structured backbone\" dependency structure.\n\nThe simulation experiments demonstrated that capturing such dependence is (of course) beneficial in a number of time series and clustering tasks.\n\nWeaknesses: It seems that there may be some complexities involved in constructing more detailed correlation structures. I didn't get too much of a sense how hard or involved this was from the paper though.\n\nIt's also not clear what the computational overheads of the proposed method are, and how this scales with increasing k.", "was_revised": true, "final_justification": "The authors have cleared up some areas that were initially unclear, provided extra simulation to demonstrate method performance in more challenging scenarios, and done a proper cost analysis of computational overheads. Overall the method is a nice contribution, and (if the revisions are made), it should be well received in the community. I have changed my score to 5 to clearly indicate support to the AC.", "rebuttal": "We thank the reviewer for the comments and suggestions. Below we respond to each point raised by the reviewer.\n\n**Q1: How difficult or involved is it to construct more detailed correlation structures within the proposed framework?**\n\nA1: Constructing more detailed correlation structures within our framework is straightforward by design. We can increase detail in two main ways: (1) by increasing the approximation order k to capture higher-order correlations, (2) by generalizing the dependency backbone from a simple chain to a tree structure to model more flexible correlations among data instances.\n\n**Q2: What are the specific computational overheads of the proposed method? How do these overheads scale as the approximation order k increases?**\n\nA2: The specific computational overheads involve three parts: (i) the overhead of neural network evaluations for re-parameterizing correlation coefficients, (ii) the overhead of sampling from the variational posterior, and (iii) entropy calculation.\n\nIn our method, to define a k-order posterior over $N$ latent variables, we need to estimate exactly $(N – 1)$ first-order, $(N – 2)$ second-order, and so on, up to $(N – k)$ k-order correlations, yielding a total of $(N-1) + (N-2) + \\cdots + (N-k) = k(2N-k-1)/2$ correlation coefficients. Each coefficient is produced by a single forward pass through our re-parameterization network $f_{\\phi}$, so this cost amount to the cost of evaluating $\\mathcal{O}(kN)$ times of the neural network $f_{\\phi}$.\n\nThe sampling and entropy calculations involve operations like inversion and determinant computation on $k \\times k$ sub-matrices, incurring a cost of $\\mathcal{O}(k^3)$ FLOPs. Considering that $k$ is typically much smaller than $N$ and the complexity of evaluating neural networks, the cost of these operations is negligible compared to the cost of neural network evaluations. \n\nTherefore, the total cost of the variational inference approximately amounts to the cost of evaluating $\\mathcal{O}(kN)$ times of the neural network $f_{\\phi}$ per epoch, which scales linearly with the approximation order k.\n\n**Q3: What are the practical limits for k? Is k=10 the reasonable maximum, or are much larger values like k=50, 100, 250 also feasible?**\n\nA3: Yes, larger values of k beyond 10 are also feasible for our method, and we have conducted experiments with approximation orders up to k=50 and k=100, as shown in Tables 1-3 below.\n\n---\n\nTable 1: Time series anomaly detection performance on SMAP dataset.\n| Methods  | Mean-field | k=1 | k=3 | k=10 | k=50 | k=100 |\n| - | - | - | - | - | - | - |\n| Runtime (s/epoch) | 1.00 | 2.44 | 8.51 | 27.60 | 142.58 | 294.53 |\n| F1     | 0.7774 | 0.8411 | 0.8552 | 0.8636 | 0.8711 | 0.8755 |\n| ELBO   | -109.2182 | -97.6057 | -95.2314 | -92.2948 | -90.4291 | -89.9577 |\n\n---\n\nTable 2: Time series forecasting performance on ETTh1 dataset with horizon 24.\n| Methods  | Mean-field | k=1 | k=3 | k=10 | k=50 | k=100 |\n| - | - | - | - | - | - | - |\n| Runtime (s/epoch) | 4.80 | 11.45 | 36.94 | 126.44 | 675.75 | 1340.22 |\n| MSE | 0.739 | 0.664 | 0.543 | 0.363 | 0.348 | 0.333 |\n| MAE | 0.716 | 0.570 | 0.505 | 0.376 | 0.362 | 0.352 |\n\n---\n\nTable 3: Constrained clustering performance on MNIST dataset.\n| Methods  | Mean-field | TreeVI | k=1 | k=3 | k=5 | k=10 | k=50 | k=100 |\n| - | - | - | - | - | - | - | - | - |\n| Runtime (s/epoch) | 0.25 | 0.52 | 0.59 | 1.84 | 3.09 | 6.22 | 29.53 | 60.59 |\n| ACC (\\%)     | 96.50 | 97.43 | 97.55 | 98.12 | 98.39 | 98.52 | 98.62 | 98.69 |\n| NMI (\\%)     | 91.37 | 93.14 | 93.44 | 93.80 | 94.17 | 94.55 | 94.63 | 94.85 |\n| ARI (\\%)     | 92.54 | 93.73 | 93.89 | 94.89 | 95.24 | 95.65 | 95.85 | 96.09 |\n\n---\n\nThe results show that performance continues to improve beyond k=10, but with diminishing marginal returns. Since the runtime scales linearly with k, the significant increase in computational cost for these very high orders may not be justified by the smaller performance improvements. Therefore, we set k to moderate values (up to 10) in our main experiments, as it provides substantial performance gains over lower-order methods while maintaining computational efficiency. Ultimately, the optimal choice for k depends on the specific application and the available computational budget."}, {"review_id": "ttjPvwBCZc", "reviewer": "Reviewer_Evn6", "summary": "This paper addresses the problem of expanding mean-field Variational Inference (VI) to account for instance-level dependencies in a computationally efficient way. By considering a particular form of Gaussian variational family, the authors are able to incorporate instance-level dependency within a predefined and tunable range into the model. Guidelines are provided to confine the iterates within a suitable region of the optimization space. Experiments on real-world data support the method's practical applicability and efficacy.", "questions": "* TreeVI appears to be the most similar work to yours. Could you elaborate on the key differences between your method and TreeVI, and highlight the improvements brought by your approach?\n* How does your method compare to concurrent methods like TreeVI in terms of effective runtime?\n* While some information about the datasets used in the experiments is provided in Appendix D.1, it is not said exactly how large these datasets actually are and how many observations they contain. For instance, SMAP is a collection of data products. Could you provide more details on this matter?\n* Lines 80-81: “the method can easily be extended to take the dimensional-level correlation into account”. Can you expand on this?", "limitations": "Yes.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Strengths**\n* The proposed method is original and its implementation seems to be quite simple. There is a clear way to dial between the expresiveness of the variational family and the model's complexity (however this trade-off is not explored in details, see Weaknesses below).\n* The method empirically outperforms the baselines in a wide variety of real-world scenarios, including anomaly detection (3 datasets), time series forecasting (4 datasets) and clustering (4 datasets).\n\n**Weaknesses**\n* There is no runtime study, Appendix D.4 is not very informative about the computational cost of the method. The trade-off between the gains yielded by more complex models and the added computational cost is not discussed. This is essential for establishing the method's superiority over existing approaches, and should be addressed by the authors.\n* The variational posteriors must be Gaussian, hence the method is not applicable to settings where such approximating families are irrelevant, e.g., directional statistics.\n* In the constrained clustering experiment, it is unclear how the constraints are incorporated into the model (even with the additional details in Appendix D.3).\n* The authors do not discuss how the proposed reparameterization effectively restricts the expressiveness of the variational family.\n\n**Additional comments**\n* Line 77: there are two product signs in the equation.\n* Table 5: typo in “Epchs”\n* Throughout the text: replace \"validness\" with validity, and \"indexes\" with indices.", "was_revised": true, "final_justification": "The paper introduces a novel, interesting, and empirically well-justified method.\nI was originally concerned that the method may be a simple extension of the already-existing TreeVI, but the authors have clarified important differences between both procedures. Responses to other reviewers go into more details about this matter, and the authors are commited to spending more time explaining the key differences in their revised manuscript.\nI was also not convinced that the performance gains were worth the added computational complexity. In their rebuttal, the authors provided a runtime study that successfully addressed this concern.\nWhile the method introduced in this paper is strictly limited to the Gaussian setting, I believe that it still appears as a relevant contribution. It is noteworthy that this method achieves state-of-the-art results on a variety of tasks.", "rebuttal": "We thank the reviewer for the comments and suggestions. Below we respond to each point raised by the reviewer.\n\n**Q1: What is the runtime study for the proposed method, and what is the trade-off between the gains from more complex models and their added computational cost?**\n\nA1: An empirical runtime study was conducted to measure the cost of increasing the approximation order k. Tables 1-3 below show the average training runtime per epoch on different tasks, demonstrating that the computational cost scales linearly with k.\n\n---\n\nTable 1: Time series anomaly detection performance on SMAP dataset.\n| Methods  | Mean-field | k=1 | k=3 | k=10 | k=50 | k=100 |\n| - | - | - | - | - | - | - |\n| Runtime (s/epoch) | 1.00 | 2.44 | 8.51 | 27.60 | 142.58 | 294.53 |\n| F1     | 0.7774 | 0.8411 | 0.8552 | 0.8636 | 0.8711 | 0.8755 |\n| ELBO   | -109.2182 | -97.6057 | -95.2314 | -92.2948 | -90.4291 | -89.9577 |\n\n---\n\nTable 2: Time series forecasting performance on ETTh1 dataset with horizon 24.\n| Methods  | Mean-field | k=1 | k=3 | k=10 | k=50 | k=100 |\n| - | - | - | - | - | - | - |\n| Runtime (s/epoch) | 4.80 | 11.45 | 36.94 | 126.44 | 675.75 | 1340.22 |\n| MSE | 0.739 | 0.664 | 0.543 | 0.363 | 0.348 | 0.333 |\n| MAE | 0.716 | 0.570 | 0.505 | 0.376 | 0.362 | 0.352 |\n\n---\n\nTable 3: Constrained clustering performance on MNIST dataset.\n| Methods  | Mean-field | TreeVI | k=1 | k=3 | k=5 | k=10 | k=50 | k=100 |\n| - | - | - | - | - | - | - | - | - |\n| Runtime (s/epoch) | 0.25 | 0.52 | 0.59 | 1.84 | 3.09 | 6.22 | 29.53 | 60.59 |\n| ACC (\\%)     | 96.50 | 97.43 | 97.55 | 98.12 | 98.39 | 98.52 | 98.62 | 98.69 |\n| NMI (\\%)     | 91.37 | 93.14 | 93.44 | 93.80 | 94.17 | 94.55 | 94.63 | 94.85 |\n| ARI (\\%)     | 92.54 | 93.73 | 93.89 | 94.89 | 95.24 | 95.65 | 95.85 | 96.09 |\n\n---\n\nThe trade-off is between this linear increase in runtime and the corresponding performance gains. While model performance consistently improves as k increases, our experiments show diminishing marginal returns, with higher values of k (such as k beyond 10) yielding progressively smaller gains. Therefore, we focused on k up to 10 in the main paper because this range demonstrates the advantages of our approach while maintaining a practical and efficient balance between performance and computational tractability.\n\n**Q2: What are the specific limitations of the method given that its variational posteriors must be Gaussian, particularly concerning its applicability to settings like directional statistics?**\n\nA2: The Gaussian posterior is a foundational and highly effective choice, used with great success in many state-of-the-art generative models like VAEs and modern diffusion models. We acknowledge that, while powerful, the Gaussian posterior is not universally applicable and is not directly suited for modeling certain data types, such as discrete data or directional statistics. The core contribution of our work is to introduce a novel and scalable method for structuring correlations, and extending this framework to more flexible distributional families to address these broader settings is a promising direction for our future research.\n\n**Q3: How are the constraints specifically incorporated into the model within the constrained clustering experiment?**\n\nA3: As detailed in Section 4.2, the constraints are not implemented as hard rules but are integrated as a probabilistic prior over the cluster assignments $\\mathbf{c}$. This soft constraint mechanism guides the model to favor clusterings that are consistent with the provided must-link and cannot-link information. Specifically, the key is to define a conditional probability distribution for the cluster assignments $\\mathbf{c}$, conditioned on a pairwise prior information matrix $\\mathbf{A}$. This distribution $p( \\mathbf{c}|\\mathbf{A})$, as shown in Line 345, replaces the uninformative prior $p(\\mathbf{c})$ used in standard generative clustering models, to assign a higher probability to clusterings that respect the constraints encoded in $\\mathbf{A}$ (If there is a must-link between instances $x_i$ and $x_j$, then $[\\mathbf{A}]_ {ij}$ is a  large positive value; if there is a cannot-link between $x_i$ and $x_j$, then $[\\mathbf{A}]_ {ij}$ is a large negative value; if there is no constraint, then $[\\mathbf{A}]_ {ij}=0$). This prior is then integrated into the overall generative model $p_{\\theta}(\\mathbf{X}, \\mathbf{Z}, \\mathbf{c} | \\mathbf{A}) = p_{\\theta}(\\mathbf{X}|\\mathbf{Z})p(\\mathbf{Z}|\\mathbf{c})p(\\mathbf{c}|\\mathbf{A})$, guiding the learning process within the variational inference framework towards a posterior that respects the given domain knowledge.\n\n**Q4: How does the proposed reparameterization method restrict the expressiveness of the variational family?**\n\nA4: Our reparameterization method trades the universal expressiveness of a general multivariate Gaussian for guaranteed validity and computational tractability. It restricts the variational family to a subset of Gaussian distributions whose precision matrices have a k-order connection structure. And the hyperparameter k allows us to explicitly control the trade-off between expressiveness and computational cost.\n\n**Q5: What are the key differences between the proposed method and TreeVI, what specific improvements does your approach offer, and what is the comparison in terms of effective runtime?**\n\nA5: The key difference between TreeVI and our method is that TreeVI is inherently limited to modeling only first-order correlations, while our method can capture higher-order correlations. This limitation of TreeVI arises from its reliance on an acyclic tree structure to construct its correlation matrix. While this is sufficient for simple pairwise relationships, attempting to model higher-order correlations inherently introduces loops into the underlying correlation structure. The construction that TreeVI uses to build a valid correlation matrix depends on the acyclic property of its backbone and is no longer valid when these loops are present. Moreover, simply modeling higher-order correlation coefficients within the framework of TreeVI does not guarantee the resulting correlation matrix to remain positive definite. So even though TreeVI can capture instance-level correlations, it cannot be easily extended to model higher-order correlations.\n\nOur approach, in contrast, is generalized to explicitly model k-order correlations. For $k=1$, our model is conceptually equivalent to TreeVI, but for $k>1$, it becomes substantially more expressive. As shown in our constrained clustering results (refer to Table 3 above), the runtime of our method is approximately k times that of TreeVI. But even with a small k (e.g., k=3), our method achieves significantly better performance across all metrics compared to TreeVI, demonstrating the benefit of modeling higher-order correlations.\n\n**Q6: More detailed information about the actual sizes of the datasets used in the experiments (e.g., the number of observations for SMAP)?**\n\nA6: The actual sizes of the datasets used in the constrained clustering task can be found in Appendix D.1. Detailed information about the other datasets is shown as follows and will be included in the appendix of our final version.\n\n| Tasks | Dataset | Dim | Size (Train, Validation, Test) | Domain  |\n| :--- | :--- | :-- | :--- | :--- |\n| **Forecasting** | ETTm1 | 7 |  (34465, 11521, 11521) | Electricity |\n| | ETTh1 | 7 |  (8545, 2881, 2881) | Electricity |\n| | Electricity | 321 |  (18317, 2633, 5261) | Electricity |\n| | Weather | 21 | (36792, 5271, 10540) | Weather |\n| | Exchange | 8 |  (5120, 665, 1422) | Exchange rate |\n| **Anomaly Detection** | SMD | 38 |(566724, 141681, 708420) | Server Machine |\n| | MSL | 55 |  (44653, 11664, 73729) | Spacecraft |\n| | SMAP | 25 | (108146, 27037, 427617) | Spacecraft |\n\n**Q7: Clarify extending the method to dimension-level correlations.**\n\nA7: While we focused on instance-level correlations in the paper for clarity, the method can be extended to model dimension-level correlations by introducing a structured covariance to the sampling procedure's noise term.\n\nIn the current framework, each dimension is sampled independently, meaning the noise vector $\\epsilon$ in the re-parameterization step is drawn from a standard Gaussian with a diagonal covariance matrix $\\mathcal{N}(0, I)$. To incorporate dimension-level correlations, we can replace this with a low-rank Gaussian distribution for the noise. By sampling $\\epsilon$ from a distribution with a low-rank covariance, we can efficiently induce correlations across the dimensions of the generated latent variables $\\mathbf{z}$. This provides a practical and computationally efficient path to a more expressive posterior that captures both instance-level and dimension-level dependencies within a unified framework."}]

Author rebuttal:

Meta-review: We thank the authors for their submission. 

The authors propose a structured variational approximation that enables structured covariance between latent variables while remaining tractable, and providing a direct tradeoff between fidelity and computational efficiency.  Reviewers voiced concerns over potential runtime issues, which the authors addressed throughout the author-reviewer discussion.  In a somewhat crowded area, the proposed method strikes a novel balance between efficiency and expressivity.  In particular, this work enables the capturing of higher order correlations, which is empirically shown to improve performance.
Decision: Accept (poster)