On the Regularization Landscape for the Linear Recommendation Models Dong Li1
Zhenming Liu2
Ruoming Jin1
Hao Zhou1
Zhi Liu3
Jing Gao3
Bin Ren2
arXiv:2609.11876v1 [cs.AI] 10 Sep 2026
1 Kent State University, USA 2 College of William and Mary, USA 3 iLambda, USA
{dli12,rjin1,hzhou6}@kent.edu {zliu,jgao}@ilambda.com {zliu,bren}@cs.wm.edu
Abstract Recently, a wide range of recommendation algorithms inspired by deep learning techniques have emerged as the performance leaders on several standard recommendation benchmarks. While these algorithms were built on different DL techniques (e.g., dropouts, autoencoder), they have similar performance and even similar cost functions. This paper studies whether the models’ comparable performance are sheer coincidence, or they can be unified under a single framework. We find that all linear performance leaders effectively add only a nuclear-norm based regularizer, or a Frobenius-norm based regularizer. The former ones possess a (surprising) rigid structure that limits the models’ predictive power but their solutions are low rank and have closed form. The latter ones are more expressive and more efficient for recommendation but their solutions are either full-rank or require executing hard-to-tune numeric procedures such as ADMM. Along this line of finding, we further propose two low-rank, closed-form solutions, derived from carefully generalizing Frobenius-norm based regularizers. The new solutions get the best of both nuclear-norm and Frobenius-norm world.
1
Introduction
Research progress on algorithms for recommendation has escalated in recent years, partially fueled by the adoption of deep learning techniques. However, recent studies have found that many new deep learning recommendation models have shown sub-par performance against the simpler linear recommendation models [Dacrema et al., 2019a, Rendle et al., 2019]. Although some studies are available to analyze linear vs non-linear models [Dacrema et al., 2019a], it remains puzzling why these seemingly different techniques all result in models with similar performance or even similar cost functions. This motivates us to ask a fundamental question: What is the barebones engine that drives the performance improvement for the recent recommendation algorithms? Specifically, do different recommendation techniques offer different “magic”, but coincidentally have similar performance, or can they be unified under a single framework, which has a potential to produce one single algorithm that gets the best of all existing works? In the latest study, Jin et al. [2021] examines the relationship between the widely used matrix factorization (MF), such as ALS [Warlop et al., 2017], and the linear
1
autoencoders (LAE) which encompasses the recent performance leaders, such as EASE [Steck, 2019] and EDLAE [Steck, 2020]. They consider two basic regularization forms (See Eq (1) and (6)) and found that the optimal (closed-form) solutions of both models recover the direction of principal components, while shrinking the corresponding singular values differently. They suggests this difference may enable LAE to utilize a larger number of latent dimensions to improve recommendation accuracy, and use this to highlight the similarity as well as difference between LAE and MF. In this paper, we go much beyond the two basic models studied in Jin et al. [2021] to analyze a large number of recent performance leaders of (linear) recommender algorithms. We found they all can be categorized into those that implement nuclear-norm regularizers, and into those that implement Frobenius-norm regularizers. We found that the former ones possess a (surprising) rigid structure that limits the models’ predictive power, and the latter ones tend to be more expressive and more effective for recommendation. In many cases, both regularizers recover the direction of principal components, while shrinking the corresponding singular values differently. Interestingly, we observe that it is not matrix factorization or LAE that determines the shrinkage structure (as Jin et al. [2021] suggested), but instead it is the forms of regularization. Thus, this paper provides a more complete and accurate characterization on how a linear recommendation model performs under different regularizations. To better understand the regularizations that can be transformed into the (weighted) nuclear-norm regularizer ∥𝑊 ∥ ∗ , we first show that Variational Linear AutoEncoders (VLAE) solves a weighted nuclear-norm regularization problem, in which the weights possess a specific combinatorial structure. We also show that this technique cannot be generalized to tackle arbitrary weight sequences. Second, it has been known that using dropout techniques is equivalent to adding a squared nuclear-norm regularizer ∥𝑊 ∥ 2∗ , and that the solution structure is strikingly similar to the regularizer that uses ∥𝑊 ∥ ∗ . we generalize the result to show that the solution structures for ∥𝑊 ∥ ∗𝑝 are highly similar for all 𝑝 ≥ 1. But when 𝑝 = 1, 2, the solution and hyper-parameters possess favorable properties so hyper-parameter search becomes easier. This also partially explains why only 𝑝 = 1, 2 have been extensively considered. Third, all nuclear-norm–based techniques possess a salient property: their estimators keep the singular vectors of the data matrix and shrink only the singular values. But this severely limits the search space and explains why models that use only nuclear-norm–based regularizers share the same performance ceiling even when hyper-parameters are extensively searched. The (weighted) Frobenius-norm regularizers ∥Λ𝑊 ∥ 2𝐹 are implemented in EASE [Steck, 2019] and EDLAE [Steck, 2020]. These models produce closed form full-rank estimators; and if the zero diagonal constraint on 𝑊 is enforced, their singular vectors will no longer coincide with those of the data, and can deliver (slightly) better performance. However, no closed form solutions for the low-rank estimator is known and the current approaches rely on ADMM or stochastic factorized gradients [Steck, 2020]. In this paper, we propose two new low-rank, closed-form estimators that deliver comparable results to the full rank models (EASE and full-rank EDLAE) as well as the ADMM based solutions [Steck, 2020]. The new closed-form solutions for low-rank models have profound implications at both practical and conceptual fronts: First, low rank solutions are often more scalable (the full rank 𝑊 can be too large to materialize) and can better serve real-world training and deployment. Second, perhaps more excitingly, these solutions get the best of nuclear (closed form and low rank) and Frobenius worlds (strong predictive power). A simple “one-liner” formula (from each solution) concisely pack all the benefits obtained by a recent long line of research and abstract out all the computation nuance (e.g., the need to tune ADMM and deal with local optimal). We believe these closed-form solutions are also powerful tools to help us compare different regularizations analytically.
2
2
Background and overview
Recommendation system background. Recommendation algorithms can be categorized into explicit ones that aim to predict unseen ratings between a user and an item and implicit ones that aim to predict actions, such as user click or add-cart [Steck, 2019, Dacrema et al., 2019b, Zhang et al., 2019]. We focus on the implicit problem because it is more economically relevant. Here, let 𝑛 be the number of items and 𝑚 be the number of users. We are given a binary matrix 𝑋 ∈ {0, 1} 𝑚×𝑛 that represents the interaction between users and items so far, i.e., 𝑋𝑖, 𝑗 = 1 iff user 𝑖 has purchased or made a rating on item 𝑗. Our goal is to produce a ˆ which we evaluate against future interactions using information retrieval metrics such real-valued matrix 𝑋, as Top-𝑘 or nDCG. Note that while the syntax of this problem resembles matrix completion (MC) [Candès and Tao, 2010], recommendation systems and MC have different evaluation criteria so MC’s results are not directly applicable here.
2.1
Nuclear-norm based regularizations
Let 𝑋 ∈ R𝑚×𝑛 be a matrix of rank at most 𝑘 with 𝑘 leading singular values being 𝜎1 (𝑋) ≥ 𝜎2 (𝑋) ≥ · · · ≥ 𝜎𝑘 (𝑋). Let 𝜔 = (𝜔1 , . . . , 𝜔 𝑘 ) ∈ (R+ ) 𝑘 . The weighted unclear norm of 𝑋 with respect to 𝜔 is defined as Í𝑘 ∥ 𝑋 ∥ 𝜔,∗ = 𝑖=1 𝜔𝑖 × 𝜎𝑖 (𝑋). We can see that this is a natural generalization of the weighted nuclear norm for low-rank matrices [Gu et al., 2014]. Also, despite its name, the weighted nuclear norm is neither convex nor differentiable unless 𝜔𝑖 ’s are sorted in descending order [Chen et al., 2013, Iglesias et al., 2020]. Nuclear-norm based regularizers perform ℓ1 -shrinkage over the estimator’s singular values, which resembles performing ℓ1 -shrinkage for coefficients in a linear model in LASSO [Tibshirani, 1996]. Therefore, Nuclearnorm regularizers also promote sparsity over the solution’s singular values (i.e., the solution is usually low rank). We note a large fraction of recent recommendation algorithms effectively add only a nuclear-norm regularizer to MF (see also Table 3 in Appendix): A1. Regularized PCA [Udell et al., 2016, Zheng et al., 2018] aims to solve min ∥ 𝑋 − 𝑃𝑄 𝑇 ∥ 2𝐹 + 𝜆∥𝑄∥ 2𝐹 + 𝜆∥𝑃∥ 2𝐹 .
(1)
𝑃,𝑄
It has been known that equation (1) is equivalent to solving min𝑋ˆ ∥ 𝑋 − 𝑋ˆ ∥ 2𝐹 + 2𝜆∥ 𝑋ˆ ∥ ∗ . To solve equation (1), one can use factored gradient descent [Bhojanapalli et al., 2016] or directly use its closed form solution [Kunin et al., 2019], which involves computation of SVD of 𝑋. A2. Matrix Factorization via dropouts. This approach use 𝑃𝑄 T (𝑃 ∈ R𝑚×𝑘 , 𝑄 ∈ R𝑛×𝑘 ) to approximate 𝑋 and uses a neural net to find 𝑃 and 𝑄. A standard dropout technique is used when we train 𝑃 and 𝑄. Cavazza et al. [2018] shows that optimization with dropout is equivalent to solving min𝑋ˆ ∥ 𝑋 − 𝑋ˆ ∥ 2𝐹 + 𝜆∥ 𝑋ˆ ∥ 2∗ , and the closed form solution is obtained by shrinking all singular values of 𝑋 by a magnitude of 𝜇, which depends on the data 𝑋 and choise of 𝜆. The next two approaches use techniques from (variational) auto-encoder. We show that they also effectively add variants of nuclear-norm regularizers although this may not be clear at the first glance. A3. Linear Regression via Denoising Linear Auto-Encoder considers the following (non-uniform) weighted ℓ2 -regularization [Bao et al., 2020]: 1
1
∥ 𝑋 − 𝑋𝑊1𝑊2 ∥ 2𝐹 + ∥𝑊1 Λ 2 ∥ 2𝐹 + ∥Λ 2 𝑊2 ∥ 2𝐹 , 3
(2)
where Λ is a diagonal matrix. Bao et al. [2020] has shown a closed-form for eq. 2 with a specific of diagonal Λ = 𝑑𝑖𝑎𝑔(𝜆1 , 𝜆2 , · · · , 𝜆 𝑘 ) when the weight is non-descending: 𝜆1 ≤ 𝜆 2 ≤ · · · ≤ 𝜆 𝑘 . It remains unclear if a closed-form exists for an arbitrary weight order. A4. Variational Linear Auto-encoders. While not explicitly studied before, it is also natural to consider linear simplification of Variational Autoencoders, such as Multi-VAE [Liang et al., 2018], which has shown to exhibit strong performance for recommendation. To optimize linear VAE, we need to find the MLE for the probabolistic model 𝑝(x | z) = N 𝑊z + µ, 𝜎 2 𝐼 and 𝑝(z | x) = N (𝑉 (x − µ), 𝐷), (3) where 𝐷 is a diagonal covariance matrix. Then we have the following observation. Lemma 1. Consider optimizing the ELBO (Evidence Lower Bound) [Kingma and Welling, 2014] for the above LVAE model (Eq.3). When the optimization is over the entire dataset and µ = 0, this optimization problem is equivalent to minimizing √ (4) L = ∥ 𝑋 − 𝑋𝑉 T𝑊 T ∥ 2𝐹 + 𝑁 ∥ 𝐷𝑊 𝑇 ∥ 2𝐹 + 𝜎 2 ||𝑋𝑉 T ∥ 2𝐹 + 𝑔(𝐷, 𝜎) where, 𝑔(𝐷, 𝜎) = −𝜎 2 𝑁 log |𝐷 | − tr(𝐷) + 𝑘 − 𝑛 log 2𝜋𝜎 2 ). Here, we set µ = 0 for simplicity and following typical practices in recommendations. The proof can be found in Appendix. Since our objective is for recommendation (not purely on recovering the low rank factors of the data), we treat the covariance matrix √ 𝐷 as hyperparameters. Thus, the term 𝑔(𝐷, 𝜎) becomes a constant, and let 𝐴 = 𝜎𝑉 𝑇 , 𝐵 = 1/𝜎𝑊 𝑇 , Λ = 𝜎 𝑁 𝐷. If we consider 𝐷 as an optimization parameter, then the optimal solution of LVAE is equivalent to that of pPCA [Tipping and Bishop, 1999]. In this case, 𝐷 = 𝜎 2 (𝚺2 /𝑁) −1 , where 𝚺2 = 𝑑𝑖𝑎𝑔(𝜎𝑖2 ) are the eigen-values of 𝜎 2𝑗 1 Í𝑛 the covariance matrix of 𝑋, and 𝜎 2 = 𝑛−𝑘 𝑗=𝑘+1 𝑁 . Further, the closed form solutions of 𝑉 (𝐴) and 𝑊 (𝐵) are characterized. However, for recommendation, the matrix 𝐷 can be considered as a hyperparameter (to be learned); in this case, the closed form solution is not studied yet. We may “clean up” equation (4) and obtain the following optimization problem: min
∥ 𝑋 − 𝑋 𝐴𝐵∥ 2𝐹 + ∥ 𝑋 𝐴∥ 2𝐹 + ∥Λ𝐵∥ 2𝐹 ,
(5)
𝐴∈R𝑛×𝑘 ,𝐵∈R 𝑘×𝑚
in which decision variables are 𝐴 and 𝐵, and the hyper-parameter is a diagonal matrix Λ ∈ R 𝑘×𝑘 . Our main results: solution structure and implications (Sec. 3). (i) We show that solving Eq. 5 (A4) is equivalent to solving ∥ 𝑋 − 𝑊 ∥ 2𝐹 + 𝜆∥𝑊 ∥ 𝜔,∗ subject to rank(𝑊) ≤ 𝑘, where 𝜔 consists of Λ’s diagonal values, sorted in ascending order. In addition, the closed form solution for 𝑊 is merely shrinking the 𝑖-th singular value of 𝑋 by a magnitude of 𝜔𝑖 . When the diagonals of Λ is already sorted (i.e., Λ1,1 ≤ · · · ≤ Λ 𝑘,𝑘 ), the problem effectively reduces to A3. When Λ is proportional to identity, the problem reduces to A1. This result has three major implications. First, A1, A3, and A4 effectively only add a variant of nuclear-norm based regularizer, and the major benefits from these algorithms are computational. Second, our result generalizes that in A3 and solves an open problem left there, i.e., the hyper-parameters Λ do not need to have sorted diagonal values, the optimization algorithm will “automatically perform the sorting”. Third, while the variational auto-encoder offers a flexibility to tailor-make the prior for each entry in the latent variable z (in eq. 3), the solution space is quite rigid due to the auto-sorting property: the 𝑖-th smallest entry in Λ will find its way to match with the 𝑖-th largest singular values in 𝑋. In other words, it is impossible to shrink 𝑋’s singular values by an arbitrary sequence via carefully choosing Λ in eq. 5. 4
(ii) We characterize the optimal solution for ∥ 𝑋 − 𝑊 ∥ 2𝐹 + 𝜆∥𝑊 ∥ ∗𝑝 for any 𝑝 ≥ 1. We shall show that regardless the choice of 𝑝, the optimal solution uniformly shrinks all 𝑋’s singular values by a constant magnitude 𝜇 (and to 0 if a singular value is already less than 𝜇). The specific 𝜇 depends on 𝜆, 𝑝, as well as the data 𝑋 unless 𝑝 = 1. It has two implications. First, 𝜇 needs to be fine-tuned to optimize test performance. Therefore, choices of 𝑝 (again) only produces computational gain. Second, when 𝑝 = 1, 𝜇 does not depend on the data so we can tune this hyper-parameter in a direct manner. When 𝑝 = 2, 𝜆 is scale-invariant, i.e., it does not need to be rescaled when all entries in 𝑋 is scaled by a constant factor. Being able to directly tune 𝜇 or having the scale invariant property helps the hyper-parameter search; when 𝑝 ≠ 1, 2, it does not offer benefit in either computation or search, which explains why we see only 𝑝 = 1, 2 in the literature. Finally, for all nuclear-norm-based approaches discussed above, the estimators always keep singular vectors of 𝑋 and shrink its singular values. Therefore, the solution space offered by nuclear-norm based regularization is quite constrained, which limits these models predictive power.
2.2
Frobenius norm based regularizations
Most algorithms below were originally motivated by the design of (denoising) auto-encoders, it has been shown that they effectively add a Frobenius-norm regularizer. See also Table 3 in Appendix. A5. EASE [Steck, 2019] aims to optimize min𝑊 ||𝑋 − 𝑋𝑊 || 2𝐹 + 𝜆 · ||𝑊 || 2𝐹 subject to the constraint that diag(𝑊) = 0. A closed form solution exists for this problem. A6. DLAE [Steck, 2020] adds a weighted Frobenius-norm regularizer so the objective becomes min ∥ 𝑋 − 𝑋𝑊 ∥ 2𝐹 + ∥Λ1/2𝑊 ∥ 2𝐹 , 𝑊
where Λ = 1−𝑝𝑝 𝑑𝑖𝑎𝑔𝑀 (𝑑𝑖𝑎𝑔(𝑋 𝑇 𝑋)). A7. EDLAE [Steck, 2020] integrates weighted Frobenius norm in DLAE with EASE’s diagonal constraint so its objective is the same as DLAE but it requires diag(𝑊) = 0. A closed form solution exists for this problem. When 𝑊 is required to be low rank, an ADMM algorithm may be used. A8. Tikhonov regularization/Low Rank Regression (LRR) [Jin et al., 2021]. Let 𝑉𝑘 be the 𝑘 leading right singular vectors of 𝑋. Jin et al. [2021] finds an estimator that solves 𝑊 = arg
min
𝑟 𝑎𝑛𝑘 (𝑊 ) ≤ 𝑘
∥ 𝑋 − 𝑋𝑊 ∥ 2𝐹 + ∥𝚪𝑊 ∥ 2𝐹 ,
(6)
1
where Γ = Λ 2 𝑉𝑘𝑇 and Λ = 𝑑𝑖𝑎𝑔(𝜆′1 , · · · , 𝜆′𝑘 ) is a hyper-parameter. Its a closed-form solution is. ∗
𝑊 = 𝑉𝑘 𝑑𝑖𝑎𝑔(
𝜎12 𝜎12 + 𝜆′1
𝜎𝑘2 ,...,
𝜎𝑘2 + 𝜆′𝑘
)𝑉𝑘𝑇
(7)
Our results. Our major goal is to design a low-rank closed form estimator whose performance is comparable to the performance leaders. We first remark that A5 and A6 produce full-rank estimators. A7 can produce either full-rank or low-rank estimator (via ADMM) and has the best performance (among all approaches we discussed). The estimator from A8 is low-rank and has a closed-form solution but it has to keep singular vectors of 𝑋 so its predictive power is also limited. Nevertheless, A8 is conceptually interesting because it uses Frobenius norm regularizers but its solution space cover the solution space offered in A4 (and thus also A1-A4). 5
Proposition 1. For any regularized instances in the form of (5) with regularization parameter Λ such that 1 𝜎𝑖 (𝑋) ≥ 𝜆 (𝑘−𝑖) for all 𝑖, there is a corresponding Tikhonov regularized instance with 𝚪 = Λ 2 𝑉𝑘T which provides the same regularization effect. The proof is in Appendix. A major implication of Prop. 1 is that we can focus on designing Frobenius-norm regularizers because it also gets the value from using nuclear-norm regularizers. Indeed, Sec. 4 will introduce two low-rank Frobenius-norm-based model with closed form solutions that have comparable performance to linear performance leaders.
3
Nuclear-norm based regularization
Rigidity of VLAE. We first analyze solution for Eq. 5 (A4). To facilitate the analysis, we also consider the following problem: 1
1
min ∥ 𝑋 − 𝑃𝑄∥ 2𝐹 + ∥Λ 2 𝑄∥ 2𝐹 + ∥𝑃 Λ 2 | 2𝐹 , or equivalently , min ∥ 𝑋 − 𝑃𝑄∥ 2𝐹 + ∥Λ𝑄∥ 2𝐹 + ∥𝑃∥ 2𝐹 , 𝑃,𝑄
(8)
𝑃,𝑄
Note that when we let 𝑄 ′ = 𝑋 𝐴∗ and 𝑃′ = 𝐵∗ (where 𝐴∗ and 𝐵∗ are an optimal solution of Equation 5), the syntax of Eq. 8 matches with that of Eq. 5. This implies that solution for Eq. 8 is a lower bound of that for Eq. 5. These two solutions coincide only when the columns in the optimal 𝑃∗ in Eq. 8 are spanned by the columns of 𝑋. We shall first find a closed-form solution (𝑃∗ , 𝑄 ∗ ) for Eq. 8, and show that indeed that the column space of 𝑃∗ is in the column space of 𝑋. Below is our major Proposition. Proposition 2. Let 𝑓 : R𝑚×𝑛 → R+ be any cost function. Let 𝑃 ∈ R𝑚×𝑘 and 𝑄 ∈ R 𝑘×𝑛 . Let Λ ∈ R 𝑘×𝑘 be a diagonal matrix such that 𝜆𝑖 = Λ𝑖𝑖 ≥ 0 (𝑖 ∈ [𝑘]). Let also 𝜔 = (𝜆 𝜋 (1) , 𝜆 𝜋 (2) , . . . , 𝜆 𝜋 (𝑘 ) ), where 𝜋 is a permutation on [𝑘] such that 𝜆 𝜋 (1) ≤ 𝜆 𝜋 (2) ≤ · · · ≤ 𝜆 𝜋 (𝑘 ) . The following two optimization problems have the same optimal values 𝑂𝑃𝑇1 : min 𝑃,𝑄
1
1
𝑓 (𝑃𝑄) + ∥Λ 2 𝑄∥ 2𝐹 + ∥𝑃Λ 2 ∥ 2𝐹 .
𝑂𝑃𝑇2 : min𝑊 𝑓 (𝑊) + 2∥𝑊 ∥ 𝜔,∗ subject to rank(𝑊) ≤ 𝑘. In addition, if (𝑃∗ , 𝑄 ∗ ) is an optimal solution for 𝑂𝑃𝑇1, then 𝑊 ∗ = 𝑃∗ 𝑄 ∗ is an optimal solution for 𝑂𝑃𝑇2. If 𝑊 ∗ is an optimal solution for 𝑂𝑃𝑇2, then there exists an optimal solution (𝑃∗ , 𝑄 ∗ ) for 𝑂𝑃𝑇1 such that 𝑊 ∗ = 𝑃∗ 𝑄 ∗ . We reiterate three points made earlier (Sec. 2). (i) Both A3 and A4 effectively add a nuclear norm regularizer. (ii) Diagonals of Λ do not need to be sorted in ascending order as stated in [Bao et al., 2020] because any permutation of the diagonals will be equivalent to 𝑂𝑃𝑇2. This also limits the search space and affects a model’s prediction power. (iii) Prop. 2 “compiles” a non-differentiable objective (𝑂𝑃𝑇2) into an equivalent differentiable one (𝑂𝑃𝑇1), which is easier to optimize. In addition, 𝑓 (·) in A3 & A4 is the reconstruction error, in which case closed form solutions exist. We next explain the intuition for proving Prop. 2 (see Appendix for the full analysis). Consider 𝑂𝑃𝑇1 and let 𝑊 = 𝑃𝑄. Our goal is to characterize the behaviors of 𝑃 and 𝑄 with the presence of the regularizers T . Because two regularizers ∥Λ 12 𝑄∥ 2 and when 𝑊 = 𝑃𝑄 is known (fixed). Let the SVD of 𝑊 be 𝑈𝑊 Σ𝑊 𝑉𝑊 𝐹 1
1
1
T , where Ω is a unitary matrix. 2 2 ∥𝑃Λ 2 ∥ 2𝐹 are symmetric, we could “guess” 𝑃 = 𝑈𝑊 Σ𝑊 Ω and 𝑄 = ΩT Σ𝑊 𝑉𝑊 Now we have
6
1
1
1
1
1
1
1
1
T 2 2 2 ∥Λ 2 𝑄∥ 2𝐹 + ∥𝑃Λ 2 ∥ 2𝐹 = ∥Λ 2 ΩΣ 2 𝑉𝑊 ΩΛ 2 ∥ 2𝐹 = 2∥Λ 2 ΩΣ𝑊 ∥ 2𝐹 . ∥ 𝐹 + ∥𝑈𝑊 Σ𝑊
Now the1 question of finding 𝑃 and 𝑄 when 𝑊 is known boils down to finding a unitary matrix Ω that minimizes 1 2 ∥Λ 2 ΩΣ𝑊 ∥ 2𝐹 , where diagonal matrices Λ and Σ𝑊 are given. Recall that 𝜆𝑖 = Λ𝑖𝑖 and let 𝜎𝑖 = (Σ𝑊 )𝑖𝑖 . Note that 𝜆𝑖 ’s could be unsorted and, and that 𝜎𝑖 ’s are sorted in descending order. If we restrict Ω to be only a permutation matrix, then we aim to find a permutation 𝜋 ∈ [𝑘] that minimizes Í 𝑖 ≤ 𝑘 𝜆 𝜋 (𝑖) 𝜎𝑖 . Using a rearrangement inequality Yue [2020], we can see that the minimal is achieved when 1
1
2 𝜆 𝜋 (1) ≤ 𝜆 𝜋 (2) ≤ · · · ≤ 𝜆 𝜋 (𝑘 ) . In this case, we indeed have minΩ a permutation ∥Λ 2 ΩΣ𝑊 ∥ 2𝐹 = ∥Σ𝑊 ∥ 𝜔,∗ , where 𝜔 = (𝜆 𝜋 (1) , . . . , 𝜆 𝜋 (𝑘 ) ).
Note that because 𝑃𝑄 = 𝑃ΩΩT 𝑄 for any unitary matrix Ω, it is always beneficial to use Ω to shuffle the rows and columns of 𝑃 and 𝑄 so that the largest 𝜎𝑖 is mapped to the smallest 𝜆𝑖 , etc. This “degree of freedom” from Ω also explains why ordering the values along Λ’s diagonal is irrelevant. Appendix shows that even when Ω is allowed to be any unitary matrix, the optimal one is still a permutation matrix. This conclusion can be viewed as a matrix version of re-arrangement inequality. Prop. 2 also leads to the following Corollary. Corollary 1. Let X ∈ R𝑚×𝑛 (𝑚 ≥ 𝑛) be a full rank matrix. Let Λ be a diagonal matrix with Λ𝑖𝑖 ≥ 0 for all 𝑖. Let 𝑃 ∈ R𝑚×𝑘 and 𝑄 ∈ R 𝑘×𝑛 . Consider the optimization problems 1
1
min ∥ 𝑋 − 𝑃𝑄∥ 2𝐹 + ∥Λ 2 𝑄∥ 2𝐹 + ∥𝑃Λ 2 ∥ 2𝐹
(9)
𝑃,𝑄
min ∥ 𝑋 − 𝑋 𝐴𝐵∥ 2𝐹 + ∥Λ𝐵∥ 2𝐹 + ∥ 𝑋 𝐴∥ 2𝐹 .
(10)
𝐴,𝐵
Let the SVD of 𝑋 be 𝑈Σ𝑉 T , and let Σ 𝑘 ∈ R 𝑘×𝑘 be a matrix comprising the 𝑘 largest singular values of 𝑋, and 𝑈 𝑘 and 𝑉𝑘 be the corresponding singular vectors. Let 𝜆 (1) ≥ 𝜆 (2) ≥ · · · ≥ 𝜆 (𝑘 ) be the sorted sequence of the diagonal values from Λ. (9) has a closed-form solution: √︃ √︃ 𝑃∗ = 𝑈 𝑘 𝑑𝑖𝑎𝑔( (𝜎1 − 𝜆 (𝑘 ) ) + ), . . . , (𝜎𝑘 − 𝜆 (1) ) + )Ω, √︃ 𝑄 ∗ = ΩT 𝑑𝑖𝑎𝑔( (𝜎1 − 𝜆 (𝑘 ) ) + ), . . . )𝑉𝑘T ,
(11)
where Ω is a unitary matrix that corresponds to the permutation 𝜋 such that 𝜆 𝜋 (1) ≤ · · · ≤ 𝜆 𝜋 (𝑘 ) . In addition, 1 1 (10) has a closed-form solution: 𝐴∗ = 𝑋 † 𝑃∗ Λ 2 and 𝐵∗ = Λ− 2 𝑄 ∗ , where 𝑋 † is the pseudo-inverse of 𝑋. Uniform solution structure for regularizing ∥𝑊 ∥ ∗𝑝 . In [Cavazza et al., 2018], it was observed that when we use a neural net 𝑋 = 𝑃𝑄 T (with learnable parameters being 𝑃 and 𝑄) to train a model and a standard dropout is used, the objective is equivalent to solving min𝑊 ∥ 𝑋 − 𝑊 ∥ 2𝐹 + 𝜆∥𝑊 ∥ 2∗ . While the regularizer ∥𝑊 ∥ 2∗ deviates from the standard one ∥𝑊 ∥ ∗ , the optimal solution here is 𝑊 = 𝑈𝑆 𝜇 (Σ)𝑉 T , where 𝑈Σ𝑉 T is SVD of 𝑋, and 𝑆 𝜇 (Σ) is a diagonal matrix such that its (𝑖, 𝑖)-th element is (Σ𝑖,𝑖 − 𝜇) + , in which 𝜇 depends on the data 𝑋 and 𝜆. In other words, the optimal solution for regularizers ∥𝑊 ∥ 2∗ and ∥𝑊 ∥ ∗ are strikingly similar. Thus, we are interested in how regularizers with different exponents are connected. Our main observation is that for any regularizer ∥𝑊 ∥ ∗𝑝 (𝑝 ≥ 1), the optimal solution has the same structure. Lemma 2. Let 𝑋 ∈ R𝑚×𝑛 , where 𝑚 ≥ 𝑛. Let the 𝑑 leading SVDs of 𝑋 be 𝑈𝑑 and 𝑉𝑑 respectively. Let 𝜎1 , . . . , 𝜎𝑛 be the singular values of 𝑋. Consider the optimization problem: 1 min ∥ 𝑋 − 𝑊 ∥ 2𝐹 + 𝜆∥𝑊 ∥ ∗𝑝 . 𝑊 2 7
(12)
Í Let 𝜇 𝑘 = 𝑘1 𝑖 ≤ 𝑘 𝜎𝑖 (𝑋), 𝜂(𝜇 𝑘 ) be the positive root of the function 𝑧 + 𝜆𝑘 𝑧 𝑝−1 − 𝑘 𝜇 𝑘 and 𝑑 be the largest value such that 𝜎𝑑 (𝑋) − 𝜆(𝜂(𝜇 𝑑 )) 𝑝−1 ≥ 0. Let 𝜇 = 𝜆(𝜂(𝜇 𝑑 )) 𝑝−1 . The optimal solution of 𝑊 is 𝑈𝑑 diag(𝜎1 − 𝜇, 𝜎2 − 𝜇, . . . , 𝜎𝑑 − 𝜇)𝑉𝑑T . Lemma 2 is a straightforward generalization of Prop. 12 in [Cavazza et al., 2018] so the contribution here is a conceptual one: it implies that regularizers ∥𝑊 ∥ ∗𝑝 with different 𝑝’s differ in how the shrinkage variable 𝜇 is obtained. Observe that 𝜇 is an important hyper-parameter that needs to be extensively tuned against data, all regularizers ∥𝑊 ∥ ∗𝑝 provide the same learning power. As noted earlier, 𝜇 is a function of 𝑋 (i.e., different 𝑋 will result in different 𝜇) unless 𝑝 = 1. In addition, 𝜆 needs to be rescaled when 𝑋 is scaled by a constant factor unless 𝑝 = 2. This implies it could be easier to tune 𝜇 when 𝑝 = 1, 2, and explains why only 𝑝 = 1, 2 have been extensively considered.
4
Low-Rank Frobenius norm based regularizations
This section presents the close-form low-rank estimators with comparable performance to the state of the art algorithms. Approximate Low rank DLAE and EDLAE. Recall that DLAE solves 1
min ||𝑋 − 𝑋𝑊 || 2𝐹 + ||Λ 2 𝑊 || 2𝐹 𝑊
𝑠.𝑡.
Λ=
𝑝 𝑑𝑀𝑎𝑡 (𝑑𝑖𝑎𝑔(𝑋 𝑇 𝑋)), 1− 𝑝
(13)
and EDLAE with the additional 𝑑𝑖𝑎𝑔(𝑊) = 0 constraint. Even though both the Nuclear norm based regularization as well as the full rank DLAE and EDLAE solutions all have closed-form solutions, such solution is unknown for the low-rank DLAE and EDLAE, whose existing solution is based on ADMMSteck [2020]. The closed form solutions will help both better understand and compare these models, and determine the hyper-parameters, which is usually difficult for the ADMM type solutions. For DLAE, since it can be considered a special form of Tikhonov regularization (Eq 6) which has a closed form Jin et al. [2021], its closed form solution, referred to LR-DLAE, is immediately available (See Line 6 in Table 3). However, for EDLAE, it has the zero diagonal constraint, which make the exact solution difficult to express. Here, we present an approximate low-rank closed-form solution of equation (13) by decomposing the optimization problem into two subproblems, which is similar to [Jin et al., 2021]: We first consider the full-rank closed form solution for EDLAE Steck [2020], which is: 𝑊 ∗ = 𝐼 − 𝐶 · 𝑑𝑀𝑎𝑡 (1 ⊘ 𝑑𝑖𝑎𝑔(𝐶)), where 𝐶 = (𝑋 𝑇 𝑋 + Λ) −1 Then, we consider two approaches to produce low-rank matrix approximation of 𝑊 ∗ : b to best approximate the performance of 𝑊 ∗ without the zero (Method 1 (LR-EDLAE-1): ) Selecting 𝑊 diagonal constraint: b = arg 𝑊 = arg
min
||𝑋𝑊 ∗ − 𝑋𝑊 || 2𝐹
min
||𝑋𝑊 ∗ − 𝑋𝑊 || 2𝐹 + ||Λ 2 (𝑊 ∗ − 𝑊)|| 2𝐹
𝑟 𝑎𝑛𝑘 (𝑊 ) ≤ 𝑘 𝑟 𝑎𝑛𝑘 (𝑊 ) ≤ 𝑘
1
(14)
𝑋 where 𝑋 = 1 . Noting, in the full rank problem equation (13), it forces the diagonal of derived matrix Λ2 (𝑊 ∗ ) to be zero. Here, we relax the zero diagonal constraint - the diagonal of low rank approximate matrix 8
b doesn’t have to be zero, which has also been discussed in [Steck, 2020]. The closed-from solution of 𝑊 equation (14) is given by: b = 𝑊 ∗ (𝑄 𝑘 𝑄 𝑇 ) 𝑊 (15) 𝑘
where 𝑄 𝑘 comes from SVD: ∗
∗
𝑌 = 𝑋𝑊 ∗ , and 𝑌 (𝑘) = 𝑃 𝑘 Σ 𝑘 𝑄 𝑇𝑘 (Method 2 (LR-EDLAE-2):) SVD approximation of 𝑊 ∗ : The alternative solution is to simply perform SVD, and which gives a low-rank estimation of 𝑊 ∗ . b is relaxed. In the next section, the Note that in both approaches, the hard constraint of zero-diagonal on 𝑊 experimental results show both approaches can provide comparable or better performance compared with the ADMM solution, and also very close to the full rank EDLAE solution.
5
Experimental Results
In this section, we experimentally study different regularizations for linear recommendation models. Our goal is to validate the effectiveness of various regularizations (all can be categorized under nuclear norm and Frobenius norm) together with their simple closed-form solutions. We aim to answer three questions: Q1. How does the closed form solution of low rank Frobenius norm perform compared with the ADMM solutions (Section 4) and how does the weighted nuclear norm regularizer for matrix factorization (Proposition 2 and Corollary 1) perform? Q2. What is the tradeoff between the number of factors (rank 𝑘) and the recommendation accuracy? Q3. How does the ordering of weights from small to large (non-descending) for adjusting the singular values (Corollary 1) affect the recommendation performance? Experimental Setup: We use three commonly used datasets for recommendation studies: MovieLens 20 Million (ML-20M) [Harper and Konstan, 2015], Netflix Prize (Netflix) [Bennett et al., 2007], and the Million Song Data (MSD)[Bertin-Mahieux et al., 2011]. We obtained these datasets and all benchmarks from authors of EASE [Steck, 2019], EDLAE [Steck, 2020], and Mult-VAE [Liang et al., 2018]. Similar to the latest study in EASE [Steck, 2019], and EDLAE [Steck, 2020], we consider the following state-of-the-art recommendation models: ALS (WMF) [Hu et al., 2008] for matrix factorization approaches, SLIM [Ning and Karypis, 2011], EASE [Steck, 2019], and EDLAE [Steck, 2020] for linear autoendoers, CDAE [Wu et al., 2016], Mult-DAE and Mult-VAE [Liang et al., 2018] for deep learning models. The experiment settings for these baseline are the same as [Liang et al., 2018, Steck, 2019, 2020]. Also we follow their practice [Liang et al., 2018, Steck, 2019, 2020] for the strong generalization by splitting the users into training, validation and tests group, and report performance metrics 𝑅𝑒𝑐𝑎𝑙𝑙@20, 𝑅𝑒𝑐𝑎𝑙𝑙@50 and 𝑛𝐷𝐶𝐺@100. Finally, we note that our code are openly available (see Appendix). Q1: Low-Rank Frobenius Norm and (Weighted) Nuclear Norm Regularization: In this experiment, we evaluate the low rank Frobenius norm regularization and the nuclear norm regularization (equation (1)) for the matrix factorization. Here EDLAE-ADMM, LRR, LR-DLAE, LR-EDLAE-1, LR-EDLAE-2, MF dropout and LVAE are listed in table 3 in Appendix. To determine the non-descending order of weights 𝜆 𝑖 for the closed-form solution in (equation (11)), we follow the practice in weighted nuclear norm regularization in [Gu et al., 2014] as well as the optimized pPCA weight [Lucas et al., 2019b]. Let 𝜆𝑖 = 𝜎𝐶𝑖 where 𝐶 is a hyperparameter, and we perform grid-search to find the optimal one. In Table 1, we can see that the weighted nuclear norm regularization (LVAE) based matrix factorization actually performs worse than the constant weighted version (Regularized PCA). And the latter shows very 9
ML 20M
Netflix
0.400
MSD
0.40
0.424 0.422
0.38
0.395
0.36
0.418 0.416
EDLAE-ADMM LR-EDLAE-2 LR-DLAE LR-EDLAE-1
0.414 0.412 0.410
2000
4000
k
6000
8000
0.390
nDCG@100
nDCG@100
nDCG@100
0.420
0.385
EDLAE-ADMM LR-EDLAE-2 LR-DLAE LR-EDLAE-1
0.380 0.375
2000
4000
(a)
k
6000
8000
0.34 0.32
EDLAE-ADMM LR-EDLAE-2 LR-DLAE LR-EDLAE-1
0.30 0.28 2000
(b)
4000
6000
k
8000
10000 12000
(c)
Figure 1: Low rank models 𝑛𝐷𝐶𝐺@100 on test data for 3 datasets. strong performance comparing against the WFM/ALS (one of the most popular implicit matrix factorization algorithm). We also observe the closed form solutions (LR-DLAE, LR-EDLAE-1 and LR-EDLAE-2) all perform very comparable with the ADMM based low rank solution and the full rank DLAE and EDLAE solutions. Q2: nDCG vs Rank 𝑘 for low-rank Frobenius norm: In this experiment, we focus on evaluating the recommendation accuracy (using nDCG) against the rank 𝑘. Specifically, we vary the rank 𝑘 from around 1𝐾 to around 10𝐾, and we tune and compare four different methods, including EDLAE-ADMM, LR-DLAE, LREDLAE-1, and LR-EDLAE-2. We have the following observations: 1) As 𝑘 increases, the recommendation accuracy also increases in general; however, most of them reaches a plateau around similar 𝐾, and for different datasets, the saturating point varies. 2) LR-DLAE performs worse than EDLAE based approaches in two out of three datasets. This partially demonstrates the benefits of zero-diagonal constraint. 3) The closed form solution of EDLAE performs comparable or even slightly better than ADMM methods as 𝐾 grows; but when 𝐾 is relatively small, ADMM method perform slightly better. But none-the-less, for most of the reasonable choices of 𝑘 when low-rank approximates full rank, the closed form solution performs comparable or better.
Table 1: The performance comparison between different regularizations. For notation, please refer table 3 for more details. Model EASE DLAE EDLAE EDLAE-ADMM Frobinius Norm LRR LR DLAE LR-EDLAE-1 LR-EDLAE-2 MF dropout Nuclear Norm Regularized PCA LVAE WMF/ALS SLIM Baseline CDAE MULT-DAE MULT-VAE # items # users # interactions
Recall@20 0.391 0.392 0.393 0.392 0.376 0.392 0.392 0.392 0.367 0.364 0.348 0.360 0.370 0.391 0.387 0.395
ML-20M Recall@50 0.521 0.527 0.523 0.524 0.511 0.527 0.523 0.523 0.501 0.501 0.474 0.498 0.495 0.523 0.524 0.537 20108 136677 10mil
nDCG@100 0.420 0.424 0.424 0.424 0.408 0.424 0.424 0.424 0.393 0.392 0.378 0.386 0.401 0.418 0.419 0.426
Recall@20 0.362 0.362 0.366 0.365 0.348 0.362 0.365 0.365 0.334 0.331 0.325 0.316 0.347 0.343 0.344 0.351
10
Netflix Recall@50 0.445 0.446 0.449 0.448 0.431 0.445 0.449 0.449 0.418 0.417 0.405 0.404 0.428 0.428 0.438 0.444 17769 463435 57mil
nDCG@100 0.393 0.395 0.398 0.396 0.380 0.395 0.398 0.398 0.365 0.365 0.357 0.351 0.379 0.376 0.380 0.386
MSD Recall@20 Recall@50 nDCG@100 0.333 0.428 0.389 0.329 0.426 0.387 0.334 0.429 0.392 0.330 0.424 0.386 0.248 0.335 0.301 0.306 0.403 0.363 0.327 0.423 0.384 0.325 0.421 0.382 0.270 0.367 0.328 0.229 0.313 0.279 0.205 0.254 0.286 0.211 0.312 0.257 no results in [Ning and Karypis, 2011] 0.188 0.283 0.237 0.266 0.363 0.313 0.266 0.364 0.316 41140 571353 34mil
Table 2: Investigating the weight ordering of Matrix Factorization Model MF/LRR weighted MF sorted
Recall@20 0.3806 0.3017
ML-20M Recall@50 0.5175 0.4507
nDCG@100 0.4102 0.3361
Recall@20 0.3484 0.2860
Netflix Recall@50 0.4320 0.3801
nDCG@100 0.3797 0.3265
Recall@20 0.2508 0.2288
MSD Recall@50 0.3390 0.3148
nDCG@100 0.3037 0.2802
Q3: Impact of Weight Ordering: Finally, we study how the ordering of weights from small to large (non-descending) for adjusting the singular values (Corollary 1) affects the recommendation performance using matrix factorization (closed-form solution in Eq. 11). Our results are in Table 2. Here, we obtain the searched optimal weight parameters from weighted Tikhonov regularization (following the approach in Jin et al. [2021]), and map it back to the parameters in the closed form solution (Proposition 1). Then we sort the parameters in the non-descending order, and then report their results in the second row of Table 2. We can see that the recommendation performance becomes significant worst. This help confirm our conjecture that the strict ordering of weight on matrix factorization and other regularizations can be an inherent limitation for those approaches.
6
Conclusion and Discussion
This work provides a complete analysis on the recently proposed linear models for recommendation systems. Despite that models leverage different deep learning techniques, they achieve similar performance. We find that this is not coincident: all the models add either a nuclear-norm-based (Lemma 1) or a Frobenius-norm based regularizer. The nuclear-norm-based approach results in estimators that keep 𝑋’s singular vectors and shrink its singular values in a quite rigid way (Proposition 2 and Lemma 2), which limit their prediction power. The Frobenius-norm models are more express (Proposition 1) and effective but their estimators are either full-rank or do not have closed form solutions. To get the best of both nuclear and Frobenius worlds, we propose two low-rank and closed-form estimators (Sec. 4) based on carefully generalizing Frobenius-norm based regularizers. These estimators have competitive performance against linear performance leaders, and thus concisely pack all the benefits obtained by a recent long line of research and abstract out all the computation nuance.
11
References Charu C. Aggarwal. Recommender Systems: The Textbook. Springer, 1st edition, 2016. ISBN 3319296574. Xuchan Bao, James Lucas, Sushant Sachdeva, and Roger B. Grosse. Regularized linear autoencoders recover the principal components, eventually. In Hugo Larochelle, Marc’Aurelio Ranzato, Raia Hadsell, Maria-Florina Balcan, and Hsuan-Tien Lin, editors, Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurIPS 2020, December 6-12, 2020, virtual, 2020. James Bennett, Charles Elkan, Bing Liu, Padhraic Smyth, and Domonkos Tikk. Kdd cup and workshop 2007. 2007. doi: 10.1145/1345448.1345459. URL https://doi.org/10.1145/1345448.1345459. Thierry Bertin-Mahieux, Daniel PW Ellis, Brian Whitman, and Paul Lamere. The million song dataset. 2011. Srinadh Bhojanapalli, Anastasios Kyrillidis, and Sujay Sanghavi. Dropping convexity for faster semi-definite optimization. In Conference on Learning Theory, pages 530–582. PMLR, 2016. Emmanuel J Candès and Terence Tao. The power of convex relaxation: Near-optimal matrix completion. IEEE Transactions on Information Theory, 56(5):2053–2080, 2010. Jacopo Cavazza, Pietro Morerio, Benjamin Haeffele, Connor Lane, Vittorio Murino, and Rene Vidal. Dropout as a low-rank regularizer for matrix factorization. In Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics. PMLR, 2018. Kun Chen, Hongbo Dong, and Kung-Sik Chan. Reduced rank regression via adaptive nuclear norm penalization. Biometrika, 100(4):901–920, 2013. Evangelia Christakopoulou and George Karypis. Hoslim: Higher-order sparse linear method for top-n recommender systems. In Advances in Knowledge Discovery and Data Mining, 2014. Maurizio Ferrari Dacrema, P. Cremonesi, and D. Jannach. Are we really making much progress? a worrying analysis of recent neural recommendation approaches. In RecSys’19, 2019a. Maurizio Ferrari Dacrema, Paolo Cremonesi, and Dietmar Jannach. Are we really making much progress? RecSys’19, 2019b. Mukund Deshpande and George Karypis. Item-based top-n recommendation algorithms. ACM Trans. Inf. Syst., 2004. Shuhang Gu, Lei Zhang, Wangmeng Zuo, and Xiangchu Feng. Weighted nuclear norm minimization with application to image denoising. In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 2862–2869, 2014. F. Maxwell Harper and Joseph A. Konstan. The movielens datasets: History and context. ACM Trans. Interact. Intell. Syst., 2015. doi: 10.1145/2827872. Y. Hu, Y. Koren, and C. Volinsky. Collaborative filtering for implicit feedback datasets. In ICDM’08, 2008. José Pedro Iglesias, Carl Olsson, and Marcus Valtonen Örnhag. Accurate optimization of weighted nuclear norm for non-rigid structure from motion. arXiv preprint arXiv:2003.10281, 2020. Ruoming Jin, Dong Li, Jing Gao, Zhi Liu, Li Chen, and Yang Zhou. Towards a better understanding of linear recommendation models. In KDD’21, 2021. URL https://arxiv.org/abs/2105.12937. Santosh Kabbur, Xia Ning, and George Karypis. Fism: Factored item similarity models for top-n recommender systems. KDD ’13, 2013. 12
Diederik P. Kingma and Max Welling. Auto-Encoding Variational Bayes. In 2nd International Conference on Learning Representations, ICLR 2014, Banff, AB, Canada, April 14-16, 2014, Conference Track Proceedings, 2014. Yehuda Koren. Factorization meets the neighborhood: A multifaceted collaborative filtering model. In KDD’08, 2008. Yehuda Koren, Robert Bell, and Chris Volinsky. Matrix factorization techniques for recommender systems. Computer, 42(8):30–37, August 2009. Daniel Kunin, Jonathan Bloom, Aleksandrina Goeva, and Cotton Seed. Loss landscapes of regularized linear autoencoders. In Proceedings of the 36th International Conference on Machine Learning, Proceedings of Machine Learning Research, pages 3560–3569. PMLR, 09–15 Jun 2019. Xiaopeng Li and James She. Collaborative variational autoencoder for recommender systems. In Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’17, page 305–314, 2017. Dawen Liang, R. G. Krishnan, M. D. Hoffman, and T. Jebara. Variational autoencoders for collaborative filtering. In WWW’18, 2018. James Lucas, George Tucker, Roger B Grosse, and Mohammad Norouzi. Don't blame the elbo! a linear vae perspective on posterior collapse. In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019a. James Lucas, George Tucker, Roger B. Grosse, and Mohammad Norouzi. Don’t blame the elbo! A linear VAE perspective on posterior collapse. In Hanna M. Wallach, Hugo Larochelle, Alina Beygelzimer, Florence d’Alché-Buc, Emily B. Fox, and Roman Garnett, editors, Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada, pages 9403–9413, 2019b. URL https://proceedings.neurips. cc/paper/2019/hash/7e3315fe390974fcf25e44a9445bd821-Abstract.html. Sahand Negahban and Martin J Wainwright. Estimation of (near) low-rank matrices with noise and high-dimensional scaling. The Annals of Statistics, pages 1069–1097, 2011. Xia Ning and George Karypis. Slim: Sparse linear methods for top-n recommender systems. ICDM ’11, 2011. Steffen Rendle, Li Zhang, and Yehuda Koren. On the difficulty of evaluating baselines: A study on recommender systems, 2019. Suvash Sedhain, Aditya Krishna Menon, Scott Sanner, and Darius Braziunas. On the effectiveness of linear models for one-class collaborative filtering. AAAI’16, 2016. Ilya Shenbin, Anton Alekseev, Elena Tutubalina, Valentin Malykh, and Sergey I. Nikolenko. Recvae: A new variational autoencoder for top-n recommendations with implicit feedback. In Proceedings of the 13th International Conference on Web Search and Data Mining, WSDM ’20, page 528–536, 2020. Harald Steck. Embarrassingly shallow autoencoders for sparse data. WWW’19, 2019. Harald Steck. Autoencoders that don’t overfit towards the identity. In NIPS, 2020. Robert Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society: Series B (Methodological), 58(1):267–288, 1996.
13
Michael E. Tipping and Chris M. Bishop. Probabilistic principal component analysis. JOURNAL OF THE ROYAL STATISTICAL SOCIETY, SERIES B, 61(3):611–622, 1999. Madeleine Udell, Corinne Horn, Reza Zadeh, and Stephen Boyd. Generalized low rank models. Found. Trends Mach. Learn., 9(1):1–118, June 2016. Romain Warlop, Alessandro Lazaric, and Jérémie Mary. Parallel higher order alternating least square for tensor recommender system. In Workshops at the Thirty-First AAAI Conference on Artificial Intelligence, 2017. Yao Wu, Christopher DuBois, Alice X. Zheng, and Martin Ester. Collaborative denoising auto-encoders for top-n recommender systems. WSDM ’16, 2016. Man-Chung Yue. A matrix generalization of the hardy-littlewood-p\’olya rearrangement inequality and its applications. arXiv preprint arXiv:2006.08144, 2020. Shuai Zhang, Lina Yao, Aixin Sun, and Yi Tay. Deep learning based recommender system: A survey and new perspectives. ACM Comput. Surv., 2019. Shuai Zheng, Chris Ding, and Feiping Nie. Regularized singular value decomposition and application to recommender system, 2018.
14
A
Related Work
There have been extensive researches on recommendation [Aggarwal, 2016]. Besides the basic user-based and item-based collaborative filtering [Deshpande and Karypis, 2004], the full rank linear autoencoder approaches include SLIM [Ning and Karypis, 2011], HOLISM [Christakopoulou and Karypis, 2014], EASE [Steck, 2019], DLAE (Denoising linear autoencoder) Steck [2020], whereas low-rank approaches include [Kabbur et al., 2013, Sedhain et al., 2016, Steck, 2020]. All the customized recommendation has been enforcing zero diagonal constraints for generalization purpose, whereas we show an approximate closed-form solution for a two-term Tikhonov regularization without the zero diagonal constraint can be as effective as these models. Matrix factorization has been been widely studied in practice, partially due to Netflix competition Koren et al. [2009]. Methods like SVD++ Koren [2008] and implicit Alternating Least Square (ALS) method Hu et al. [2008] (also weighted matrix factorization) have been very influential. [Jin et al., 2021] shows the relationship between linear autoencoders and matrix factorization, and pointed out a potential advantage of linear autoencoders. In this work, we take a step further to reveal a deeper relationship between Tikhonov regularized linear autoencoders and a few other regularizations including matrix factorization, and show the potential limitation of the class of regularization. We also utilize the linear variational autoencoders (LVAE) to study how the deep VAE based recommendation approaches [Li and She, 2017, Liang et al., 2018, Shenbin et al., 2020] relate to linear autoencoders and matrix factorization. Outside recommendation, there have been a few recent studies on regularization landscapes of linear (variational) autoencoders [Kunin et al., 2019, Bao et al., 2020, Lucas et al., 2019a]. They do not provide the general weighted ℓ2 regularization and thus did not find the inherent limitation on the regularization (for MF). Our LVAE inspired regularization is also never studied before. Nuclear norm regularizers can recover low-rank matrices in the vector regression setting [Negahban and Wainwright, 2011]. Its weighted generalization can be applied in the area of image processing [Gu et al., 2014]. Because weighted nuclear-norm is usually not convex or differentiable, finding optimal solutions is difficult except for a few special cases [Chen et al., 2013].
B
Proofs
B.1
Proof of Lemma 1
The Linear Variational AutoEncoder (LVAE) is defined in the same way as [Lucas et al., 2019b]:
𝑝(𝑥 | 𝑧) = N 𝑊 𝑧 + µ, 𝜎 2 𝐼
(16)
𝑞(𝑧 | 𝑥) = N (𝑉 (𝑥 − µ), 𝐷) For simplification, we set 𝜇 = 0 in following context. And the ELBO of LVAE is known as: L 𝑥 = −𝐾 𝐿 (𝑞(𝑧|𝑥)|| 𝑝(𝑧)) + E𝑞 (𝑧 | 𝑥 ) [log 𝑝(𝑥|𝑧)] 𝐾 𝐿(𝑞(𝑧|𝑥)|| 𝑝(𝑧)) = − log |𝐷 | + 𝑥 𝑇 𝑉 𝑇 𝑉𝑥 + 𝑡𝑟 (𝐷) − 𝑘 𝑛 1 E𝑞 (𝑧 | 𝑥 ) [log 𝑝(𝑥|𝑧)] = − 2 𝑡𝑟 (𝑊 𝐷𝑊 𝑇 ) + 𝑥 𝑇 𝑉 𝑇 𝑊 𝑇 𝑊𝑉𝑥 − 2𝑥 𝑇 𝑊𝑉𝑥 + 𝑥 𝑇 𝑥 − log 2𝜋𝜎 2 2 2𝜎 Again, the (maximizing) ELBO can be written as: 15
(17)
𝑛 1 − log |𝐷 | + 𝑥 𝑇 𝑉 𝑇 𝑉𝑥 + 𝑡𝑟 (𝐷) − 𝑘 − log 2𝜋𝜎 2 2 2 1 𝑇 𝑇 𝑇 𝑇 𝑇 𝑇 − 𝑡𝑟 (𝑊 𝐷𝑊 ) + 𝑥 𝑉 𝑊 𝑊𝑉𝑥 − 2𝑥 𝑊𝑉𝑥 + 𝑥 𝑥 2𝜎 2 1 1 √ 2 2 = − ||𝑉𝑥|| 22 − ||𝑊 𝐷 || + ||𝑥 − 𝑊𝑉𝑥|| 𝐹 2 + 𝑓 (𝐷, 𝜎) 2 2𝜎 2
L𝑥 = −
(18)
where 𝑓 (𝐷, 𝜎) = 21 log |𝐷 | − 12 𝑡𝑟 (𝐷) + 𝑘2 − 𝑛2 log 2𝜋𝜎 2 , 𝑥 ∈ R𝑛 and 𝑧 ∈ R 𝑘 . For whole data, it is equivalent to minimize: √ L = ||𝑋 − 𝑊𝑉 𝑋 || 2𝐹 + 𝑁 ||𝑊 𝐷 || 2𝐹 + 𝜎 2 ||𝑉 𝑋 || 2𝐹 + 𝑔(𝐷, 𝜎) √ = ||𝑋 𝑇 − 𝑋 𝑇 𝑉 𝑇 𝑊 𝑇 || 2𝐹 + 𝑁 || 𝐷𝑊 𝑇 || 2𝐹 + 𝜎 2 ||𝑋 𝑇 𝑉 𝑇 || 2𝐹 + 𝑔(𝐷, 𝜎)
(19)
where 𝑔(𝐷, 𝜎) = −𝜎 2 𝑁 log |𝐷 | − 𝑡𝑟 (𝐷) + 𝑘 − 𝑛 log 2𝜋𝜎 2 ).
B.2
Proof of Proposition 1
Note that when 𝜎𝑖 ≤ 𝜆 (𝑘−𝑖) , the new singular value shrinks to zero, and can be removed. Basically, for any 𝜆 (1) ≥ · · · ≥ 𝜆 (𝑘 ) , we can build the corresponding Tikhonov regularized instance by setting 𝜎𝑖2
= ′
𝜎𝑖2 + 𝜆 𝑖
𝜎𝑖3 𝜎𝑖 − 𝜆 (𝑘−𝑖) , i.e.,𝜆′𝑘 = − 𝜎𝑖2 . 𝜎𝑖 𝜎𝑖 − 𝜆 (𝑘−𝑖)
(20)
Discussion of Proposition 1: Further, the same observation holds true for the regularization (2), and the weighted-nuclear norm regularization in when the weights are in the non-ascending order. This observation suggests a potentially limitation of the earlier regularization as they will always try to maintain the larger singular values: when a singular value is large, the shrinkage will be small.Such regularization has shown to work well in the areas such as image processing Gu et al. [2014]. But it has not been studied or confirmed if it will work for the recommendation. In Section 5, we report our experimental study which shows such regularization could be too restrictive for recommendation.
B.3
Proof of Proposition 2
By slightly abusing the notation, we shall let 𝑂𝑃𝑇1 (𝑂𝑃𝑇2) be the value of the optimal solution for 𝑂𝑃𝑇1 (𝑂𝑃𝑇2). We need to show that 𝑂𝑃𝑇1 = 𝑂𝑃𝑇2. We need two directions. 𝑂𝑃𝑇2 ≥ 𝑂𝑃𝑇1: Let 𝑊 ∗ be an optimal solution for 𝑂𝑃𝑇2. Let the SVD of 𝑊 ∗ be 𝑈 ∗ Σ∗ (𝑉 ∗ ) T . Recall that 𝑊 ∗ needs to satisfy the rank constraint rank(𝑊 ∗ ) ≤ 𝑘 so 𝑈 ∗ ∈ R𝑚×𝑘 , Σ∗ ∈ R 𝑘×𝑘 , and 𝑉 ∗ ∈ R𝑛×𝑘 . Let 𝜋 be a permutation on [𝑘] such that 𝜆 𝜋 (1) ≤ 𝜆 𝜋 (2) ≤ · · · ≤ 𝜆 𝜋 (𝑘 ) . Let also Ω be the corresponding permutation matrix. Specifically, Ω ∈ {0, 1} 𝑘×𝑘 and there is exactly one entry in each row of Ω is 1: 1 if 𝑗 = 𝜋(𝑖). Ω𝑖, 𝑗 = 0 otherwise. For example, consider a case in which 𝜆1 > 𝜆 2 > · · · > 𝜆 𝑘 . Then we set 𝜋 = (𝑘, 𝑘 − 1, . . . , 1), and correspondingly,
16
0 © 0 Ω= « 1 1
... 0 1 ª ... 1 0 ® ®. ... ® ... 0 0 ¬
1
Next, let 𝑃 = 𝑈 ∗ (Σ∗ ) 2 Ω and 𝑄 = ΩT (Σ∗ ) 2 (𝑉 ∗ ) T . We have 𝑊 ∗ = 𝑃𝑄 and 𝑓 (𝑊 ∗ ) = 𝑓 (𝑃𝑄). In addition, ∑︁ 1 1 1 1 𝜆 𝜋 (𝑖) 𝜎𝑖 = 2∥𝑊 ∗ ∥ 𝜔,∗ , ∥Λ 2 𝑄∥ 2𝐹 + ∥𝑃Λ 2 ∥ 2𝐹 = 2∥Λ 2 Ω(Σ∗ ) 2 ∥ 2𝐹 = 2 𝑖≤𝑘
where 𝜎𝑖 is the 𝑖-th largest singular value of 𝑊 ∗ . In other words, we have found a (𝑃, 𝑄) pair such that 1
1
𝑓 (𝑃𝑄) + ∥Λ 2 𝑄∥ 2𝐹 + ∥𝑃Λ 2 ∥ 2𝐹 = 𝑓 (𝑊 ∗ ) + 2∥𝑊 ∗ ∥ 𝜔,∗ = 𝑂𝑃𝑇2, which shows that 𝑂𝑃𝑇1 ≤ 𝑂𝑃𝑇2. 𝑂𝑃𝑇2 ≤ 𝑂𝑃𝑇1. Let 𝑃∗ and 𝑄 ∗ be an optimal solution for 𝑂𝑃𝑇1. Let the singular values of 𝑃∗ be 𝜎1 (𝑃∗ ) ≥ 𝜎2 (𝑃∗ ) ≥ · · · ≥ 𝜎𝑘 (𝑃∗ ) and those of 𝑄 ∗ be 𝜎1 (𝑄 ∗ ) ≥ 𝜎2 (𝑄 ∗ ) ≥ · · · ≥ 𝜎𝑘 (𝑄 ∗ ). Let also 𝜎1∗ ≥ · · · ≥ 𝜎𝑘∗ be the singular values of 𝑃∗ 𝑄 ∗ . 1
1
We shall find a lower bound of ∥Λ 2 𝑄∥ 2𝐹 + ∥𝑃Λ 2 ∥ 2𝐹 expressed in terms of 𝜎𝑖∗ ’s. In fact, we shall show that 1
1
∥Λ 2 𝑄∥ 2𝐹 + ∥𝑃Λ 2 ∥ 2𝐹 ≥ 2∥𝑃∗ 𝑄 ∗ ∥ 𝜔,∗ .
(21)
One can see that if Eq. 21 were true, we have 1
1
𝑂𝑃𝑇2 ≤ 𝑓 (𝑃∗ 𝑄 ∗ ) + 2∥𝑃∗ 𝑄 ∗ ∥ 𝜔,∗ ≤ 𝑓 (𝑃∗ 𝑄 ∗ ) + ∥Λ 2 𝑄 ∗ ∥ 2𝐹 + ∥𝑃∗ Λ 2 ∥ 2𝐹 = 𝑂𝑃𝑇1. Thus, it remains to prove Eq. 21. Let 𝜆 (1) ≥ 𝜆 (2) ≥ · · · ≥ 𝜆 (𝑘 ) be a sorted sequence of 𝜆𝑖 ’s i.e., 𝜆 (𝑘 ) = 𝜆 𝜋 (1) , 𝜆 (𝑘−1) = 𝜆 𝜋 (2) , . . . , 𝜆 (1) = 𝜆 𝜋 (𝑘 ) . Í𝑘 Í𝑘 1 1 First, we show that ∥𝑃Λ 2 ∥ 2𝐹 ≥ 𝑖=1 𝜆 (𝑘−𝑖+1) × 𝜎𝑖2 (𝑃∗ ) and ∥Λ 2 𝑄∥ 2𝐹 ≥ 𝑖=1 𝜆 (𝑘−𝑖+1) × 𝜎𝑖2 (𝑄 ∗ ). We need the following Lemma (see e.g., Theorem 2 in Yue [2020]): Lemma 3. Let 𝐴 and 𝐵 be two positive definite matrices in R 𝑘×𝑘 . Then it holds that 𝑘 ∑︁
1
1
𝜎𝑖 ( 𝐴)𝜎𝑘−𝑖+1 (𝐵) ≤ tr(𝐵 2 𝐴𝐵 2 ).
(22)
𝑖=1
Let the SVD of 𝑃∗ be 𝑈 𝑃∗ Σ 𝑃∗ 𝑉𝑃T∗ and that of 𝑄 ∗ be 𝑈𝑄∗ Σ𝑄∗ 𝑉𝑄T ∗ . We have 1
1
1
∥𝑃∗ Λ 2 ∥ 2𝐹 = ∥𝑈 𝑃∗ Σ 𝑃∗ 𝑉𝑃T∗ Λ 2 ∥ 2𝐹 = ∥Σ 𝑃∗ 𝑉𝑃T∗ Λ 2 ∥ 2𝐹 = tr(Σ 𝑃∗ 𝑉𝑃T∗ Λ𝑉𝑃∗ Σ 𝑃∗ ).
(23)
We now apply Lemma 3 by setting 𝐴 = 𝑉𝑃T∗ Λ𝑉𝑃∗ and 𝐵 = Σ2𝑃∗ , and obtain that ∗
∥𝑃 Λ
1 2
∥ 2𝐹 = tr(Σ 𝑃∗ 𝑉𝑃T∗ Λ𝑉𝑃∗ Σ 𝑃∗ ) ≥
𝑘 ∑︁ 𝑖=1
17
𝜆 (𝑘+1−𝑖) × 𝜎𝑖2 (𝑃∗ ).
(24)
1
We may similarly prove that ∥Λ 2 𝑄 ∗ ∥ 2𝐹 ≥ 1
2 ∗ 𝑖=1 𝜆 (𝑘−𝑖+1) × 𝜎𝑖 (𝑄 ). Therefore,
Í𝑘 1
∥Λ 2 𝑄 ∗ ∥ 2𝐹 + ∥𝑃∗ Λ 2 ∥ 2𝐹 ≥
𝑘 ∑︁
𝜆 (𝑘+1−𝑖) × (𝜎𝑖2 (𝑃∗ ) + 𝜎𝑖2 (𝑄 ∗ ))
(25)
𝑖=1 1
1
(25) provides a lower bound of ∥Λ 2 𝑄 ∗ ∥ 2𝐹 + ∥𝑃∗ Λ 2 ∥ 2𝐹 in terms of 𝜎𝑖 (𝑃∗ ) and 𝜎𝑖 (𝑄 ∗ ). We next aim to express the lower bound in terms of 𝜎𝑖∗ ’s (singular values of 𝑃∗ 𝑄 ∗ ) directly. 1
1
1
1
The following program gives a lower bound for ∥Λ 2 𝑄 ∗ ∥ 2𝐹 + ∥𝑃∗ Λ 2 ∥ 2𝐹 : ∥Λ 2 𝑄 ∗ ∥ 2𝐹 + ∥𝑃∗ Λ 2 ∥ 2𝐹
min :
∗
subject to 𝑊 = 𝑃 𝑄
(26)
∗
𝜎𝑖 (𝑊) = 𝜎𝑖∗
for 𝑖 ≤ 𝑘.
T . Also, let 𝑃˜ = 𝑈 T 𝑃 ∗ and 𝑄 ˜ = 𝑄 ∗𝑉𝑊 . Noting that the columns in 𝑃∗ are Write the SVD of 𝑊 be 𝑈𝑊 Σ𝑊 𝑉𝑊 𝑊 ∗ ˜ and in the column space of 𝑊 and the rows in 𝑄 are in the row space of 𝑊, we have (i) 𝜎𝑖 (𝑃∗ ) = 𝜎𝑖 ( 𝑃) 1 1 1 1 ∗ 2 ∗ 2 2 2 ∗ ˜ ˜ ˜ 𝜎𝑖 (𝑄 ) = 𝜎𝑖 (𝑄) for 𝑖 ≤ 𝑘, and (ii) ∥Λ 2 𝑄 ∥ 𝐹 + ∥𝑃 Λ 2 ∥ 𝐹 = ∥Λ 2 𝑄∥ 𝐹 + ∥ 𝑃Λ 2 ∥ 𝐹 .
Therefore, (26) can be equivalently written as 1
1
˜ 2 + ∥ 𝑃Λ ˜ 2 ∥ 2𝐹 ∥Λ 2 𝑄∥ 𝐹 Σ𝑊 = 𝑃˜𝑄˜
min : subject to
(Σ𝑊 )𝑖,𝑖 = 𝜎𝑖∗
(27)
for 𝑖 ≤ 𝑘.
Now 𝑃˜𝑄˜ is positive definite. Using a similar technique developed in [Bao et al., 2020] (Theorem 1), one 1 2 can see that 𝑃˜ = 𝑄˜ T . See also Lemma 4. This implies that 𝑃˜ = Σ𝑊 Ω for some unitary matrix Ω and √︁ ∗ ∗ ∗ ˜ ˜ 𝜎𝑖 (𝑃 ) = 𝜎𝑖 ( 𝑃) = 𝜎𝑖 (𝑄 ) = 𝜎𝑖 (𝑄) = 𝜎 for 𝑖 ≤ 𝑘. Together with (25), we have 𝑖
1
1
∥Λ 2 𝑄 ∗ ∥ 2𝐹 + ∥𝑃∗ Λ 2 ∥ 2𝐹 ≥ 2
∑︁
𝜆 (𝑘+1−𝑖) × 𝜎𝑖2 (𝑃∗ ) = 2
𝑖≤𝑘
B.4
∑︁
𝜆 (𝑘+1−𝑖) × 𝜎𝑖∗ = 2∥𝑃∗ 𝑄 ∗ ∥ 𝜔,∗ .
𝑖≤𝑘
Proof of Corollary 1
We first find an optimal solution for (9). Let the SVD of 𝑋 be 𝑋 = 𝑈𝑋 Σ𝑋𝑉𝑋T , where 𝑈𝑋 ∈ R𝑚×𝑛 , Σ𝑋 ∈ R𝑛×𝑛 , and 𝑉𝑋 ∈ R𝑛×𝑛 . Let 𝑈¯ 𝑋 be an arbitrary basis for the subspace that is orthogonal to 𝑋’s column space so 𝑈¯ 𝑋 ∈ R𝑚× (𝑚−𝑛) and [𝑈𝑋 , 𝑈¯ 𝑋 ] form a basis for R𝑚 . We have 1
1
∥ 𝑋 − 𝑃𝑄∥ 2𝐹 + ∥𝑃Λ 2 ∥ 2𝐹 + ∥Λ 2 𝑄∥ 2𝐹 T T T 2 2 1 1 𝑈𝑋 𝑈𝑋 𝑈𝑋 2 𝑋𝑉𝑋 − ¯ T 𝑃𝑄𝑉𝑋 + + ∥Λ 2 𝑄𝑉𝑋 ∥ 2𝐹 . = ¯ T 𝑃Λ 𝑈 𝑈¯ 𝑋T 𝑈𝑋 𝑋 𝐹 𝐹 Let 𝑃˜ =
𝑈𝑋T 𝑈¯ 𝑋T
𝑃 and 𝑄˜ = 𝑄𝑉𝑋 . Then our objective becomes
min ˜ 𝑄˜ 𝑃,
Σ𝑋 0 (𝑚−𝑛) ×𝑛
2
− 𝑃˜𝑄˜ 𝐹
18
˜ 2. ˜ 12 ∥ 2𝐹 + ∥Λ 12 𝑄∥ + ∥ 𝑃Λ 𝐹
(28)
Let 𝑊˜ = 𝑃˜𝑄˜ and the singular values of 𝑊˜ be 𝜎1∗ ≥ 𝜎2∗ ≥ · · · ≥ 𝜎𝑘∗ . Let also Σ̃ =
Σ𝑋 0 (𝑚−𝑛) ×𝑛
.
Recall also that 𝜎𝑖 is the 𝑖-th largest singular value of 𝑋. We next show that
Σ𝑋 0
2
˜ 2 ≥ = ∥ Σ̃ − 𝑃˜𝑄∥ 𝐹
− 𝑃˜𝑄˜ 𝐹
𝑘 ∑︁
(𝜎𝑖 − 𝜎𝑖∗ ) 2 +
𝑖=1
𝑛 ∑︁
𝜎𝑖2 .
𝑖=𝑘+1
Note first that ˜ ∥ Σ̃ − 𝑊˜ ∥ 2𝐹 = ∥ Σ̃∥ 2𝐹 + ∥𝑊˜ ∥ 2𝐹 − 2⟨Σ̃, 𝑊⟩.
(29)
Next, we have ([Zheng et al., 2018]): ˜ = |tr( Σ̃𝑊˜ T ) ∥ ≤ |tr( Σ̃Σ𝑊˜ )| = |⟨Σ̃, 𝑊⟩|
𝑘 ∑︁
𝜎𝑖 𝜎𝑖∗ .
𝑖=1
˜ is maximized when Therefore, ⟨Σ̃, 𝑊⟩ 𝑊˜ 𝑖, 𝑗 =
𝜎𝑖∗ if 𝑖 = 𝑗 ≤ 𝑘 0 Otherwise.
When we plug in this optimized 𝑊˜ to Eq. 29, we get ∥ Σ̃ − 𝑊˜ ∥ 2𝐹 ≥
𝑘 ∑︁
(𝜎𝑖 − 𝜎𝑖∗ ) 2 +
𝑖=1
𝑛 ∑︁
𝜎𝑖2 .
𝑖=𝑘+1
Next, from Proposition 2, we have ˜ 2 ≥ ˜ 21 ∥ 2𝐹 + ∥Λ 12 𝑄∥ ∥ 𝑃𝑉 𝐹
𝑘 ∑︁
𝜆 (𝑘−𝑖+1) 𝜎𝑖∗ .
𝑖=1
Therefore, we can find a lower bound for Eq. 5 in terms of 𝜎𝑖∗ ’s:
L (𝜎1∗ , . . . , 𝜎𝑘∗ ) =
𝑘 ∑︁
(𝜎𝑖 − 𝜎𝑖∗ ) 2 + 2
𝑖=1
𝑘 ∑︁
𝜆 (𝑘−𝑖+1) 𝜎𝑖∗ +
𝑖=1
𝑚 ∑︁
𝜎𝑖2
(𝜎1∗ ≥ · · · ≥ 𝜎𝑘∗ ≥ 0).
(30)
𝑖=𝑘+1
We next find a minimal value of L (by treating 𝜎𝑖∗ ’s as decision variables). This will give us a lower bound (and is independent of 𝜎𝑖∗ ) on our optimization problem. We then show that this lower bound can be achieved ˜ This means such 𝑊˜ is optimal. by carefully constructing 𝑊˜ (as well as 𝑃˜ and 𝑄). Specifically, we need to find an optimal solution for the following program:
minimize 𝜎1∗ ,..., 𝜎𝑘∗ subject to:
L (𝜎1∗ , . . . , 𝜎𝑘∗ ) 𝜎𝑖∗ ≥ 0 𝜎1∗ ≤ 𝜎2∗ ≤ · · · ≤ 𝜎𝑘∗ 19
(31) (Ordering constraint)
We shall first find an optimal solution for
L (𝜎1∗ , . . . , 𝜎𝑘∗ )
minimize 𝜎1∗ ,..., 𝜎𝑘∗
(32)
𝜎𝑖∗ ≥ 0
subject to:
Note here, the ordering constraint is removed so the optimal value for (32) should be no more than that for (31). We shall see that the optimal solution for (31) also satisfies the ordering constraint so indeed optimal solutions for (31) and (32) are the same. The problem (32) boils down to finding min (𝜎𝑖 − 𝜎𝑖∗ ) 2 + 2 ∗
𝑘 ∑︁
𝜎𝑖 ≥0
𝜆 (𝑘−𝑖+1) 𝜎𝑖∗ .
𝑖=1
We note that 𝜎𝑖∗ ’s do not interact with each other so we can optimize each 𝜎𝑖∗ ’s independently. We get 𝜎𝑖∗ = (𝜎𝑖 − 𝜆 (𝑘−𝑖+1) ) + . We can check that 𝜎1∗ ≥ · · · ≥ 𝜎𝑘∗ . Therefore, the optimal value for (31) is 𝑘 ∑︁
+ 2
(𝜎𝑖 − (𝜎𝑖 − 𝜆 (𝑘−𝑖+1) ) ) + 2
𝑖=1
𝑘 ∑︁
𝑛 ∑︁
+
𝜆 (𝑘−𝑖+1) (𝜎𝑖 − 𝜆 (𝑘−𝑖+1) ) +
𝜎𝑖2 .
𝑖=𝑘+1
𝑖=1
This is also a lower bound for (9). One can check that when we set 𝑃 and 𝑄 as
√︃ √︃ 𝑃∗ = 𝑈 𝑘 𝑑𝑖𝑎𝑔( (𝜎1 − 𝜆 (𝑘 ) ) + ), . . . , (𝜎𝑘 − 𝜆 (1) ) + )Ω, √︃ √︃ 𝑄 ∗ = ΩT 𝑑𝑖𝑎𝑔( (𝜎1 − 𝜆 (𝑘 ) ) + ), . . . , (𝜎𝑘 − 𝜆 (1) ) + )𝑉𝑘T ,
(33)
the lower bound is achieved so (33) gives an optimal solution. Here, 𝑈 𝑘 and 𝑉𝑘 are leading left and right singular vectors of 𝑋. Now we move to analyze (10). Our goal is to reduce (10) to (9). Let 1
1
𝑃 = 𝑋 𝐴Λ− 2
𝑄 = Λ 2 𝐵.
Then (10) becomes minimize 𝑃,𝑄 subject to
1
1
∥ 𝑋 − 𝑃𝑄∥ 2𝐹 + ∥𝑃Λ 2 ∥ 2𝐹 + ∥Λ 2 𝑄∥ 2𝐹 𝑃 = 𝑋 𝐴Λ
− 12
1 2
𝑄=Λ 𝐵
(34)
(Constraint P)
(Constraint Q).
Here, 𝑋 and Λ are given, whereas 𝑃, 𝑄, 𝐴, and 𝐵 are decision variables. The (Constraint P) says that each column of 𝑃 needs to be in a column space of 𝑋 (it is a necessary and sufficient condition for 𝐴 to exist). The (Constraint Q) simply says 𝑄 and 𝐵 are linearly related and does not have tangible impact to the optimization problem. 20
But we note that when we put aside the constraints, an optimal (𝑃, 𝑄) is specified by (33). The columns of the optimal 𝑃 indeed is in the column space of 𝑋. So (𝑃, 𝑄) is also an optimal solution for (34). We may find the corresponding 𝐴 and 𝐵: 1
𝐴∗ = 𝑋 † 𝑃 ∗ Λ 2
B.5
1
and 𝐵∗ = Λ− 2 𝑄 ∗ ,
Symmetric lemma
˜ 𝑄˜ ∈ R 𝑘×𝑘 be full rank, Λ be a diagonal matrix, and Σ𝑊 be a diagonal matrix so that Lemma 4. Let 𝑃, ∗ (Σ𝑊 )𝑖,𝑖 = 𝜎𝑖 , where 𝜎𝑖∗ ’s are sorted in descending order. Consider the optimization problem: 1 ˜ 2 + ∥ 𝑃Λ ˜ 21 ∥ 2𝐹 ∥Λ 2 𝑄∥ 𝐹 Σ𝑊 = 𝑃˜𝑄˜
min : subject to
(Σ𝑊 )𝑖,𝑖 = 𝜎𝑖∗
(35)
for 𝑖 ≤ 𝑘.
There is an optimal solution such that 𝑃˜ = 𝑄˜ T ˜ The program (36) is equivalent to ˜ 12 and 𝑄ˆ = Λ− 21 𝑄. Proof. Let 𝑃ˆ = 𝑃Λ
ˆ 2 + ∥ 𝑃∥ ˆ 2𝐹 ∥Λ𝑄∥ 𝐹 Σ𝑊 = 𝑃ˆ𝑄ˆ
min : subject to
(Σ𝑊 )𝑖,𝑖 = 𝜎𝑖∗
(36) for 𝑖 ≤ 𝑘.
Let the SVD of 𝑄ˆ be 𝑈𝑄ˆ Σ𝑄ˆ 𝑉 Tˆ so 𝑄ˆ −1 = 𝑉𝑄ˆ Σ −1 𝑈 Tˆ . We can also see that 𝑃ˆ = Σ𝑊 𝑄ˆ −1 . Therefore, the 𝑄 𝑄ˆ 𝑄 objective term becomes −1 T 2 2 −1 2 ∥Λ𝑈𝑄ˆ Σ𝑄ˆ 𝑉𝑄Tˆ ∥ 2𝐹 + ∥Σ𝑊 𝑉𝑄ˆ Σ𝑄 ˆ 𝑈𝑄ˆ ∥ 𝐹 = ∥Λ𝑈𝑄ˆ Σ𝑄ˆ ∥ 𝐹 + ∥Σ𝑊 𝑉𝑄ˆ Σ𝑄ˆ ∥ 𝐹 .
Let us consider the stationary points 𝑈𝑄ˆ and 𝑉𝑄ˆ when Σ𝑄ˆ is fixed. We can see that they need to be permutation matrices to minimize both terms in the objective (using the rearrangement inequality again). Therefore, we 1 can see 𝑄ˆ = Σ1 Σ𝑄ˆ Σ2 for two permutation matrices Σ1 and Σ2 . This implies that 𝑄˜ = Λ 2 Σ1 Σ𝑄ˆ Σ2 , i.e., each row (column) of 𝑄˜ has exactly one non-zero entry. We may similarly show that each row (column) of 𝑃˜ has exactly one non-zero entry. In addition, the locations of non-zero entries of 𝑃˜ and 𝑄˜ T are identical because 𝑃˜𝑄˜ is a diagonal matrix. We may thus write 𝑃˜ = Σ (1) Σ 𝑃˜ Σ (2)
T 𝑄˜ = ΣT(2) Σ ( 𝑄) ˜ Σ (1) ,
˜ where 𝜏 is a permutation on [𝑘]. The set of (possibly ˜ and (Σ ( 𝑄) where (Σ 𝑃˜ )𝑖,𝑖 = 𝜎𝑖 ( 𝑃) ˜ )𝑖,𝑖 = 𝜎𝜏 (𝑖) ( 𝑄), ˜ Thus, we can see that there exists a permutation 𝜋¯ ˜ 𝜏 (𝑖) (𝑄). unsorted) singular values for 𝑃˜𝑄˜ thus is 𝜎𝑖 ( 𝑃)𝜎 such that ∑︁ 1 1 ˜ 2 + ∥ 𝑃Λ ˜ 2 ∥2 = ∥Λ 2 𝑄∥ 𝜎𝑖2 (𝑃∗ )𝜆 𝜋¯ (𝑖) + 𝜎𝜏2 (𝑖) (𝑄 ∗ )𝜆 𝜋¯ (𝑖) 𝐹 𝐹 𝑖≤𝑘
≥
∑︁
2𝜎𝑖 (𝑃∗ )𝜎𝜏∗ (𝑖) (𝑄)𝜆 𝜋¯ (𝑖)
𝑖≤𝑘
≥ 2∥𝑃∗ 𝑄 ∗ ∥ 𝜔,∗ . One can see that we can set 𝑃˜ = 𝑄˜ T to make all inequality becomes equality so there is an optimal solution such that 𝑃˜ = 𝑄˜ T . □ 21
C
Experimental Details
C.1
DLAE hyperparameters tuning
This section presents the hyperparameter tuning process on the validation data over three (ML-20M, Netflix, MSD) datasets for the full rank DLAE formula, which was introduced by Steck [2020] yet not investigated: min ||𝑋 − 𝑋𝑊 || 2𝐹 + ||Λ1/2 · 𝑊 || 2𝐹 𝑊
𝑝 𝑑𝑀𝑎𝑡 (𝑑𝑖𝑎𝑔(𝑋 𝑇 𝑋)) 1− 𝑝 𝑊ˆ = (𝑋 𝑇 𝑋 + Λ) −1 𝑋 𝑇 𝑋 Λ=
In practical, 𝑙 2 regularization is also imposed: 𝑊ˆ = (𝑋 𝑇 𝑋 + Λ + 𝜆) −1 𝑋 𝑇 𝑋 The tables 4 to 6 show the results of 𝑛𝐷𝐶𝐺@100 over three datasets respectively. And the optimal parameters are highlighted.
C.2
Matrix Factorization with Dropout Hyperparameters Tuning
Cavazza et al. [2018] shows that optimization with dropout (allowing rank optimizing) is equivalent to solving a matrix approximation problem with nuclear norm: min ||𝑋 − 𝑃𝑄 𝑇 || 2𝐹 + 𝑑 𝑃,𝑄,𝑑
min ||𝑋 − 𝑌 || 2𝐹 + 𝑌
𝑑 1 − 𝑝 ∑︁ · ||𝑃 𝑘 || 22 · ||𝑄 𝑘 || 22 𝑝 𝑘=1
1− 𝑝 ||𝑌 || 2∗ 𝑝
and the solution is given by: SVD
𝑋 = 𝑈Σ𝑉 𝑇 𝑌 ∗ = 𝑃∗ · (𝑄 ∗ ) 𝑇 = 𝑈 · 𝑆 𝜇 (Σ) · 𝑉 𝑇 𝑆 𝜇 (𝜎) = max(𝜎 − 𝜇, 0) 𝑑¯ ∑︁ 1− 𝑝 𝜇= 𝜎𝑖 (𝑋) 𝑝 + (1 − 𝑝) 𝑑¯ 𝑖=1
where 𝑑¯ denotes the largest integer such that: 𝑑¯ ∑︁ 1− 𝑝 𝜎𝑖 (𝑋) 𝜎𝑑¯ (𝑋) > 𝑝 + (1 − 𝑝) 𝑑¯ 𝑖=1
Hence, there is only one parameter 𝑝 to tuning. We present the tuning process on the validation set below, see tables 7 to 9. Optimal parameters as well as induced rank 𝑑¯ are highlighted.
C.3
Resources
Our code are mainly implemented in Numpy 1.19, Pytorch 1.7.1 on CUDA 11.0. Our experiments are performed on nodes with two sockets, each containing a 24-core Intel(R) Xeon(R) Platinum 8268 CPU @ 2.90GHz and 4 GeForce RTX 3090 24GB memory GPU. 22
Table 3: Investigating the closed/analytic solutions of linear models. 𝑑𝑀𝑎𝑡 (·) denotes a diagonal matrix, 𝑑𝑖𝑎𝑔(𝑋) is the vector on the diagonal of 𝑋. Model
regularization
solution
min ||𝑋 − 𝑋𝑊 || 2𝐹 + 𝜆 · ||𝑊 || 2𝐹 𝑊
1. EASE(full rank) [Steck, 2019]
𝐶 = (𝑋 𝑇 𝑋 + 𝜆𝐼) −1 𝑊 = 𝐼 − 𝐶 · 𝑑𝑀𝑎𝑡 (𝑑𝑖𝑎𝑔(1 ⊘ 𝐶))
𝑑𝑖𝑎𝑔(𝑊) = 0
𝑠.𝑡.
min ||𝑋 − 𝑋𝑊 || 2𝐹 + ||Λ1/2 · 𝑊 || 2𝐹
Frobenius norm
𝑊
2. DLAE(full rank) [Steck, 2020]
Λ=
𝑊 = (𝑋 𝑇 𝑋 + Λ) −1 𝑋 𝑇 𝑋
𝑝 𝑑𝑀𝑎𝑡 (𝑑𝑖𝑎𝑔(𝑋 𝑇 𝑋)) 1− 𝑝
min ||𝑋 − 𝑋𝑊 || 2𝐹 + ||Λ1/2 · 𝑊 || 2𝐹 𝑊
𝐶 = (𝑋 𝑇 𝑋 + Λ) −1
𝑝 𝑑𝑀𝑎𝑡 (𝑑𝑖𝑎𝑔(𝑋 𝑇 𝑋)) 1− 𝑝 𝑠.𝑡. 𝑑𝑖𝑎𝑔(𝑊) = 0
3. EDLAE(full rank) [Steck, 2020]
Λ=
𝑊 = 𝐼 − 𝐶 · 𝑑𝑀𝑎𝑡 (𝑑𝑖𝑎𝑔(1 ⊘ 𝐶))
min ||𝑋 − 𝑋 𝐴𝐵𝑇 || 2𝐹 + ||Λ1/2 · 𝐴𝐵𝑇 || 2𝐹 4. EDLAE-ADMM [Steck, 2020]
𝐴,𝐵
ADMM update 𝐴, 𝐵
𝑑𝑖𝑎𝑔(𝑊) = 0
𝑠.𝑡.
∗
min
5. LRR [Jin et al., 2021]
𝑟 𝑎𝑛𝑘 (𝑊 ) ≤ 𝑘
SVD
𝑌 = 𝑋𝑊 ∗ = 𝑈Σ𝑉 b = (𝑋 𝑇 𝑋 + Γ𝑇 Γ) −1 𝑋 𝑇 𝑋 (𝑉𝑘 𝑉 𝑇 ) 𝑊
||𝑋 − 𝑋𝑊 || 2𝐹 + ||Γ𝑊 || 2𝐹
𝑘
min
𝑟 𝑎𝑛𝑘 (𝑊 ) ≤ 𝑘
6. LR-DLAE(this paper)
Λ=
𝑊 ∗ = (𝑋 𝑇 𝑋 + Λ) −1 𝑋 𝑇 𝑋
||𝑋 − 𝑋𝑊 || 2𝐹 + ||Λ1/2 · 𝑊 || 2𝐹
∗
min
𝑟 𝑎𝑛𝑘 (𝑊 ) ≤ 𝑘
𝑘
𝐶 = (𝑋 𝑇 𝑋 + Λ) −1
||𝑋 − 𝑋𝑊 || 2𝐹 + ||Λ1/2 · 𝑊 || 2𝐹
𝑊 ∗ = 𝐼 − 𝐶 · 𝑑𝑀𝑎𝑡 (𝑑𝑖𝑎𝑔(1 ⊘ 𝐶))
𝑝 𝑑𝑀𝑎𝑡 (𝑑𝑖𝑎𝑔(𝑋 𝑇 𝑋)) 1− 𝑝 𝑠.𝑡. 𝑑𝑖𝑎𝑔(𝑊) = 0
7. LR-EDLAE-1(this paper)
SVD
𝑌 = 𝑋𝑊 ∗ = 𝑈Σ𝑉 𝑇 b = 𝑊 ∗ (𝑉𝑘 𝑉 𝑇 ) 𝑊
𝑝 𝑑𝑀𝑎𝑡 (𝑑𝑖𝑎𝑔(𝑋 𝑇 𝑋)) 1− 𝑝
Λ=
∗
SVD
𝑌 = 𝑋𝑊 ∗ = 𝑈Σ𝑉 𝑇 b = 𝑊 ∗ (𝑉𝑘 𝑉 𝑇 ) 𝑊 𝑘
min
𝑟 𝑎𝑛𝑘 (𝑊 ) ≤ 𝑘
𝐶 = (𝑋 𝑇 𝑋 + Λ) −1
||𝑋 − 𝑋𝑊 || 2𝐹 + ||Λ1/2 · 𝑊 || 2𝐹
𝑊 ∗ = 𝐼 − 𝐶 · 𝑑𝑀𝑎𝑡 (𝑑𝑖𝑎𝑔(1 ⊘ 𝐶))
𝑝 𝑑𝑀𝑎𝑡 (𝑑𝑖𝑎𝑔(𝑋 𝑇 𝑋)) 1− 𝑝 s.t. 𝑑𝑖𝑎𝑔(𝑊) = 0
8. LR-EDLAE-2(this paper)
Λ=
SVD
𝑊 ∗ = 𝑈Σ𝑉 𝑇 b = 𝑈𝑘 Σ𝑘 𝑉 𝑇 𝑊 𝑘
𝑃∗ = 𝑈𝑘
min ||𝑋 − 𝑃𝑄 𝑇 || 2𝐹 + 𝜆 · (||𝑃|| 2𝐹 + ||𝑄|| 2𝐹 )
𝑄 ∗ = 𝑉𝑘 Ω √︁ Ω = (𝜎𝑖 − 𝜆)+
𝑃,𝑄
9. Regularized PCA [Zheng et al., 2018]
SVD
𝑋 = 𝑈Σ𝑉 𝑇
Nuclear Norm
min ||𝑋 − 𝑃𝑄 𝑇 || 2𝐹 + 𝑑 10. MF dropout [Cavazza et al., 2018]
𝑃,𝑄,𝑑
min ||𝑋 − 𝑌 || 2𝐹 + 𝑌
𝑑 1 − 𝑝 ∑︁ · ||𝑃 𝑘 || 22 · ||𝑄 𝑘 || 22 𝑝 𝑘=1
SVD
𝑋 = 𝑈Σ𝑉 𝑇 𝑌 ∗ = 𝑃∗ · (𝑄 ∗ ) 𝑇
1− 𝑝 ||𝑌 || 2∗ 𝑝
= 𝑈 · 𝑆 𝜇 (Σ) · 𝑉 𝑇
1
11. LAE [Bao et al., 2020]
1
𝑊1∗ = 𝑃(𝐼 − Λ𝑆 −2 ) 2 𝑈 𝑇
1
min ∥ 𝑋 − 𝑋𝑊1𝑊2 ∥ 2𝐹 + ∥𝑊1 Λ 2 ∥ 2𝐹 + ∥Λ 2 𝑊2 ∥ 2𝐹 ,
1 𝑊2∗ = 𝑈 (𝐼 − Λ𝑆 −2 ) 2 𝑃𝑇
𝑊1 ,𝑊2
SVD
min ||𝑋 − 𝑃𝑄|| 2𝐹 + ||Λ1/2 𝑄|| 2𝐹 + ||𝑃Λ1/2 || 2𝐹
𝑋 = 𝑈Σ𝑉 𝑇
𝑃,𝑄
12. LVAE(this paper)
min ||𝑋 − 𝑋 𝐴𝐵|| 2𝐹 + ||Λ𝐵|| 2𝐹 + ||𝑋 𝐴|| 2𝐹 𝐴,𝐵
min
𝑟 𝑎𝑛𝑘 (𝑊 ) ≤ 𝑘
||𝑋 − 𝑊 || 2𝐹 + 2||𝑊 || 𝑤,∗
23
∗
√︁ √︁ 𝑃 = 𝑈 𝑘 · 𝑑𝑖𝑎𝑔( 𝜎1 − 𝜆 (𝑘 ) , . . . , 𝜎1 − 𝜆 (1) ) · Ω √︁ √︁ ∗ 𝑇 𝑄 = Ω · 𝑑𝑖𝑎𝑔( 𝜎1 − 𝜆 (𝑘 ) , . . . , 𝜎1 − 𝜆 (1) ) · 𝑉𝑘𝑇 1
𝐴∗ = 𝑋 † 𝑃 ∗ Λ 2
1
𝐵∗ = Λ − 2 𝑄 ∗
Table 4: ml-20m, DLAE full rank, parameter tuning on validation dataset by 𝑛𝐷𝐶𝐺@100 𝜆 0.1 0.2 0.3 0.4 0.5
p
800 0.42024 0.43132 0.43203 0.43001 0.42754
900 0.42063 0.43139 0.43211 0.43001 0.42745
1100 0.42102 0.4314 0.43206 0.42996 0.42718
1000 0.42073 0.43154 0.43214 0.42995 0.42729
1200 0.42131 0.43147 0.43203 0.42984 0.42715
1300 0.4212 0.43136 0.43196 0.42978 0.42704
Table 5: netflix, DLAE full rank, parameter tuning on validation dataset by 𝑛𝐷𝐶𝐺@100
p
0.2 0.25 0.3 0.35 0.4 0.45 0.5
800 0.3904 0.39247 0.39359 0.39402 0.39399 0.39346 0.39249
900 0.3904 0.39252 0.39359 0.39403 0.39393 0.3935 0.39241
1000 0.39027 0.39248 0.39366 0.394 0.39395 0.39343 0.39247
𝜆 1100 0.39024 0.39249 0.39358 0.39405 0.39393 0.39344 0.39241
1200 0.3902 0.3925 0.39362 0.39399 0.39389 0.39338 0.39242
1300 0.3903 0.3925 0.39368 0.39403 0.39388 0.3933 0.3923
1400 0.39018 0.39256 0.39369 0.39397 0.39387 0.39329 0.39224
Table 6: msd, DLAE full rank, parameter tuning on validation dataset by 𝑛𝐷𝐶𝐺@100 𝜆
p
0.3 0.4 0.5 0.6
10 0.38514 0.38596 0.38556 0.38382
20 0.38515 0.38599 0.38555 0.3838
30 0.38517 0.38602 0.38553 0.38381
40 0.38505 0.386 0.38557 0.38374
50 0.38492 0.38597 0.38553 0.38373
60 0.38474 0.38592 0.38549 0.38366
Table 7: ml-20m, matrix factorization with dropout, hyper parameter tuning by 𝑛𝐷𝐶𝐺@100 on validation dataset and its induced rank . p induced rank 𝑑¯ nDCG@100
0.9 10 0.29723
0.99 200 0.39369
0.995 385 0.40045
0.996 467 0.40046
0.997 602 0.39925
Table 8: netflix, matrix factorization with dropout, hyper parameter tuning by 𝑛𝐷𝐶𝐺@100 on validation dataset and its induced rank . p induced rank 𝑑¯ nDCG@100
0.9 9 0.26026
0.99 209 0.35462
0.996 524 0.36453
0.997 653 0.36495
0.998 883 0.36406
Table 9: msd, matrix factorization with dropout, hyper parameter tuning by 𝑛𝐷𝐶𝐺@100 on validation dataset and its induced rank . p induced rank 𝑑¯ nDCG@100
0.99 249 0.18986
0.999 2054 0.28532
24
0.9995 3783 0.307
0.9999 11380 0.32634
0.99995 19308 0.30995