On PAC-Bayes Bounds for Linear Autoencoders

As authors didn't bother to actually answer my questions, in particular about the actual usage for RS systems, I lower my score.

Question 3: I wonder if your bound could be helpful for optimistic exploration strategies when the user (row) is sampled from a uniform distribution and the item is selected by the algorithm. This might be challenging since independence is broken, and you'd likely need additional assumptions (such as rank-1 matrix only) - which, while commonly used in state-of-the-art approaches, may be overly simplistic.

Answer: Thanks for sharing this idea with us! It will be interesting to investigate how to incorporate feedback/responses and item (response) dependency into our bound, which can be non-trivial. We found that one latest paper [Tasdighi24] seems to have an interesting solution on this, where the authors introduce a reinforcement learning algorithm that leverages PAC-Bayes analysis to model epistemic uncertainty in the critic component. This approach enables the algorithm to dynamically adapt its exploration strategy, thereby enhancing performance in continuous control tasks.

Question 4: In all EASE-based system implementations I'm aware of, the zero-diagonal constraint is used. While the paper studies this to some extent, can your analysis provide new insights beyond confirming it as a useful bias?

Answer: Our bound applies to $W$ with or without zero-diagonal constraint. One effect of the zero-diagonal constraint is that it can reduce the PAC-Bayes bound, thereby improving the model's generalizability.

In details, we have shown in Section 4.3 that if the constraint of $W$ is changed from $\text{diag}(W) \sim N(0, \sigma^2I)$ to zero diagonal constraint $\text{diag}(W) = 0$, then the covariance matrix of each column $\Sigma_r^{1/2}(W^* - W)_{*i}$ will be reduced from $\sigma^2\Sigma\_{r}$ to $\sigma^2(\Sigma_r - (\Sigma\_r^{1/2})\_{*i}(\Sigma_r^{1/2})\_{*i}^T)$, i.e., $\sigma^2\Sigma\_{r} \succeq \sigma^2(\Sigma\_r - (\Sigma_r^{1/2})\_{*i}(\Sigma\_r^{1/2})\_{*i}^T)$. The reduction of the covariance matrix of the distribution of $W$ leads to a smaller PAC-Bayes bound, as shown in the proof of Theorem 4 in Appendix A. This suggests that the model has better generalizability.

Question 5: You mention producing 'non-vacuous' bounds. Could you elaborate on what constitutes such a bound? Specifically, if your definition implies practical usefulness, could you explain how these bounds could be used to improve recommender systems?

A PAC-Bayes bound is considered non-vacuous if the gap between the bound and the test error is small. To the best of our knowledge, there is no universal criterion for how small the gap must be to classify a bound as non-vacuous. [Dziugaite17] showed in their experiments that PAC-Bayes bounds within $10$ times the test error can be considered non-vacuous, which is the criterion we adopted.

If a theoretical bound is too vacuous, for example, $1000$ times the test error, then it will predict almost nothing in practice. A discussion of opinions on vacuous bounds can be found in Section 3.1 of [Alquier21].

A non-vacuous bound benefits a recommender system by explaining experimental results. For example, if experiments show that a model has a large test error, the bound can help determine whether this is due to overfitting or a biased test set. However, it is important to note that the bound does not guide model training; it simply justify whether the model, once trained by whatever training method, should perform well on unseen data or not. We have also added a discussion on the empirical implication and potential application of our bounds in Section 6 (Conclusion and Discussions).

References:

[Alquier21] Pierre Alquier. User-friendly introduction to pac-bayes bounds. arXiv preprint arXiv:2110.11216, 2021.

[Tasdighi24] Probabilistic Actor-Critic: Learning to Explore with PAC-Bayes Uncertainty, Bahareh Tasdighi, Nicklas Werge, Yi-Shan Wu, Melih Kandemir, arXiv preprint arXiv:2402.03055, 2024.

We sincerely apologize for the delayed response, as we have tried to submit a few rounds of revisions to address the suggestions and concerns from you and other reviewers before the revision submission deadline.

Question 1 Since the initial error appears in a paper published at AAAI, why not contact the Program Committee from that edition? I believe AAAI should take responsibility for errors made during their review process.

Answer: Thanks for your suggestion. As Reviewer nWHo has pointed out, simply claiming it as being an error is a bit unfair to the original paper as their analysis is a bit informal. The overall result/bound is still correct, but the convergence analysis is missing some necessary/sufficient conditions for the distribution $\theta$. In the revision, we have further clarified this lack of "additional" conditions for the convergence analysis. We have made this point clear in the introduction (Blue text in revision).

Note that one of our main contributions is to provide the first PAC-Bayes bound for the multivariate linear regressions which in special cases reduces to the AAAI paper bound; we also present a sufficient condition by specifying the distributions of $\theta$ (corresponding to $W$ in our bound) under which the bound converges or diverges. We have presented this result and generalization in new Section 3.

