PAC-Bayes Bounds for Multivariate Linear Regression and Linear Autoencoders
Authors: Ruixin Guo
Ruoming Jin
Xinyu Li
Yang Zhou
Venue: NeurIPS 2025 poster ()
Keywords: PAC-Bayes bound
multivariate linear regression
linear autoencoder
recommender system
Primary area: theory

Abstract: Linear Autoencoders (LAEs) have shown strong performance in state-of-the-art recommender systems. However, this success remains largely empirical, with limited theoretical understanding.  In this paper, we investigate the generalizability -- a theoretical measure of model performance in statistical learning -- of multivariate linear regression and LAEs. We first propose a PAC-Bayes bound for multivariate linear regression, extending the earlier bound for single-output linear regression by Shalaeva et al., and establish sufficient conditions for its convergence. We then show that LAEs, when evaluated under a relaxed mean squared error, can be interpreted as constrained multivariate linear regression models on bounded data, to which our bound adapts. Furthermore, we develop theoretical methods to improve the computational efficiency of optimizing the LAE bound, enabling its practical evaluation on large models and real-world datasets. Experimental results demonstrate that our bound is tight and correlates well with practical ranking metrics such as Recall@K and NDCG@K.

PAC-Bayes Bounds for Multivariate Linear Regression and Linear Autoencoders
Ruixin Guo Kent State University rguo5@kent.
Ruoming Jin Kent State University rjinickent... @kent edu
Xinyu Li Kent State University x1i740kent edu
Yang Zhou Auburn University
yangzhou@auburn.co edu
Abstract
Linear Autoencoders (LAEs) have shown strong performance in state-of-the-art recommender systems. However, this success remains largely empirical, with limited theoretical understanding. theoretical measure of model performance in statistical learning - of multivariate linear regression and L.F... We first propose a PAC-Bayes bound for multivariate linear regression, extending the earlier bound for single-output linear regression by Shalaeva et al. [45, and establish sufficient conditions for its convergence. We then show that LAEs, when evaluated under relaxed mean squared error, interpreted as constrained multivariate linear regression models on bounded data, to which our bound adapts. Furthermore, develop theoretical methods improve the computational efficiency of optimizing the LAE bound, enabling its practical evaluation large models and real-world datasets. Experimental results demonstrate that our bound is tight and correlates well with practical ranking metrics such Recall @K K and NDCG@K.
this paper, we investigate the generalizability
Introduction
In recent years, simple linear recommendation models have consistently demonstrated impressive performance, often rivaling deep learning models 122 24 35). In particular, linear autoencoders (LAFs) such as EASE [48] and EDLAE
(49 have shown a surprising edge over classical linear methods like ALS 233 Despite their empirical success and widespread adoption, the theoretical understanding of why LAES perform so well remains limited. Moreover, much of recommender system research has focused heavily on empirical comparisons, where weak baselines and unreliable sampled metrics often render evaluations biased and difficult to reproduce [12 II]. A solid theoretical foundation is therefore urgently needed explain and justify the true performance of recommendation models beyond purely empirical assessments.
Statistical learning theory [11 provides such foundation by estimating model's theoretical perfor- mance over the underlying data distribution. classic result in this area is the uniform convergence PAC bound established by Vapnik and Chervonenkis [52]. Although this bound guarantees conver- gence, it characterizes the worst-case generalization gap and is typically vacuous for large models such as neural networks, thus failing to reflect true model performance in practice (39 [3] To address this limitation, a variant PAC framework incorporates information-theoree
techniques to bound the expected generalization gap, which often yields tighter results [20]. One notable example is the PAC-Bayes bound, first introduced by McAllester [37]. Recently, Dziugaite and Roy [13 empirically demonstrated that PAC-Bayes bounds can remain non-vacuous even for large neural networks, suggesting that PAC-Bayes theory provides more accurate and practically meaningful characterization of model performance.
While statistical learning has been extensively developed for a wide range of machine learning and deep learning models [1 6 53 its application recommendation systems remains largely
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

underexplored, with only few exceptions 44 been directly applied to LAEs. In this work, we aim to advance the theoretical understanding of LAE performance through PAC-Bayes theory.
46, 16]. To the best of our knowledge, it has not yet
When evaluated under mean squared error, LAFS are closely related multivariate linear regression models. Prior works have developed several PAC-Bayes bounds for linear regression. Notably, Alquier et al. [5] proposed PAC-Bayes bound for general machine learning models with unbounded loss; Germain et al. 17 adapted this bound to linear regression models with mean squared loss, but their result does not converge; and Shalaeva et al. [45] improved Germain' bound deriving a strictly tighter and convergent version. We aim to extend existing PAC-Bayes bounds for linear regression to the setting of LAEs. Several challenges arise in doing so:
Multivariate Data: Existing PAC-Bayes bounds for linear regression primarily focus on the single- output setting. Extending them to the multivariate (multi-output) case is nontrivial, as the outputs can exhibit statistical dependencies. To derive a valid bound in this setting, more general assumptions on the data distribution must first be established to capture potential output dependence.
LAEs differ from multivariate linear regression in several important aspects. They typically operate bounded data, so the standard Gaussian data assumption used in linear regression does not apply. Moreover, LAEs impose unique structural constraints, such as the zero-diagonal constraint on the weight matrix and the hold-out constraint on data. These characteristics must be rigorously defined and integrated into the theoretical framework.
spe-siecific Characteristics:
Computational Inefficiency: Optimizing PAC-Bayes bounds is typically computationally expensive, as estimating the distance between prior and posterior in high-dimensional spaces is complex [4]] Since LAE models often contain hundreds of millions of parameters and are computed on large real-world datasets (Table 2). developing computationally efficient methods is crucial for practical bound evaluation.
This paper addresses the aforementioned challenges and makes the following key contributions:
(Section 3) We generalize Shalaeva's 's Gaussian data assumption 455 to the multivariate setting (Assumption 3. and propose a corresponding PAC-Bayes bound for multivariate linear regression (Theorem 3 2 extending Shalaeva's single-output bound 45 We further establish sufficient
conditions (Theorem 3.3) that guarantee convergence for both bounds.
(Section 4 We propose relaxed mean squared error for evaluating LAE models and show that, under this loss, LAFS can be viewed as constrained multivariate linear regression models on bounded data. Building on this, we adapt our PAC-Bayes bound to LAEs by replacing the Gaussian data assumption with bounded one (Assumption 4.1 and incorporating LAE-specific constraints: the zero-diagonal constraint on weights and the hold-out constraint on data (Section 4.2).
(Section We develop theoretical methods to improve the computational efficiency of optimizing the bound. Following Dziugaite and Roy 13, we restrict the prior and posterior distributions be Gaussian (Assumption
[.J), leading to a closed-form expression for the tightest bound (Theorem 5.2). We further establish a practical upper bound with reduced complexity to mitigate the computational cost introduced by the zero-diagonal constraint (Theorem 5.4). 5.44).
(Section 6 We evaluate the bound for LAES on real-world datasets. Experimental results demon- strate that our bound is tight and correlates well with practical evaluation metrics such as Recall K and NDCG@N K, suggesting that it effectively reflects the actual model performance of LAEs.
All proofs of the theorems, lemmas and propositions in the main paper are provided in Appendix A Related Works are in Appendix D Conclusions and Discussions are in Appendix
2 Preliminaries
Notation: We denote S [x,, ym as the dataset, where 2; E R" and ys E RP for all i, and m is the number of samples. For each sample (I,, 2; is the input and y; is the target. Let fw : Rn > R be the linear regression model parameterized by W E Rox" The model prediction is fw (x.) = Wx, and the mismatch between target and prediction is measured by the squared Frobenius norm loss 11 - Wr
Denote X = T2,,, Lm] E Rx as the input matrix and Y = [vn V2. Vm E Rexm as the target matrix. The empirical risk, representing the average loss on the observed dataset S, is then defined Homp (W) 1 -n |y Wx,ll? || W | To evaluate performance

---

unseen data, we assume that each (I; V/) is i.i.d. sampled from an unknown distribution D, and define the true risk as (W" = B(a,))> [lly |
To construct PAC-Bayes bound, we treat fw as stochastic model by considering W random variable, Denote т as the prior distribution over W and p as the posterior distribution. я represents our initial belief about W before observing any data, whereas p represents our updated belief after incorporating information from the dataset 5 (88]
Multivariate Linear Regression [25]: From the definition above, the linear regression equation can be written as Y WX + E, where E [e1, C22 ..: Cm] E Rexm is the error matrix.
Usually the first dimension of every x; is set (i.e., X1, is a vector of all 1s) to represent the bias term. We say the linear regression is multivariate (or multi-output) if 1.
Multivariate linear regression typically assumes that the error vectors e, and ej are independent for i # j. but allows dependencies among the elements within each ei. This leads to our statistical assumption stated Assumption 3.1
Alquier's Bound [5] adapted to Linear Regression: Alquier's bound is general bound that can be applied any model with unbounded loss. When adapted the linear regression model using our notation, Alquier's bound can be stated as follows: Given for any À >>> 0 and 5 >>> 0,
P Vp, (R < Ew [R (W)) + [(ollT)+ + + v (2 2 21-8 (I)
where Vr, m) = In Rw Bs (W) and D(plI -w [In p(aw) denotes the Kullback-Leine respect to The loss is typically assumed following a light-tailed distribution, such as sub- uub-Gaussian or sub-exponential, to ensure that ,,, (A, m) is bounded [7 18 19].
=(dW)
(KL) Divergence where the measure is absolutely continuous with
Note that Alquier's 's bound holds simultaneously for all posteriors If replacing v with any single such as Gaussian posterior or a Gibbs posterior, the bound still holds.
Shalaeva's Bound [55] Shalaeva's 's bound is an application of Alquier's bound to the single-output linear regression model, i.e., the case p = 1. further assumes D is Gaussian: Given constants Tn and Ge, for any draw (I. ~ D, 1. x ~ (N o²/), and 2. there exists W E R Rixa such that y W e, where e (N o is Gaussian noise. Under this assumption, D is fixed in L (A, while π remains unspecified. Germain et al. [17] showed that, if π is Gaussian, Vn, (A, m) is bounded, since the true risk R! (W) with W ~ π is sub- gamma; however, their bound is independent of m and does not guarantee convergence. Shalaeva et al. [45] improve this by showing that V ",, m) fact has a strictly tighter upper bound: For any π,
exp(Auw)
22A22
V (A, m) = In EW (1
< In EW-a exp
w
Alwe pm/2 m/2
where Ur = - - W +o² This bound depends on mn and can be used to establish convergence.
Convergence of Shalaeva's Bound: The convergence analysis in Shalaeva et al.'s paper 45) is
presented informally. Here we formally state their results as follows: Since limm-+so exp (Au,,), ), for any À > 0, the convergence of V Vr,D(M (A, m) follows from
A"wy ym/2
m/?
exp(Auw) A )m/2 m/2
exp(Auw)
lim V (n,D(A, m) lim InE X--50
In Ee lim
0 (3)
(1
(1 An )m//
m/2
Upon careful examination of their analysis, we found that additional conditions are required to guarantee (3). which were not discussed in their original paper. Specifically, swapping lim and is valid only under certain conditions. For example, by the dominated convergence theorem 44] the condition can be (cxp<</[c[[Xxrr convergence is not guaranteed. These issues are discussed in Section 3.2 and Appendix B
< xc. If the choice of (X, π) does not satisfy this condition,
Collaborative Filtering Recommenders for Implicit Feedback: In collaborative filtering, an implicit feedback dataset is typically represented as a binary user-item interaction matrix H E (0. 1) |) with items and m users (Section 1.3.1.1, D). Fach H, denotes an interaction: Hy means that user has interacted with item i, while Hag means observed interaction.
3

---

Suppose H is a test set. To evaluate a model, we typically hold out fraction 1 - p (p E (0. 1)) of 1s in H (Section 7.4.2. M: also [22 48 (8)) Formally, we use a binary mask matrix A {0, 1 m to perform the hold-out operation. A denotes a random matrix, where each Aij is independently drawn from a Bernoulli distribution conditioned on Ha:: P(A:j = 1 Hij 1) p. P(A.j
0 H3 = 1) = 1 - and P(Ais 1) = H, = 0) =1
Let A be a realization of A. Define the input matrix as Hinpat = A C H and the target matrix as A) @ H, where @ denotes the Hadamard (element- product and 1 [1] For any j.j such that Hi 1, we say Hi is held out if Aij = O, which yields H = 0 and Mig 1. Consequently, Hiinput retains p fraction of the 1s in H, while the remaining 1 p fraction are held-out and moved to H/uage
Hurget
collaborative filtering model takes Hinp" as imput and generates a prediction matrix H E R Model performance is typically evaluated by how well (1 - A) C H approximates H The masked prediction (1 -) @ Hreed indicates that only the entries 1/pied that coincide with the held-out interactions contribute to the evaluation. This approximation quality is typically measured using ranking-basee metrics such as Recall@ K or NDCG@ K (Section 7.5.3, 7.5.4, D.
LAE Models and EASE [48]: : LAE models are a class of collaborative filtering models. They are typically represented by a square matrix W E R=x= and trained by solving arg minw /| W H|] where I/ denotes the training set. The model takes // as input and generates a prediction W H. which aims to reconstruct H itself. W is commonly constrained by zero diagonal (i.e., diag(W) o to
prevent overfitting towards the identity matrix I [48]]]]]]]]]] 49 50]. Some studies relax this constraint by allowing a diagonal with bounded norm instead
38
EASE is one of the most popular LAE models, obtained by solving
arg min H - H 2- + w/ s.t. diag(W) =
(4)
W
where 7 is the regularization parameter. Let Wo be the solution of (4), then Wo has closed from:
Let P = (HII" + then (Wo)ss = if =j and (Wo) = PPP./ if $
3
PAC-Bayes Bound for Multivariate Linear Regression
3.1 The Statistical Assumption and the Bound
We first generalize Shalaeva et al.'s Gaussian data assumption 45 to the multivariate data with dependent outputs and potentially degenerate covariance. Based on this assumption, we derive our bound and show that Shalaeva's 's bound is a special case of ours.
Assumption 3.1. Let Ma E R", Σx Rexn be positive semi-definite, and Eg E R be positive- definite. Suppose yz ~ D satisfies: 1. x ~ NN Σz); 2. there exists W* e R2x= such that y W e, where (0, E.); in other words, y|ı: ~NWWx, Ee).
The positive semi-definite assumption of Σ, allows it to be singular, implying a degenerate Gaussian distribution whose support lies on a lover-dimensional manifold embedded R". This includes the standard multivariate linear regression setting in which the first element of x is and the remaining n elements are Gaussian. In this case, the first row and first column of Σx, are
Under Assumption 3.1
for any model W E RRxx,, the prediction error y - Wx = (W* W)x + follows the Gaussian distribution N (P E ). where
Pw [[WW w)x + e] (W W)E[=] E(e] (W* W)pz
Σn = [[(W" - w) - m2) (W - (z - -z) -
(W" W) (W*-W WW)T + M
Note that Σn is positive definite due to the positive definiteness of E.. Let Eny S AS be its eigenvalue decomposition where S is orthogonal and A = diag(m)) .... 7pp with 1): > 0 for all i. Both S and depend on W. The PAC-Bayes bound for multivariate linear regression is then
stated as follows:
'The same symbol will be used elsewhere in this paper to represent all-ones matrices of different sizes. 'The notation diag is defined as follows: If W E 1=*x" then diag(W)& E R² denotes the vector consisting of the diagonal elements of W. If w E R", then diag(w) gax" denotes the diagonal matrix whose diagonal entries the elements of w.

---

Theorem 3.2. Denote b = // pm. Given for any X > 0 and 5 > 0.
P Vp,, St < Sw-pltt" (W)] D( (( lu + vr V..
21-2 (5)
where
exp E
Amminn
V Vx m) = In Bw exp (I(1(Z()) [ + + ww)
2
S In En exp
2 11
(1 + 2 /m)
m
The bound of Theorem 3.2 is general case Shalaeva's bound. can be reduced to Shalaeva's bound by taking p Mr O, Ex o and Σ. of for some constants ozy de.
3.2 Convergence Analysis
This section presents the convergence analysis of Theorem 3.2 We provide sufficient condition based on the dominated convergence theorem that guarantees convergence, thereby completing and rigorously formalizing the convergence analysis Shalaeva's bound [45]. This condition stated as follows:
Theorem 3.3. If X and satisfies Ew [exp (AI(E, MaP W) oo, then
limm-+oo RR (A, m) 0.
Here are some examples of the combinations (A, T) that satisfy the condition of Theorem 3.3
Example 3.4. Let be distribution with bounded support, then for any À >> O, the condition holds, because there exists constant G > with [[W]] F G such that
Ew-w
eep (AIl(E> + pap1) 1/2 w & Es (exp (>I(EE + mane " |W" -w1B3)]
exp (AI(EE= + + 1/1 ||= + [[ exp (A1l(E Mape 1/2 (IIW"|| G)²)
Ew-x
Example 3.5. Let π be Gaussian distribution parameterized by Uo E R and σ > 0, such that each Wi is independently drawn from (V o²). Let Σx + P QTQQQ be its eigenvalue decomposition, where A = diag(V), V2...., Vn) and is the largest eigenvalue, then
Auy (e, -ub)...²
exp
122ho
Bw [exp (A]((E. p.p 1/2 -w)
(1 2X0 1/2
In this case, the sufficient condition holds for any X (0, 2
Applying Theorem 3.3 Ew [exp (oo W W
Shalaeva's bound, then 3 []]] < 2o. Moreover, since
guaranteed if and क satisfies
Ew [exp (oo [
w|l) < Ew [exp W - w1 + +=o = Ew (ww)()
sufficient condition is therefore Ew-- (x <0.
4 PAC-Bayes Bound for LAEs
This section presents a PAC-Bayes bound for LAEs. The model W can be obtained using any training method, such as EASE, EDLAE or ELSA; our bound only focuses on analyzing its test performance and is independent of the training procedure.
4.1 Adapting to Bounded Data Assumption
Most real-world recommendation datasets are bounded rather than Gaussian. For example, the user-item interaction matrix introduced in Section 2 binary. we assume that each user vector H.; is i.i.d. sampled from an underlying distribution, this distribution must have bounded support. Therefore, adapt the PAC-Bayes bound for multivariate linear regression to recommendation datasets, we replace the Gaussian assumption on D (Assumption 3.1) with the following bounded- support assumption:
5

---

Assumption 4.1. Suppose D is characterized by three finite cross-correlation matrices Exe ((xy)) [2.r Eny ((rr))=(ay"" and Eny ((z.))=luyy)) where Σg. is positive definite.
This assumption is indeed general. It holds for all D with bounded support, and also holds for certain D with unbounded support such as Gaussian, since Assumption [3.1 jippies Assumption 4.1 Consequently. it enables the derivation of more general bound, although in this work we focus on its application bounded-data settings.
We first note that the true risk under Assumption 4.1 can be expressed in an explicit form:
Lemma 4.2. Given any W, the true risk can be expressed
R (W) ||WE / -1 - 2
+ tr(_pn)
(6)
Then, in I.. we have the upper bound Vn.D(A. m) S In E,
e ( which is obtained by
removing -eemp (W) due to its non-positivity. This upper bound, originally used by Germain et al. (Appendix A.4, (17), does not ensure convergence since it is independent of m. but it simplifies computation. By plugging (6) we get
Proposition 4.3. Denote B = --ET Eny Emp/2 and C e> ((((EEv)) 1-revv)-1E 11 1 Then ! holds, with
L Vr, (A, m)) apper-bounded by
V/ (x, S E. ( In E.
4.2 Applying the PAC-Bayes Bound to LAES
Recall from Section 2 that the LAE model is defined by a squared matrix W E Raxx, and evaluated by comparing the closeness between the masked prediction (1 1 A) @ (W o H)) and the target metrics such as Recall@ K and
(1-A) C H. This closeness is typically measured by ranking- NDCG @ K. which are discrete and difficult to analyze statistically. To simplify the analysis, use
mean squared error (MSE; see Section 7.5. 1, M instead. The classic MSE is defined as the mean squared Frobenius norm of all held-out interactions:
(1 - A) o H (1 A) O (W(A @ H))
m
where the mask 1 - A on prediction distinguishes it from multivariate linear regression. We relax the MSE by removing this mask, allowing all predicted interactions to participate in the evaluation rather than only the held-out ones:
A)
11 W(A H)
m
Let A O H be the input and Y (1- A) の H be the target. The is in X.j represent items observed by user j. while the Is Y.₁ represent items that are hidden but potentially of interest to the user. The classic MSE only evaluates on held-out items, i.c., those with (X,, 0, 0. = 1). In contrast, our relaxed MSE (8) also accounts for items with with (X,, = 0, Y. = 0). indicating that items unobserved and unlikely to interest the user should not be recommended; and those with (Xe, = 1, Yú = 0). indicating that the model should avoid recommending items already observed by the user (see Section 7.3.4, M.. Consequently, (8) can be viewed special case of the empirical
risk of linear regression Re (W) m Y WX| P. with output dimension p = n. and under the following LAE-specific constraints:
Hold-out constraint on X and Y: For any i,j. X., and Yi, are either or 1, but cannot both be
2. Zero saro-diagonal constraint on W: diag(W) = 0 (Optional).
Since both X and Y are derived from H and A, the true risk can be defined by introducing statistical assumptions on and respectively. For H, we assume that each user vector H., is i.i.d. sampled from a multivariate Bermoulli distribution M, and denote ΣnA = [ as the cross-correlation matrix of M. For A, denote B the distribution from which each column A.j is independently drawn. Note that B depends on M, and this dependence encodes the hold-out mechanism: for 8 ~ B
and h ~ M, P(o. 1 1) p. P(o, 0bh, 1) 1 p and P(o, 0 hi 0) 1. In the
true risk 1e (W) E(D [lly W x: [1%.]. by plugging in x = d o h. y (1 8) @ h. get

---

(W) = BB M [[[(1 6) @ -
(9)
By Lemma 4.2 further obtain the following result:
Lemma 4.4.
9 can be written in the same form as
(6) by plugging in
Egg= - p- En hd + p(1 - (T @ EnN), Eny (1 -1- PI + p(1 P) @ Enhd, Egy pp(1 p) - TO Eph)
Furthermore, Egz is positive definite if Σia. is positive definite.
If W is subject to zero-diagonal constraint, such as models trained from EASE, EDLAE or ELSA, then an additional condition (iag(W) = 0 is applied to both the empirical risk (8) and the true risk
The PAC-Bayes bound for LAFs formed by plugging (8) and (9) into ① Since Assumption
holds for this bound, it directly leads to 77)
Practical Methods for Computing the PAC-Bayes Bound for LAES
In the PAC-Bayes bound for LAEs proposed in Section 4 the choice of a and p is so far unspecified, and the computation of the bound has not yet been addressed. This section develops theoretical methods improve the computational efficiency, enabling the bound to be evaluated on large models and datasets.
Our goal is to optimize the right hand side of (T): Given 6, find x, p, A that minimize
Ew [ (W)) D(PllT) In 8 + 4 MA,((A, m)
(10)
part
part2
with (W) given by { and Rinee (W) given by (9). However, solving for X, n. p simultancously is generally intractable
We therefore consider a weaker problem: Given À and æ,, we optimize (10) with respect to p. We compute the optimal p for different choices of À and æ,, and select the combination yielding the tightest bound.
We discuss the computation of part 1 of (10) in Section 5.1 and part 2 in Section 5.2
5.1
Closed-form Solution for the Optimal p under Gaussian Constraint
Given and X, we search for the optimal p by
1 arg min Ew (W)) + D(oll=)
(11)
Not all choices of and p make (T) easy solve. For example, given any n, optimal p is (
the Gibbs posterior, defined as p(dW) B. [ generally a complex distribution without closed-form density function or parameterization, making
However, the Gibbs posterior
analytically intractable practice.
To obtain a tractable and efficient solution, we instead optimize (TT) under the constraint that π, p restricted to specific distribution families. Notably, Dziugaite and Roy [13] proposed a practical way to compute PAC-Bayes bounds for deep neural networks by assuming x and p to be entry-wise Gaussian distributions, which allows the KL--Divergence D( [ [3] and formally state it as follows:
be computed analytically. We
follow Dziugaite and Roy's s assumption
Assumption 5.1. Let A, _Raxn with B O (entry-wise non-negative), and denote (A, B) as the entry-wise Gaussian distribution such that W ~ [AA, B) means each Wij is independently drawn from fromA((A,, (VAAj B.)). Assume p = (24, S) and π = (i,,, J²J), where U, S gaxn with S > 0 (entry-wise positive), J E (1 11nnm is the all-ones matrix, and > 0.
Applying the constraint diag(W) 0 to and implies setting diag(l() 0, diag(S) 0. diag(14) = and diag(o?J) = 0, since a deterministically zero random variable has zero mean and zero variance.
Since neural network models are typically non-linear and do not admit closed-form solutions for the optimal Dziugaite and Roy [13] solved for the optimal p using stochastic gradient descent
7

---

with gradients estimated via Monte Carlo sampling, which requires a trade-off between sample size and computational cost. Surprisingly, due the linearity of to LAE models, we find that their assumption adapted to LAFs admits closed-form solution for the optimal p. as shown in Theorem 5.2 This fosed-ffomm solution enables efficient, direct computation of p. avoiding sampling, iteration or trade-off procedures.
Theorem 5.2. (a) Under Assumption 5.1.
the closed-form solution of the optimal p of 77 is given by
== LYxT + 2h02
1
L [xxx]_ xx
2x0=²)
1
S. EX,77 for 1,j E (1, 2,...,n) n}
m
(b) If we add the constraint diag(W) = to both p and n, then the optimal becomes
1
Sez
Su = for i,J (1, ..., nit and i +
22 m
X
X7.
+
jn
u
1 Y
1
Uo ddag(())
L xxx + 2x
m
22o
where
2 2-aig
1 yx
xxx
2 1)
| diag
L xxx
-1
2x0
and Ø denotes element-wise division.
5.2 Reducing the Complexity of (n, m) under the Zero-Diagonal Constraint
Under Assumptions 4.1 and 5.1, the closed-fore expression of >,, m) is too complex derive explicitly, making direct computation infeasible. We therefore compute the upper bound given in 77 instead. This section shows that, under the zero-diagonal constraint, this computation has a high complexity of O(n*), which we reduce to O(n") ) by establishing a simpler upper bound.
We first consider the case without the constraint diag(W) (W) 0. Since assume -N o²J) in Assumption
5.1 WT ~ (N(()) I²T), thus (W,, 1//2 + B..)T = 2// + B ~
N (S (Mo) BB. (..)). In this case, the E, (W)) term in 7 can be further expressed
as follows:
Proposition 5.3. Let A OZ and STA be the eigenvalue decomposition where S is orthogonal and A = diag(1) 171 T2,,... mm) with mn 2... 2m
Denote j >> (Uo)). + BT"
A(b:)
12 n exp
E
(R)
where D' SS2
(12)
i=1j=1 (1- 22n;) 1/2
The computational complexity of 12) is O(n²), mainly due to the eigenvalue decomposition of A. Now we discuss the case that diag(W) W) = is applied. Denote n' as the distribution with the constraint diag(W) (W 0, that is, for W n', W;; = 0 for all i. Then n' =No(oo, [²(J- -)) where diag((4g) O, and W N/ ((4)) o2 - r')) where I' is matrix with I3 1 and other entries being Therefore, (W,. 1/2 + B,,) )T ~ N 2>1// (U)).. + BTT, o2
).. (- )
Denote 10) = O² (E1/2) ) ( (E1/) 1/2,T (T)), then A is singular and positive semi-definite. Let (0)) = S(1T A 5() be the eigenvalue decomposition where s(i) is orthogonal and (0)
diag(r) (1/ 112 (()) with m/ ... n 0. Then
(1) (6)
exp
Eg [ (W))
C E11=1 (1 (1-2X1
1-22n
1/2 where 6(0) = (A())))))
(13)
(2)
'We slightly abuse the notation of S and A. In Theorem 3.2 here they denote the decomposition of A.
they denote the decomposition of Vn whereas
8

---

The issue with 13) is its high computational complexity: We need to compute the eigenvalue decomposition for each A in order to obtain SK1) and A 1 Since each eigenvalue decomposition
costs O(n"), the computation of (13) costs O(n²), which is impractical for large n.
Since the direct computation of T (Ar""(W) is difficult, we instead compute upper bound, as established by the following theorem:
Theorem 5.4. Given (12) and 13). for any X (0: 2m1 where n/ is the largest eigenvalue of I ET [2)
E
Theorem 5.4 holds for any Uo including the special case where diag(Uo) = 0 for both n' and This theorem allows us to compute Er (e))e((W) instead of Ex (xll" (w) thereby reducing the complexity from 3) to Onn)).
5.3 The Final Bound and the Algorithm for its Computation
The final step in computing the bound is to select À that yields the tightest bound. Following 4 we search over finite grid 1 = {A1, λ2,... .... Λr} where I. is the number of elements in A and ) > 0 for i (1,2, ...; ., L}; details are provided in Appendix C Applying this grid search to (1. we obtain the final bound: Given æ, for any À E A and ó > O, with probability at least 5, the following bound holds for any p:
Ew R (RTe"(")) (W) S
D(ell=)+ +
L + In By l [
(14)
We now summarize computation 14 under the LAE setting in Algorithm T By default, the algorithm assumes that the dir-oniaa constraint W is applied. For unconstrained W the algorithm can be adapted by switching the solution for p from Theorem 5.2 (b) to Theorem 5.2 (a),
and by computing D( pPll || using the unconstrained case (25) instead of (42).
Note that this algorithm requires Σna Ba-M [hh") as input, which depends on M. In practice, M may be an oracle distribution, making Eaa unknown and inaccessible. However, there practical scenarios where a non-oracle M is available. For example, if there exists larger and fixed dataset H E {0, 1 nxxm' (m' m) such that the columns of are sampled without replacement from the columns of # we can take M as the population distribution over the columns of #h This is similar to the matrix completion setting of Srebro et al. [47], where all entries of a fixed matrix treated as underground truth and observed entries sampled from them. Under this assumption. EAh m H whole (H whwoleey? is known and accessible. We provide further discussion Appendix E.
Algorithm 1 Computing the PAC-Bayes bound for LAEs
Input: Enh, p. 8, G. A {A1 λ2, ..., AL}. X, Y, and an LAE model W (with tiag(W) 0).
Compute Exz; Eny: Enz with Ess, by Lemma 4.4
Set N(W, V(W,o21) (i.c., let Uo = W such that W is the mean prior of тT).
Let G (} be
store the results.
for each A; in A:
Compute p (U,, S) with v, A, by Theorem 5.2 (b).
Compute D [ with p, by (42) in Appendix F
Compute E [ Remp (W) with p. X.Y by 40 A Appendix
Compute E
Rimee [W) with p. En Ezgs by by 41 in Appendix F and let it be the left
hand side of denoted LH,
Compute E ARm AR"
(W) with n, Exr, Eny, Eggs X by 12
Compute the right hand side of 14 denoted as RH, with Ew RaMP (W))) D(pll ((P))))))) (
Append (LH2, RH;) to G.
Output: the pair (LH", RH*) in G, where RH" minigist (RH.).

---

Experiments
In this section, we conduct experiments compute the PAC-Bayes bound for LAES using Algorithm on real-world datasets, evaluate the tightness of the bound, and empirically assess its correlation with practical ranking metries such as Recall K and NDCG@K.
We adopt the strong generalization evaluation setting, which divides the entire dataset into a training set and a test set with disjoint users [48] 38 Let the entire dataset be H Wh E {0, 1 We split it into training set H & {0, 1) ] -m) and test set H {0, 1 jaxm by setting = 0.3m'. The test set Histt further split into an input matrix X and a target matrix Y, with hold-out fraction 1 2. The LAE model W is obtained by solving the EASE objective (4) on the training set H (The EASE model can also be replaced by other LAE models such EDLAE ELSA. Our bound only focuses on evaluation and is independent of the training method.). We set ry in 4 to values of 50, 100, 200, 500, 1000, 2000 and 5000 to generate seven different LAE models and evaluate them accordingly. Other inputs of the algorithm are set as follows: 0.01, o = 0.001,
A = (1, 4, 16, 32, 64, 128, 25,, 512).
Our experiments run on a machine with 500 GB RAM and an Nvidia A100 GPU. The GPU has 80 GB RAM. We use three datasets: MovieLens 20M (ML 20M), Netflix and MSD, with their details shown in Table in Appendix
The results are presented in Table 1 On the left side, each pair (LH. RH) is the output of Algorithm I where LH and RH represent the left-hand side and right- side of (I4), respectively. Detailed values of the components of RH are provided in Table B in Appendix FT The results demonstrate that our bound is tight: in all cases, RH is within times L.H, contrast to Dziugaite and Roy's non-vacuous bound for deep neural networks [33] where RH can reach up to 10 times LH in their experiments.
The right side of Table I reports Recall @ 50 and NDCG@ 100 for each model on test sets, with both metrics referenced from the EASE paper (48). Across all datasets, models with smaller LH and RH generally achieve higher Recall 50 and NDCG@ 100. This negative correlation reflects the expected relationship: LH and RH are lase-d where lower values are better; Recall @ 50 and NDCG@ 100 are ranking metrics, where higher values are better. Although minor deviations exist - for example, on MSD, the best LH/RH occur at y 1000, while the best Recall 50/NDGG@ 100 occur at y 500 the overall trend demonstrates strong alignment between LH/RH and Recall @ 50/NDDG 100, suggesting that our bound effectively reflects the practical performance of LAE models. Further discussion on this alignment is provided in Appendix E.4
Table 1: Experiment results
PAC- BBy-Bayes Bound for LAEs
Ranking Performance
Models
ML 20M
Netflix
MSD
ML 20M
Netflix
MSD
LH
61.66
87.22
15.96
Recal
0.3434
0.2567
0.3454
7 = 50
RH
128.66
178.11
32.60
NDCG@ 100
0.4342
0.3766
0.3187
LH
60.75
86.54
15.85
Recall 50
0.3453
0.2580
0.3472
7 = 100
RH
125.90
176.25
32.26
NDCG@ 100
0.4373
0.3785
0.3205
LH 60.06
85.96
15.76
Recall@50
0.3471
0.2592
0.3486
200
RH
123.67
174.55
31.94
NDCG@ 100
0.4402
0.3804
0.3220
LH 59.46
85.35
15.66
Recall@50
0.3489
0.2605
0,3490
7 500
RH
121.41
172.64
31.62
NDCG@ 100
0.4439
0.3826
0.3225
LH
59.19
85.00
15.64
Recall@50
0.3502
0,2612
0,3475
7 1000
RH
120.17
171.44
31.50
NDCG@ 100
0.4464
0.3840
0.3210
LH 59.09
84.72
15.68
Recall@50
0.3510
0.2619
0.3434
7 2000
RH
119.34
170.45
31.52
NDCG@ 100
0.4487
0.3854
0.3171
LH 59.19
84,48
15.83
Recall@:0
0.3506
0.2625
0.3340
7 5000
RH
118.91
169.47
31.77
NDCG@ 100
0.4509
0.3871
0.3079
10

---

Acknowledgments and Disclosure of Funding
This research was partially supported by NSF grant IIS 2142675. RG acknowledges the NeurlIPS 2025 Financial Aid Award.
We thank all anonymous reviewers of the current and previous submissions of this manuseript for their valuable comments. In particular, the use of the relaxed MSE as a testing metric for LAES was suggested by Reviewer nWHo, and the experimental comparison between our bound and practical ranking metrics was suggested by Reviewer TXNp.
References
Charu C Aggarwal. Recommender systems: the textbook. Springer, 2016.
[2] Pierre Alquier. Bayesian methods for low-rank matrix estimation: short survey and theoretical study. In International Conference on Algorithmic Learning Theory, pages 309-323. Springer, 2013.
[3] Pierre Alquier and Gérard Biau. Sparse single-index model. Journal of Machine Learning Research, 14(1):243-28)) 280, 2013.
[4] Pierre Alquier et al. User-friendly introduction to PAC-Bayes bounds. Foundations and Trends in Machine Learning, 17(2):174-33 303, 2024.
Pierre Alquier, James Ridgway, and Nicolas Chopin. On the properties of variational approximations of Gibbs posteriors. Journal of Machine Learning Research, 17(236); 1-1,,, 2016.
[6] Peter Bartlett, Andrea Montanari, and Alexander Rakhlin. Deep learning: a statistical viewpoint. Acta momerica, 30:87-200 2017-201,220 2021.
[7] Emmanuel J Candès and Terence Tao. The power of convex relaxation: Near-optimal matrix completion. IEEE transactions Information Theory, 56(5):2053-2000 2080. 2010.
Ioar Casado, Luis A Ortega, Aritz Pérez, and Andrés R Masegosa. PAC-Bayes- Chermoff bounds for unbounded losses. Advances in Neural Information Procexsing Systems, volume 37, pages 24350-24374, 2024.
[9] Olivier Catoni. A PAC-Bayesial approach to adaptive classification. preprint LPMA, 840());6, 2003.
[10] Badr-Eddine reconstruction guarantees for VAEs. In International Conference Artificial Intelligence and Statistics, pages 3066-0000077 3079, PMLR, 2022.
Yuyang Shi, Arnaud Doucet, and Benjamin Guedj. On PAC- Bay-Bayesian
[11] Paolo Cremonesi Dietmar Jannach. Progress recommender systems research: Crisis? what crisis? Al Magazine, 42(3):43-55 54, 2021.
[12] Maurizio Ferrari Dacrema, Paolo Cremonesi, and Dietmar Jannach. Are we really making much progress? worrying analysis of recent neural recommendation approaches. In Proceedings of the 13th ACM Conference on Recommender Systems, pages 101-109, 2019.
[13] Gintare Karolina Dziugaite and Daniel M Roy. Computing nonvacuous generalization bounds for deep (stochastic) neural networks with many more parameters than training data. In Proceedings of Conference on Uncertainty in Artificial Intelligence, 2017.
[14] Xinlong Feng and Zhinan Zhang. The rank of a random matrix. Applied mathematics and computation, (55(1)688966)) 6999, 2007.
[15] Gerald Folland. Real analysis: modern techniques and their applications. John Wiley & Sons, 1999.
[16] Rina Foygel, Ohad Shamir, Nati Srebro, and Russ R Salakhutdinov. Learning with the weighted trace-norm under arbitrary sampling distributions. Advances Neural Information Processing Systems, 24. 2011.
[17] Pascal Germain, Francis Bach, Alexandre Lacoste, and Simon Lacoste-Julicn. PAC-Bayesie Baycsia theory meets Bayesian inference. In Advances in Neural Information Processing Systems, volume 29, 2016.
[18] Maxime Haddouche and Benjamin Guedj. PAC- DAC-Bayes generalisation bounds for heavy-tailed losses through supermartingales. Transactions Machine Learning Research, 2023.

---

[19] Maxime Haddouche, Benjamin Guedj. Omar Rivasplata, and John Taylor, PAC-Bayes unleashed: Generalisation bounds with unbounded losses. Entropy. 23(10):1330, 2021.
[20] Fredrik Hellström, Giuseppe Durisi, Benjamin Guedj, Maxim Raginsky, et al. Generalization bounds: Perspectives from information theory and PAC-Bayes. Foundations and Trends in Machine Learning. 18(1):1-225) 2025.
[21] Roger A Hom and Charles R Johnson. Topics in matrix analysis. Cambridge university press, 1991.
[22] Roger A Horn and Charles R Johnson. Matrix analysis. Cambridge university press, 2012.
[23] Yifan Hu, Yehuda Koren, and Chris Volinsky. Collaborative filtering for implicit feedback datasets. In IEEE International Conference on Data Mining, pages 263-272. IEEE, 2008.
[24] Ruoming Jin, Dong Li, Jing Gao, Zhi Liu, Li Chen, and Yang Zhou. Towards a better understanding of linear models for recommendation. In Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining, pages 776-785 2021.
[25] Richard A Johnson and Dean W Wichern. Applied maltivariate statistical analysis. Prentice Hall Upper Saddle River, NJ, 2007.
[26] Yehuda Koren. Factorization meets the neighborhood: multifaceted collaborative filtering model. In Proceedings of the 14th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 426-434, 2008.
[27] Ychuda Koren, Robert Bell, and Chris Volinsky, Matrix factorization techniques for recommender systems. Computer, 22-3), 2009.
[28] John Langford and Robert Schapire. Tutorial on practical prediction theory for classification. Journal of Machine Learning Research, 6(3), 2005.
[29] John Langford and Matthias Seeger. Boundy for averaging classifiers. School of Computer Science, Carnegie Mellon University, 2001.
[30] Antoine Ledent and Rodrigo Alves. Generalization analysis of deep nonlinear matrix completion. In International Conference on Machine Learning. PMLR, 2024.
[31] Antoine Ledent, Rodrigo Alves, Yunwen Lei, and Marius Kloft. Fine- Fine-grained generalization analysis of inductive matrix completion. In Advances in Neural Information Processing Systems, volume 34, pages 25540-25552, 2021.
[32] Dawen Liang, Rahul G Krishnan, Matthew D Hoffman, and Tony Jebara. Variational autoencoders for collahorative filtering. In Proceedings of the World Wide Web Conference, pages 689-698, 2018.
[33] David G Luenberger and Yinyu Ye. Linear and nonlinear programming. 3rd Edition. Springer, 2008.
[34] The Tien Mai. From bilinear regressioon inductive matrix completion: quasi-Bayesian analysis. Entropy, 25()33333, 2023.
[35] Kelong Mao, Jieming Zhu, Jinpeng Wang, Quanyu Dai. Zhenhua Dong, Xi Xiao, and Xiuqiang He. Simplex: A simple and strong baseline for collaborative filtering. In Proceedings of the 30th ACM International Conference on Information & Knowledge Management, pages 1243-1252, 2021.
[36] Arakaparampil M Mathai and Serge B Provost. Quadratic forms in random variables: theory and applications. Marcel Dekker, 1992.
[37] David A McAllester. Some PAC-Bayesian theorems. In Proceedings the 11th Annual Conference on Computational Learning Theory, pages 230-234, 1998.
[38] Jaewan Moon, Hye-young Kim, and Jongwuk Lee. It's enough: Relaxing diagonal constraints in linear autoencoders for recommendation. In Proceedings of the 46th International ACM SIGIR Conference on Research and Development in Information Retrieval, pages 1639-1648, 2023.
[39] Vaishnavh Nagarajan and Zico Kolter. Uniform convergence may be unable to explain generalization in deep lcarning. In Advances in Neural Information Processing Systems, volume 32, 2019.
[40] Xia Ning and George Karypis. SLIM: Sparse linear methods for top-N recommender systems. In International Conference Data Mining, pages 497-506. IEEE, 2011.
12

---

[41] Luca Oneto, Sandro Ridella. and Davide Anguita. Do really need new theory to understand over- parameterization? Neurocomputing. 543: 126227, 2023.
[42] Benjamin Recht. A simpler approach to matrix completion. Journal of Machine Learning Research, 12(12), 2011.
[43] Borja Gálve, Ragnar Thobaben, and Mikael Skoglund. More PAC-Bayes bounds: From bounded losses, losses with general tail behaviors, to anytime validity. Journal of Machine Learning Research, 25(110): 43, 2024.
[44] Walter Rudin. Principles of mathematical analysis. McGraw-hill New York, 1976.
[45] Vera Shalaeva, Alireza Fakhrizadeh Esfahani, Pascal Germain, and Mihaly Petreczky. Improved PAC- Bayesian bounds for linear regression. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 34, pages 16600-566.220 2020.
[46] Ohad Shamir and Shai Shalev-Shwartz. Matrix completion with the trace norm: Learning, bounding, and transducing. Journal of Machine Learning Research, 15(1):3401-3423, 2014.
[47] Nathan Srebro, Noga Alon, and Tommi Jaakkola. Generalization error bounds for collaborative prediction with low-rank matrices. Advances Neural Information Processing Systems, volume 17, 2004.
[48] Harald Steck. Embarrassingly shallow autoencoders for sparse data. In The World Wide Web Conference, pages 3251-3257, 2019.
[49] Harald Steck. Autoencoders that don't 't overfit towards the identity. In Advances Neural Information Processing Systems, volume 33, pages 19598-19608 2020.
[50] Vojtëch Vančura, Rodrigo Alves, Petr Kasalický, and Pavel Kordík. Scalable linear shallow autoencoder for collaborative filtering. In Proceedings of the 16th ACM Conference on Recommender Systems, pages 604-609, 2022.
[51] Vladimir Vapnik. The nature of statistical learning theory. Springer science & business media, 1999.
[52] Vladimir Vapnik and Alexey Chervonenkis. On the uniform convergence of relative frequencies of events to their prohabilities. Theory of Probability & Applications, 16(2):264-280, 1971.
[53] Chiyuan Zhang, Samy Bengio, Moritz Hardt, Benjamin Recht, and Oriol Vinyals. Understanding deep learning requires rethinking generalization. International Conference on Learning Representations, 2017.
13

---

A Proofs of Theorems, Lemmas, and Propositions
Proof of Theorem 3.2
Given W and (x, y) ~ D, denote v = Wx. x. then V ~ N(/ w. Let Q E R such that En = QQ". Such Q exists since we can take Q 1 STA but we do not assume to be unique. Let (0, I), then we can write e = QE + Mw- Thus,
Re (W) = E( [lly Wxlli] = E. [[Q + Mw [1 = E, [[Qe + Mw " + pw)]
= Q Qe + p Qc + +=TQ""m Pw + 1 Swrlll tr(Q" Q) + Mw
=t(QQ")== + w Pw tr(2ww p Hw
(15)
Also, we can express the random variable ||0112: in quadratic form (Representation 3.1a.1, [36]:
||ul = v²u (Qc Mw) )w) 3[Qcc+ (Ww+
(Qe + (ww )2 ">11/", E, 1/2 (w)
( + +- 5--1/2 kw)
kw) TT "Ew (E + (/w
(ww)
-2= + / kw) 1T AS(2 + 1111 (ww)
= -SE + ---11/2 pw )" A(SE 1/2 Q2 + -5.11/2 1/2 plw)
(16)
Denote e' = SE -/,1/2Qe,, then é' ~ (0, 1). because Ele') 1/ Q[E[]] 0 and
Cov(c] = EE]e E = SE 1/2 QE 10 x 1/2 = I
As b 1/2 ,w,
write
= - -c' (d A(e' +b Σ m(e b,)²
Hence, each c, + b, is independently from (b; 1). and each (<) + b,)² is independently from the non-central chi- chi-squared distribution of noncentrality parameter by and with degree of freedom. Thus the MGF (c) + b,)² is
exp 1-27
Mic(+a,)a(t) E te [et
(17)
(1
2t)
Given i..d. samples {(r,, My) 1 from D, let vj V/
W Waz. Then U1, U2,... Un are i.idd from
N( Ew), and
m 11 Wx,1l 2,,11
1 m
(W"
=
m 1=1
1=1
Hence the MGF of (W") is
M (t) ESnDm
(1))
ESnDm
exp
115113
m
j=1
L exp
m |1011/ ))
L
m
EgwDm
Esndm
exp
Ce + bz)
m
tbf
tyy (e( + bz)2 }1)
exp
Bee(+o.)a
exp
372
1/2
exp
Σ
209p
(18)
(1
20ms/mm m/2
14

---

By and 18
15
we can expand MP, Wx,D(A, m) as
Vn, Vz,D(A, m) = In BEE Bs Dm | ( (w) (w)) (
(W -RR
Amopn
exp Di=l m-22A
In Eww
exp (A (tr(Ew) m Mw))
(19)
IP: (1 + 2Am m mm/
Use the inequality that for any x 0 and k 0, =FF we have
1)th
and the fact tr(Sww Ei=1 !/i
Σ
-mmb, my m 2Am
exp
In Ew
exp (A (Ir(Iw) + 1 (w
Il (1 2Amss //mps/mm)))))) m/2 mm/m
- Ammmmm n m+22A
exp
Ei=1
< In Bw
exp (> (tr(2w)))) + + )w))
[T:=1 exp
mAni m+22mm
mb² 1/d
mdri m 2Am:
In Ew exp
Api Hw
A(n=
m 2ATR
2X²1 - >mmbf
In Ew exp
Aplly Pow
T
m
2Am,
2X 2X 1/
22>(>>
1p²)
< In Ew exp
Xp Mw
In Ew exp
m
m71
The last equality above is because
byn = Ab Wh 1-1/2 ASE Mw po
Since
Σ 11
tr((... x.T = S5 11 = tr(
we have
(2>((E)) mi
22 17
In In exp
In Ewar exp
Proof of Theorem 3.3
By 19
given any x > let (fm)met be sequence of functions where
exp E
Ambdfn m=+2An
Sm [m) exp (A (ir(Ew) + u (w))
11 (1 + 22m:/m)
for m > 0, and
fo(W exp
exp (A ) + p Hw))
Note that each fi is a non-negative function.
*Since _ < In(2 + 1) for any x -1, replacing with z, and taking exponential both sides, get 3FI
eath < { + 1)..
15

---

Since W ~ т and W E Rexp, let E = Rexp, then fm is measurable function on E, x is Borel probability measure on E. and we can express Ex [mm] as a Lebesgue integral (Section 10.1, [[5]:
f fm du
Ex [fml
Now prove the following three conditions:
(a) (W) S fu(W) for any m 2 and W E.
For any given W, and b; are fixed for all i, where f > and bf 2 0. Thus for any m 2 O, numerator of fm satisfies exp [E=
-Ambfy m+22Am
S 1, and the denominator of fm satisfies
IT (1+ + 2mm/mm)""/2 2 Note that the numerator monotonically decreasing with m and the denominator is monotonically increasing with m..
(b) fm pointwisely - For any W.
exp E
Ambb 2
lim Sm (w) exp (> (ir(Ew) + I ww) lim
(=1 (1 +1+ 2 ///mmmm /m) 12/2
> 1 m--00 22m
lim Ammb, oz
exp
exp (> (tr(Ziw) ) + p (w))
I1
lim (1 2dT /
exp (-> -- b9m))
exp (A (ir(Ew, ) + p Mwr (w))
= 1
II;-11 exp (211))
The last inequality uses the facts that Σ; b,
1 pw and Σ!-1 ?7 (r((ww)
(c) Je d= [So] Dc.
E [fo]
En exp () (tr(Ew) ] + p Hw)) ))
B exp (> [ur((W* (W" -W) 2.) +1(("* -W)p=1111
1
B exp
( W..' W) ir(E)) Σ (w*
-Ww*-Waante - Wi (W w
A Σ (W* w)in im [Ex + p=p22] (W* W)th + (r(EEe)
En exp
||(E2m + poll 1/2 (W* - (w)--W>(? + (m(Ees)
E, exp
1/2 (W* w))l)/) o + voll") +
exp (Atr(E.)) By exp
(a (x|l(22 + + 1/2 (W--w)ll oc is our
The last inequality holds because Exx exp assumption and exp (Atr(E.)) is a constant.
Since the conditions (a), (b) and (c) hold, by the dominated convergence theorem (Theorem 11.32. we have
44
lim / fm di m-+00 /
[.] =1 { lim fm dm ta== dd =
lim ES [mm]
m m0
Since In is continuous on (0, 00), we can interchange lim and In. Therefore,
lim V (A, m) 000
lim In Ex [fm]] In lim E, 1f In 1=1=0 m=+00
=00
16

---

Proof of Lemmat 4.2
R (()) E [llv Woll.a
Bllly: W
W R. 1x [ 2W Elyrod Ely
W 22 W 2w fn (Eiry)+ +i (Day)it
(Wi, 21/2
(W
2(W (. (0 + (Eng)it
(W, 11/2
(Ezz)
T T5-11/2 1/2 (W 5:1/
(Eary) YT, 12 (Fery) 1 + (F
2z
Σ
WW, in 23
(Eay !! 1i 1/2 - 1/2 20
(Ezy) )aal( + (Eny)i
||E!// Fm
Pイ 1
- || -n 1/2 113- + tr (Evy)
zz
Since assume Σgz is positive definite, Σ, exists.
Proof of Proposition
4.3
By Lemma
4.2
E (n))
-y -. || +
= F
<NllE!!-- 1err://--__-_ ST (M((EEg)) I. |لا / C. Ex (NI
= Ex
(20)
Proof of Lemma 4.4
Since 6.6 h and y (1 8) O h. we have
Σ E [x=-| E [8 h) o () h)²]
Exy E [ayy " = E. [6 @ ())
Eny = E [Vv"]] = E [((1-5) 5) 50 1)
o)
We first prove Egx. For i.j {1, 2, .....] n) with i + j. (Enz); = Elō,ö Note that R,6 5, a Bernoulli random variable (as its value can either be 0 or 1), 5, depends on h, and o, depends on hy. We have
E = P (6,66 = 1) =
= 1,44 1,,, P =P = 116, = 11, = -1
(h, = 1,h; = 1)
-6, = 1|h, = P6, 1|h, = 1) P 1,,= 1)
= [hayl = =p²(Enn);
(21)
For any i, (Enz)i E(((,,,,))). Using the property that Bernoulli random variable X has
E[x²] = E[X],
h(,)),,) = P (6 1) P(. = 1, = 1) P (6, = 1|,, = 1) P (hy = 1) = pE[h:]
= pEEhh -p( P(Enn)s
(22)
Combining 21 and 22 we get
Ege =p²2 M'nnn + pp1 (1 O En,)
(23)
17

---

Since (Ey )ij = E - 6:)(1 6, h and (Eng) )a B[((1 - 5,)h,)9. replacing p with 1
in 23 we get Eyg = (1 - E + p( p) ④ InN)-
Since (Erylig ( 8
Ep) and = 6.) 0 (Note that 6,(1 b.g 0 regardless of whether ò; is 0 or 1.), we have Eay - p(1 - p) I Ens)-
Note that in (33 if Σna is positive definite, then O Ena is also positive definite (Theorem 7.5.3 (b),
22 which implies that Eee is positive definite.
Proof of Theorem 5.2
(a) Since Ep [P] = U and E, E2|WTTT W| = 4Tu + diag (2 2k=1 SKI. Σi= SK2, ..... .k. Skn).
m
Ep R (W))
E W
Ep W Wx..l
m71
m
m
Bp((Y,, WX..) WX.))
Y.Y. - 2Y B. Xa1 + (WTT X.:
m
n
Σ -
X + x diag Σ Ski, Σ
S k2
Skn
Xx1
m
(24)
And D(pl(T) can be expressed as
1 D(p](T)= = 2
12)2
Σ
1 n
SI
Sety o2
1121 4oll
(25)
o²
k= 1=1
Denote o, A) = [W)] D pl xF.. ) our optimization problem becomes
min mig f(U, 0,
(26)
u.s
The optimal U/ and S are obtained by solving au fu, S/M, o, 0 and à su, SVo, o, A)
0.
First we show the partial derivatives of the 3D
term:
i x D(el(*)=
(4y - (44)))))
a 1 D(Pl(*)
-
1
Aa
as.
X
2X Sig oz
Then we show the partial derivatives of the Ep [Hemp (W)) term. By 24p. for any i,
a E, (e
θ 1
x diag Σ SkI
Ska
Σ Skn
Xx1
as,
aS. m
0 1
X
X X3.
aS m
m
0
6 1 m
Eunxi
E
-2 u X.1 + XT
-2YaX
WMg m
m
-Y
0
(UnxX.1))
(-2YnX + x(((x.x)(x.1)
(-Y + U4 Xx..)
Wrap up the above results, get
o5. Su, o,
_
X,.X7
1
1
1
(27)
m
2X
S
18

---

0
2 (-YiaX3+ + +U.xXX]))
(Ui (2o)is)
s(u 81,,, o, 0,
(28)
m
do?
Therefore, by
27
the solution of 0 fu,,,,,Uu o, A) 0 is that
as
Sij
for any i,j E (1,
(29)
22 X
+
By
28
we have
2 (-YX +uxx²)+ 222 (4 -20- T
f [U, SKu, o.
(30)
au
do
Thus the solution of 0 au
fu, SVL4, 0, is
1 L =
1 2o
_xx
L
u
(31)
2X
Now we show that f(4, SK,. or A) is a convex function, such that the solutions of S in 29 and 2 in (31 are the global minimizer of (26). By (27) and (28) we have
f²
if = =i=kkj- k =
f a
x X7
+
20
if = k,j =
2A(S
2 X X7.
if k,j # otherwise
0S
otherwise
Denote R2n² where for i 2, ..,, n and j 2, ..., nJ, M(i-1)n+t Uij and Mn +a
8 Sig. Let H E IR 222 x2n be the Hessian matrix where (Hş)iz Ou,āv;" Then can write
1 xo= T,2 and B is a n n² diagonal matrix with A 0 where A = 2 (XX") @ In +
H =
(( (1-1)+) 2x Here means Kronecker product.
The Kronecker product has property that, let {Aii== 1..., m] be the eigenvalues of P mx and (nylī ...., n) be the eigenvalues of Q E Ruxn then {Aauuli = 1,... m, 1, are the eigenvalues of P @ Q (Theorem 4.2.12, 22I). Since XX" is positive semi-definite and I, is Thus A is positive definite. Since all elements of S are positive. B is positive definite. Therefore, Hy is positive definite matrix for any 2/ and S. which means f(,, Si a, X) a convex function. Thus, the solutions of (9) and in 31 give the global minimum.
positive definite, (XXX7) @ I, positive
(b) Applying the constraint diag(W) = 0 to and implies taking diag(U) = 0, diag(S) = 0, diag(Ud) 0, and diag(o"J) In this case, (66 becomes a constrained optimization problem.
min fu4, a, X) u.S
s.t. diag(U) 0, diag(S)
(32)
Then we remove the constraint diag(S) by defining f as a function of only the off-diagonal elements of S. S- (Si i,j {1, ..., n), i + j}. then (32) is equivalent to
min f(u, S A, s.t. diag(U) = u.s
(33)
To solve 33
construct the Lagrangian function
"uu,s,x'/4, a, A) f(4, S 14, a, A) + xf diag(()
where r E R", and solve
ØL
dr [diag(U))
(34)
0L à
( or A) + diag(x) 0
(35)
au au
aL. 0
S Mo. a, A) = for E n)
(36)
as, as,
19

---

Since (36 is the i + case of 27 the optimal S is obtained by 29 with +
The optimal U is obtained by solving (35) and 34 By 35
2 (-YXT + +XXX) + Ao
(U 40) diag(v)
=0
u Lyx" + 2A0 2o0z Uo diag(z) )
1
(
-xxx
+ 2\
(37)
m
Then solve to satisfy 33
-1
1
1
diag(U)) diag
YXT
2202
+
diag
diag(r)
XXX
+
m
270
m
22
yx
1
(=xx"
1
-17
[xxx
-1
diag
2x0
+
x O diag 2
+
=
0
m
m
2Ao
m
we get
2.diag
L LYXT +
2o0241
1
_xx"
1
Ødiag
(=XXT
2
2A
35 and (36) gives the global minimum of the problem (33) Let Ht be the Hessian matrix of the Lagrangian L. It is easy to verify that, by removing the rows and columns of H, corresponding to S S22,...S..n to obtain H' E (22n"=)))(2n=== -n), we have
Now show that the solution of (34),
H2 H This shows that H is positive definite for any LI and S Hence, by the second-order sufficiency conditions (Section 11.5, (33]), any solution S .x)
satisfying (34), (35) and (36 is strict local minimum. Since the solution is unique, it is also strict global minimum.
Proof of Proposition 5.3
Denote v, = (W,, [W,>' I//2 + B,.)", and write e = A + μ² where e ~ N)0, Since A = STS. using the quadratic form shown in [T6), we have
||l - c (1) + +2) (A2/2 + / A5A -1 11))
= -Se + SS-112 Hd) A(SE+ A 1A-1/2
= (Se+ + A + b²) Σ ny + BB
j=1
It is easy show that each Syc i..dd from (0 1) for all j. thus each Sy.c b} is independently from N N(B5, 1). Since each 0; is independent, by 20 we have
E, (w)) CE A DE
I[a, E//2 BB.
C
Ex
1
=
E Any (S).ee833))
1=1
i=1 j=
A(B3)
cxp
1-2Am;
(1 2
The last equality above follows from 17
Proof of Theorem 5.4
Let P.Q E R xx be two symmetric matrices, write P > Q P Q is positive semi-definite, and write P > Q i Q is positive definite.
By Corollary 7.7.4 (c) of (22) if P Z Q. then (P(P) (Q) for any j, where 11; (P) and m (Q) denote the jth largest eigenvalues P and Q. respectively. Since A ( = o²(221/2) (>11211 > 0 for any i, we have 77 2 ni for any i.j.
(8)
20

---

Since 6(0) = (A - p1' we have
(8)
= 10 "1 = ( pz
Therefore, 13 can be expressed as
M(n)) (s5))) (S))))
exp
exp
1 Ex'
1-2An
1-2An5
C
i=1j=1
[1 11-2A1
2Am;
(1)
1/2
==1j=1
(1-2Xr))
1/2
(
exp
exp ( ((1)) (())))(()))00)
1-22n
(1 1 11-2201 1-2211²))
1/2
IT== 2Xn²z") (10
j=1
where ( = diag
1-2An/ (TT) 1-22n 1 1
Similarly,
12 can be expressed as
exp ST
E
C
T==1 (1 22m1;)
where Ā diag 2Am
2Anz
1-22mm
Now show that 5Th ÄS = (S()) (0) g(+) for any i. By Corollary 7.7.4 (a) of [22), If > 0 and
Q > 0, then P = Q if and only Q 2 P-1. Since we assume 0 < < 27p it follows that
1 1- 2-2 > 0 and 1 2An > 0 for any i. Thus, all diagonal elements of A and A are positive, implying that (0) > 0 and STH AA 0.
Since (( -1 = (T (1 - 2AA ()) s(i) = -244(0 and (5*AS) = - have
(15)))AA5(T) -1 = (S* + I 21A I 2AA = A A
Thus, STAS = (SO) (t) holds, implying that STA 2 (5)) any n¹. Therefore,
5(t) 4 for
Es (W)) C
1
exp (A(x') (S0))500
exp (T
1 En [CA="((ms]
I1m (1- 21-22nn;)) 1/2
II ===1 (1 2\11))
C
B Further Discussion on the Convergence of Shalaeva's Bound
Another convergence result by Shalaeva el al [45] is that: For any fixed p, such that D(pl T <
oo, take A m²/d for some constant d > 2, then the right hand side of Shalaeva's bound (See Section 2 converges to the left hand side m - oc.
2x
This result is based on the convergence of the upper bound Ew exp i
By taking X =
//d/ we have Ew exp
2A ?
i:: = Ew exp (2m2/d-1 ). In this case, the following
convergence statement
21

---

lim x [(ollT] + 1n + v. ( m)
s
lim m 1//d D( [ m=-00
x))
+ lim m 1/d In Be exp (2m (2m 2/ 1,
m
holds if
lim 1/d
m-1/d In Ew exp 22 22 2om2/d-112 (
(38)
Shalaeva et al. 445 provided only one condition d > 2 to ensure
38 and they did not discuss the choice of However, find that > 2 alone is not sufficient to guarantee convergence, since (38) does not hold for all т. A few examples are given below:
Example B.I. I is distribution with bounded support, then (38 holds. This is because there exists constant > 0 such that ||W]|2 < G. We can show that for any À > O, Ew exp (Ao <o:
Ew-r exp (Aui Ew-= [exp (A(o2) w112 + -2)
< E [exp (N(o ||]ll22 +
BT [exp (A ( G)² + + = exp (I 112 G)2² oz))) < 00
Thus, when d > 2.
lim m 11dd In Ew exp 12m 2/d (2m 2/1-11, lim 1/d In EEw ~a exp (2v 0
Example B.2. Let क be a Gaussian distribution. We show that 38 does not hold.
We first show that Ew exp (Yu = oc for any X > 0. Denote = W W E R xm where W n. then a Gaussian random vector. Thus
Ew exp (Ae = Ew-= exp (A(o2) W11 + +o?))) 2 E eexp (A(o W W1
)) ]
= Ew /exp (Ao wl1l = E.
cxp
oo
2 Ew /exp (Ao=
Here w1 E is the first element of w, which is Gaussian random variable. Suppose w1 ~ NN o?), then
(-(w) (wi 1001-00 14)
En
[exp (Ao
exp (No
exp
dwi
V2TO
2
=
(Wi - 11)2
exp
dwy 00
1TT
202
because Ao [w] 11)
00 as W1 + 0c. Hence, Ew exp (Av = oo for any 0.
2x
If (38) holds, then by the definition of limit, for any E > 0, there exists a finite integer M such
that for any m > M, m -1/d In Ew exp (2m2/d-1 ww < e. The negation of this statement that, there exists e such that for any finite integer M, there exists > M satisfying
m m-/// In Ew exp (2m2/d-1,²? b2 2 e. Let 1, M be any finite integer, and = M + 1. Then
En cxp (2m2/d-1,22 negation is true, and (38) does not hold.
oo and c S 1ld In Ew exp (2m2/d-1 2 co. So the
In summary, the validity of 388 depends on म, and sufficient condition for convergence requires that x be appropriately specified.
C Allowing Multiple Trails on
Finding the optimal À that yields the tightest bound is nontrivial. As suggested in Section 2.4 of 4 we approximate the optimal À by searching
finite grid A (21, ) where each
A and L denotes the cardinality of A.
P VA E A,V,, Ew-p(R (T (W)] < Bw-p n0 (W)) 1 X (( in 2/f v...(hm))
4
22

---

This is because
P Yy E V., Bw-p! (W]) < Ew- (It" (W) D(pl(=) + lu L v vn.
=1 =1-P - P
EAEA,3 Bw [1" (W)) Bw [1| (W)] + x
1
D (t) +
In 4
+ (n.
L
1
=1-P
3p, ( Bar Bw R _ (W))
x
\1=1
L
21
Σ
3p, Ew F. (W]] Ew R (W))
D(ll( vnn.(ds,,)))
X
i=1
21
8
The grid search is a standard method for optimizing À and is also used by Dziugaite and Roy [S] Note that Ả cannot be directly optimized via gradient descent while keeping P and æT, as this would make the optimal Ả depend on the dataset S, which is random variable 43). Consequently, the optimal Ả would itself become random variable, contradicting the assumption that it is fixed.
D Related Works
In statistical learning, generalization bounds commonly provide an upper bound on the generalization gap, R. (W)- Remp (W) (or, in the two-sided case, [R)) - Reme (W)l. thereby estimating the true risk Rume (W) for any given model W. A generalization bound called a Probably Approximately Correct (PAC) bound if it guarantees, with probability least - 6, that the gap does exceed a small e. It is PAC-Bayesian if. in addition, the model W is treated as a random variable and KL-divergence term, D(p || ). is introduced allow Bayesian inference between the prior and posterior of W. Classic PAC bounds typically estimate the worst-case generalization gap, supw { {Rt (W) Remp (W))
(1I) which often becomes loose as the number of parameters in W increases (39 41). In contrast, PAC-Bayes bounds estimate the expected generalization gap, Ew Switiiee (W) Demp (W)). which is typically much tighter than the worst-case bound
Early PAC-Bayes bounds, including those by McAllester 37), Catoni 9. Langford and Seeger 29
primarily focus binary or bounded losses. In recent years, PAC PAC-Bayes bounds have been extended to unbounded losses, as in Alquier et al. [5], Haddouche et al. [19]. Haddouche and Guedj [8]. Rodríguez-Gállezz al. [43. Casado et al. S These works typically assume that the loss follows light-tailed or heavy-tailed distribution. However, such assumptions usually rely on oracle knowledge of the underlying data distribution; consequently, these bounds often oracle bounds and can only be computed when the true data distribution is known.
Linear regression typically employs the squared loss, which is unbounded. Several PAC-Bayes bounds targeting the generalization gap have been proposed for linear regression. One major line of work adapts the unbounded squared loss the framework of Alquier et al. [5). such as Germain et al. [77 and Shalaeva et al. [45], which we follow here. Other studies, such as Haddouche et al. [19], consider the E loss for linear regression instead of the standard &g (squared) loss. To the best of our knowledge, all these bounds focus on the single-output setting and have not yet been extended to the multivariate case.
Another class of generalization bounds for linear regression targets the excess risk H pit"(W"*) where W' denotes the minimizer of the true risk. The advantage of this formulation is that, when the tail distribution the loss satisfies the Bernstein assumption (Definition 4.1, [4],, such bounds (including both PAC and PAC-Bayes types) typically achieve 1/m convergence rate, which faster than the 1/ m rate commonly observed for PAC- WAC-Bayes bounds targeting the generalization gap (Section 4.2, (4). Several bounds for single-output and multivariate linear regression follow this setting, including those by Alquier and Bieu (3]. Alquier (2). and Mai (34. However, these bounds do not depend on the empirical risk, making them less practical for real-world applications.
Collaborative filtering can be viewed as special case matrix completion problem, i.e., predicting missing values in matrix (Section 1.3.1.2, [[). In recent years, LAEs have become popular
23

---

model for collaborative filtering due to their simplicity and effectiveness. Unlike other models, LAEs have the distinctive property that their training objective, such as in 4 resembles constrained linear regression problem where the input and target matrices are the same, highlighting their close relationship with linear regression. The earliest LAE model can be traced back to SLIM [4), followed by models such as EASE [48], EDLAE [49] and ELSA 50). All of these models introduce a zero-diagonal constraint on the weight matrix W, preventing items from learning themselves and thereby avoiding overfitting toward the identity. A recent study by Moon et al. (38) shows that this zero-diagonal constraint can be slightly relaxed to a diagonal with small bounded norm, potentially improving performance.
While the generalizability of LAE models remains unexplored, related work has investigated general- izability in matrix completion for collaborative filtering, including Srebro et al. [47], Shamir et al. 46). and Foygel al. [16] Other studies focus on the generalizability of general matrix completion, such as Candès and Tao 77, Recht 422, Ledent et al. [1]. and Ledent and Alves [30] Our work differs from these studies in that we analyze LAEs mainly from the perspective of linear regression rather than matrix completion. In addition, Variational Autoencoders (VAEs) are another type of collaborative filtering model [32], and PAC-Bayes bounds have been developed for VAEs [10].
should be noted that VAEs and LAEs are fundamentally different models, so the PAC-Bayes bounds for VAEs are not directly comparable to our bounds, as discussed in Appendix E.3
Therefore, to the best of our knowledge, we propose the first PAC-Bayes bound for multivariate linear regression targeting the generalization gap, and the first PAC-Bayes bound for LAEs.
E
Conclusions and Discussions
This paper studies the generalizability of multivariate linear regression and LAEs. We first propose a PAC-Bayer bound for multivariate linear regression under a Gaussian data assumption, extending Shalaeva's bound for single-output linear regression, and establish sufficient condition that guaran- tees convergence. Next, we build the connection between multivariate linear regression and LAE models by introducing a relaxed MSE as an evaluation metric, under which LAE models can be interpreted as multivariate linear regression on bounded data, subject to zero-diagonal constraint on weights and a hold-out constraint on the input and target data. This connection allows us to adapt our bound for multivariate linear regression to LAEs.
In practice, LAEs are typically large models evaluated on large datasets, which makes computing the tightest bound inefficient. To address this, we develop theoretical methods to improve computational efficiency. Specifically, by restricting both the prior and posterior to be Gaussian, we obtain an efficient sub-optimal bound in dosed-form. We then address the computational cost imposed by the zero-diagonal constraint by establishing and computing an upper bound with reduced complexity. Experimental results demonstrate that our bound is tight and correlates strongly with practical metrics such as Recall @ K and NDCG@K K, suggesting that it effectively reflects the real-world performance of LAE models.
Below are the discussions of our work.
E.1 Limitations
One limitation of our work lies in Algorithm I which takes Σna = (a-M(Hhi") as input and requires to be known, which is only possible when M is non-oracle. Whether M is non-oracle depends on how the dataset is modeled statistically. If the dataset H is drawn from a meta-dataset H that is known and of fixed size, we may assume M to be the population distribution of this meta-dataset, thereby making M non-oracle. However, this assumption is rather restrictive: In most real-world scenarios, such a meta-dataset may not exist, as data are continuously collected or expand over time. In this case, M is typically modeled as an unknown oracle distribution to account for unseen data. Therefore, whether M is non-oracle ultimately depends on how the dataset is modeled, and our work shows that the bound is at least computable in restricted scenarios where datasets can be represented by a non-oracle M.
To avoid introducing oracle M, one may also consider reconstructing the bound by applying
linear regression to empirical PAC-Bayes bounds such as those of McAllester [37], and Langford and Seeger (29), rather than to the Alquier's oracle bound [5]. However, this approach is not feasible
24

---

because empirical PAC-Bayes bounds are typically derived under the assumption of bounded loss, whereas the loss in linear regression is unbounded since W itself is unbounded in | W=;ll?-
Moreover, PAC-Bayes bounds for unbounded losses are inherently oracle, as they are derived from tail-distribution assumptions that are themselves oracle in nature. These distributions are usually characterized by bounded exponential moment Ex-~ [e>(Ex=><(X)) the Hoeffding assumption [5] and the bounded Cumulant Generating Function (CGF) assumption 43 19] Assuming these quantities are bounded implicitly presumes access to oracle information
for some À R, in
about M, which is rarely realistic in practice.
Another limitation is that the computational methods introduced in Section 5 are primarily designed for full rank or nearly full-rank L.FSs. By 'nearly-full rank', we refer to matrices rank n formed by applying the zero-diagonal constraint to a full-rank X matrix. Some LAE models, however, impose low-rank constraints [49] on W, where a W of rank &o (k < n) can be decomposed as W UV with U and vl*= Assumption 5.1 does not hold in this case, since requires W to be a random Gaussian matrix, implying full rank [14]. Consequently, results that rely on this assumption, including Theorems 55.2 and [5.4] are not applicable to low-rank LAEs. One possible approach to adapting our bound for low-rank W is to impose distributional assumptions on U and V so that any realization of W UV is always low-rank; however, this contradicts Assumption 5.T
E.2
Comparison with Mai's Excess Risk Bounds
3
Mai proposed generalization bounds for bilinear regression (Theorem 2, (34 and matrix completion (Theorem 3, (34) based excess risk. Their bilinear regression setting is defined follows: Given two input matrices E R and X E R=x", let Y E Roxm be the target matrix following a distribution conditioned where E E IR is a random noise matrix whose entries Eij has zero mean and finite variance. If taking r = p and 7 T, it reduces the multivariate linear regression presented in Section 2
Z and X, then there exists W R such that Y ZW X + E,
Given a model W RPxm, let the true risk be defined as R. (W)- = p. Eyz,x(1 ZW X 11:.) then the excess risk can be expressed as
1
(ttue()) - R (W*) = IZWX - ZW
pm
If we further suppose that W is a random variable following distribution أ. and define the expected excess risk Ew-pLR (W)] Rime (W*), then by Jensen's inequality we have
1
E [Tttr (W)]- 1
(Bw->(ZWXi, - ZW* XII
(39)
pm
Mai first derived a generalization bound for bilinear regression by establishing an upper bound on They then adapted this bound to the matrix completion setting, which
[|Ew~>[ZW] X] - ZW* XI2 - assumes that (A < pm) pairs in {[ZZW* (WWx")y, X).. Y(y) are observed, where the trained model W is used to recover the remaining pairs. Further, their bound converges linearly with respect to k.
It is well known that if the loss satisfies Bernstein assumption (Definition 4. [4), can construct an upper bound on the excess risk (or the expected excess risk) that converges at a linear rate of 1 m.. By (39). Mai's bound is indeed upper bound on a lower bound of the excess risk, which makes its linear convergence rate reasonable.
In contrast, our bounds are based on the generalization gap T (W)= namp (W), rather than the excess risk. PAC-Bayes bounds on the generalization gap typically converge at a slower rate of 1//m (Section 4.2, (4). The Bernstein assumption, which enables linear convergence for excess risk bounds, is not applicable in the setting of the generalization gap.
E.3
Comparison with PAC-Bayes Bounds for Variational Autoencoders 10
Like LAEs, Variational Autoencoders (VAEs) are another class of autoencoders that have been applied to collaborative filtering (32. Recently, PAC-Bayes bounds for VAEs were proposed by Abdella-- et al. [[0].
Although LAES and VAEs are both autoencoders, they differ fundamentally in model architecture, learning objectives, and loss formulations, making our bound for LAFS not directly comparable to the -Byes bound for VAEs.
25

---

In detail, the architecture of VAES 10 is:
encoder Bo (x|px) latent code z
decoder po (x]2)
input x with the loss defined as (x [og po (x|22)) And the architecture of LAFs is
prediction
target
encnder and decnder W. prediction Wx
input x
target y
with the loss defined as the relaxed MSE lyy W
Wx 1 in our framework, see Section 4.2
The differences between VAEs and LAEs are summarized as follows:
Difference in Architecture: A VAE consists of an encoder and a decoder and explicitly involves the latent code z. In an LAE,, the model W serves as both the encoder and the decoder, with the latent code being implicit. While one could decompose W AB and treat B as the encoder and as the decoder, such a decomposition is uncommon in recommender systems. Most LAE recommender models do not decompose W [40 48 49
50]; instead, they focus on studying W directly by introducing constraints such as a zero diagonal. As result, LAEs and VAEs are generally regarded as two distinct model classes.
identical. In contrast, for LAEs, the input x and target y are the same during training but differ during evaluation due to the hold-out constraint. Our bound is defined for evaluation, where x and y are indeed different. This fundamental difference in framework.
Difference in Input and Target: VAEs require the input and target
configuration makes LAFs incompatible with the VAE
Difference in Loss: The VAE loss L (=/xx) [log pa (x [z)] aims to minimize the mismatch between the encoder go (z(x) and the decoder Pa (x) while the LAE loss l Wx }. aims to minimize the mismatch between the target y and the prediction Wx. These losses reflect fundamentally different learning objectives, so the two models cannot share the same loss function. Moreover, the LAE loss aligns more closely with a multivariate linear regression loss than with the VAE loss, which is why our bound for multivariate linear regression can be naturally extended to LAE models.
Therefore, due the fundamental differences between VAEs and LAEs, our bound for LAEs is not directly comparable to the PAC-Bayes bound for VAEs proposed by Chérief-Abdellatif et al. [0]
E.4 Relationship between MSE and Ranking Metries
While minimizing MSE is not theoretically consistent with achieving optimal ranking, empirical studies have observed strong correlation between reductions in MSE and improvements in ranking metrics [26]. Moreover, ranking metrics are typically set-based, discrete, and non-differentiable, making them difficult to optimize directly using gradient-bare become standard practice collaborative filtering to employ regression-styyee or reconstruction-basee surrogate losses most commonly the squared error as training objectives (Section 2.6, (]), which resemble the MSE used during evaluation. This approach underlies many successful algorithms, including SLIM 40), EASE [48] EDLAE [4). ELSA (50), and matrix factorization models (33 27).
methods. Consequently, it has
which consistently demonstrate strong empirical performance on ranking tasks despite optimizing non-ranking loss.
Since LAES typically use a linear regression loss for training and ranking metrics for evaluation and model performance is ultimately measured by the latter derive our generalization bound solely with respect to the evaluation metric. As discussed in Section 4.2) because ranking metrics are difficult to analyze statistically, we instead adopt MSE as the evaluation metric, which follows the form of a linear regression loss but differs from the training loss in both definition and purpose.
E.5 Broader Impacts
This work advances the theoretical foundations of machine learning by introducing the first PAC- Bayes bound for multivariate linear regression targeting generalization gap. extending beyond single-output regression to handle multiple dependent variables simultaneously. This establishes new generalization guarantees for structured prediction, multi-task learning, and recommendation systems.
26

---

Additionally, we identify and correct limitation in an existing PAC-Bayes proof for single-output linear regression, further strengthening the theoretical foundation of regression analysis.
Building on this, we apply our bound to LAFS in recommendation systems, delivering their first rigorous generalization analysis. Our approach accounts for key structural constraints, such the zero-diagonal weight requirement, ensuring applicability to models like EASE and EDLAE.
Beyond theory, our work has direct practical implications for model evaluation and selection. Our bound provides a post-training diagnostic tool for assessing the generalization of any LAE model, regardless of its training process. While not directly guiding training hyperparameter tuning, smaller PAC-Bayes bound suggests better generalization on unseen data. Empirical results confirm that our bound remains within a reasonable multiple of the test error, offering reliable probabilistic estimates of true risk independent of training error.
Our work focuses on theoretical generalization analysis and poses no immediate ethical risks. How- ever, recommendation systems shape content exposure and user behavior in domains like e-commerce and social media. Strengthening generalization theory alongside other recommendation criteria may help mitigate bias, enhance fairness, and improve trust Al-driven systems.
Other Supplemental Materials
F.1 Details of Algorithm
Here we provide details on the computation ( (W)), in Algorithm | which are not fully described in the main paper.
and D ||x)
Given À and x ( o²I). f). the optimal P = (u, S) that minimizes the right hand side of 14
is obtained by Theorem 5.2 Once is obtained, we can compute (.wpp[Rom (W)) by (24). which
can be simplified as
_ || -x m
n-1 diag(S m
Ew-p! 10 [W)]
(40)
The n - term in (40) due to the zero-diagonal constraint, which enforces diag(S) 0. Without this constraint, the term becomes instead of n 1.
Similarly, by (6). Ew [Atrr (W)) can be expressed as
E Ew-pl mim (W) = | E u
13 + (n 1)Il diag(S diag(Si.)11/' 1/25_/1111
ur( (ng) -
(41)
The derivation of (41 analogous 24 by substituting Y with 2T 1/1 and X with 2-1/ 1/2
The D(pl[=) x) term under zero-diagonal constraint is obtained from 25 by removing the diagonal elements of
1 D(p[(T)=
Sat 02
11 1
(n² (22-n)(210) m2)
Σ
(In
(42)
2
1=1,4
F22
Dataset Description
The following table shows the details of the datasets used in the experiments.
Table 2: Dataset description
Dataset
ML 20M
Netflix
MSD
#users (m)
138493
480189
1017982
#items (n)
26744
17770
40000
#interactions
2000263
100480507
33687193
27

---

F.3 Details of the Results in T
The following table presents the detailed values of the components of the RH terms in Table [] illustrating how the results were obtained. This information may be helpful for reproducing the experiments.
Table 3: Details of the terms of each RH in Table
PAC-Bayes Bound for LAFS
Models
MI. 20M
Netflix
MSD
512
512
512
Ew u Nomp (W))
66.99
90.87
16.58
Y 50
D ((plTT)
0.28
0.18
0.0019
In ET
(W)
31571.14
44659.37
8196.30
512
512
512
Ew- [Reaa (W)]
65.14
89.68
16.34
7 = 1010
D ((pl TT)
0.27
0.17
0.0018
In E8, (w))
311.2353
44313.39
8141.72
512
512
512
Ew Pomp (W))
63.59
88.57
16.12
7 = 200
D P(plT ))
0.26
0.17
0.0018
E (W)
30753.19
44014.86
8092.62
512
512
512
Ew Reup (W))
61.93
87.26
15.89
7 = 500
D p In EE, an (W)
0.23
0.17
0.0016
30444.39
43703.53
8044.64
512
512
512
Ew~pp Homp (W)
60.96
86.42
15.79
7 1000
D p|| In E2 eAn" (W)
0,23
0.16
0.0016
30310.47
43522.96
8033.10
512
512
512
Ew [Rop (W)
60.23
85.71
15.78
7 2000
D ppl] In E,, An (W)
0.22
0.15
0.0015
30255,43
43382.46
8052.90
512
512
512
Ew pemp (W))
59.70
84.97
15.88
7 5000
D p
0.20
0.14
0.0014
E,
en (W)
30308.30
43255.35
8128.97
28

---

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:
We have accurately stated our contributions in Abstract and Introduction.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in 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 perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results
be expected to generalize other settings.
is fine to include aspirational goals as motivation as long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification:
We have discussed the limitations of our work in Appendix
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in 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 the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low images are taken in low lighting. Or speech-to-text system might 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 address problems of privacy and faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
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]

---

We provide full set of assumptions (Assumption 3.1 and 5.1. and all proofs are presented in Appendix A We have aimed to make the proofs both complete and easy to follow.
Justification:
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 properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: We have disclosed all the information needed to reproduce the main experi- mental results in Section and Section G
Guidelines:
The answer NA means that the paper does not include experiments. If the paper includes experiments, a No answer 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.
the contribution is a 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 a novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, 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 often one good way to accomplish this, but reproducibility can also be provided via detailed 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, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require 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 (e.g., a large language model), then there should either be way access this model for reproducing the results way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reprodueibility may be tricky in some cases, in which authors are welcome to 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 (e.g., to registered users), but it should be possible for other researchers have some path reproducing verifying the results.
5. Open access to data and code
30

---

The authors should provide scripts to reproduce all experimental results for the new reproducible, they Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification: We have provided the code and the links to the open datasets (MovieLens 20M, Netflix, MSD) in the supplemental material.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.c cc/
for more details.
While we encourage the release of code and data, we understand that this might be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.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
/ https/cccppblligggudeeeeeeeeeeee Coaee///ddeeebmmmcia ssi onPolicy) for more details.
The authors should provide instructions data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
proposed method and baselines. only subset of experiments 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 in supplemental material (appended to the paper) 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:
We have specified all experimental settings and details Section 6
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can provided either with the code, in appendix, or as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [No]
Justification: Our models are obtained from closed-form solution, which is deterministic and does not involve random emrrs..
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
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 with given experimental conditions).
31

---

The method for calculating the error bars should be explained (closed form formula, call 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 or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma enror 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 or figures symmetric error bars that would yield results that are range (e.g. negative enror rates).
If error bars reported in tables plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables the text.
8. Experiments compute resources
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?
Question:
Answer: [Yes]
Justification: We have provided information on the computer resources in Section O and supplemental material.
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 as 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 it 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:///rrr....cccoooo./c/tttaaa
Answer: [Yes]
Justification: We have reviewed NeurlPS Code of Ethics and ensured that our research conforms to it in every respect.
Guidelines:
The answer NA means that the authors have not reviewed the NeurlPS 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 sure to preserve anonymity (e.g., there a special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification:
We have discussed the broader impacts in Appendix
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
32

---

Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), faireess considerations (e.g., deployment 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 particular applications, let al.one deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate 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 is 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 (c.g., gated release models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of MI).
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 work is mainly theoretical and is unlikely to involve risks potential misuse in application.
Guidelines:
The answer NA means that the paper poses no 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 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: We use open datasets (MovieLens 20M, Netflix, MSD) in accordance with their licenses, and we do not use existing code or models.
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.
The authors should state which version of the asset used and, 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 source should be provided.
33

---

If assets are released, the license, copyright information, and terms of use in the thcode. has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
package should be provided. For popular datasets, paperswi
latasets.paperssi
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 not available online, the authors are encouraged to reach out to the asset's s creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification:
We have included a README file as documentation along with the code.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moodl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing 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 details about compensation (if any)?
Answer: [NA]
Justification: Our paper 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 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 NeurIPS Code of Ethics, workers involved in 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 equivalent approval/review based the requirements of your country institution) were obtained?
Answer: [NA]
Justification: Our paper does not involve crowdsoureing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending 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.
34

---

We recognize that the procedures for this may vary significantly between institutions and locations, and we 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 is an 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: We use LLMs for writing purposes and for generating code unrelated the core methodology. All code related to the core methodology is written by the authors.
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 our LLM for what should or should not be described.
35


Reviews:
[{"review_id": "DPeUMbyxbd", "reviewer": "Reviewer_npXJ", "summary": "The authors develop PAC-Bayes bounds for\n1) multivariate linear regression problems (that is, when the output variable Y is a vector),\n2) linear autoencoders (LAE).", "questions": "1) Please compare the results to the existing works on multivariate regression and LAE/VAE's.\n2) Explain what are the improvements over these results.", "limitations": "Yes.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 1, "originality": 1, "strengths_and_weaknesses": "General comments:\n- the paper seems technically correct, and is overall well written,\n- but while authors know the basic PAC-Bayes techniques, they ignore the existing literature on PAC-Bayes bounds for multivariate regression and for VAEs.\n\n1) Multivariate regression:\nCovered as a special case of https://www.mdpi.com/1099-4300/25/2/333 (take Z= Identity),\nAlso included in https://arxiv.org/abs/2206.08619\nOne observation: the rates in the above references are in 1/n (n = sample size), while all rates derived in this submission cannot be better than 1/sqrt(n).\n\n2) LAEs:\nThey should compare with existing bound on VAEs:\nhttps://proceedings.mlr.press/v151/cherief-abdellatif22a.html\nhttps://proceedings.neurips.cc/paper_files/paper/2023/hash/b29500824d22ee9bbd25e4cd97c49b55-Abstract-Conference.html\nIt's not OK not to cite these papers, and the results stated here cannot be assessed if there is no rigorous comparison with these existing works.\n\nMinor comments:\n- the claim that Alquier et al.'s 2016 bound cannot be used in practice because Psi depends on the unknown data-generating process is a little weird, as the paper provides many examples where we can upper bound Psi by distribution independent quantities.", "was_revised": true, "final_justification": "The authors included a fair discussion with the existing papers on the topic. As a consequence, the contributions of the paper are more clearly visible now. I think this is a reasonably good contribution to the literature on generalization bounds for VAE/LAEs and as such, I think it can be published in NeurIPS.", "rebuttal": "Thank you for your thoughtful and constructive feedback. We especially appreciate your comments regarding the existing literature on PAC-Bayes bounds in the context of multivariate regression and variational autoencoders (VAEs). These works are indeed highly relevant, and we will ensure they are appropriately cited and discussed in our revised manuscript.\n\nIn the rebuttal below, we clarify how our contributions relate to these prior efforts. While there are conceptual connections, the prior works address different problem settings and rely on assumptions or formulations that are not directly applicable to our framework. Through this comparison, we highlight the novel technical aspects and the distinct contributions of our approach; and we hope this adequately addresses your concerns regarding the positioning and originality of our work within the broader landscape of existing research.\n\n> Multivariate regression: Covered as a special case of [2] (take Z= Identity), Also included in [3]. One observation: the rates in the above references are in 1/n (n = sample size), while all rates derived in this submission cannot be better than 1/sqrt(n).\n\nThank you for highlighting this relevant literature. Mai's bound [2, 3] indeed address the multivariate linear regression setting and use Bayesian techniques in deriving the bound. \n\nThe focus and formulation of their bounds differ fundamentally from ours: Mai's bound is for the excess risk $R^{\\text{true}}(W) - R^{\\text{true}}(W^\\*)$, the difference of true risks between an arbitrary model $W$ and the optimal model $W^\\*$ -- such bounds are often referred to as ``non Bayesian'' PAC bounds and typically achieve a $1/n$ rate under the Bernstein assumptions [1, Section 4.2; 10]. Our bound is for the generalization gap $R^{\\text{true}}(W) - R^{\\text{emp}}(W)$, the difference between the true and empirical risk of the same model $W$ -- such bounds are standard PAC-Bayes bounds and typically have a $1/\\sqrt{n}$ rate [1, Section 4.2]. These two types of bounds estimate different quantities and are therefore not directly comparable.\n\nIn detail, we explain how Mai's bound satisfies the Bernstein assumption. By [1, Definition 4.1], let $l(W)$ denote the loss function with the weights $W$ and $S$ be the dataset, we say the *Bernstein assumption* is satisfied if there exists a constant $K$ such that \n\\begin{equation*}\n\\mathbb{E}\\_{S}[(l(W) - l(W^\\*))^2] \\le K(R^{\\text{true}}(W) - R^{\\text{true}}(W^\\*))\n\\end{equation*}\n\nIn the proof of Lemma 3 of Mai's bound [3], they define $l(M) = \\sum\\_{i=1}^n (Y\\_i - \\prod\\_C(XM)\\_{\\mathcal{I}\\_i})^2$ and $l(M^\\*) = \\sum\\_{i=1}^n (Y\\_i - (XM^\\*)\\_{\\mathcal{I}\\_i})^2$, and showed that under Assumption 1 and 2 in [3], there exists constant $C\\_1$ such that \n\\begin{equation*}\n\\mathbb{E}\\_{S}[(l(W) - l(W^\\*))^2] \\le nC\\_1(R^{\\text{true}}(W) - R^{\\text{true}}(W^\\*))\n\\end{equation*}\nwhich implies the Bernstein assumption with $K = nC\\_1$.\n\nHowever, the Bernstein assumption is not applicable to our bound because it only applies to the excess risk $R^{\\text{true}}(W) - R^{\\text{true}}(W^\\*)$, whereas we focus on the generalization gap $R^{\\text{true}}(W) - R^{\\text{emp}}(W)$. Therefore, the rate of our bound cannot be improved in the same way as Mai's bound.\n\n> LAEs: They should compare with existing bound on VAEs [4]. It's not OK not to cite these papers, and the results stated here cannot be assessed if there is no rigorous comparison with these existing works.\n\nThank you for pointing out this related line of work. We will certainly cite [4] in our revised manuscript and expand our literature review to include a discussion of its relevance. However, we would like to clarify that the PAC-Bayesian bounds developed for VAEs in [4] are not directly comparable to our bounds for LAE models, due to fundamental differences in model architecture, learning objectives, and loss formulations.\n\nTo make these distinctions precise, we summarize the key differences below:\n\nThe architecture of VAE [4] is: \n\\begin{equation*}\n\\text{input}\\\\;\\mathbf{x} \\xrightarrow[]{\\text{encoder}\\\\; q\\_{\\phi}(\\mathbf{z}|\\mathbf{x})} \\text{latent code}\\\\;\\mathbf{z} \\xrightarrow[]{\\text{decoder}\\\\; p\\_{\\theta}(\\mathbf{x}|\\mathbf{z})} \\text{prediction} \\quad\\approx\\quad \\text{target} \\\\;\\mathbf{x}\n\\end{equation*}\nwith the loss defined as $-\\mathbb{E}\\_{q\\_{\\phi}(\\mathbf{z}|\\mathbf{x})}[\\log p\\_{\\theta}(\\mathbf{x}|\\mathbf{z})]$.\n\nThe architecture of LAE is \n\\begin{equation*}\n\\text{input} \\\\;\\mathbf{x} \\xrightarrow[]{\\text{encoder and decoder}\\\\; W} \\text{prediction}\\\\; W\\mathbf{x} \\quad\\approx\\quad \\text{target}\\\\;\\mathbf{y}\n\\end{equation*}\nwith the loss defined as $\\|\\|\\mathbf{y} - W\\mathbf{x}\\|\\|\\_F^2$ in our framework.\n\nThe differences between VAE and LAE are as follows: \n\n**Difference in Architecture**: VAE consists of an encoder and a decoder and explicitly involves the latent code $\\mathbf{z}$. In LAE, the model $W$ serves as both the encoder and the decoder, with the latent code being implicit. While one could decompose $W = AB$ and treat $B$ as the encoder and $A$ as the decoder, such a decomposition is uncommon in recommender systems. Most LAE recommender models do not decompose $W$ [6, 8, 9]; instead, they focus on studying $W$ directly by introducing constraints such as a zero diagonal. As a result, LAE and VAE are generally regarded as two distinct models.\n\n**Difference in Input and Target**: VAE requires the input and target to be the same. In contrast, the input $\\mathbf{x}$ and target $\\mathbf{y}$ in LAE are different. This difference stems from the data splitting strategy described in Lines 201 – 208, known as *weak generalization*, a widely used evaluation method in recommender systems [5, 6, 7]. This fundamental difference in input–target configuration makes LAE not adaptable to VAE.\n\n**Difference in Loss**: The VAE loss $-\\mathbb{E}\\_{q_{\\phi}(\\mathbf{z}|\\mathbf{x})}[\\log p\\_{\\theta}(\\mathbf{x}|\\mathbf{z})]$ aims to minimize the mismatch between the encoder $q\\_{\\phi}(\\mathbf{z}|\\mathbf{x})$ and the decoder $p\\_{\\theta}(\\mathbf{x}|\\mathbf{z})$, while the LAE loss $\\|\\|\\mathbf{y} - W\\mathbf{x}\\|\\|\\_F^2$ aims to minimize the mismatch between the target $\\mathbf{y}$ and the prediction $W\\mathbf{x}$. These different losses reflect fundamentally different learning objectives, meaning the two models cannot share the same loss function. Moreover, the LAE loss is less similar to the VAE loss and is instead more closely related to a multivariate linear regression loss, which is why our bound for multivariate linear regression can be naturally extended to LAE models.\n\nTo sum, VAEs and LAEs serve fundamentally different modeling purposes. Their loss functions and assumptions are not aligned. Consequently, their PAC-Bayesian analyses address incompatible questions. We will include a detailed discussion of [4] in the revised manuscript, along with a clear articulation of these distinctions.\n\n> The claim that Alquier et al.'s 2016 bound cannot be used in practice because Psi depends on the unknown data-generating process is a little weird, as the paper provides many examples where we can upper bound Psi by distribution independent quantities.\n\nThank you for this thoughtful question. We appreciate the opportunity to clarify our claim regarding the practical applicability of Alquier et al.'s bound.\n\nAs shown in Line 264 of our manuscript, the upper bound of $\\Psi$ used in our analysis is $ \\ln \\mathbb{E}\\_{\\pi}\\left[ e^{\\lambda R^{\\text{true}}(W)} \\right]$,\nwhere $R^{\\text{true}}(W)$ denotes the true risk and depends on the data distribution $\\mathcal{D}$. Therefore, it cannot be computed unless $\\mathcal{D}$ is known, i.e., a non-oracle distribution.\n\nIn our experimental setting, we work with a fixed dataset and treat the empirical distribution as the population distribution. Under this assumption, $\\mathcal{D}$ is non-oracle, and the bound is computable. Therefore, Alquier’s bound is applicable in such *controlled or retrospective scenarios*, where the entire dataset is observed and finite.\n\nHowever, in more general and realistic settings, such as online or streaming learning, data is collected continuously, and its underlying distribution remains unknown. In these scenarios, it is more appropriate to model the data as being drawn from an unknown oracle distribution $\\mathcal{D}$, rendering $\\Psi$ uncomputable and the bound inapplicable.\n\nTo be clear, we do *not* claim that Alquier’s bound cannot be used in practice. Rather, as discussed in Lines 302--309 and 638--646, we argue that its applicability is limited to specific cases where the data distribution can be assumed known or accurately approximated. We will revise the manuscript to better reflect this nuance and avoid any potential overstatement of our critique.\n\n**References:**\n\n[1] Alquier et al. User-friendly introduction to PAC-Bayes bounds. Foundations and Trends in Machine Learning, 2024.\n\n[2] The Tien Mai. From bilinear regression to inductive matrix completion: a quasi-Bayesian analysis. Entropy, 2023.\n\n[3] The Tien Mai, Pierre Alquier. Optimal quasi-Bayesian reduced rank regression with incomplete response. arXiv:2206.08619, 2022.\n\n[4] Chérief-Abdellatif et al. On PAC-Bayesian reconstruction guarantees for VAEs. AISTATS, 2022.\n\n[5] Jaewan Moon et al. It's enough: Relaxing diagonal constraints in linear autoencoders for recommendation. SIGIR, 2023.\n\n[6] Harald Steck. Embarrassingly shallow autoencoders for sparse data. WWW, 2019.\n\n[7] Xiang Wang et al. Neural graph collaborative filtering. SIGIR, 2019.\n\n[8] Harald Steck. Autoencoders that don’t overfit towards the identity. NIPS, 2020.\n\n[9] Vojtech Vancura et al. Scalable linear shallow autoencoder for collaborative filtering. RecSys, 2022.\n\n[10] Peter Bartlett et al. Convexity, classification, and risk bounds. Journal of the American Statistical Association, 2006."}, {"review_id": "1ww7TNzCor", "reviewer": "Reviewer_7BK6", "summary": "This work develops PAC-Bayes generalization bound for multivariate linear regression and extends that bound to the Linear AutoEncoders (LAEs) as a constrained form of multivariate linear regression under the bounded data assumption. To make the theory practical, an algorithm for computing the bound is developed. The bound is further applied to Linear Autoencoder (LAE) recommender systems using the EASE model for empirical validation. Experiments conducted on standard recommendation datasets confirm the bound's effectiveness in a practical setting.", "questions": "Related to weakness 2: Why do we see correlation in the used loss  and the ranking metrics, given that optimal loss may not be the best outcome for ranking?", "limitations": "yes", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "**Strengths**\n\n- Bounds for a more general case: The paper provides a nice theoretical extension of PAC-Bayes bounds to multivariate linear regression case, while also formalizing and clarifying the the established convergence proof conditions from Shalaeva et al.\n- Extension to Linear AutoEncoder: The paper does a good job of linking its theory to Linear Autoencoders (LAEs) in recommender systems. The derived multivariate bound is adjusted to account for the specific case of LAEs, like the zero-diagonal constraint on the weight matrix and how the input/target data is structured in implicit feedback scenarios (Section 4.2).\n- The paper also provides an algorithm to compute the bound in practice. A neat closed-form solution for the optimal posterior is also provided under the Gaussian data assumption (Theorem 5.2). The theoretical results are backed with real-world tests in Table 1, which highlight the tightness of the bound.\n\n**Weaknesses**\n\n- Limited applicability to low-rank scenario: The practical computation methods, particularly the closed-form solutions in Section 5, are derived under the assumption that the weight matrix W follows a multivariate Gaussian distribution, which implies it is full-rank. This framework does not directly apply to low-rank structure: W = UV^T. Such low-rank models are common in recommendation systems, e.g. in collaborative filtering and matrix completion. Adapting the theory to impose distributions on the factors U and V would break the Gaussian assumption for W, rendering the current computational methods inapplicable. This limits the scope of the paper's impact.\n- Mismatch Between Theoretical Loss and Recall/NDCG metrics: The commonly used metrics in recommendation systems (Recall/NDCG) are set-based ranking metrics (e.g. computed as Recall@K and NDCG@k for a set of k predictions). It is not clear to me why the chosen loss ($(||Y - WX||_F^2)$) in the analysis is correlated to these metrics. A model that is optimal under squared error is not guaranteed to be optimal for ranking.", "was_revised": true, "final_justification": "The authors have provided clarifications related to squared loss vs. ranking metrics (Weakness 2 in the review). I will maintain my current rating.", "rebuttal": "Thank you for your insightful feedback and comments. We hope the rebuttal below fully addresses your concerns, and we will incorporate this important discussion into our paper.\n\n> Weakness 2: Mismatch Between Theoretical Loss and Recall/NDCG metrics: The commonly used metrics in recommendation systems (Recall/NDCG) are set-based ranking metrics (e.g. computed as Recall@K and NDCG@k for a set of k predictions). It is not clear to me why the chosen loss ($\\|\\|Y - WX\\|\\|\\_F^2$) in the analysis is correlated to these metrics. A model that is optimal under squared error is not guaranteed to be optimal for ranking.\n\n> Related to weakness 2: Why do we see correlation in the used loss and the ranking metrics, given that optimal loss may not be the best outcome for ranking?\n\nThank you for raising this important point regarding the potential mismatch between our theoretical loss function and the evaluation metrics used in recommender systems. We hope the explanation below clarifies the rationale and empirical justification for our use of the squared loss $\\|\\|Y - WX\\|\\|\\_F^2$, and its observed correlation with ranking-based metrics such as Recall@K and NDCG@K.\n\nFirst, it is well understood that ranking metrics, particularly Recall@K and NDCG@K in implicit feedback settings, are set-based, discrete, and non-differentiable. As such, they are difficult to optimize directly using gradient-based methods. Consequently, it has become standard practice in collaborative filtering and recommender system research to employ regression-style or reconstruction-based surrogate losses, most notably the squared error $\\|\\|Y - WX\\|\\|_F^2$, as training objectives. This approach underlies many successful algorithms, including SLIM [1], EASE [2], EDLAE [3], and matrix factorization models, which consistently demonstrate strong empirical performance on ranking tasks despite optimizing a non-ranking loss.\n\nSecond, while minimizing squared loss does not theoretically guarantee optimal ranking, numerous empirical studies (including our own) observe a strong correlation between improvements in squared error and gains in Recall/NDCG. This is especially evident in *latent linear models*, where $WX$ serves as a real-valued scoring function over items. Since ranking metrics (e.g., NDCG) depend only on the relative order of scores, and squared loss penalizes deviations in score magnitudes, optimizing squared error often leads to improved ranking performance, even if this relationship is not theoretically guaranteed.\n\nThird, we believe this correlation is further reinforced in our setting by the *LAE-specific constraints* imposed on both the data and the model, as described in Section 4.2. These constraints, such as the train/test split structure and zero-diagonal regularization, are applied consistently in our bounds and in the evaluation using ranking metrics. This constraint alignment encodes *shared information* about what constitutes a good model, promoting consistency between the evaluation signals. That is, although the loss and ranking metrics target different objectives, the shared modeling constraints lead both to emphasize similar aspects of model quality.\n\nThat said, these metrics do not share all information. For example, squared loss evaluates absolute prediction accuracy, while Recall and NDCG focus on the relative ranking of positive items. As a result, the correlation is not perfect. This explains the observed mismatch: the model that minimizes squared loss may not be the same as the one that maximizes Recall@K or NDCG@K, and vice versa.\n\nThis effect is visible in our experiments. For example, in Table 1 on the MSD dataset, the model achieving the best Recall@50 and NDCG@100 occurs at $\\gamma = 500$, while the lowest LH and RH values are attained at $\\gamma = 1000$. This highlights that the loss-optimal and ranking-optimal models may differ. Nevertheless, we observe a general trend: lower LH/RH tends to correspond to better ranking performance, even if the relationship is not strictly monotonic.\n\nFinally, we emphasize that our theoretical focus is on deriving a PAC-Bayesian generalization bound for the squared error—an analytically tractable and widely used surrogate loss. While deriving bounds directly for Recall or NDCG would be desirable, such metrics are discontinuous, non-convex, and intractable for theoretical analysis. Our work thus offers principled insight into a practical surrogate loss function, bridging the gap between theoretical guarantees and empirical evaluation.\n\nIn summary, although squared loss and ranking metrics optimize different objectives, their empirical alignment is well-documented and widely accepted in the literature. Our analysis provides theoretical justification for a surrogate loss function that is tractable, interpretable, and highly relevant to real-world recommendation tasks.\n\n**References:**\n\n[1] Xia Ning and George Karypis. Slim: Sparse linear methods for top-N recommender systems. ICDM, 2011.\n\n[2] Harald Steck. Embarrassingly shallow autoencoders for sparse data. WWW, 2019.\n\n[3] Harald Steck. Autoencoders that don’t overfit towards the identity. NIPS, 2020."}, {"review_id": "SIlCXoqRea", "reviewer": "Reviewer_dPHG", "summary": "This paper develops a PAC-Bayes generalization bound tailored for multivariate linear regression. The authors then derive a sufficient conditions to guarantee its theoretical convergence under certain assumption (stated in Thm 3.3), e.g. Gaussian data assumption. They then adapt this framework to handle bounded data distributions and then specifically apply it to analyze Linear Autoencoders (LAEs). After that, under Gaussian data assumption, the authors presents a practical computational method for this LAE-specific bound. At the end, experimental results on validating the tightness of the proved bounds are provided.", "questions": "My question will be mainly about the assumptions on input data. I think it will help greatly if the paper can provide more detailed justifications on those assumptions, e.g., maybe they are not common, but they are important and indeed appearing in XX by evidence/literature XX.\n\nMy next question might be helpful (hopefully) for enriching the paper or future direction, if possible. Can the value of the PAC-Bayes bound after a bunch of update steps during the training process help adaptively improve the training process? Is it possible?", "limitations": "yes", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strengths:\n(1) The paper provides the first PAC-Bayes bound for LAEs.\n(2) The paper derives a sufficient condition for convergence in the case of Gaussian data.\n(3) The paper also fixes the issues in Shalaeva's bound and proves the claimed bound in Shalaeva's work.\n\nWeaknesses:\n(1) The calculation of PAC-Bayes bound requires the data distribution to be known.\n(2) A couple of the results requires the Gaussian data assumption.", "was_revised": true, "final_justification": "I am satisfied with the authors' responses to our questions. Since my initial rating was already positive and I still find the data assumption to be not general enough, I will maintain my score.", "rebuttal": "Thank you for your insightful and constructive comments. We hope our rebuttal addresses your concerns.\n\n> My question will be mainly about the assumptions on input data. I think it will help greatly if the paper can provide more detailed justifications on those assumptions, e.g., maybe they are not common, but they are important and indeed appearing in XX by evidence/literature XX.\n\nThanks for highlighting this important question. The assumptions made on the input and target matrices $X$ and $Y$, as described in Lines 108–120 and 201–208, follow a standard evaluation protocol in collaborative filtering known as \\textit{weak generalization} [1, 2]. Specifically, the user–item interaction matrix $R$ is binarized and then randomly split into disjoint subsets to create the input matrix $X$ and the held-out target matrix $Y$, where only the 1s (positive interactions) are split.  \n\nThis setup is widely adopted in recommender system literature to evaluate generalization performance on unseen interactions and to avoid data leakage between training and testing. We follow this protocol as used in prior work such as [1, 2, 3]. Thank you for raising this point, we will include these references and a clearer explanation of the weak generalization setup in the revised manuscript.\n\n> My next question might be helpful (hopefully) for enriching the paper or future direction, if possible. Can the value of the PAC-Bayes bound after a bunch of update steps during the training process help adaptively improve the training process? Is it possible?\n\nThank you very much for this insightful and forward-looking question. We sincerely appreciate your suggestion, as it touches on a promising direction that goes beyond the current scope of our work and could meaningfully enrich future research.\n\nTo clarify, in the current version of our paper, the PAC-Bayesian bound is used solely for *evaluation* purposes. As described in Lines 201–202, the LAE bound measures the generalization gap between the true and empirical risks on a held-out test set. The model $W$ used in this bound can originate from any training method and is not required to minimize the bound itself.\n\nThat said, we observe in Table 1 that the bound correlates strongly with practical performance metrics such as Recall@K and NDCG@K. This empirical observation suggests that the bound could potentially serve as a surrogate objective during training. Minimizing the bound might lead to models with stronger ranking performance -- thus bridging the gap between theory and practical evaluation.\n\nMoreover, as established in Theorem 5.2, our bound admits a closed-form expression, which enables efficient computation and avoids the need for complex iterative updates. This adds to its appeal as a candidate training objective.\n\nAgain, we are grateful for your suggestion, and we agree that incorporating the PAC-Bayesian bound into the training loop presents an exciting future direction. We will highlight this possibility more explicitly in our revision and hope to explore it further in future work.\n\n**References:**\n\n[1] Jaewan Moon et al. It's enough: Relaxing diagonal constraints in linear autoencoders for recommendation. SIGIR, 2023.\n\n[2] Harald Steck. Embarrassingly shallow autoencoders for sparse data. WWW, 2019.\n\n[3] Xiang Wang et al. Neural graph collaborative filtering. SIGIR, 2019."}, {"review_id": "ZCJ2nDW3ud", "reviewer": "Reviewer_z6h4", "summary": "The paper extends PAC-Bayes theory to multivariate linear regression, deriving a novel generalization bound that remains tight for full-rank models under unbounded loss. It adapts this bound to Linear Autoencoders by framing them as constrained multivariate regression on bounded data, addressing challenges like zero-diagonal weight constraints and input-target dependencies. Experiments demonstrate tightness and correlation with recommendation metrics like Recall@K and NDCG@K.\n\nStrengths:\n\nFirst PAC-Bayes bound for multivariate regression and LAEs, filling a theoretical gap.\nRigorous convergence analysis with sufficient conditions.\nComputationally efficient methods for bound optimization.\nStrong empirical validation on real-world datasets.\n\nLimitations:\nThe author has elaborately described his framework and has carried out numerous theoretical studies and experiments to verify its effectiveness.\nDue to the large number of formulas and symbols, some parts may appear rather messy (for example, the lower half of the second page). It is hoped that the author can adjust the format and structure of the paper to make it more aesthetically pleasing and easy to read.", "questions": "The author has elaborately described his framework and has carried out numerous theoretical studies and experiments to verify its effectiveness.\nDue to the large number of formulas and symbols, some parts may appear rather messy (for example, the lower half of the second page). It is hoped that the author can adjust the format and structure of the paper to make it more aesthetically pleasing and easy to read.", "limitations": "The author has elaborately described his framework and has carried out numerous theoretical studies and experiments to verify its effectiveness.\nDue to the large number of formulas and symbols, some parts may appear rather messy (for example, the lower half of the second page). It is hoped that the author can adjust the format and structure of the paper to make it more aesthetically pleasing and easy to read.", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "This is the first PAC-Bayes bound for multivariate regression and LAEs, and fill a theoretical gap.\nThe rigorous convergence analysis with sufficient conditions shows the effectiveness of the framework.\nThe author has proved it is a computationally efficient methods for bound optimization.\nThe framework has a strong empirical validation on real-world datasets.", "was_revised": true, "final_justification": "We are satisfied that the author's response has resolved our issues. Since our initial assessment was already positive, we have decided to maintain this score.", "rebuttal": "Thank you very much for your insightful reviews and constructive suggestions.\n\nWe will incorporate your suggestions in the revised version of the paper, including adjustments to the format and structure, particularly in the lower half of the second page as noted. We believe these improvements will enhance the clarity and presentation of our work."}, {"review_id": "qANySKUdya", "reviewer": "Reviewer_fRvV", "summary": "This paper provides a novel theoretical understanding of Linear AutoEncoders (LAEs) through PAC-Bayes theory. First, the authors extend an existing result of Shalaeva et al. to the case of multivariate linear regression and explicit sufficient conditions to ensure convergence of such bounds. Then, they adapt their results to LAEs by focusing on the particular case of multivariate linear regression on bounded data. From such theoretical results, they derive a novel learning algorithm for LAEs and conclude their work with some experiments.", "questions": "- L. 208 what is $X_{*,j}$ in practice? The star suggest this is a theoretical quantity in some sense, is this true? \n- Lemma 4.2 why highlighting the specific form of $R_{true}$ given this quantity cannot be computable in practice?\n- Theorem 3.2: can you compare your result PAC-Bayes bounds existing for linear regression here (instead of the appendix)? For instance, with either [1,2] or Shalaeva's work to understand precisely what is original in your bound?\n- Theorem 3.2 is obtained from the general Alquier's bound, which is often achieved by subgaussian or subexponential assumptions, which is a choice you do not make here. In recent works PAC-Bayes bounds for heavy-tailed losses (i.e. admitting finite variances) have been obtained via supermartingales argument see [3,4]). Couldn't you reach tighter results using these bounds? This would avoid the exponential moment which may ba hard to control and (potentially) computationally instable\n- l. 87 the assumption on $\\mathcal{D}$ is not valid throughout the whole paper, you should detail it\n- l.196, isn't it the exact same statement than Theorem 3.2? What is new here?\n- l. 214, 'And the PAC-Bayes bound for LAEs is formed by applying the above two constraints to (8)' what does this mean?  \n- l. 252 You affirm that 'neural network models typically do not admit closed-form solutions for the optimal ρ for (11).' This is actually not true if you choose to assimilate a neural net to its vectors of weights. (11) is minimized by Gibbs posteriors (see e.g. [5, Section 5.1])\n- l. 266-267 Are such derivation useful in the main text?\n- It seems that \n- Proposition 5.3 Again, why highlighting the specific form of the exponential moment (a function of $R_{true}$) given this quantity cannot be computable in practice?\n- l. 291 you say this exponential moment is much easier to compute. How can you make this bound tractable as it is supposed to be unknown on real-life experiments?\n\n\nIn conclusion, while some theoretical novelties seems to appear in this work (PAC-Bayes bounds for the specific case of multilinear regression without light-tailed assumptions), additional discussions and comparisons with modern methods in PAC-Bayes is severely lacking and make the paper hard to parse.  As I am not an expert in LAE, I may miss the particular relevance of the proposed results but, in my humble opinion, the paper has not been truly pedagogical about it. I believe then it needs to be enriched and rewritten before acceptance. \n\n\nReferences\n[1] Germain et al., PAC-Bayesian Theory meets Bayesian inference, 2016\n\n[2] Haddouche et al. PAC-Bayes Unleashed: Generalisation Bounds with Unbounded Losses, 2021\n\n[3] Haddouche et al. PAC-Bayes Generalisation Bounds for Heavy-Tailed Losses through Supermartingales; 2023\n\n[4] Chugg et al. A unified recipe for deriving (time-uniform) PAC-Bayes bounds, 2024\n\n[5] Catoni, A PAC-Bayesian approach to adaptive classification, 2003", "limitations": "Literature review is missing to understand the precise novelty of such results. Also, a strength of Theorem 3.2, which is to avoid assumptions such as subgaussianity, is not highlighted.", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 1, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Strengths: \n- The question of PAC-Bayes for LAEs is, to my knowledge, new.\n- Authors avoid classical assumptions in PAC-Bayes such as subgaussianity for their specific setting\n\nWeaknesses\n- The paper is hard to read and there is little discussion about the results.\n- Lit review is incomplete in the main text?\n- The motivation behind the various lemma and theorem is not always clear (e.g Lemmas 4.2,4.3 Proposition 5.3)", "was_revised": true, "final_justification": "Authors answered to all my concerns about the writing of the paper. Furthermore, they seem to have answered in a convincing manner (to me) to the questions of other reviewers, while plugging those change implies a significant rewriting, it seems it can be done in a reasonable amount of time.", "rebuttal": "Thank you very much for your detailed review and thoughtful feedback regarding the clarity of our paper. We hope the rebuttal below addresses all of your concerns. Due to character limits, we may not be able to address some questions in full detail. Please do not hesitate to reach out if any points remain unclear or if further concerns arise; we sincerely welcome the opportunity to clarify and improve our work.\n\n**Questions**:\n\n> Q1: L. 208 what is $X\\_{*j}$ in practice? The star suggest this is a theoretical quantity in some sense, is this true?\n\nThanks for pointing out this symbol. $X\\_{*j}$ denotes the $j$th column of $X$, which indeed has practical meaning, as discussed in Lines 108-120. $X\\_{\\*j}$ represents the items user $j$ have observed, while $Y\\_{*j}$ represents the items treated as unseen for user $j$. An LAE model $W$ that minimizes $\\|\\|Y\\_{\\*j} - WX\\_{\\*j}\\|\\|\\_F^2$ aims to learn latent relationships between the observed and unseen items. This modeling of $X\\_{\\*j}$ and $Y\\_{\\*j}$ follows a widely adopted in collaborative filtering evaluation approach called *weak generalization* [1, 2, 3]. Our Constraint 2 in Section 4.2 is formulated based on this weak generalization setting.\n\n> Q2: Lemma 4.2 why highlighting the specific form of $R^{\\text{true}}$ given this quantity cannot be computable in practice?\n\nThanks for your question. Our original claim is that (7) is computable in practice, but only when the data distribution $\\mathcal{D}$ is *non-oracle*. Since $\\Sigma_{xx}, \\Sigma_{xy}$ and $\\Sigma_{yy}$ depend on $\\mathcal{D}$, a non-oracle $\\mathcal{D}$ ensures these quantities are known, thereby making (7) computable. \n\nHowever, a non-oracle $\\mathcal{D}$ is quite restrictive for modeling real-world data -- it can only represent datasets that are fixed and finite, with all training and testing samples drawn from it, which are difficult (though not impossible) to find in practice. Many real-world datasets are continuously collected or grow over time and are often modeled as samples drawn from an unknown oracle distribution $\\mathcal{D}$. Our contributions show that the bound is at least computable in restricted scenarios where datasets can be modeled with non-oracle distributions.\n\n> Q3: Theorem 3.2: can you compare your result PAC-Bayes bounds existing for linear regression here (instead of the appendix)? For instance, with either [1,2] or Shalaeva's work to understand precisely what is original in your bound?\n\nThank you for this helpful suggestion. We will revise the main text to include a more explicit comparison with Shalaeva et al's work.\n\nTo clarify the relationship: Shalaeva's bound (1) applies to *single-output* linear regression (i.e., $p = 1$ in our notation), while our bound (4) extends this result to the more general *multi-output* case where $p \\geq 1$. As noted in Line 148, our result strictly generalizes theirs: under appropriate parameter choices, our bound reduces exactly to Shalaeva’s.\n\nMore importantly, our formulation captures scenarios where multiple output dimensions may be dependent, cases not addressed by prior single-output analyses. This extension is nontrivial, as it requires a redefinition of key variance terms and a careful treatment of the loss geometry in higher-dimensional output spaces.\n\n> Q4: Theorem 3.2 is obtained from the general Alquier's bound, which is often achieved by subgaussian or subexponential assumptions, which is a choice you do not make here. In recent works PAC-Bayes bounds for heavy-tailed losses (i.e. admitting finite variances) have been obtained via supermartingales argument see [3,4]). Couldn't you reach tighter results using these bounds?\n\n> Limitations: A strength of Theorem 3.2, which is to avoid assumptions such as subgaussianity, is not highlighted.\n\nThank you for this insightful question and for pointing us to these important recent developments.\n\nIn Alquier’s general framework, sub-Gaussian or sub-exponential assumptions are commonly used to ensure that the exponential moment $\\Psi\\_{\\pi, \\mathcal{D}}(\\lambda, m) = \\ln \\mathbb{E}\\_{\\pi} \\mathbb{E}\\_{\\mathcal{D}} \\left[ e^{\\lambda (R^{\\text{true}}(W) - R^{\\text{emp}}(W))} \\right]$ is bounded, but does not guarantee the convergence, that $\\Psi\\_{\\pi, \\mathcal{D}}(\\lambda, m) \\to 0$ as $m \\to \\infty$. For instance, Germain et al.[4, A.4] showed that in linear regression, $\\Psi$ is sub-gamma, but Shalaeva et al. [5] later demonstrated that this bound does not converge and proposed an alternative upper bound that is both strictly tighter and convergent.\n\nOur bound builds on the work of Shalaeva et al, using a more specific assumptions on $(\\lambda, \\pi)$ rather than the general sub-gamma assumption on $\\Psi$. As discussed in Section 3.2, certain combinations of $(\\lambda, \\pi)$ yield a both sub-gamma and convergent upper bound for $\\Psi$, indicating that our assumptions on $(\\lambda, \\pi)$ are stronger than the standard sub-gamma assumption and lead to a tighter bound.\n\nRegarding the supermartingale-based bounds, we note that they also ensure convergence at a $1/\\sqrt{m})$ rate (see [6], Page 8), which is the same as our bound (See Appendix B in our paper).\n\n> Q5: l. 87 the assumption on $\\mathcal{D}$ is not valid throughout the whole paper, you should detail it.\n\nThanks for pointing out this issue. We apologize for the confusion caused by the abuse of notation $\\mathcal{D}$. $\\mathcal{D}$ is redefined in Assumption 3.1 and again in Assumption 4.1, and both definitions differ from the one introduced in Line 87. We will make this abuse of notation explicit in the revision.\n\n> Q6: l.196, isn't it the exact same statement than Theorem 3.2? What is new here?\n\nPlease see Line 195. The definitions of $R^{\\text{emp}}$ and $R^{\\text{true}}$ differ here: they are derived from Assumption 4.1 (bounded data), rather than Assumption 3.1 (Gaussian data). Apart from this difference, everything else in (8) remains the same as in (4).\n\n> Q7: l. 214, 'And the PAC-Bayes bound for LAEs is formed by applying the above two constraints to (8)' what does this mean?\n\nWe apply the two constraints to the loss $\\|\\|y - Wx\\|\\|_F^2$ in (8). Constraint 1 applies to $W$, and Constraint 2 applies to $(x, y)$. Together, these constraints adapt the original loss to the LAE setting.\n\n> Q8: l. 252 You affirm that 'neural network models typically do not admit closed-form solutions for the optimal $\\rho$ for (11).' This is actually not true if you choose to assimilate a neural net to its vectors of weights. (11) is minimized by Gibbs posteriors (see e.g. [5, Section 5.1])\n\nThank you for your thoughtful comment and for pointing this out. While it is true that the bound in (11) is minimized by a Gibbs posterior, the Gibbs posterior is typically a complex distribution without a closed-form density function or parameterization. Therefore, although the Gibbs posterior is well-defined in theory, it is often not analytically tractable in practice.\n\nIn contrast, when choosing a Gaussian posterior (as we do in our work), we can derive an explicit closed-form solution for the optimal parameters (mean and variance) under linear model assumptions, such as in the case of linear autoencoders. While this choice may lead to a looser bound compared to the Gibbs posterior, it offers significant computational advantages, especially for large-scale datasets, by avoiding iterative approximation techniques. This reflects a trade-off between bound tightness and computational efficiency.\n\n> Q9: l. 266-267 Are such derivation useful in the main text?\n\nThe derivation in Lines 266 - 267 defines the constants $B$ and $C$, which are then used in forming (13) and (14).\n\n> Q10: Proposition 5.3 Again, why highlighting the specific form of the exponential moment (a function of $R^{\\text{true}}$) given this quantity cannot be computable in practice?\n\nPlease see our answer to Q2. The quantity in (13) of Proposition 5.3 is computable when $\\mathcal{D}$ is non-oracle. Non-oracle data distributions have been used in earlier generalization bounds, for example, the bound by Srebro et al., which we cite in Line 309.\n\n> Q11: l. 291 you say this exponential moment is much easier to compute. How can you make this bound tractable as it is supposed to be unknown on real-life experiments?\n\nPlease see our answers to Q2 and Q10, in real-life experiments, if the dataset remains fixed and finite over time, and all training and testing samples are drawn from it, then it can be modeled by its population distribution, which is non-oracle and makes the bound computable.\n\n**Weaknesses:**\n\n> Lit review is incomplete in the main text?\n\nWe will consider moving Literature Review from Appendix to the main text if space permits in the revision.\n\n> The motivation behind the various lemma and theorem is not always clear (e.g Lemmas 4.2,4.3 Proposition 5.3)\n\nThank you for raising this point.\n\nLemmas 4.2 and 4.3 provide the explicit expansion of (8) under Assumption 4.1. These expansions are essential for implementing Algorithm 1, which computes the objective in (8) efficiently. Thus, the lemmas directly support the algorithmic component of our work.\n\nProposition 5.3 gives a closed-form expression for the quantity $\\mathbb{E}\\_{\\pi}\\left[e^{\\lambda R^{\\textnormal{true}}(W)}\\right]$, which is used in the computational complexity analysis in Lines 276--282. It also serves as a crucial intermediate step in deriving the low-complexity upper bound in Theorem 5.4.\n\n**References:**\n\n[1] Jaewan Moon et al. It's enough: Relaxing diagonal constraints in linear autoencoders for recommendation. 2023.\n\n[2] Harald Steck. Embarrassingly shallow autoencoders for sparse data. 2019.\n\n[3] Xiang Wang et al. Neural graph collaborative filtering. 2019.\n\n[4] Germain et al. PAC-Bayesian Theory meets Bayesian inference. 2016.\n\n[5] Shalaevaet al. Improved PAC-Bayesian bounds for linear regression. 2020.\n\n[6] Haddouche et al. PAC-Bayes Generalisation Bounds for Heavy-Tailed Losses through Supermartingales. 2023."}]

Author rebuttal:

Meta-review: ### Summary of scientific claims and findings

The paper develops the first PAC-Bayes generalization bound for multivariate linear regression targeting the generalization gap, with sufficient convergence conditions under Gaussian data assumptions. It extends this framework to Linear Autoencoders (LAEs), introducing LAE-specific constraints, and proposes a computationally efficient evaluation method. Empirical validation on large-scale recommender datasets shows the bound is tight and correlates with practical metrics such as Recall and NDCG.

### Strengths

* **Novelty**: First PAC-Bayes bound for both multivariate regression and LAEs, addressing a theoretical gap.
* **Technical Soundness**: Rigorous and technically sound analysis, including convergence conditions.
* **Practical Relevance**: Efficient computational method.
* **Empirics**: Empirical results show tightness and correlation with performance.

### Weaknesses and missing elements

* Practical computation restricted by Gaussian/full-rank assumptions, limiting applicability.
* Discrepancy between theoretical squared-loss bound and ranking-based evaluation metrics.
* Some presentation and clarity issues (notation reuse, dense exposition, formatting).
* Initial literature review was incomplete (e.g., earlier PAC-Bayes results for regression and VAEs).

### Key reasons for acceptance

1. Novel and technically rigorous contribution filling an acknowledged gap in PAC-Bayes theory.
2. Relevance and applicability to recommender systems, with practical evaluation.
3. Rebuttal and discussion stage convincingly addressed literature positioning, clarity, and empirical concerns.

### Summary of discussion and rebuttal

Before rebuttal:

* `npXJ`: Criticized lack of references to prior PAC-Bayes work on multivariate regression and VAEs; recommended *Reject*.
* `7BK6`: Positive on theory and empirical validation but raised concerns about Gaussian assumptions and mismatch with ranking metrics.
* `dPHG`: Positive but questioned assumptions on input data and potential use of the bound in training.
* `fRvV`: Raised clarity and novelty positioning concerns; found paper difficult to parse.
* `z6h4`: Strongly positive; main concern was formatting/structure.

### Authors’ rebuttal

* Clarified distinctions between their bound and Mai (2022, 2023) on excess risk, and between LAEs and VAEs.
* Acknowledged limitations of Gaussian/full-rank assumptions, suggesting possible low-rank extensions.
* Justified squared-loss use as a surrogate, citing literature and empirical correlation with ranking metrics.
* Clarified notation, assumptions, and motivation for key lemmas and theorems.
* Agreed to improve literature review and presentation.

### Post-rebuttal positions

* `npXJ`: Upgraded to positive, calling the contribution “reasonably good” and suitable for NeurIPS.
* `7BK6`: Maintained borderline accept, satisfied with clarification but noting scope limits.
* `dPHG`: Maintained positive, satisfied with responses.
* `fRvV`: Updated to positive, stating changes requested would secure acceptance.
* `z6h4`: Maintained positive, satisfied with rebuttal.

### Weighing these points

All reviewers converged toward acceptance after rebuttal, with initial concerns on literature positioning, clarity, and evaluation satisfactorily addressed. While some limitations remain, these are well-documented, acknowledged by authors, and do not undermine the core contributions.
Decision: Accept (poster)