Year: 2021

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 deep learning 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 competitive linear models effectively add only a nuclear-norm regularizer, or a Frobenius-norm 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 so they have better performance. However, the performance leaders EASE and EDLAE have only full-rank closed-form solutions. Their low-rank counterparts can be solved only by hard-to-tune numeric procedures such as Alternating Directions Method of Multipliers (ADMM). Along this line of findings, we further propose two lowrank, closed-form solutions, derived from carefully generalizing Frobenius-norm regularization techniques.

heading: 1 INTRODUCTION
text: Research progress on algorithms for recommendation has escalated in recent years, partially fueled by the adoption of deep learning (DL) techniques. However, recent studies have found that many new deep learning recommendation models have shown sub-par performance against simpler linear recommendation models (Dacrema et al., 2019; Rendle et al., 2019). Although some studies are available to analyze linear vs non-linear models (Dacrema et al., 2019), it remains puzzling why these seemingly different techniques all result in models with similar performance or even similar cost functions. In the latest study, Jin et al. (2021) examine the relationship between the widely used matrix factorization (MF), such as ALS (Hu et al., 2008), and the linear autoencoder (LAE) which encompasses the recent performance leaders, such as EASE (Steck, 2019) and EDLAE (Steck, 2020). They considered two basic regularization forms (see Eqs. (1) and (6)) and found that the optimal (closed-form) solutions of both models recover the directions of principal components, while shrinking the corresponding singular values differently. They suggest that the difference may enable LAE to utilize a larger number of latent dimensions to improve recommendation accuracy. This finding highlights the similarity and 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) recommendation algorithms, and aim to provide an in-depth understanding of various regularized objective functions and towards unifying them through their closed-form solutions. In return, the closed-forms serve as the barebones engine to help reveal what drives the performance improvement for the recent recommendation algorithms.
We examine three open and closely related questions: (i) Characterization of models: Recent recommendation models are built upon a diverse set of techniques, including dropouts (Srivastava et al., 2014), variational methods (Kingma & Welling, 2014), and matrix denoising (Tipping & Bishop, 1999b). We aim to identify a unified regularization framework to (re)-interpret different models and understand how they are related. We are specifically interested in how models based on variational autoencoder (VAE, a natural generalization of recent approaches) relate to other models (Liang et al., 2018). (ii) Weighted generalization of regularizers: A key idea in a recent performance leader EDLAE (Steck, 2020) is to utilize of weighted/non-uniform regularizers on the
parameters’ Frobenius norm. Applying dropouts(Cavazza et al., 2018) is shown to be equivalent to re-weighting the exponents in the regularizers, such as designing the weighted sum of regularizers based on other norms or tuning the exponent weights. We are also interested in determining what circumstance closed-form solutions still exist. (iii) Low-rank closed-form solutions for EDLAE. While most linear recommendation models are shown to have closed-form solutions, the low-rank version of recent performance leader EDLAE (Steck, 2020) remains an exception. Fitting this model requires using ADMM (Steck, 2020). Although its performance may not be better than its full-rank counterpart, low-rank solutions are easier to interpret, use less storage, and can be more scalable with respect to the number of items. In addition, a closed-form solution disentangles key performance drivers from nuances (e.g., need to tune learning rate or deal with local optimal), and can help reveal the key driven factors comparing with other closed-form solutions. More importantly, such solutions are significantly easier to implement and be tested using a generic matrix computation platform (without specialized recommendation library). Can we approximate low-rank EDLAE with closed-form solutions?
Our investigation led to the following discovery and resolves the above questions:
Regularizer dichotomy (Section 2): We found that all of the leading (linear) recommendation models can be categorized into those that implement possibly weighted nuclear-norm regularizers, or those that implement Frobenius-norm regularizers. Specifically, we characterize the Variational Linear Autoencoders (VLAE) as a form of the weighted nuclear-norm regularization problem, in which the weights possess a specific combinatorial structure. 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 form of regularization. Thus, our paper provides a more complete and accurate characterization of a linear recommendation model’s performance under different regularizations.
Rigidity of nuclear-norm regularizers (Section 3): With the dichotomy result, we next aim to understand whether the weighted sum idea for Frobenius-norm regularizers (Steck, 2020; Jin et al., 2021) can be generalized to nuclear-norm regularizers, and whether tuning exponent weights can be beneficial. First, while VLAE is equivalent to nuclear-norm regularization and is easier to optimize (weighted nuclear norms are non-smooth but VLAE’s objectives are smooth), its closedform solution possesses a (surprisingly) rigid structure, i.e., the weighted regularization will lead to an auto-sorting singular value reduction and the larger single values tend to receive smaller reduction (non-ascending reduction). Second, it has been shown that the solution structure for a model with a squared nuclear-norm regularizer ∥W∥2∗ (i.e., dropout’s equivalence) is strikingly similar to that for using ∥W∥∗-regularizer (Gu et al., 2014). We generalize the result to show that the solution structures for ∥W∥p∗ are highly similar for all p ≥ 1 (Regularization invariant/rigidity with respect to p). But when p = 1, 2, the solution and hyperparameters possess favorable properties which make hyperparameter searches easier. This also partially explains why only p = 1, 2 have been extensively considered. These rigidity properties severely limit the search space and explains why models that use only nuclear-norm regularizers share the same performance ceiling even when hyperparameters are extensively searched.
Closed-form solution for low-rank EDLAE (Section 4): The (weighted) Frobenius-norm regularizers ∥ΛW∥2F (Λ is the hyperparameter diagonal matrix) 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 W is enforced, their singular vectors will no longer coincide with those of the data, and will deliver (slightly) better performance. However, no closed-form solution for the low-rank estimator is known and existing approaches rely on ADMM or Alternating Least Square (Steck, 2020). In this paper, we propose two 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), and thus resolve the aforementioned third open question.