Question 2: I'm curious about how similar or different this is to the analysis of random projections and LinUCB. While we're working in a Bayesian framework, the mathematical tools used are quite similar, and some research has already bridged the gap between these approaches (such as the work in https://www.jmlr.org/papers/volume17/14-087/14-087.pdf)

Answer: Thanks for providing the reference for the very interesting paper by Russo & Van Roy on information theoretical bound of Thompson sampling. Based on our understanding, both the PAC-Bayes bound and Russo & Van Roy's bound are derived using tools from information theory.

The PAC-Bayes bound, use Alquier's bound [Alquier16] for example, is

$

P\left(\mathbb{E}\_{\theta \sim \rho}[R^{\text{true}}(\theta) - R^{\text{emp}}(\theta)] <  \frac{1}{\lambda}\left[D(\rho||\pi) + \ln \frac{1}{\delta} + \Psi\_{\pi, \mathcal{D}, l}(\lambda, m)\right]\right) \ge 1 - \delta

$

where $R^{\text{emp}}(\theta)$ is the empirical risk, $R^{\text{true}}(\theta)$ is the true risk, $\pi$ is the prior distribution of $\theta$, $\rho$ is the posterior distribution of $\theta$.

The regret bound of Thompson sampling (Proposition 4, [Russo16]) is

$

\mathbb{E}[R(Y_{A^*}) - R(Y_A)] \le \sqrt{\frac{1}{2}I(A^*, (A; Y))}

$

where $A^\*$ is a random variable that maximizes the reward $R$, while $A$ is any other random variable.

We notice some differences between the two bounds:

1. In the PAC-bayes bound, $R^{\text{true}}(\theta)$ is the expectation of the random variable $R^{\text{emp}}(\theta)$. Thus the distance between $R^{\text{true}}(\theta)$ and $R^{\text{emp}}(\theta)$ can be measured by concentration inequalities, on which the PAC-bayes bound is based. In the regret bound, $R(Y\_{A^\*})$ is the upper bound of the random variable $R(Y_A)$. The regret bound is deterministic (it is fixed and does not rely on probability), while the PAC-Bayes bound is probabilistic.

2. The regret bound uses the Pinsker's inequality as its information bound (Step (b) in the proof of Proposition 4, relating to Fact 9 [Russo16]).

$

\mathbb{E}_{x\sim P}[g(x)] - \mathbb{E}_{x\sim Q}[g(x)] \le \sqrt{\frac{1}{2}D(P||Q)}

$

while the PAC-Bayes bound uses Donsker and Varadhan’s variation (Lemma 2.2, [Alquier21]):

$

\mathbb{E}_{x\sim P}[g(x)] - \log\mathbb{E}_{x\sim Q}[e^{g(x)}] \le D(P||Q)

$

Could you please provide more details regarding the analysis of random projection or LinUCB? We have reviewed the original papers on random projection [Bingham01] and LinUCB [Li10], but we did not find much (mathematical) analysis in them, which has left us uncertain. Please let us know if we have referenced the wrong papers. Thank you for your understanding.
 

References:

[Alquier16] Pierre Alquier, James Ridgway, and Nicolas Chopin. On the properties of variational approximations of gibbs posteriors. Journal of Machine Learning Research, 17(236):1–41, 2016

[Alquier21] Pierre Alquier. User-friendly introduction to pac-bayes bounds. arXiv preprint arXiv:2110.11216, 2021.

[Li10] Lihong Li, et al. A contextual-bandit approach to personalized news article recommendation. Proceedings of the 19th international conference on World wide web. 2010.

[Russo16] Daniel Russo, and Benjamin Van Roy. An information-theoretic analysis of thompson sampling. Journal of Machine Learning Research, 2016.

[Bingham01] Ella Bingham, and Heikki Mannila. Random projection in dimensionality reduction: applications to image and text data. SIGKDD. 2001.

Weakness 1: The paper's positioning relative to previous work in recommender systems is inadequate and imprecise. For instance, citing Rendle 2022 for Matrix Factorization/ALS is an unusual choice, as Rendle is primarily known for his 2010 work on Factorization Machines.

Answer: Thanks for pointing this out. Indeed, Rendle's 2022 work is less well-known and is a recent effort to revitalize ALS. Based on your suggestion, we have removed it; and we add a more complete references of matrix factorization and matrix completion in the revised Related Works (Appendix D).

Weakness 2: The practical applicability of these findings to real-world recommender systems is unclear. While the experimental section successfully demonstrates the theoretical bounds, it falls short of the standards expected in recommender systems research, lacking comprehensive evaluation on metrics and scenarios that matter in practice.

Answer: We apologize for the confusion. In Section 6 (Conclusion and Discussions), we have added a detailed discussion on the empirical implications and potential applications of our bounds in recommendation systems. We acknowledge that PAC-Bayesian bounds are derived for surrogate losses (e.g., sum of squared errors), but they can be indirectly related to recommendation evaluation metrics by decomposing the bound into two components: (1) a generalization bound on the surrogate loss, and (2) the empirical correlation between the surrogate loss and top-$k$ metrics.

This connection implies that PAC-Bayes bounds provide guarantees on the generalization of surrogate losses, which is a necessary condition for achieving strong downstream performance. While surrogate losses may not directly align with top-$k$ metrics such as Recall@$k$ or NDCG@$k$, demonstrating low generalization error on the surrogate loss establishes a strong theoretical foundation for the model’s ability to perform well on these evaluation metrics.

For instance, when evaluating a model on finite test sets using metrics like test error, NDCG@$k$, or Recall@$k$, these metrics often fluctuate depending on the specific data in the test set. If a model performs poorly on these metrics, PAC-Bayes bounds can help identify whether the poor results are likely due to model uncertainty (indicated by a large bound) or other factors; for example: 1) A large PAC-Bayes bound could indicate high model uncertainty, suggesting insufficient training or data sparsity. 2) A small PAC-Bayes bound coupled with poor performance might point to issues like sub-optimal surrogate metrics, data distribution shifts, or model design.

Thus, PAC-Bayes bound can help diagnose and improve model performance more effectively.

Weakness 3: The analysis primarily focuses on the non-regularized model, which has limited practical relevance as real-world systems invariably use regularization.

Answers: Actually, the bound primarily focuses on performance evaluation and does not directly guide the training process. See our reply https://openreview.net/forum?id=XYG98d5bCI&noteId=D4nDlcgnZB. The model $W$ used in the bound can come from any training process (whether optimized via EASE, EDLAE, or even obtained from a random initialization. Note that a randomly initialized $W$ is still a LAE model, though its parameters are unoptimized), with or without regularization. As long as $W$ is an $n \times n$ matrix, it can be considered as a LAE model, and its performance can be evaluated by the PAC-Bayes bound by comparing the test error and expected test error (true risk).

The relationship between the regularization hyperparameter $\gamma$ and the bound is indirect: Different values of $\gamma$ results in different $W$ (The relationship between $\gamma$ and $W$ can be complicated), and $W$ is a variable in the bound; but $\gamma$ itself is not a variable in the bound. Thus, the bound does not provide guidance on how to choose $\gamma$ to optimize it.

Broadly speaking, PAC-Bayes bound is based on the framework of statistical learning theory [Vapnik99]. In this theory, the regularizer $\gamma\|\|W\|\|_F^2$, with some constant $\gamma$, is equivalent to the constraint $\|\|W\|\|_F^2 < c$ for some $c$, which is part of the structural risk minimization principle [Vapnik91]. The structural risk minimization principle states that as $\gamma$ increases, $c$ decreases. A smaller $c$ leads to a smaller parameter space of $W$ (Note $\|\|W\|\|_F^2 < c$ implies $\|\|W - 0\|\|_F < \sqrt{c}$, meaning the parameter space is a ball of radius $\sqrt{c}$ under the Frobenius norm). Suppose the optimal $W^*$ (the one that minimizes the bound) lies in somewhere of the space. If $c$ is too large, the space of $W$ will be too large, and searching for the optimal $W^*$ will be hard (use any optimization algorithm). If $c$ is too small, the space of $W$ may exclude $W^*$, so you are not able to find it no matter how you search in the space.

However, there is no theoretical guideline for choosing $c$ in statistical learning theory (including in the PAC-Bayes bound). Therefore, $c$ is typically determined empirically, with the general recommendation that it should not be too large or too small.
 

References:

[Vapnik91] Vladimir Vapnik. Principles of risk minimization for learning theory. Advances in neural information processing systems, 4, 1991.

[Vapnik99] Vladimir, Vapnik. The Nature of Statistical Learning Theory. Springer, 1999.

Thanks for the clarifications, I appreciate the updates to the paper I still feel than "tight" would be a better wording than "non-vacuous" since I would reserve the later one to scenarios where the prediction made by the bound is actually used to improve a system at a lower cost than cross validation. But if the wording is accepted in the PAC community, I have no strong opposition. I update my score accordingly.

Dear Reviewer Yb5B,

We appreciate you updating the score and apologize for our late response and any inconvenience it may have caused.

Once again, we sincerely appreciate the time and effort you dedicated for reviewing our paper.

Best Regards,

Authors

Thank you for your insightful suggestions and questions.

Question 1.1: Repeat the analysis in the more general context of multi-output linear regression.

Answer: We have applied our analysis from LAE to the general multivariate (multi-output) linear regression model and propose a new bound (Theorem 1). It can be shown that Shalaeva's bound is a special case of our bound. Besides, our bound is based on a more general statistical assumption that the multivariate Gaussian distribution of $x$ can be dependent or degenerate, while in Shalaeva's bound it is assumed to be i.i.d..

Weakness 4: The bounds without the diagonal constraint are meaningless.

Answer: We completely integrate the constraint $\text{diag}(W) = 0$ to the bound, see Eq (5). The issue is that, the upper bound of part 2 of Eq (5) is difficult to calculate under the constraint $\text{diag}(W) = 0$ since it costs $O(n^4)$. We address this issue by finding an easy-to-calculate upper bound. This upper bound is established by replacing $\text{diag}(W)$ with i.i.d. zero mean Gaussian random variables, is validated by Theorem 4 (This theorem is newly added!), and costs $O(n^3)$.

Question 4: Incorporate the approximation of $\mathbb{E}[r^\top r]$ to the main bound.

Answer: We tackle this problem by providing two separate error bounds for estimating $\mathbb{E}[r^\top r]$ (which are not incorporated to the main bound): a theoretical bound based on Hoeffding's inequality and an empirical bound based on central limit theorem, see Appendix C. The theoretical bound is vacuous and thus not applicable in practice, while the empirical one is much tighter.

Weakness 1.2, Question 1.2, Question 3: The reason why taking $\lambda = m^{1/n}$ and causing $\lambda$ depending on $n$ is not mentioned. Also, why not taking $\lambda=\sqrt{m}$.

Answer: We have replaced $n$ with $d$. $\lambda$ can be independent of $n$. The setting $\lambda = m^{1/n}$ is from Shalaeva's paper and we realize it can cause confusion.