heading: 2 BACKGROUND AND REGULARIZER DICHOTOMY
text: This section first explains the background, and then give an overview of the dichotomy results. Some key theorems are deferred to Section 3 and Appendices.
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., 2019; Zhang et al., 2019). We focus on the implicit
Table 1: Investigating the closed/analytic solutions of linear models. dMat(·) denotes a diagonal matrix, diag(X) is the vector on the diagonal of X . Λ is the (hyperparameter) diagonal matrix as a coefficient of the regularization term. W ∗ (or P ∗, Q∗, etc) is the optimal solution for corresponding
case, except for cases 9-12, where Ŵ is the low-rank closed-form solution.
Model Regularization Solution
Nuclear norm
1. Regularized PCA (Zheng et al., 2018) min P,Q
||X − PQ||2F + λ · (||P ||2F + ||Q||2F ) X SVD = UΣV T Ω =
√ (σi − λ)+
P ∗ = Uk Q ∗ = ΩV Tk
2. MF dropout (Cavazza et al., 2018) min P,Q,d
||X − PQ||2F + d 1− p p
· d∑
k=1
||Pk||22 · ||QTk ||22
min Y ||X − Y ||2F + 1− p p ||Y ||2∗
X SVD = UΣV T Y ∗ = P ∗ ·Q∗ = U · Sµ(Σ) · V T
3. WLAE (Bao et al., 2020) min W1,W2
∥X −XW1W2∥2F + ∥W1Λ 1 2 ∥2F + ∥Λ 1 2W2∥2F ,
W ∗1 = V (I − ΛS−2) 1 2PT W ∗2 = P (I − ΛS−2) 1 2V T
4. VLAE (this paper)
min P,Q
||X − PQ||2F + ||Λ1/2Q||2F + ||PΛ1/2||2F
min A,B
||X −XAB||2F + ||ΛB||2F + ||XA||2F
min rank(W )≤k
||X −W ||2F + 2||W ||w,∗
X SVD = UΣV T P ∗ = Uk · diag( √ σ1 − λ(k), . . . , √ σk − λ(1)) · Ω
Q∗ = ΩT · diag( √ σ1 − λ(k), . . . , √ σk − λ(1)) · V Tk A∗ = X†P ∗Λ 1 2 B∗ = Λ− 1 2Q∗
Frobenius norm
5. EASE (full rank) (Steck, 2019) min W ||X −XW ||2F + λ · ||W ||2F s.t. diag(W ) = 0
C = (XTX + λI)−1 W ∗ = I − C · dMat(diag(1⊘ C))
6. DLAE(full rank) (Steck, 2020) min W
||X −XW ||2F + ||Λ1/2 ·W ||2F
Λ = p
1− p dMat(diag(XTX))
W ∗ = (XTX + Λ)−1XTX
7. EDLAE(full rank) (Steck, 2020)
min W
||X −XW ||2F + ||Λ1/2 ·W ||2F
Λ = p
1− p dMat(diag(XTX))
s.t. diag(W ) = 0
C = (XTX + Λ)−1 W ∗ = I − C · dMat(diag(1⊘ C))
8. LRR (Jin et al., 2021) min rank(W )≤k
||X −XW ||2F + ||ΓW ||2F Y ∗ = XW ∗ SVD = UΣV
Ŵ = (XTX + ΓTΓ)−1XTX(VkV T k )
9. EDLAE-ADMM (Steck, 2020) min A,B
||X −X(AB − dMat(diag(AB)))||2F
+||Λ1/2 · (AB − dMat(diag(AB))||2F ADMM update A,B
10. LR-DLAE (this paper) min rank(W )≤k ||X −XW ||2F + ||Λ1/2 ·W ||2F
Λ = p
1− p dMat(diag(XTX))
W ∗ = (XTX + Λ)−1XTX Y ∗ = XW ∗ SVD = UΣV T Ŵ = W ∗(VkV T k )
11. LR-EDLAE-1 (this paper)
min rank(W ′)≤k
||X −XW ||2F + ||Λ1/2 ·W ||2F
W = W ′ − dMat(diag(W ′)) Λ = p
1− p dMat(diag(XTX))
C = (XTX + Λ)−1 W ∗ = I − C · dMat(diag(1⊘ C)) Y ∗ = XW ∗ SVD = UΣV T Ŵ = W ∗(VkV T k )
12. LR-EDLAE-2 (this paper)
min rank(W ′)≤k
||X −XW ||2F + ||Λ1/2 ·W ||2F
W = W ′ − dMat(diag(W ′)) Λ = p
1− p dMat(diag(XTX))
C = (XTX + Λ)−1 W ∗ = I − C · dMat(diag(1⊘ C)) W ∗ SVD = UΣV T Ŵ = UkΣkV T k
problem because it is more economically relevant. Here, let n be the number of items and m be the number of users. We are given a binary matrix X ∈ {0, 1}m×n that represents the interaction between users and items so far, i.e., Xi,j = 1 iff user i has purchased or made a rating on item j. Our goal is to produce a real-valued matrix X̂ , which we evaluate against future interactions using information retrieval (Top-k) metrics such as Recall or nDCG.
The problem is closely related to matrix completion (MC) because Xi,j = 0 can be viewed as “missing a data point”. But MC’s evaluation metric is mean-squared error and is different from ours (Candès & Tao, 2010). The connection between two problems results in models with similar objectives (Zheng et al., 2018). A technique developed for one problem often finds its counterpart for the other. Nevertheless, because of the difference in evaluation, efficacy of an algorithm for one problem does not imply its performance guarantee for the other. Thus, the non-overlapping component between two problems remains substantial. We also remark that our structural results on weighted nuclear-norm is new and applicable to MC.

heading: 2.1 DICHOTOMY OF THE MODELS
text: This section explains that (i) All linear recommendation models can be categorizeed into those that use nuclear-norm regularizers and those that use Frobenius norm (see also Table 1), and (ii) The form of regularization, instead of whether the problem shall be cast as matrix factorization or LAE, determines the shrinkage structure (Proposition 1).
Nuclear-norm regularizers. Let X ∈ Rm×n be a matrix of rank at most k with k leading singular values being σ1(X) ≥ σ2(X) ≥ · · · ≥ σk(X). Let ω = (ω1, . . . , ωk) ∈ (R+)k. The weighted nuclear norm of X with respect to ω is defined as ∥X∥ω,∗ = ∑k i=1 ωi · σi(X). This is a natural generalization of the weighted nuclear norm for low-rank matrices (Gu et al., 2014). Despite its name, the weighted nuclear norm is neither convex nor differentiable unless ωi’s are sorted in descending order (Chen et al., 2013; Iglesias et al., 2020).
Nuclear-norm 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, nuclear-norm regularizers also promote sparsity over the solution’s singular values (i.e., the solution is usually low rank). Algorithms below effectively add only a nuclear-norm regularizer to MF.
A1. Regularized PCA (Udell et al., 2016; Zheng et al., 2018) aims to solve
min P,Q
∥X − PQ∥2F + λ∥P∥2F + λ∥Q∥2F (1)
where λ is the hyperparameter. It has been known that Eq. (1) is equivalent to minX̂ ∥X − X̂∥2F + 2λ∥X̂∥∗. To solve Eq. (1), one can use factored gradient descent (Bhojanapalli et al., 2016) or directly derive its closed-form solution (Kunin et al., 2019), which involves computation of SVD of X (see case 1 in Table 1).
A2. Matrix Factorization with dropouts. This approach uses PQ (P ∈ Rm×k, Q ∈ Rk×n) to approximate X . A standard dropout technique is utilized when we train P and Q. Cavazza et al. (2018) show that optimization with dropout is equivalent to solving minX̂ ∥X − X̂∥2F + λ∥X̂∥2∗, and the closed-form solution is obtained by shrinking all singular values of X by the same magnitude of µ, which depends on the data X and the choice of hyperparameter λ. See details in Appendix E.3.
The next two approaches leverage techniques from (variational) autoencoder. We show that they also effectively add variants of nuclear-norm regularizers although this may not be clear at the first glance.
A3. Weighted Linear Autoencoder (WLAE). Consider the following (non-uniform) weighted ℓ2regularization (Bao et al., 2020):
∥X −XW1W2∥2F + ∥W1Λ 1 2 ∥2F + ∥Λ 1 2W2∥2F , (2)
where Λ is a hyperparameter diagonal matrix and W1, W2 denote the encoder and the decoder network respectively. Bao et al. (2020) have shown a closed-form solution for Eq. (2) with a specific of diagonal matrix Λ = dMat(λ1, λ2, · · · , λk) when the weight is non-descending: λ1 ≤ λ2 ≤ · · · ≤ λk. It remains unclear if a closed-form solution exists for an arbitrary weight order.
A4. Variational Linear Autoencoders (VLAE). 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 VLAE, we need to find the maximum likelihood estimation (MLE) for the probabilistic model
p(x | z) = N ( Wz+ µ, σ2I ) and p(z | x) = N (V (x− µ), D) (3)
where V , W denote the encoder and the decoder networks, σ is a parameter, and D is a diagonal covariance matrix (for the data X ∈ Rm×n). We have the following observation. Lemma 1. Consider optimizing the ELBO (Evidence Lower Bound) (Kingma & Welling, 2014) for the above LVAE model (Eq. (3)). When µ = 0, this optimization problem is equivalent to minimizing
L = ∥X −XV TWT∥2F +m∥ √ DWT ∥2F + σ2||XV T∥2F + g(D,σ) (4)
where g(D,σ) = −σ2m ( log |D| − tr(D) + k − n log(2πσ2) ) .
Here, we set µ = 0 for simplicity and following standard practices in recommendations. The proof can be found in Appendix D.1. We may further “clear up” Eq. (4) and obtain the following optimization problem (see Appendix D.3):
min A∈Rn×k,B∈Rk×n
∥X −XAB∥2F + ∥XA∥2F + ∥ΛB∥2F , (5)
in which decision variables are A and B, and the hyperparameter is a diagonal matrix Λ ∈ Rk×k. In Section 3, we show that solving Eq. (5) (A4) is equivalent to solving ∥X −W∥2F + λ∥W∥ω,∗ subject to rank(W ) ≤ k, where ω consists of Λ’s diagonal values, sorted in non-descending order. Furthermore, we characterize the optimal solution for ∥X − W∥2F + λ∥W∥ p ∗ for any p ≥ 1. We shall show that regardless the choice of p, the optimal solution uniformly shrinks all X’s singular values by a constant magnitude µ (and to 0 if a singular value is already less than µ). Finally, for all nuclear-norm-based approaches discussed above, the estimators always keep singular vectors of X and shrink its singular values. Therefore, the solution space offered by nuclear-norm regularizers is quite constrained, which limits their predictive power.
Frobenius-norm regularizers. Most algorithms below were originally motivated by the design of (denoising) autoencoders. It has been shown that they effectively add a Frobenius-norm regularizer.
A5. EASE (Steck, 2019) aims to optimize minW ||X −XW ||2F + λ · ||W ||2F subject to the constraint that diag(W ) = 0. A closed-form solution exists for this problem (case 5 in Table 1).
A6. DLAE (Steck, 2020) adds a weighted Frobenius-norm regularizer so the objective becomes minW ||X −XW ||2F + ||Λ1/2 ·W ||2F , where Λ = p 1−pdMat(diag(X
TX)), dMat(·) denotes a diagonal matrix, diag(X) is the vector on the diagonal of X . The closed-form solution also exists (case 6 in Table 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(W ) = 0. A closed-form solution exists for this problem (case 7 in Table 1). For the low-rank EDLAE (case 9 in Table 1), an ADMM algorithm is used.
A8. Tikhonov regularization/Low Rank Regression (LRR) (Jin et al., 2021). Let Vk be the k leading right singular vectors of X . Jin et al. (2021) find an estimator that solves
W = arg min rank(W )≤k
∥X −XW∥2F + ∥ΓW∥2F , (6)
where Γ = Λ 1 2V Tk and Λ = diag(λ ′ 1, · · · , λ′k) is a hyperparameter. The closed-form solution is
W ∗ = Vk · dMat( σ21
σ21 + λ ′ 1
, . . . , σ2k
σ2k + λ ′ k
)V Tk (7)
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 X 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).
Proposition 1. For any regularized instances in the form of Eq. (5) with regularization parameter Λ such that σi(X) ≥ λ(k−i) for all i, there is a corresponding Tikhonov regularized instance with Γ = Λ
1 2V Tk that provides the same regularization effect.
The proof is in Appendix D.2. A major implication of Proposition 1 is that we can focus on designing Frobenius-norm regularizers because it also gets the value from using nuclear-norm regularizers.
As mentioned earlier, no closed-form solution exists for the low-rank EDLAE (case 9 in Table 1) so existing approaches rely on ADMM or Alternating Least Square algorithms (Steck, 2020). In Section 4, we will introduce two low-rank models based on Frobenius-norm regularizers that have both competitive performance and closed-form solutions.

heading: 3 RIGIDITY OF NUCLEAR-NORM REGULARIZERS
text: This section presents two results. (i) Eq. (5) (A4) is equivalent to a model with weighted nuclear-norm, where the weights are diagonal values of Λ arranged in non-descending order. Because weighted nuclear norms are non-differentiable in general so Eq. (5) compiles non-differentiable objectives into differentiable ones, which are easier to optimize. Also, the auto-sorting property restricts Eq. (5) from expressing an arbitrary weight sequence in ∥W∥ω,∗, which limits its predictive power. (ii) We give a closed-form solution for models with regularizer ∥W∥p∗ for all p ≥ 1. The solutions share the same structure so tuning p will not improve a model’s predictive power.
Rigidity of closed-form solutions. We first analyze Eq. (5) (A4). We start with the problem:
min P,Q
∥X−PQ∥2F+∥Λ 1 2Q∥2F+∥P Λ 1 2 ||2F , or equivalently, min P,Q ∥X−PQ∥2F+∥ΛQ∥2F+∥P∥2F , (8)
Note that when we let P ′ = XA∗ and Q′ = B∗, where A∗ and B∗ are an optimal solution of Eq. (5), the syntax of Eq. (8) matches with that of Eq. (5). This implies that the solution for Eq. (8) is a lower bound of that for Eq. (5). These two solutions coincide only when the columns in the optimal P ∗ in Eq. (8) are spanned by the columns of X .
We next find a closed-form solution (P ∗, Q∗) for Eq. (8), and show that indeed the column space of P ∗ is in the column space of X . Proposition 2. Let f : Rm×n → R+ be any cost function. Let P ∈ Rm×k and Q ∈ Rk×n. Let Λ ∈ Rk×k be a diagonal matrix such that λi = Λii ≥ 0 (i ∈ [k]). Let also ω = (λπ(1), λπ(2), . . . , λπ(k)), where π is a permutation on [k] such that λπ(1) ≤ λπ(2) ≤ · · · ≤ λπ(k). The following two optimization problems have the same optimal values
OPT1 : minP,Q f(PQ) + ∥Λ 1 2Q∥2F + ∥PΛ 1 2 ∥2F .
OPT2 : minW f(W ) + 2∥W∥ω,∗ subject to rank(W ) ≤ k.
In addition, if (P ∗, Q∗) is an optimal solution for OPT1, then W ∗ = P ∗Q∗ is an optimal solution for OPT2. If W ∗ is an optimal solution for OPT2, then there exists an optimal solution (P ∗, Q∗) for OPT1 such that W ∗ = P ∗Q∗.
See Appendix D.4 and Appendix D.5 for the full analysis. We note that diagonals of Λ do not have to be sorted in ascending order as stated in (Bao et al., 2020) because any permutation of the diagonals will be equivalent to OPT2. In addition, Proposition 2 is applicable to A3 & A4. f(·) in A3 & A4 is the reconstruction error, in which case closed-form solutions exist: Corollary 1. Let X ∈ Rm×n (m ≥ n) be a full rank matrix. Let Λ be a diagonal matrix with Λii ≥ 0 for all i. Let P ∈ Rm×k and Q ∈ Rk×n. Consider the optimization problems
min P,Q
∥X − PQ∥2F + ∥Λ 1 2Q∥2F + ∥PΛ 1 2 ∥2F (9)
min A,B
∥X −XAB∥2F + ∥ΛB∥2F + ∥XA∥2F . (10)
Let the SVD of X be UΣV T, and let Σk ∈ Rk×k be a matrix comprising the k largest singular values of X , and Uk and Vk be the corresponding singular vectors. Let λ(1) ≥ λ(2) ≥ · · · ≥ λ(k) be the sorted sequence of the diagonal values from Λ. Eq. (9) has a closed-form solution:
P ∗ = Ukdiag( √ (σ1 − λ(k))+), . . . , √ (σk − λ(1))+)Ω,
Q∗ = ΩTdiag( √ (σ1 − λ(k))+), . . . , √ (σk − λ(1))+)V Tk , (11)
where Ω is a unitary matrix that corresponds to the permutation π such that λπ(1) ≤ · · · ≤ λπ(k). In addition, Eq. (10) has a closed-form solution: A∗ = X†P ∗Λ 1 2 and B∗ = Λ− 1 2Q∗, where X† is the pseudo-inverse of X .
Rigidity of ∥W∥p∗ In (Cavazza et al., 2018), it was observed that when we use a neural net X = PQ (with learnable parameters being P and Q) to train a model and a standard dropout is used, the
objective is equivalent to solving minW ∥X−W∥2F +λ∥W∥2∗. While the regularizer ∥W∥2∗ deviates from the standard one ∥W∥∗, the optimal solution here is W = USµ(Σ)V T, where UΣV T is SVD of X , and Sµ(Σ) is a diagonal matrix such that its (i, i)-th element is (Σi,i−µ)+, in which µ depends on the data X and λ. In other words, the optimal solutions for regularizers ∥W∥2∗ and ∥W∥∗ are strikingly similar. Thus, we are interested in how regularizers with different exponents are connected. Our main observation is that for any regularizer ∥W∥p∗ (p ≥ 1), the optimal solution has the same structure. Lemma 2. Let X ∈ Rm×n, where m ≥ n. Let the d leading SVDs of X be Ud and Vd respectively. Let σ1, . . . , σn be the singular values of X . Consider the optimization problem:
min W
1 2 ∥X −W∥2F + λ∥W∥p∗. (12)
Let µk = 1k ∑ i≤k σi(X), η(µk) be the positive root of the function z + λkz p−1 − kµk and d be the largest value such that σd(X)− λ(η(µd))p−1 ≥ 0. Let µ = λ(η(µd))p−1. The optimal solution of W is Ud · dMat(σ1 − µ, σ2 − µ, . . . , σd − µ) · V Td .
Lemma 2 implies that regularizers ∥W∥p∗ with different p’s differ in how the shrinkage variable µ is obtained. Observing that µ is a hyperparameter that requires extensive tuning, all regularizers ∥W∥p∗ provide the same learning power. As noted earlier, µ is a function of the data X unless p = 1; λ needs to be rescaled when X is scaled by a constant factor unless p = 2. This implies it could be easier to tune µ when p = 1, 2, and explains why only p = 1, 2 have been extensively considered.