It is true that $\lambda=\sqrt{m}$ (i.e., $d = 2$) guarantees the convergence of $m^{-1/d}\ln \mathbb{E}_{\theta \sim \pi} \exp\left(2m^{2/d - 1} v\_{\theta}^2\right)$ as long as $\mathbb{E}\_{\theta \sim \pi}\exp\left(v\_{\theta}^2\right) < \infty$. We keep $d > 2$ because we want to show that the divergence $\lim\_{m\rightarrow \infty} \ln \mathbb{E}\_{\theta \sim \pi} \exp\left(2m^{2/d - 1} v\_{\theta}^2\right) = \infty$ cannot be addressed by taking any $d > 2$ if $\mathbb{E}\_{\theta \sim \pi}\exp\left(t v\_{\theta}^2\right) = \infty$ for any $t > 0$, see Section 3.2 (4).

Weakness 2, Question 1.3: Theorem 7 (Srebro's bound) is always vacuous for predicting high dimensional vectors, since it is originally used for low rank weight matrix on element-wise prediction. It should be removed altogether.

Answer: We realize that our bound makes an unfair comparison with Srebro's bound, as our bound assumes a nearly full-rank weight matrix (if there is no constraint $\text{diag}(W) = 0$, the Gaussian random matrix $W$ is of full rank [Feng07]) for vector-wise prediction. We have removed Theorem 7.

Weakness 6: It is unsuitable to call Shalaeva's convergence analysis an error.

Answer: We have revised our description of Shalaeva's convergence analysis, noting that it is not rigorous due to missing conditions.
 

Reference:

[Feng07] Xinlong Feng and Zhinan Zhang. The rank of a random matrix. Applied mathematics and computation, 185(1):689–694, 2007.

Dear Reviewer nWHo,

Thanks for reminding us! We have realized that the main task of LAE is not to predict the entire vector, but only the hidden items (the items that not shown in training stage but shown in test/evaluation stage) in a vector. Our bound needs to be revised to adapt to your settings before the experiments can be conducted. We have developed a new theoretical bound under your settings as well as its practical method for calculation, as shown below:

Changes in the Main Settings (Section 4.1):

We first rewrite the empirical risk function from $\frac{1}{m}\|\|R - RW\|\|_F^2$ to $\frac{1}{m}\|\|R^T - WR^T\|\|_F^2$ by taking the transpose and redefine $W^T$ as $W$, so that it aligns with out settings in Section 3. Assume $R\in \\{0, 1\\}^{m\times n}$ and $W \in \mathbb{R}^{n\times n}$.

Now we generate $X \in \\{0, 1\\}^{n\times m}$ and $Y \in \\{0, 1\\}^{n \times m}$ in this way: For any $i \in \\{1, 2, ..., n\\}$ and $j \in \\{1, 2, ..., m\\}$, if $(R^T)\_{ij} = 0$, we set $X\_{ij} = Y\_{ij} = 0$; if $(R^T)\_{ij} = 1$, we set either $X\_{ij} = 1, Y\_{ij} = 0$ or $X\_{ij} = 0, Y\_{ij} = 1$, each case with $1/2$ probability. In this way, we randomly split the $1$s in $(R^T)\_{\*j}$ into two groups $X\_{\*j}$ and $Y\_{\*j}$, and the number of $1$s in each group is roughly equal. Note that $X\_{\*j}$ corresponds to the $r\_{\text{train}}$ and $Y\_{\*j}$ corresponds to $r\_{\text{test}}$ in your settings.

Suppose we have a LAE model $W$ (The $W$ can be obtained through any training method, and we do not focus on the specific approach used, whether it's solved by EASE, EDLAE, or even a randomly initialized one. The constraint $\text{diag}(W) = 0$ can be applied but is not necessary.), its performance is evaluated by the empirical risk

$

R^{\text{emp}}(W) = \frac{1}{m}\|\|Y - WX\|\|\_F^2 \quad\quad\quad\quad \text{(R1)}

$

Assume each pair $(X\_{\*j}, Y\_{\*j})$ is i.i.d. sampled from an $2n$ dimensional distribution $\mathcal{D}$, then we can define the true risk as

$

R^{\text{true}}(W) = \mathbb{E}\_{(x, y)\sim \mathcal{D}}\left[\|\|y - Wx\|\|\_F^2\right] \quad\quad\quad \text{(R2)}

$

The new bound is obtained by replacing the old bound Eq (5) with the $R^{\text{true}}(W)$ in Eq (R2) and $R^{\text{emp}}(W)$ in Eq (R1).

Statistical assumptions for $\mathcal{D}$:

Just like we use Gaussian assumption for $\mathcal{D}$ in Assumption 1 in our paper, we need additional assumptions for the $\mathcal{D}$ in order to build our statistical model. Note that now the support of $\mathcal{D}$ is $\\{0, 1\\}^{2n}$, obviously not Gaussian. Thus we re-assume $\mathcal{D}$ as follows:
 

• Assumption 3: Suppose $\mathcal{D}$ is characterized by 3 finite cross-correlation matrices $\Sigma\_{xx} = \mathbb{E}\_{(x, y) \sim \mathcal{D}}[xx^T], \Sigma\_{xy} = \mathbb{E}\_{(x, y) \sim \mathcal{D}}[xy^T]$ and $\Sigma_{yy} = \mathbb{E}\_{(x, y) \sim \mathcal{D}}[yy^T]$, and $\Sigma\_{xx}$ is positive definite.
   

Under Assumption 3, given any $W$, we can expand $R^{\text{true}}(W)$ in the following way:

$

&\quad R^{\text{true}}(W) = \mathbb{E}\left[\|\|y - Wx\|\|\_F^2\right] = \sum_{i=1}^n\mathbb{E}[\|\|y\_i - W_{i\*}x\|\|_F^2] = \sum\_{i=1}^n W\_{i\*}\mathbb{E}[xx^T]W\_{i\*}^T - 2W\_{i\*}\mathbb{E}[y\_ix] + \mathbb{E}[y\_i^2] \\\\

&= \sum_{i=1}^n W_{i*}\Sigma_{xx}W_{i*}^T - 2W_{i*}(\Sigma_{xy})_{*i} + (\Sigma_{yy})_{ii} = \sum_{i=1}^n (W_{i*}\Sigma_{xx}^{1/2})(W_{i*}\Sigma_{xx}^{1/2})^T - 2(W_{i*}\Sigma_{xx}^{1/2})\Sigma_{xx}^{-1/2}(\Sigma_{xy})_{*i} + (\Sigma_{yy})_{ii} \\ &= \sum_{i=1}^n (W_{i*}\Sigma_{xx}^{1/2} - (\Sigma_{xy})_{*i}^T\Sigma_{xx}^{-1/2})(W_{i*}\Sigma_{xx}^{1/2} - (\Sigma_{xy})_{*i}^T\Sigma_{xx}^{-1/2})^T - (\Sigma_{xy})_{*i}^T\Sigma_{xx}^{-1}(\Sigma_{xy})_{*i} + (\Sigma_{yy})_{ii} \\ &= \sum_{i=1}^n ||W_{i*}\Sigma_{xx}^{1/2} - (\Sigma_{xy})_{*i}^T\Sigma_{xx}^{-1/2}||_F^2 - ||\Sigma_{xx}^{-1/2}(\Sigma_{xy})_{*i}||_F^2 + (\Sigma_{yy})_{ii} \\ &= ||W\Sigma_{xx}^{1/2} - \Sigma_{xy}^T\Sigma_{xx}^{-1/2}||_F^2 - ||\Sigma_{xy}^T\Sigma_{xx}^{-1/2}||_F^2 + \text{tr}(\Sigma_{yy}) \quad\quad\quad\quad \text{(R3)}

$

Experiments:

We obtain $W$ through minimizing the EASE function:

$

\min_{W}\|\|X - WX\|\|_F^2 + \gamma \|\|W\|\|_F^2 \quad\quad \text{s.t. } \text{diag}(W) = 0  \quad\quad\quad\quad \text{(R5)}

$

But note that this is not the only method to obtain $W$. We just use it as an example.

We also need to approximate $\Sigma\_{xx}, \Sigma\_{xy}$ and $\Sigma\_{yy}$, which is carried out using the same way in Section 4.4. Suppose there exists a test set $(X', Y')$ with $X' \in \\{0, 1\\}^{n \times m'}$ and $Y' \in \\{0, 1\\}^{n \times m'}$ and $(X\_{\*j}, Y\_{\*j})$ is i.i.d. from $\mathcal{D}$, Denote $X\cup X' = [X; X'] \in \\{0, 1\\}^{n\times (m + m')}$ and $Y\cup Y' = [Y; Y'] \in \\{0, 1\\}^{n\times (m + m')}$, then we have

$

\Sigma\_{xx} &\approx \hat{\Sigma}\_{xx} = \frac{1}{m + m'} [X \cup X'][X \cup X']^T \\\\
\Sigma\_{xy} &\approx \hat{\Sigma}\_{xy} = \frac{1}{m + m'} [X \cup X'][Y \cup Y']^T \\\\
\Sigma\_{yy} &\approx \hat{\Sigma}\_{yy} = \frac{1}{m + m'} [Y \cup Y'][Y \cup Y']^T

$

The error of the approximation above can be measured with the bounds in Appendix C of our paper. The $B$ and $c$ in Eq (R4) can be calculated with the approximated $\Sigma\_{xx}, \Sigma\_{xy}$ and $\Sigma\_{yy}$.

Denote ${R^{\text{emp}}}'(W) = \frac{1}{m'}\|\|Y' - WX'\|\|_F^2$, i.e., $R^{\text{emp}}(W)$ is the test error on the training set while ${R^{\text{emp}}}'(W)$ is the test error on the test set. Both $R^{\text{emp}}(W)$ and ${R^{\text{emp}}}'(W)$ are unbiased estimators of $R^{\text{true}}(W)$ so we can apply PAC-Bayes bound on them. The evaluation of the bound is implemented by comparing the gap between the PAC-Bayes bound and ${R^{\text{emp}}}'(W)$ (or $R^{\text{emp}}(W)$).

Remarks:

1. The Assumption 3 is very general since it not only includes the case when $\mathcal{D}$ is of binary support, like $\\{0, 1\\}^{2n}$, but also any bounded support, like $[a, b]^{2n}$.

2. Since we use EASE function Eq (R5) to obtain $W$, the training error can be expressed as $\frac{1}{m}\|\|X - XW\|\|_F^2$. However, it is hard to build a bound between the training error $\frac{1}{m}\|\|X - XW\|\|_F^2$ and the test error $\frac{1}{m}\|\|Y - XW\|\|_F^2$ since there are of different metrics and there seems to be no direct statistical connections between them.

3. It is also unnecessary to build a bound between the training error $\frac{1}{m}\|\|X - XW\|\|_F^2$ and the test error $\frac{1}{m}\|\|Y - XW\|\|_F^2$. As we mentioned in the answer of Question 5 of Reviewer Yb5B, the $W$ can be obtained by whatever training method, and PAC-bayes bound only evaluates the model's test performance. Since the PAC-Bayes bound can be independent of the training process, it is meaningless to build relationship between training error and testing error with a PAC-bayes bound, especially when training error and test error are under different metrics.

Please let us know if this solution is acceptable or you have any suggestions. Also, please give us one more day to carry out and post experiment results, as we need to rewrite programs and redo experiments. Thanks again for your patience!

To calculate the bound, we use Assumption 2 from our paper: Let $W$ be a Gaussian random matrix, with its prior distribution $\pi$ being $\bar{\mathcal{N}}(\mathcal{U}_0, \sigma^2J)$ and its posterior distribution $\rho$ being $\bar{\mathcal{N}}(\mathcal{U}, \mathcal{S})$.

Changes in calculating the optimal $\rho$ (Section 4.2):

The optimal $\rho$ obtained by Theorem 3 is revised as follows:

Without constraint $\text{diag}(W) = 0$, the optimal $\rho$ is given by

$

\mathcal{S}\_{ij} = \frac{1}{\frac{2\lambda}{m}X\_{i*}X_{i\*}^T + \frac{1}{\sigma^2}} \text{ for } i, j \in \\{1, 2, ..., n\\}

$

$

\mathcal{U} = \left(\frac{1}{m}XX^T + \frac{1}{2\lambda\sigma^2}I\right)^{-1}\left(\frac{1}{m}YX^T + \frac{1}{2\lambda\sigma^2}\mathcal{U}_0\right)

$

with constraint $\text{diag}(W) = 0$, the optimal $\rho$ is given by

$

\mathcal{S}_{ij} = \frac{1}{\frac{2\lambda}{m}X\_{i\*}X\_{i*}^T + \frac{1}{\sigma^2}} \text{ for } i, j \in \\{1, 2, ..., n\\}, i\ne j \quad\quad \mathcal{S}\_{ii} = 0 \text{ for } i \in \\{1, 2, ..., n\\}

$

$

\mathcal{U} = \left(\frac{1}{m}XX^T + \frac{1}{2\lambda\sigma^2}I\right)^{-1}\left(\frac{1}{m}YX^T + \frac{1}{2\lambda\sigma^2}\mathcal{U}\_0 - \frac{1}{2}\text{Diag}(x)\right)

$

where

$

x = 2 \cdot\textnormal{diag}\left[\left(\frac{1}{m}XX^T + \frac{1}{2\lambda\sigma^2}I\right)^{-1}\left(\frac{1}{m}YX^T + \frac{1}{2\lambda\sigma^2}\mathcal{U}\_0\right)\right] \oslash \textnormal{diag}\left[\left(\frac{1}{m}XX^T + \frac{1}{2\lambda\sigma^2}I\right)^{-1}\right]

$

The proof is similar to the original one.

Changes in calculating $\Psi$ (Section 4.3):

Since by Eq (R3), $R^{\text{true}}(W) = \|\|W\Sigma\_{xx}^{1/2} - \Sigma\_{xy}^T\Sigma_{xx}^{-1/2}\|\|\_F^2 - \|\|\Sigma\_{xy}^T\Sigma\_{xx}^{-1/2}\|\|\_F^2 + \text{tr}(\Sigma\_{yy})$, denote $B = \Sigma\_{xy}^T\Sigma_{xx}^{-1/2}$ and $c = - \|\|\Sigma_{xy}^T\Sigma\_{xx}^{-1/2}\|\|\_F^2 + \text{tr}(\Sigma\_{yy})$, then

$

R^{\text{true}}(W) = \sum\_{i=1}^n\|\|W\_{i\*}\Sigma\_{xx}^{1/2} - B\_{i\*}\|\|\_F^2 + c \quad\quad\quad\quad \text{(R4)}

$

If $W \sim \bar{\mathcal{N}}(\mathcal{U}\_0, \sigma^2J)$ with $\text{diag}(W) = 0$ constraint, then $W\_{\*i}^T \sim \mathcal{N}((\mathcal{U}\_0)\_{i\*}^T, \sigma^2(I - I^i))$ and $(W\_{i\*}\Sigma_{xx}^{1/2} - B\_{i\*})^T \sim \mathcal{N}(\Sigma_{xx}^{1/2}(\mathcal{U}\_0)\_{i\*}^T - B\_{i\*}^T , \sigma^2(\Sigma\_{xx} - [(\Sigma\_{xx})\_{i\*}^{1/2}]^T[(\Sigma\_{xx})\_{i\*}^{1/2}]))$.

Therefore, if in Section 4.3 we redefine $\mu^i = \Sigma\_{xx}^{1/2}(\mathcal{U}\_0)\_{i\*}^T - B\_{i\*}^T$ and $A = \sigma^2\Sigma\_{xx}$, then the upper bound of $\Psi$ given by Theorem 4 can be calculated in the same way.

Let $A = \sigma^2\Sigma\_{xx} = S^T\Lambda S$ be the eigenvalue decomposition where $\Lambda = \text{diag}(\eta_1, \eta_2, ..., \eta_n)$. The upper bound of $\Psi$ now becomes

$

&\Psi_{\pi, \mathcal{D}}(\lambda, m) \le \ln \mathbb{E}\left[e^{R^{\text{true}}(W)}\right] = \ln \mathbb{E}\left[\exp\left(||W\Sigma_{xx}^{1/2} - B||_F^2 + c\right)\right] = \ln \mathbb{E}\left[\exp\left(\sum_{i=1}^n||W_{i*}\Sigma_{xx}^{1/2} - B_{i*}||_F^2 + c\right)\right] \\ &= c + \ln \mathbb{E}\left[\exp\left(\sum_{i=1}^n||W_{i*}\Sigma_{xx}^{1/2} - B_{i*}||_F^2\right)\right] \le c + \ln \prod_{i=1}^n \frac{\exp\left(\frac{\frac{1}{\sigma^2}\lambda\eta_i||S_{i*}(\mathcal{U}_0^T - \Sigma_{xx}^{-1/2}B^T)||_F^2}{1 - 2\lambda \eta_i}\right)}{(1 - 2\lambda \eta_i)^{n/2}}

$

If without the constraint $\text{diag}(W) = 0$, the last $\le$ above becomes $=$.

Dear Reviewer nWHo,

Thanks for your patience! Below are our experimental results on MovieLens 20M and Netflix datasets.
 

MovieLens 20M:

Netflix:

$\gamma$ is the hyperparameter of EASE, and we set $\gamma = 50, 100, 200, 400$ to train different EASE models $W$ for evaluation. The experimental results show that our bound is tight, not exceeding $3\times$ of $R^{\text{emp}}(W)$ or ${R^{\text{emp}}}'(W)$.

The settings used in the experiment are mostly the same as in our paper: We split the $85\\%$ samples from the entire dataset into the training set $R$ and the rest $15\\%$ into the test set $R'$. Then we generate $X, Y$ from $R$ and $X', Y'$ from $R'$. The only difference is that we have transformed the range of ratings in the datasets from $[0, 5]$ to $[0, 1]$ by setting all non-zero entries to $1$. This change results in smaller test errors compared to the experimental results in our paper.

We use an $80-20$ splitting strategy ($80\\%$ of the 1s in $R \cup R'$ randomly assigned to $X \cup X'$ and $20\\%$ to $Y \cup Y'$) for both MovieLens 20M and Netflix. It needs to be mentioned that both MovieLens 20M and Netflix are sparse datasets, and MovieLens 20M is sparser than Netflix. We observed an issue with MovieLens 20M dataset: using a $50-50$ splitting strategy makes $X \cup X'$ too sparse, causing $\hat{\Sigma}\_{xx} = \frac{1}{m + m'} [X \cup X'][X \cup X']^T$ to be singular. Since $\Sigma\_{xx}^{-1/2}$ in Eq (R3) is approximated by $\hat{\Sigma}\_{xx}^{-1/2}$, a singular $\hat{\Sigma}\_{xx}$ will result in an undefined $\hat{\Sigma}\_{xx}^{-1/2}$, which makes the bound uncomputable. Thus one limitation when calculating the theoretical bound in practice is that $X \cup X'$ must be sufficiently dense to ensure that $\hat{\Sigma}_{xx}$ is non-singular.

It should also be noticed that $\hat{\Sigma}\_{xx}^{-1/2}$ being singular does not imply $\Sigma_{xx}^{-1/2}$ is singular. While $\Sigma\_{xx}^{-1/2}$ is unknown to us, if its non-singularity holds true, then given more data, we can obtain a non-singular approximation $\hat{\Sigma}\_{xx}^{-1/2}$ for $\Sigma\_{xx}^{-1/2}$. This suggests that the singularity of $\hat{\Sigma}\_{xx}$ is likely due to insufficient data for a non-singular approximation. Our experiments on MovieLens 20M support this. We found that if $X'$ is not included in calculating $\hat{\Sigma}\_{xx}$, i.e., simply use $\hat{\Sigma}\_{xx} = \frac{1}{m}XX^T$, then a $95 - 5$ splitting fails. If including $X'$, then an $80 - 20$ splitting works, while a $70 - 30$ splitting fails. This implies that adding more data can help address the singularity issue in the approximation $\hat{\Sigma}\_{xx}$.

In summary, if the approximation $\hat{\Sigma}\_{xx} = \frac{1}{m}XX^T$ is singular, we suggest two ways for addressing this issue: (1) Split more 1s to $X$, ensuring that $X$ is sufficiently dense to guarantee the non-singularity of $\frac{1}{m}XX^T$. (2) Introduce additional data $X'$ (sampled from the same distribution as $X$) and set $\hat{\Sigma}\_{xx}$ as $\frac{1}{m+ m'}[X \cup X'][X \cup X']^T$.

If you are satisfied with these updated results, we would be happy to include them in the paper. We kindly ask for your feedback or any further suggestions.

Question 3, Weakness 3: How the bound explains the fluctuation in performance caused by hyperparameter tuning, and how it informs model selection, regularization strategies, or error expectations.

Answer: The bound does not directly guide the training process of a model; instead, it provides information about the model's test performance. Specifically, once we obtain a model $W$ from any training method (whether optimized via EASE, EDLAE, or even obtained from a random initialization. Note that a randomly initialized $W$ is still a LAE model, though its parameters are unoptimized.), as long as $W$ an $n\times n$ matrix, it can be treated as a LAE model, and we can evaluate its test performance with the bound.

The changes of hyperparameters affect the bound indirectly: Different hyperparameter results in different $W$, and $W$ is a variable of the bound. There is no direct relationship between the hyperparameters and the bound, as the hyperparameters are not included as variables in the bound.

The bound does not provide guidance on hyperparameter tuning or regularization strategies, as these relate to the training process. The bound informs model selection in a way that a model $W$ with a smaller PAC-bayes bound is better, since it has better generalizability, i.e., likely to make more accurate prediction on unseen data. Regarding error expectations, the bound provides a probabilistic relationship between test error and the expected test error (true risk); it is independent of training error.

Weakness 1.1: Given that there may be space available, consider moving the related work section from the appendix (if present) to the main text.

Thanks for your suggestions. We have significantly revised the paper and there is no additional space for related works due to the 10-page limit. We have also expanded the related works by adding papers related to the matrix completion and trace norm.

Weakness 1.2: A more in-depth discussion on how the proposed bounds compare to existing work would make the paper more comprehensive and inclusive.

To the best of our knowledge, the PAC-Bayes bounds for linear regression are primarily covered in two papers. [Germain16] proposed the first bound for multiple linear regression (single dependent variable), but the bound does not guarantee convergence. [Shalaeva20] improved upon [Germain16]'s bound by ensuring convergence, although the conditions required for convergence are not discussed. Both [Germain16] and [Shalaeva20] focus on single-output linear regression. In the revision, we have proposed the first PAC-Bayes bound for multivariate linear regression, where the bound in [Shalaeva20] is a special case of our general bound. In addition, we have provided sufficient conditions that guarantee convergence of both our bound and Shalaeva's bound. Finally, we have applied our bound for study LAE on the zero-diagonal case. We have highlighted and clarified the above problems and contributions/differences in the revised Section 1 (Introduction).
 

References:

[Germain16] Pascal Germain, Francis Bach, Alexandre Lacoste, and Simon Lacoste-Julien. Pac-bayesian theory meets bayesian inference. Advances in Neural Information Processing Systems, 29, 2016.

[Shalaeva20] 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, pp. 5660–5667, 2020.

We would like to first thank you for your insightful reviews. Our responses to your questions and concerns are listed below.

Question 1, Weakness 2: It is unclear whether the proposed PAC-Bayes bounds can be generalized to all linear autoencoder models. The paper should clarify whether these bounds are specific to the conditions outlined or if they have broader applicability across different linear autoencoder architectures. For example, is it possible to calculate a general PAC-Bayes bound for any linear autoencoder using Theorem 6? If there are limitations or constraints, it would be valuable to elaborate on them to clarify the scope of the theoretical results.

Answer: Our PAC-bayes bound can be generalized to any LAE model $W$, including those with the zero diagonal constraint $\text{diag}(W) = 0$, as discussed in our reply https://openreview.net/forum?id=XYG98d5bCI&noteId=D4nDlcgnZB . In fact, any $n \times n$ matrix $W$, regardless of how it is obtained, is an LAE model, and our bound can be applied to it. Several prominent LAE-based recommender system models like EASE, EDLAE and ELSA adopt zero diagonal constraint, to which our bound is also applicable.

One limitation of our bound is that some LAE models may impose specific constraints on $W$, making the bound difficult to compute. For example, some LAE models require $W$ to be of low rank. While the bound can still be formed for a low rank $W$, its calculation can be difficult (We may end up with a theoretical bound that cannot be computed). This is because, the calculation of the bound relies Assumption 2 in our paper, where we assume $W$ to be a random Gaussian matrix ($W$ does not have to be Gaussian, but assuming it to be Gaussian makes the bound easy to calculate, while other distributions may not). A random Gaussian matrix $W$ is of full rank [Feng07] and cannot be reduced to the low rank case.

Another limitation is that our bound uses square loss (Frobenius norm) as the metric for the error between target and prediction, since the square loss is easy for statistical analysis. Many real-world recommender systems, however, use other metrics like NDCG@K, Recall@K, and our bound cannot be directly be applied to these metrics. Moreover, these metrics are challenging to analyze from a statistical perspective.

Question 2.1: The methodology used to calculate the PAC-Bayes bound in Table 2 requires further elaboration. Did you use a fixed $\epsilon$ in your calculations? If so, what was the reasoning behind this choice, and how might different $\epsilon$ values affect the results?

Answer: Thanks for your suggestions. The calculation of $\epsilon$ is related to Theorem 7, which we have removed. Due to Reviewer nWHo's suggestion and we also agreed, there is some unfairness in directly comparing Srebro's bound with the new bound as the settings and goals of these bounds are different: the former applies to low-rank weight matrices in element-wise prediction (explicitly recommendation setting), while our bound applies to nearly full-rank weight matrices in vector-wise prediction (implicit recommendation setting). Given this, we have removed Theorem 7.

Question 2.2: How sensitive are the practical training and test errors to hyperparameter choices? If there are significant fluctuations, a discussion on how to interpret these variations (e.g., how fluctuations influence conclusion about 2x or 4x test error) would be very helpful.

Answer: We added experiments where we set $\gamma$ (the only hyperparameter of EASE, see Eq (2)) to 50, 200, and 400. The results show that our bound remains within 3 times the test error on MovieLens 20M, Netflix, and MSD, and 4 times the test error on Yelp2018, despite fluctuations. Thus, this demonstrates our results are not sensitive due to the hyperparameter choices.

(Update) We have some new experimental results on a different LAE setting, which can be found here: https://openreview.net/forum?id=XYG98d5bCI&noteId=T7jUjy5Ige
 

Reference:

[Feng07] Xinlong Feng and Zhinan Zhang. The rank of a random matrix. Applied mathematics and computation, 185(1):689–694, 2007.