heading: 4 LOW-RANK FROBENIUS-NORM REGULARIZATIONS
text: This section presents two closed-form low-rank estimators with competitive performance.
Approximate low-rank DLAE and EDLAE. Recall that (full-rank) DLAE solves:
min W
||X −XW ||2F + ||Λ 1 2W ||2F s.t. Λ =
p
1− p dMat(diag(XTX)) (13)
and EDLAE with the additional diag(W ) = 0 constraint. While both the regularizations use nuclearnorm and the full-rank DLAE/EDLAE solutions all have closed-form solutions, such solution is unknown for the low-rank DLAE and EDLAE, which still require ADMM (Steck, 2020). The closed-form solutions will help both better understand and compare these models, and determine the hyperparameters such as learning rates, which is usually difficult for the ADMM type solutions.
Low-rank DLAE is a special case of Tikhonov regularization (Eq. (6)) so it has a closed-form solution, namely LR-DLAE (Jin et al., 2021). See case 10 in Table 1. However, for EDLAE, it enforces the zero-diagonal constraint, which makes the exact solution difficult to express. Recall that that low-rank EDLAE aims to solve :
W = arg min rank(W )≤k
||X −X(W − dMat(diag(W ))||2F + ||Λ 1 2 (W − dMat(diag(W ))||2F (14)
Here, we present an approximate low-rank closed-form solution of Eq. (14) by decomposing the optimization problem into two subproblems, which is similar to (Jin et al., 2021). For detailed analysis and process, please refer to Appendix B. We first consider the full-rank closed-form solution for EDLAE (Eq. (13) with zero-diagonal constraint), which is given by:
W ∗ = I − C · dMat(1⊘ diag(C)), where C = (XTX + Λ)−1
Then, we consider two (closed-form) approaches to produce low-rank matrix approximation of W ∗:
Method 1 (LR-EDLAE-1): Selecting Ŵ to best approximate the performance of W ∗ without the zero-diagonal constraint:
Ŵ = arg min rank(W )≤k
||XW ∗ −XW ||2F
= arg min rank(W )≤k
||XW ∗ −XW ||2F + ||Λ 1 2 (W ∗ −W )||2F
(15)
where X = [ X
Λ 1 2
] . Noting, in the full-rank problem Eq. (13), it forces the diagonal of derived matrix
(W ∗) to be zero. Here, we relax the zero-diagonal constraint. This relaxation corresponds to lower
bounding an approximate low-rank EDLAE objective function when rank k is sufficiently large. (Details in Appendix B). The closed-from solution of Eq. (15) is given by:
Ŵ = W ∗(QkQ T k ) (16)
where Qk takes the first k rows of matrix Q: XW ∗ SVD = PΣQT .
Method 2 (LR-EDLAE-2): SVD approximation of W ∗: The alternative solution is to simply perform SVD: W ∗ SVD= UΣV T , and which gives a low-rank estimation of W ∗: Ŵ = UkΣkV Tk . This solution corresponds to upper bounding the minimization objective in Eq. (15) (Details in Appendix B).
Both approaches relax the hard zoer-diagonal constraint on Ŵ . In the next section, the 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.

heading: 5 EXPERIMENTAL RESULTS
text: This section experiments different regularizations for linear recommendation models, with a goal to validate the efficacy of various regularizations (all can be categorized under nuclear or Frobenius norms) together with their closed-form solutions. We answer three questions: Q1. How do the closed-form solutions 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 k) and the recommendation accuracy? Q3. How does the ordering of weights from small to large (nondescending) for adjusting the singular values (Corollary 1) affect the recommendation performance?
Experimental Setup: We use three commonly used datasets for recommendation studies: MovieLens 20 Mil. (ML-20M) (Harper & Konstan, 2015), Netflix Prize (Netflix) (Bennett et al., 2007), and the Million Song Data (MSD)(Bertin-Mahieux et al., 2011). We obtained the 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 & 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 Recall@20, Recall@50 and nDCG@100. Finally, we note that our code are openly available (see Appendix E.1).
Q1: Low-Rank Frobenius-Norm and (Weighted) Nuclear-Norm Regularization: We evaluate the low-rank Frobenius-norm-based regularization and the nuclear-norm-based regularization (Eq. (1)) for the matrix factorization. Here EDLAE-ADMM, LRR, LR-DLAE, LR-EDLAE-1, LR-EDLAE-2, MF dropout and VLAE are listed in Table 1. To determine the non-descending order of weights λi for the closed-form solution in Eq. (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 λi = Cσi where C is a hyperparameter, and we perform grid-search to find the optimal one.
Table 2 shows that the weighted nuclear-norm regularization (VLAE) based matrix factorization performs worse than the constant weighted version (Regularized PCA). The latter demonstrates strong performance comparing against the WFM/ALS (one of the most popular implicit matrix factorization algorithms). The closed-form solutions (LR-DLAE, LR-EDLAE-1 and LR-EDLAE-2) all perform comparable with the ADMM-based low-rank solution and the full-rank DLAE and EDLAE solutions.
Q2: nDCG vs Rank k for low-rank Frobenius norm: In this experiment, we focus on evaluating the recommendation accuracy (using nDCG) against the rank k. Specifically, we vary k from around 1K to around 10K, and we tune and compare four different methods, including EDLAE-ADMM, LR-DLAE, LR-EDLAE-1, and LR-EDLAE-2. We have the following observations from Fig. 1: 1) As k increases, the recommendation accuracy also increases in general; however, most of them reach a plateau around similar K, 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 K grows; but when K is relatively small, ADMM method performs slightly better. But none-the-less, for most of the reasonable choices of k when low rank approximates full rank, the closed-form solution performs comparable or better.
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 3. 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 3. 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.

heading: 6 CONCLUSION AND DISCUSSION
text: 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 X’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 expressive ( 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 ( Section 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.

heading: A RELATED WORK
text: There have been extensive researches on recommendation (Aggarwal, 2016). Besides the basic user-based and item-based collaborative filtering (Deshpande & Karypis, 2004), the full-rank linear autoencoder approaches include SLIM (Ning & Karypis, 2011), HOLISM (Christakopoulou & Karypis, 2014), EASE (Steck, 2019), DLAE (Steck, 2020), whereas low-rank approaches include (Kabbur et al., 2013; Sedhain et al., 2016; Steck, 2020). All the customized recommendation models have 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 variational linear autoencoders (VLAE) to study how the deep VAE based recommendation approaches (Li & 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 VLAE inspired regularization is also never studied before.
Nuclear-norm regularizers can recover low-rank matrices in the vector regression setting (Negahban & 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).

heading: B APPROXIMATE LOW RANK EDLAE

heading: B.1 LOW-RANK CLOSED-FORM SOLUTION FOR DLAE
text: Recall that full-rank DLAE problem:
min W
||X −XW ||2F + ||Λ 1 2W ||2F s.t. Λ =
p
1− p dMat(diag(XTX)) (17)
And now consider the following basic low-rank DLAE model:
arg min rank(W )≤k
||X −XW ||2F + ||Λ 1 2W ||2F s.t. Λ =
p
1− p dMat(diag(XTX)) (18)
We leverage the strategy used in (Jin et al., 2021) to derive the low-rank closed-form solution. We formalize the low-rank problem as a standard regression problem:
Y = [ X 0 ] X = [ X Λ 1 2 ] (19)
Equiped with this, the above regression problem can be formulated as:
min rank(W )≤k
||Y −XW ||2F
= min rank(W )≤k
||Y −XW ∗||2F + ||XW ∗ −XW ||2F ,
where W ∗ = argmin ||Y −XW ∗||2F , is the optimal solution for the full-rank regression problem Eq. (17). The primary loss function ||Y −XW ||2F then can be decomposed into two parts: ||Y − XW ∗||2F (full-rank regression problem), and ||XW ∗ −XW ||2F (low-rank approximation problem).
Since W ∗ is the optimum, the vector ( Y −XW ∗) is orthogonal to XW ∗ −XW = X(W ∗ −W ) (van Wieringen, 2020), the above equation holds.
Thus the primary low-rank regression problem Eq. (18) can be broken into two subproblems:
(Subproblem 1:) full-rank closed-from solution:
W ∗ = argmin ||Y −XW ∗||2F = (XTX + Λ)−1XTX
(Subproblem 2:) low-rank approximation solution:
Ŵ = arg min rank(W )≤k
||XW ∗ −XW ||2F
= arg min rank(W )≤k
||XW ∗ −XW ||2F + ||Λ 1 2 (W ∗ −W )||2F
Denote Y ∗ = XW ∗, the best rank k approximation of Y ∗ in Frobenius norm can be derived from SVD (Eckart & Young, 1936):
Y ∗ SVD = PΣQT
Y ∗ (k) = PkΣkQ T k
where Qk takes the first k rows of matrix Q, as well as P and Σ. Since, Q are orthogonal matrices, we notice that:
Y ∗ (k) = PkΣkQ T k = PΣQ T (QkQ T k )
= XW ∗(QkQ T k ) = XW
Imediately, we have
Ŵ = W ∗(QkQ T k )
The complete low-rank closed-form solution for Eq. (18) is:
Ŵ = (XTX + Λ)−1XTX(QkQ T k )

heading: B.2 APPROXIMATE LOW-RANK CLOSED-FORM SOLUTION FOR EDLAE (LR-EDLAE-1)
text: Recall that that low-rank EDLAE aims to solve
W = arg min rank(W )≤k
||X −X(W − dMat(diag(W ))||2F + ||Λ 1 2 (W − dMat(diag(W ))||2F (20)
When the rank k of the low-rank model is sufficiently large (for instance k ≈ 1000), its objective function can be approximate as (Steck, 2020):
W = arg min rank(W )≤k
||X −XW ||2F + ||Λ 1 2W ||2F
s.t. diag(W ) = 0 (21)
Note that for most of our experimental studies, the rank k in the low-rank model is typically fairly large (over 1K), and thus the above approximation is reasonable. Compared to Eq. (18), low-rank EDLAE imposes only an additional zero-diagonal constraint. Our goal is to mimic the strategy developed in Appendix B.1 to build a two-staged algorithm to produce a low-rank approximation. We have the following key Lemma. Lemma 3. Let W ∗ be a (full rank) solution to
min W∗ ∥Y −XW ∗∥2F diag(W ∗) = 0.
where,
Y = [ X 0 ] X = [ X Λ 1 2 ] Let W be any matrix such that diag(W ) = 0. We have
⟨Y −XW ∗, XW ∗ −XW ⟩ = 0.
Proof. First, observe that W ∗ can be found by solving a sequence of independent linear regressions, each of which correpsonds to a column of Y . Let X(−1) be the matrix that removes the i-th column of X . Then the i-th column of W ∗, namely W ∗i can be found using two steps:
• Step 1. Ŵi = (X T (−i)X(−1)) −1X T (−i)Y .
• Step 2. Construct W ∗i by inserting a zero between the (i− 1)-st and i-th entries of Ŵi.
This is because enforcing W ∗i,i = 0 is the same as removing the i-th column in X to fit Yi.
Next, we move to show that for any i,
⟨Y i −XW ∗i , XW ∗i −XWi⟩ = 0.
Observing that W ∗i is the optimal solution for ∥Y i − XW ∗i ∥22, Y i − XW ∗i is orthogonal to the column space of X(−i). On the other hand, XW ∗i −XWi = X(W ∗i −Wi) and (W ∗i −Wi)i = 0. This implies that XW ∗i −XWi is in the column space of X(−i). Therefore, Y i−XW ∗i is orthogonal to XW ∗i −XWi, which directly implies the lemma.
Using Lemma 3, we see that
min rank(W )≤k
||Y −XW ||2F
s.t. diag(W ) = 0
= min rank(W )≤k
||Y −XW ∗||2F + ||XW ∗ −XW ||2F ,
s.t. diag(W ∗) = 0 diag(W ) = 0
≥ min rank(W )≤k ||Y −XW ∗||2F + ||XW ∗ −XW ||2F ,
s.t. diag(W ∗) = 0
(22)
The last inequality holds because one constraint diag(W ) = 0 is removed so the objective will only improve. With the above analysis, we are ready to describe our algorithm.
(Subproblem 1:) full-rank closed-from solution:
W ∗ = argmin ||Y −XW ∗||2F = argmin
W ||X −XW ∗||2F + ||Λ 1 2W ∗||2F
s.t. diag(W ∗) = 0
(23)
The closed-form solution is:
C = (XTX + Λ)−1
W ∗ = I − C · dMat(1⊘ diag(C))
(Subproblem 2:) low-rank matrix approximation:
Ŵ = arg min rank(W )≤k
||XW ∗ −XW ||2F
= arg min rank(W )≤k
||XW ∗ −XW ||2F + ||Λ 1 2 (W ∗ −W )||2F
(24)
Note, in the full-rank subproblem Eq. (23), it forces the diagonal of derived matrix (W ∗) to be zero. Here, we relax the zero-diagonal constraint - the diagonal of low rank approximate matrix Ŵ doesn’t have to be strictly zero.
Given this, the closed-from solution (LR-EDLAE-1) of Eq. (24) is given by:
Ŵ = W ∗(QkQ T k )
where Qk comes from: XW ∗ SVD = PΣQT .
Intuition on the efficacy of the proposed approximation. We believe that because the full-rank W ∗ performs well for recommendation accuracy and better than the low-rank solutions (Steck, 2020; Jin et al., 2021), when W ’s estimation is closer to W ∗ with respect to the final recovery of Y , XW ≈ XW ∗, it has a comparable performance against XW ∗. In other words, it may not be essential for the diagonal of low-rank approximate matrix Ŵ to be zero.

heading: B.3 SVD OF W ∗ (LR-EDLAE-2)
text: Furthermore, if we consider the inequality
min rank(W )≤k ||XW ∗ −XW ||2F ≤ min rank(W )≤k ||X||2F ||W ∗ −W ||2F ,
or a tighter one
min rank(W )≤k ||XW ∗ −XW ||2F ≤ min rank(W )≤k ||X||22||W ∗ −W ||2F
Both suggest if we simply perform SVD on W ∗, we can obtain an upper bound error between XW ∗
and XW .
Simply, the other closed-form solution Ŵ (LR-EDLAE-2) is:
W ∗ SVD = UΣV T
Ŵ = UkΣkV T k

heading: C TIME COMPLEXITY ANALYSIS
text: The running time of all linear models discussed in this paper is determined by one or more of the following factors: (i) Computation of matrix multiplication and inversion. The best known theoretical bound is O(n2.373)(Alman & Williams, 2021) but to the authors’ knowledge the algorithm is not implemented in any mainstream libraries. (ii) Computation of k leading singular vectors. The running
time is O(kn2) for an n× n matrix. But recent trends believe that the constant terms in the big-O notation matter so the asymptotic analysis becomes less fashionable (Stathopoulos & McCombs, 2010; Wu et al., 2017). (iii) Convergence rate of a gradient-based algorithm. Except for special cases, when an objective is not convex (e.g., when W is decomposed into W = PQ), the convergence rate of an algorithm is remarkably difficult to analyze and often has substantial discrepancies with its practical performance (Bhojanapalli et al., 2016).
We remark that all three factors exhibit significant gaps between theory and practice, i.e., the best theoretical algorithms are not implemented, constants in big-O matter, and accurate characterization of convergence rate is difficult. Therefore:
In theory, (i) the running time for full-rank Frobenius-norm models (cases 5. EASE, 6. DLAE, and 7. EDLAE in Table 1) is O(n2.373) because they involve merely matrix inversion and multiplication operations; (ii) the running time for most low-rank models (cases 3. WLAE, 8. LRR, 10. LR-DLAE, 11. LR-EDLAE-1, and 12. LR-EDLAE-2 in Table 1) are O(n2.373 + kn2) because they require computation of k leading singular vectors or eigenvectors; (iii) the nuclear-norm models require the multiplication with X ∈ Rm×n. So the running time for cases 1 (Regularized PCA), 2 (MF dropout), and 4 (VLAE) are O(kn2 +mnk), O(m2n+mn2), O(km2 + kn2) respectively; and (iv) The running time ADMM-based model (case 9. EDLAE-ADMM in Table 1) is O(tn2.373), where t is the number of iteration needed to converge. We summarize the results in Table 4
In practice, solving low-rank SVD and computing matrix multiplication require similar time especially after carefully using HPC techniques. So our proposed algorithms LR-EDLAE-1 and LR-EDLAE-2 (Section 4) are neither slower nor faster than existing ones that do not use gradient algorithms (e.g. full-rank EDLAE, etc.). The clock time for gradient-based (e.g. ADMM) algorithms is usually worse than that in theory because determination of hyper-parameters such as learning rate is usually a guesswork. A rule of thumb is that LR-EDLAE-1 and LR-EDLAE-2 are faster for small k and slower (than e.g., ADMM-based methods) for large k.
Finally, we would like to address again on the advantage of having any closed-form solutions: they are easy to implement and use, and can be quickly tested using any standard matrix computation platform without the specialized recommendation library. They are more robust to error and bug, and can serve as the competitive baseline for recommendation evaluation tasks.

heading: D PROOFS

heading: D.1 PROOF OF LEMMA 1
text: Consider X ∈ Rm×n. Suppose latent variables z ∈ Rk and generate data x ∈ Rn. The Variational Linear Autoencoder (VLAE) is defined in the same way as (Lucas et al., 2019b):
p(x | z) = N ( Wz + µ, σ2I ) q(z | x) = N (V (x− µ), D)
(25)
For simplification, we set µ = 0 in following context. And the ELBO of VLAE is known as:
Lx =−KL(q(z|x)||p(z)) + Eq(z|x)[log p(x|z)] KL(q(z|x)||p(z)) =− log |D|+ xTV TV x+ tr(D)− k
Eq(z|x)[log p(x|z)] =− 1
2σ2
( tr(WDWT ) + xTV TWTWV x− 2xTWV x+ xTx ) − n
2 log 2πσ2
(26)
Again, the (maximizing) ELBO can be written as:
Lx = − 1
2
( − log |D|+ xTV TV x+ tr(D)− k ) − n
2 log 2πσ2
− 1 2σ2
( tr(WDWT ) + xTV TWTWV x− 2xTWV x+ xTx ) = −1
2 ||V x||22 −
1
2σ2
( ||W √ D||2F + ||x−WV x||22 ) + f(D,σ)
(27)
where f(D,σ) = 12 log |D| − 1 2 tr(D) + k 2 − n 2 log 2πσ 2, x ∈ Rn and z ∈ Rk. For whole data, it is equivalent to minimize:
L = ||XT −WVXT ||2F +m||W √ D||2F + σ2||V XT ||2F + g(D,σ)
= ||X −XV TWT ||2F +m|| √ DWT ||2F + σ2||XV T ||2F + g(D,σ)
(28)
where g(D,σ) = −σ2m ( log |D| − tr(D) + k − n log 2πσ2).

heading: D.2 PROOF OF PROPOSITION 1
text: Note that when σi ≤ λ(k−i), the new singular value shrinks to zero, and can be removed. Basically, for any λ(1) ≥ · · · ≥ λ(k), we can build the corresponding Tikhonov regularized instance by setting
σ2i σ2i + λ ′ i = σi − λ(k−i) σi , i.e.,λ′k = σ3i σi − λ(k−i) − σ2i . (29)
Discussion of Proposition 1: Further, the same observation holds true for the regularization (Eq. (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.
D.3 “CLEAR UP” EQ. (4)
Since our objective is for recommendation (not purely on recovering the low rank factors of the data), we treat the covariance matrix D as hyperparameters. Thus, the term g(D,σ) becomes a constant, and let A = σV T , B = 1/σWT , Λ = σ √ mD.
If we consider D as an optimization parameter, then the optimal solution of VLAE is equivalent to that of pPCA (Tipping & Bishop, 1999a). In this case, D = σ2(Σ2/m)−1, where Σ2 = diag(σ2i ) are
the eigen-values of the covariance matrix of X , and σ2 = 1n−k ∑n j=k+1 σ2j m . Further, the closed-form solutions of V (A) and W (B) are characterized. However, for recommendation, the matrix D can be considered as a hyperparameter (to be tuned); in this case, the closed-form solution is not studied yet.
It is easy to see we can “clear up” Eq. (4) and obtain the following optimization problem :
min A∈Rn×k,B∈Rk×n
∥X −XAB∥2F + ∥XA∥2F + ∥ΛB∥2F , (30)
in which decision variables are A and B, and the hyperparameter is a diagonal matrix Λ ∈ Rk×k.
D.4 INTUITION FOR PROVING PROPOSITION 2
Here, We explain the intuition for proving Appendix D.5 Proposition 2 (see Appendix D.5 for the full analysis). Consider OPT1 and let W = PQ. Our goal is to characterize the behaviors of P and Q with the presence of the regularizers when W = PQ is known (fixed). Let the SVD of W be UWΣWV TW . Because two regularizers ∥Λ 1 2Q∥2F and ∥PΛ 1 2 ∥2F are symmetric, we could “guess” P = UWΣ 1 2 WΩ and Q = Ω TΣ 1 2 WV T W , where Ω is a unitary matrix. Now we have
∥Λ 12Q∥2F + ∥PΛ 1 2 ∥2F = ∥Λ 1 2ΩΣ 1 2V TW ∥2F + ∥UWΣ
1 2
WΩΛ 1 2 ∥2F = 2∥Λ 1 2ΩΣ
1 2
W ∥ 2 F .
Now the question of finding P and Q when W is known boils down to finding a unitary matrix Ω that minimizes ∥Λ 12ΩΣ 1 2
W ∥2F , where diagonal matrices Λ and ΣW are given. Recall that λi = Λii and let σi = (ΣW )ii. Note that λi’s could be unsorted and, and that σi’s are sorted in descending order.
If we restrict Ω to be only a permutation matrix, then we aim to find a permutation π ∈ [k] that minimizes ∑ i≤k λπ(i)σi. Using a rearrangement inequality (Yue, 2020), we can see that the minimal is achieved when λπ(1) ≤ λπ(2) ≤ · · · ≤ λπ(k). In this case, we indeed have minΩ a permutation ∥Λ 1 2ΩΣ 1 2 W ∥2F = ∥ΣW ∥ω,∗, where ω = (λπ(1), . . . , λπ(k)).
Note that because PQ = PΩΩTQ for any unitary matrix Ω, it is always beneficial to use Ω to shuffle the rows and columns of P and Q so that the largest σi is mapped to the smallest λi, etc. This “degree of freedom” from Ω also explains why ordering the values along Λ’s diagonal is irrelevant.
Appendix D.5 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.

heading: D.5 PROOF OF PROPOSITION 2
text: By slightly abusing the notation, we shall let OPT1 (OPT2) be the value of the optimal solution for OPT1 (OPT2). We need to show that OPT1 = OPT2. We need two directions.
OPT2 ≥ OPT1: Let W ∗ be an optimal solution for OPT2. Let the SVD of W ∗ be U∗Σ∗(V ∗)T. Recall that W ∗ needs to satisfy the rank constraint rank(W ∗) ≤ k so U∗ ∈ Rm×k, Σ∗ ∈ Rk×k, and V ∗ ∈ Rn×k. Let π be a permutation on [k] such that λπ(1) ≤ λπ(2) ≤ · · · ≤ λπ(k). Let also Ω be the corresponding permutation matrix. Specifically, Ω ∈ {0, 1}k×k and there is exactly one entry in each row of Ω is 1:
Ωi,j = { 1 if j = π(i). 0 otherwise.
For example, consider a case in which λ1 > λ2 > · · · > λk. Then we set π = (k, k − 1, . . . , 1), and correspondingly,
Ω =  0 ... 0 10 ... 1 0... 1 ... 0 0  . Next, let P = U∗(Σ∗) 1 2Ω and Q = ΩT(Σ∗) 1 2 (V ∗)T. We have W ∗ = PQ and f(W ∗) = f(PQ). In addition,
∥Λ 12Q∥2F + ∥PΛ 1 2 ∥2F = 2∥Λ 1 2Ω(Σ∗) 1 2 ∥2F = 2 ∑ i≤k λπ(i)σi = 2∥W ∗∥ω,∗,
where σi is the i-th largest singular value of W ∗. In other words, we have found a (P,Q) pair such that
f(PQ) + ∥Λ 12Q∥2F + ∥PΛ 1 2 ∥2F = f(W ∗) + 2∥W ∗∥ω,∗ = OPT2,
which shows that OPT1 ≤ OPT2.
OPT2 ≤ OPT1. Let P ∗ and Q∗ be an optimal solution for OPT1. Let the singular values of P ∗ be σ1(P ∗) ≥ σ2(P ∗) ≥ · · · ≥ σk(P ∗) and those of Q∗ be σ1(Q∗) ≥ σ2(Q∗) ≥ · · · ≥ σk(Q∗). Let also σ∗1 ≥ · · · ≥ σ∗k be the singular values of P ∗Q∗.
We shall find a lower bound of ∥Λ 12Q∥2F + ∥PΛ 1 2 ∥2F expressed in terms of σ∗i ’s. In fact, we shall show that ∥Λ 12Q∥2F + ∥PΛ 1 2 ∥2F ≥ 2∥P ∗Q∗∥ω,∗. (31)
One can see that if Eq. (31) were true, we have
OPT2 ≤ f(P ∗Q∗) + 2∥P ∗Q∗∥ω,∗ ≤ f(P ∗Q∗) + ∥Λ 1 2Q∗∥2F + ∥P ∗Λ 1 2 ∥2F = OPT1.
Thus, it remains to prove Eq. (31). Let λ(1) ≥ λ(2) ≥ · · · ≥ λ(k) be a sorted sequence of λi’s i.e., λ(k) = λπ(1), λ(k−1) = λπ(2), . . . , λ(1) = λπ(k).
First, we show that ∥PΛ 12 ∥2F ≥ ∑k i=1 λ(k−i+1)×σ2i (P ∗) and ∥Λ 1 2Q∥2F ≥ ∑k i=1 λ(k−i+1)×σ2i (Q∗). We need the following Lemma (see e.g., Theorem 2 in (Yue, 2020)):
Lemma 4. Let A and B be two positive definite matrices in Rk×k. Then it holds that
k∑ i=1 σi(A)σk−i+1(B) ≤ tr(B 1 2AB 1 2 ). (32)
Let the SVD of P ∗ be UP∗ΣP∗V TP∗ and that of Q ∗ be UQ∗ΣQ∗V TQ∗ . We have
∥P ∗Λ 12 ∥2F = ∥UP∗ΣP∗V TP∗Λ 1 2 ∥2F = ∥ΣP∗V TP∗Λ 1 2 ∥2F = tr(ΣP∗V TP∗ΛVP∗ΣP∗). (33)
We now apply Lemma 4 by setting A = V TP∗ΛVP∗ and B = Σ 2 P∗ , and obtain that
∥P ∗Λ 12 ∥2F = tr(ΣP∗V TP∗ΛVP∗ΣP∗) ≥ k∑
i=1
λ(k+1−i) × σ2i (P ∗). (34)
We may similarly prove that ∥Λ 12Q∗∥2F ≥ ∑k i=1 λ(k−i+1) × σ2i (Q∗). Therefore,
∥Λ 12Q∗∥2F + ∥P ∗Λ 1 2 ∥2F ≥ k∑ i=1 λ(k+1−i) × (σ2i (P ∗) + σ2i (Q∗)) (35)
provides a lower bound of ∥Λ 12Q∗∥2F + ∥P ∗Λ 1 2 ∥2F in terms of σi(P ∗) and σi(Q∗). We next aim to express the lower bound in terms of σ∗i ’s (singular values of P ∗Q∗) directly.
The following program gives a lower bound for ∥Λ 12Q∗∥2F + ∥P ∗Λ 1 2 ∥2F :
min : ∥Λ 12Q∗∥2F + ∥P ∗Λ 1 2 ∥2F (36)
subject to W = P ∗Q∗
σi(W ) = σ ∗ i for i ≤ k.
Write the SVD of W be UWΣWV TW . Also, let P̃ = U T WP ∗ and Q̃ = Q∗VW . Noting that the columns in P ∗ are in the column space of W and the rows in Q∗ are in the row space of W , we have (i) σi(P ∗) = σi(P̃ ) and σi(Q∗) = σi(Q̃) for i ≤ k, and (ii) ∥Λ 1 2Q∗∥2F + ∥P ∗Λ 1 2 ∥2F = ∥Λ 12 Q̃∥2F + ∥P̃Λ 1 2 ∥2F .
Therefore, Eq. (36) can be equivalently written as
min : ∥Λ 12 Q̃∥2F + ∥P̃Λ 1 2 ∥2F (37)
subject to ΣW = P̃ Q̃ (ΣW )i,i = σ ∗ i for i ≤ k.
Now P̃ Q̃ is positive definite. Using a similar technique developed in (Bao et al., 2020) (Theorem 1), one can see that P̃ = Q̃T. See also Lemma 5. This implies that P̃ = Σ 1 2
WΩ for some unitary matrix Ω and σi(P ∗) = σi(P̃ ) = σi(Q∗) = σi(Q̃) = √ σ∗i for i ≤ k. Together with Eq. (35), we have
∥Λ 12Q∗∥2F + ∥P ∗Λ 1 2 ∥2F ≥ 2 ∑ i≤k λ(k+1−i) × σ2i (P ∗) = 2 ∑ i≤k λ(k+1−i) × σ∗i = 2∥P ∗Q∗∥ω,∗.

heading: D.6 PROOF OF COROLLARY 1
text: We first find an optimal solution for Eq. (9). Let the SVD of X be X = UXΣXV TX , where UX ∈ Rm×n, ΣX ∈ Rn×n, and VX ∈ Rn×n. Let ŪX be an arbitrary basis for the subspace that is orthogonal to X’s column space so ŪX ∈ Rm×(m−n) and [UX , ŪX ] form a basis for Rm. We have
∥X − PQ∥2F + ∥PΛ 1 2 ∥2F + ∥Λ 1 2Q∥2F
= ∥∥∥∥( UTXŪTX ) XVX − ( UTX ŪTX ) PQVX ∥∥∥∥2 F + ∥∥∥∥( UTXŪTX ) PΛ 1 2 ∥∥∥∥2 F + ∥Λ 12QVX∥2F .
Let P̃ = (
UTX ŪTX
) P and Q̃ = QVX . Then our objective becomes
min P̃ ,Q̃ ∥∥∥∥( ΣX0(m−n)×n ) − P̃ Q̃ ∥∥∥∥2 F + ∥P̃Λ 12 ∥2F + ∥Λ 1 2 Q̃∥2F . (38)
Let W̃ = P̃ Q̃ and the singular values of W̃ be σ∗1 ≥ σ∗2 ≥ · · · ≥ σ∗k. Let also Σ̃ = (
ΣX 0(m−n)×n
) .
Recall also that σi is the i-th largest singular value of X . We next show that∥∥∥∥( ΣX0 ) − P̃ Q̃ ∥∥∥∥2 F = ∥Σ̃− P̃ Q̃∥2F ≥ k∑
i=1
(σi − σ∗i )2 + n∑
i=k+1
σ2i .
Note first that
∥Σ̃− W̃∥2F = ∥Σ̃∥2F + ∥W̃∥2F − 2⟨Σ̃, W̃ ⟩. (39) Next, we have ((Zheng et al., 2018)):
|⟨Σ̃, W̃ ⟩| = |tr(Σ̃W̃T)∥ ≤ |tr(Σ̃ΣW̃ )| = k∑
i=1
σiσ ∗ i .
Therefore, ⟨Σ̃, W̃ ⟩ is maximized when
W̃i,j = { σ∗i if i = j ≤ k 0 Otherwise.
When we plug in this optimized W̃ to Eq. (39), we get
∥Σ̃− W̃∥2F ≥ k∑
i=1
(σi − σ∗i )2 + n∑
i=k+1
σ2i .
Next, from Proposition 2, we have
∥P̃ V 12 ∥2F + ∥Λ 1 2 Q̃∥2F ≥ k∑ i=1 λ(k−i+1)σ ∗ i .
Therefore, we can find a lower bound for Eq. (5) in terms of σ∗i ’s:
L(σ∗1 , . . . , σ∗k) = k∑
i=1
(σi − σ∗i )2 + 2 k∑
i=1
λ(k−i+1)σ ∗ i + m∑ i=k+1 σ2i (σ ∗ 1 ≥ · · · ≥ σ∗k ≥ 0). (40)
We next find a minimal value of L (by treating σ∗i ’s as decision variables). This will give us a lower bound (and is independent of σ∗i ) on our optimization problem. We then show that this lower bound can be achieved by carefully constructing W̃ (as well as P̃ and Q̃). This means such W̃ is optimal.
Specifically, we need to find an optimal solution for the following program:
minimizeσ∗1 ,...,σ∗k L(σ ∗ 1 , . . . , σ ∗ k) (41)
subject to: σ∗i ≥ 0 σ∗1 ≤ σ∗2 ≤ · · · ≤ σ∗k (Ordering constraint)
We shall first find an optimal solution for
minimizeσ∗1 ,...,σ∗k L(σ ∗ 1 , . . . , σ ∗ k) (42)
subject to: σ∗i ≥ 0
Note here, the ordering constraint is removed so the optimal value for Eq. (42) should be no more than that for Eq. (41). We shall see that the optimal solution for Eq. (41) also satisfies the ordering constraint so indeed optimal solutions for Eq. (41) and Eq. (42) are the same.
The problem Eq. (42) boils down to finding
min σ∗i ≥0
(σi − σ∗i )2 + 2 k∑
i=1
λ(k−i+1)σ ∗ i .
We note that σ∗i ’s do not interact with each other so we can optimize each σ ∗ i ’s independently. We get σ∗i = (σi − λ(k−i+1))+. We can check that σ∗1 ≥ · · · ≥ σ∗k. Therefore, the optimal value for Eq. (41) is
k∑ i=1 (σi − (σi − λ(k−i+1))+)2 + 2 k∑ i=1 λ(k−i+1)(σi − λ(k−i+1))+ + n∑ i=k+1 σ2i .
This is also a lower bound for Eq. (9). One can check that when we set P and Q as
P ∗ = Ukdiag( √ (σ1 − λ(k))+), . . . , √ (σk − λ(1))+)Ω, (43)
Q∗ = ΩTdiag( √ (σ1 − λ(k))+), . . . , √ (σk − λ(1))+)V Tk ,
the lower bound is achieved so Eq. (43) gives an optimal solution. Here, Uk and Vk are leading left and right singular vectors of X .
Now we move to analyze Eq. (10). Our goal is to reduce Eq. (10) to Eq. (9). Let
P = XAΛ− 1 2 Q = Λ 1 2B.
Then Eq. (10) becomes
minimizeP,Q ∥X − PQ∥2F + ∥PΛ 1 2 ∥2F + ∥Λ 1 2Q∥2F (44)
subject to P = XAΛ− 1 2 (Constraint P)
Q = Λ 1 2B (Constraint Q).
Here, X and Λ are given, whereas P , Q, A, and B are decision variables. The (Constraint P) says that each column of P needs to be in a column space of X (it is a necessary and sufficient condition for A to exist). The (Constraint Q) simply says Q and B are linearly related and does not have tangible impact to the optimization problem.
But we note that when we put aside the constraints, an optimal (P,Q) is specified by Eq. (43). The columns of the optimal P indeed is in the column space of X . So (P,Q) is also an optimal solution for Eq. (44). We may find the corresponding A and B:
A∗ = X†P ∗Λ 1 2 and B∗ = Λ− 1 2Q∗,

heading: D.7 SYMMETRIC LEMMA
text: Lemma 5. Let P̃ , Q̃ ∈ Rk×k be full rank, Λ be a diagonal matrix, and ΣW be a diagonal matrix so that (ΣW )i,i = σ∗i , where σ ∗ i ’s are sorted in descending order. Consider the optimization problem:
min : ∥Λ 12 Q̃∥2F + ∥P̃Λ 1 2 ∥2F (45)
subject to ΣW = P̃ Q̃ (ΣW )i,i = σ ∗ i for i ≤ k.
There is an optimal solution such that P̃ = Q̃T
Proof. Let P̂ = P̃Λ 1 2 and Q̂ = Λ− 1 2 Q̃. The problem Eq. (45) is equivalent to
min : ∥ΛQ̂∥2F + ∥P̂∥2F (46) subject to ΣW = P̂ Q̂
(ΣW )i,i = σ ∗ i for i ≤ k.
Let the SVD of Q̂ be UQ̂ΣQ̂V T Q̂ so Q̂−1 = VQ̂Σ −1 Q̂ UT Q̂ . We can also see that P̂ = ΣW Q̂−1. Therefore, the objective term becomes
∥ΛUQ̂ΣQ̂V T Q̂ ∥2F + ∥ΣWVQ̂Σ −1 Q̂ UT Q̂ ∥2F = ∥ΛUQ̂ΣQ̂∥ 2 F + ∥ΣWVQ̂Σ −1 Q̂ ∥2F .
Let us consider the stationary points UQ̂ and VQ̂ when ΣQ̂ 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 can see Q̂ = Σ1ΣQ̂Σ2 for two permutation matrices Σ1 and Σ2. This implies that Q̃ = Λ 1 2Σ1ΣQ̂Σ2, i.e., each row (column) of Q̃ has exactly one non-zero entry. We may similarly
show that each row (column) of P̃ has exactly one non-zero entry. In addition, the locations of non-zero entries of P̃ and Q̃T are identical because P̃ Q̃ is a diagonal matrix. We may thus write
P̃ = Σ(1)ΣP̃Σ(2) Q̃ = Σ T (2)Σ(Q̃)Σ T (1),
where (ΣP̃ )i,i = σi(P̃ ) and (Σ(Q̃))i,i = στ(i)(Q̃), where τ is a permutation on [k]. The set of (possibly unsorted) singular values for P̃ Q̃ thus is σi(P̃ )στ(i)(Q̃). Thus, we can see that there exists a permutation π̄ such that
∥Λ 12 Q̃∥2F + ∥P̃Λ 1 2 ∥2F = ∑ i≤k ( σ2i (P ∗)λπ̄(i) + σ 2 τ(i)(Q ∗)λπ̄(i) )
≥ ∑ i≤k 2σi(P ∗)σ∗τ(i)(Q)λπ̄(i) ≥ 2∥P ∗Q∗∥ω,∗.
One can see that we can set P̃ = Q̃T to make all inequality becomes equality so there is an optimal solution such that P̃ = Q̃T.

heading: E EXPERIMENTAL DETAILS

heading: E.1 ANONYMOUS CODE
text: https://anonymous.4open.science/r/ICLR-2022-Anonymous-Demo-Code-B9FF/ README.md

heading: E.2 DLAE AND HYPERPARAMETER TUNING
text: This section presents the hyperparameter tuning process on the validation data over three (ML-20M, Netflix, MSD) datasets for the full-rank DLAE formula (case 6 in Table 1), which was introduced by Steck (2020) yet not investigated:
min W
||X −XW ||2F + ||Λ1/2 ·W ||2F
Λ = p
1− p dMat(diag(XTX))
Ŵ = (XTX + Λ)−1XTX
In practical, l2 regularization (hyperparameter λ) is also imposed:
Ŵ = (XTX + Λ+ λ)−1XTX
The Tables 5 to 7 show the results of nDCG@100 over three datasets respectively. And the optimal parameters are highlighted.

heading: E.3 MATRIX FACTORIZATION WITH DROPOUT AND HYPERPARAMETER TUNING
text: Cavazza et al. (2018) show that optimization with dropout (allowing rank d to be optimized) is equivalent to solving a matrix approximation problem with nuclear norm:
min P,Q,d ||X − PQ||2F + d 1− p p
· d∑
k=1
||Pk||22 · ||QTk ||22
min Y ||X − Y ||2F + 1− p p ||Y ||2∗
(47)
and the solution is given by:
X SVD = UΣV T
Y ∗ = P ∗ ·Q∗
= U · Sµ(Σ) · V T
Sµ(σ) = max(σ − µ, 0)
µ = 1− p
p+ (1− p)d̄ d̄∑ i=1 σi(X)
where d̄ denotes the largest integer such that:
σd̄(X) > 1− p
p+ (1− p)d̄ d̄∑ i=1 σi(X)
Hence, there is only one parameter p to tuning. We present the tuning process on the validation set below, see Tables 8 to 10. Optimal parameter p as well as induced rank d̄ is highlighted.

heading: E.4 RESOURCES
text: 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.

Review summaries:
This paper unifies under one theoretical framework the different state-of-the-art regularization approaches for linear collaborative-filtering recommendation models which leverage the user-item interaction data matrix.   The authors classify these algorithms among two families: 1) the nuclear-norm-based regularizers and 2) the Frobenius-norm-based regularizers. Both have their advantages and drawbacks.  While nuclear-norm-based regularizer's solutions are low-rank (and therefore more scalable) and have a closed form, their performance is limited.  On the other hand, Frobenius-norm-based regularizer's solutions are full-rank or difficult to train. Hence, focusing on generalizing Frobenius-norm-based regularizers, the authors bridge the gap between the two families by providing two new low-rank, closed-forms solutions - LR-EDLAE-1 and LR-EDLAE-2 - peforming similarly to full-rank EASE and EDLAE state-of-the-art models.   To that purpose, the authors start from the analysis work from Jin et al. (2021) about the relationship between Matrix Factorization approaches and linear autoencoders-based ones.  Besides EASE and EDLAE algorithms, they also analyze VLAE, DLAE and LRR. They show that VLAE solves a weighted nuclear-norm regularixation problem. They also generalize the result on the equivalence between dropout technique and the adding of a squared nuclear-norm regularizer to show that solution structures for ||W||^p_* are similar for all p>= 1.

This paper categorizes various linear models based on the applied regularization: nuclear norm vs Frobenius norm regularization. It is claimed that this categorization is key to understanding the prediction accuracies of the various linear model-classes. The paper also makes connections between nuclear norm regularization and Frobenius-norm regularization, and also discusses L_p norm regularization of the singular values. The paper also describes 2 methods that yield a closed-form solution of a low-rank model with Frobenius norm regularization: while it is not motivated/explained why these methods should work from a conceptual/theoretical perspective, the experiments show that they indeed work surprisingly well. 

This paper analyzed the linear models for recommendation from their regularization terms, more specifically, via nuclear norm and Frobenius norm. Based on the analyzes, this paper proposed two alternatives for linear recommendation models with closed form solutions. Experiments on three large scale datasets showed that the proposed alternatives are comparable to the state-of-the-art baseline methods.