Year: 2022

Abstract:
Recently, deep linear and nonlinear matrix factorizations gain increasing attention in the area of machine learning. Existing deep nonlinear matrix factorization methods can only exploit partial nonlinearity of the data and are not effective in handling matrices of which the number of rows is comparable to the number of columns. On the other hand, there is still a gap between deep learning and tensor decomposition. This paper presents a framework of multi-mode deep matrix and tensor factorizations to explore and exploit the full nonlinearity of the data in matrices and tensors. We use the factorization methods to solve matrix and tensor completion problems and prove that our methods have tighter generalization error bounds than conventional matrix and tensor factorization methods. The experiments on synthetic data and real datasets showed that the proposed methods have much higher recovery accuracy than many baselines.

heading: 1 INTRODUCTION
text: Low-rank matrices and tensors are pervasive in sciences and engineering. Low-rank matrix completion (LRMC) (Srebro & Shraibman, 2005; Candès & Recht, 2009; Recht, 2011; Hu et al., 2013; Hardt, 2014; Chen et al., 2014; Shamir & Shalev-Shwartz, 2014; Sun & Luo, 2016; Fan & Chow, 2017; Fan et al., 2020; Kümmerle & Sigl, 2018; Fan et al., 2019) and low-rank tensor completion (LRTC) (Gandy et al., 2011; Acar et al., 2011; Liu et al., 2012; Kressner et al., 2014; Yuan & Zhang, 2016; Foster & Risteski, 2019), as shown in Figure 1, aim to recover the missing entries of a lowrank matrix or tensor. LRMC and LRTC are very usefully in data preprocessing, image and video inpainting, and collaborative filtering (Fan & Cheng, 2018; Liu et al., 2012). For instance, in collaborative filtering (e.g. the recommendation problem of Netflix), the rating matrix given by users on items is often highly incomplete because each user can only rate a few items and each item is often rated by a few users. If the unknown entries of the rating matrix are predicted, recommendation can then be made, accordingly, for users or items. Since rating matrices often have potentially low-rank structures, we can use LRMC or even LRTC to recover the missing entries for recommendation.
Notation Throughout the paper, ‘x’, ‘x’, ‘X’, and ‘X ’ denote scalar, columns vector, matrix, and tensor respectively. ‘X(j)’ (or ‘Xj’) denotes a matrix with index j. ‖ · ‖F denotes the Frobenius norm of matrix or tensor. ‖ · ‖∗ denotes the nuclear norm of matrix, i.e. sum of the singular values. ‖ · ‖∞ denotes the maximum absolute element in a matrix or tensor. ‘C ×j A’ denotes the jmode product of a tensor C with a matrix A. PΩ(X) denotes the projection onto the set Ω, i.e., [PΩ(X)]ij = [X]ij if (i, j) ∈ Ω and [PΩ(X)]ij = 0 otherwise. |Ω| denotes the cardinality of Ω. g(·) denotes an activation function and hl denotes the number of units in layer l of a neural network.
LRMC In the past decade, a lot of algorithms with theoretical guarantee or/and empirical success have been proposed for LRMC. For example, Candès & Recht (2009) proved that a low-rank matrix can be recovered exactly with high probability from a few entries sampled uniformly at random, via nuclear norm minimization. Nie et al. (2012) proposed to minimize the Schatten-p quasi-norm for LRMC. Nuclear norm and Schatten-p quasi-norm minimizations are based on singular value decomposition and hence have high computational cost when the size of the matrix is very large. In addition, they do not incorporate the rank prior that are possibly available in practice, though could be inaccurate. In contrast, low-rank factorization based methods (Srebro & Shraibman, 2005; Srebro & Salakhutdinov, 2010; Wen et al., 2012; Hardt, 2014; Sun & Luo, 2016; Jin et al., 2016; Shang et al., 2016) are scalable to big matrices and can provide high recovery accuracy if the factorization size is properly determined. In fact, regularized matrix factorizations are closely related to nuclear norm and Schatten-p quasi-norm minimizations. For example, it is known (Srebro et al., 2005; Rennie & Srebro, 2005) that ‖X‖∗ = minAB=X ‖A‖F ‖B‖F = minAB=X 12 ( ‖A‖2F + ‖B‖2F ) . For LRMC, we can solve
minimize A,B
1 2 ‖PΩ(Y −AB)‖2F + λ 2
( ‖A‖2F + ‖B‖2F ) , (1)
where A ∈ Rm×d, B ∈ Rd×n, and 0 < d ≤ min{m,n}. Interestingly, Gunasekar et al. (2017) found and proved that gradient descent for (1) with λ = 0 and d = m = n converges to the minimum nuclear norm solution, provided that the learning rate is small enough and the initializations are close enough to the origin. Fan et al. (2019) proposed a class of rank regularizers called factor group sparse regularizer (FGSR) as the variational forms for Schatten-p quasi-norms, e.g.,
min X=AB
1 q ‖A‖q2,q + α‖B‖2,1 =(1 + 1/q)αq/(q+1)‖X‖ q/(q+1) Sq/(q+1) , (2)
where α > 0 and q ∈ {1, 12 , 1 4 , . . .}. When q = 1, one gets Schatten-1/2 quasi-norm. Later, Giampouras et al. (2020) extended the variational form to arbitrary p ∈ (0, 1].
LRTC A few LRMC methods have been directly extended to LRTC (Acar et al., 2011; Liu et al., 2012; Kressner et al., 2014). One may categorize LRTC algorithms according to the different tensor factorization models such as CP (CANDECOMP/PARAFAC) decomposition (Kolda & Bader, 2009; Jain & Oh, 2014), Tucker decomposition (Xu et al., 2013; Xie et al., 2018), tensor singular value decomposition (Kilmer et al., 2013; Zhang & Aeron, 2016), and tensor ring decomposition (Zhao et al., 2016; Huang et al., 2020). Take the CP decomposition based LRTC as an example, one can just learn a CP decomposition from the partial observations (Acar et al., 2011). Yuan & Zhang (2016) provided the sample complexity for the exact recovery of LRTC via nuclear norm minimization, though the optimization is intractable because computing the nuclear norm of a tensor is NP-hard (Hillar & Lim, 2013). Barak & Moitra (2016) and Foster & Risteski (2019) exploit the sum-ofsquares hierarchy for LRTC and provided theoretical guarantees for the recovery. Razin et al. (2021) studied the implicit regularization in tensor factorization and proved that gradient descent with small learning rate and near-zero initialization gives rise to a bias towards low tensor rank solutions. Liu et al. (2019) introduced convolutional neural network to tensor completion. The method can model the possibly complex interactions inside tensors while preserving the desired low-rank structure.
To complete a third-order (w.l.g) tensor based on the Tucker decomposition, we may solve
minimize X ,A1,A2,A3,C
1 2 ‖X − C ×1 A1 ×2 A2 ×3 A3‖2F , subject to PΩ(X ) = PΩ(Y), (3)
where we can also consider regularization or constraint on the factors A1, A2, A3, and C. For instance, Xu et al. (2013) reformulated (3) to all-mode matricization and solved the problem by alternating minimization. Xie et al. (2018) proposed to minimize the number of nonzero elements of C and the ranks of the matricizations of X , in addition to (3).
Deep Matrix Factorization Recently, deep matrix factorization (DMF) (Fan & Cheng, 2018; Arora et al., 2019) gains increasing attention in machine learning. There are two types of DMF methods: linear DMF (Trigeorgis et al., 2016; Zhao et al., 2017; Arora et al., 2019) and nonlinear DMF (Xue et al., 2017; Wang et al., 2017; Fan & Cheng, 2018). A general formulation is
X ≈ g1(A1g2(A2 · · · gL−1(AL−1AL) · · · )), (4)
where {gi}L−1l=1 are activation functions. When all gi are linear, (4) is linear DMF. Otherwise, (4) is nonlinear DMF. Arora et al. (2019) found that adding depth to a linear matrix factorization can enhance an implicit tendency towards low-rank and can provide higher recovery accuracy than depth-2 factorization, though the implicit regularization in linear DMF may not be captured using simple mathematical norms. In (Fan & Cheng, 2018), the numerical results showed that nonlinear DMF is able to recover (with high accuracy) the missing entries of a full-rank matrix of which the columns are generated by a nonlinear low-dimensional latent variable model.
Contributions of this paper This paper first provides theoretical analysis for why and when nonlinear DMF outperforms linear DMF in matrix completion. Second, this paper proposes a new matrix factorization method called two-mode nonlinear DMF. The method can exploit the full nonlinearity of the data while classical nonlinear DMF uses only partial nonlinearity of the data. Third, the paper presents a new tensor decomposition method called multi-mode nonlinear deep tensor factorization, which can be regarded as a high-order generalization of the two-mode nonlinear DMF. Finally, the paper shows that the multi-mode deep matrix and tensor factorization methods have tighter generalization error bounds than conventional factorization methods in matrix and tensor completion. Extensive numerical results verified the effectiveness of the proposed methods.

heading: 2 WHY AND WHEN DOES NONLINEAR DMF OUTPERFORM MF?
text: Nonlinear DMF methods can outperform linear MF methods empirically in many real tasks such as collaborative filtering (Xue et al., 2017; Fan & Cheng, 2018), image inpainting (Fan & Cheng, 2018), and disease modeling (Wang et al., 2017). The fundamental reason is that the data matrices in these tasks indeed have some nonlinear structures that cannot be exploited by linear MF methods. Nevertheless, in nonlinear DMF, the theoretical guarantee is very limited. In contrast, linear MF or even linear DMF have been well studied and guaranteed with sample complexity (Hardt, 2014; Sun & Luo, 2016; Jin et al., 2016), generalization bound (Srebro & Shraibman, 2005; Fan et al., 2019), or convergence (Gunasekar et al., 2017; Arora et al., 2019). In this section, we will provide some theoretical guarantee for nonlinear DMF in matrix completion. Assumption 1. Suppose Y = X + E, where Y ∈ Rn1×n2 . The columns of X are generated by
x = f(z),
where f : Rd → Rn1 denotes a nonlinear mapping, d < min(n1, n2), and z ∈ Rd denotes a latent variable. The entries of the noise matrix E are drawn from N (0, σ2).
Suppose we observed a few entries of Y randomly (sampling without replacement) and want to recover X from the incomplete Y . Consider the following method (Fan & Cheng, 2018):
minimize W1,...,WL,Z
1 2 ‖PΩ(Y −WLg(WL−1 · · · g(W1Z) · · · ))‖2F + λw 2 L∑ l=1 ‖Wl‖2F + λz 2 ‖Z‖2F , (5)
where Z ∈ Rd×n2 , Wl ∈ Rhl×hl−1 , l = 1, . . . , L, hL = n1, h0 = d, and h−1 = n2. λw and λz are regularization parameters. Here we have let all activation functions be the same without loss of generality and omitted the bias terms of the neural network for simplicity. According to the universal approximation theorems (Cybenko, 1989; Pinkus, 1999; Lu et al., 2017; Hanin & Sellke, 2017), f can be well approximated by the neural network provided that L is sufficiently large or at least one of h1, . . . , hL−1 is sufficiently large. The following theorem provides the excess risk bound for (5).
Theorem 1. Suppose Z and {Wl}Ll=1 are given by (5). Let X̂ = WLg(WL−1 · · · g(W1Z) · · · ). Suppose ‖Z‖F ∏L l=1 ‖Wl‖F ≤ κ, max(‖Y ‖∞, ‖X̂‖∞) ≤ ξ, the Lipschitz constant of g is η. Let N = n1n2. Then with probability at least 1− 2N−1, there exists a numerical constant c such that
1√ N ‖X − X̂‖F ≤ 1√ |Ω| ‖PΩ(Y − X̂)‖F + 1√ N ‖E‖F + cξ
(∑L l=0 hlhl−1 log(η L−1ξ−1κ)
|Ω|
)1/4 .
In Theorem 1, the last term of the RHS of the inequality is equivalent to
ψ1 = cξ ( |Ω|−1(n2d+ hL−1n1 + L−1∑ l=1 hlhl−1) log(η L−1ξ−1κ) )1/4 . (6)
Most major activation functions (e.g. ReLU and sigmoid) are, at worst, 1-Lipschitz with respect to `2 norm, which indicates that ηL−1 is not large. When n2 is large enough, ψ1 is linear with the latent dimension d of the columns of X when using (5). In practice, it is difficult to know the true d in advance. But the Frobenius norm regularizations on Z and W1 are able to reduce the nuclear norm or even the rank of W1Z theoretically, which means the d in formulation (5) and Theorem 1 can be larger than the d (ground truth) in Assumption 1.
Note that classical MF (depth-1) is a special case of (5) when L = 1 and g is linear. However, because of the nonlinearity in the data generating model of Assumption 1, the rank of X can be much higher than the intrinsic dimension of the data, i.e. r := rank(X) d. For example, if f is a q-order polynomial function, r can be as large as ( d+q q ) . For classical MF, the last term of the RHS of the inequality in Theorem 1 reduces to
ψ2 = cξ ( |Ω|−1(n1r + rn2) log(ξ−1κ′) )1/4 . (7)
Here κ′ may be less than the κ in (6) but they are in the log operator. Suppose all the Frobenius norms are the same (the worst case, may not happen because the norms of the factors in a deep factorization are usually smaller than the norms of the factors in a shallow factorization), denoted by β. Then κ = βL+1 and κ′ = β2, but r d. In the experiments, we observed that κ ≈ (κ′)1.3 when L = 2 and κ ≈ (κ′)1.7 when L = 3, which are much better than the worst cases. Now comparing (7) with (6), we conclude that the nonlinear DMF provides a tighter generalization bound than the classical MF when n2 is much larger than n1 and ∑L−1 l=1 hlhl−1 is not too large. The second condition is realistic because we usually let h1, . . . , hL−1 < n1 or we can just assume h1, . . . , hL−1 < r.
Consider the case that n2 is not large compared to n1, e.g. n2 = n1. In (6), suppose d = r/2, L = 3, h1 = 2d, h2 = 3d, and η = 1, we have ψ1 = cξ ( |Ω|−1(2n2r + 2r2) log(κ) )1/4 and
ψ2 = cξ ( |Ω|−1(2n2r) log(κ′) )1/4 , where ψ1 > ψ2. It means the upper bound given by the classical MF is tighter than that given by the nonlinear DMF. In general, the superiority of the nonlinear DMF over the classical MF increases when n2/n1 increases, which is consistent with the numerical results shown in Figure 5 of the appendix.

heading: 3 TWO-MODE NONLINEAR DEEP MATRIX FACTORIZATION
text: Data matrices in many problems are square or nearly square, which restricts the capability of nonlinear DMF. For instance, the difference between the width and height of an image is often small. The covariance matrices and graph adjacency matrices are square matrices. In recommendation system, the number of users and the number of items may be of the same order of magnitude. A covariance matrix or graph adjacency matrix may be constructed naturally by data with nonlinear latent structures. Similarly, one may imagine that the rating given by a user on an item is an interaction of the user’s features and item’s features and the features may have some nonlinear low-dimensional latent structures. The nonlinear DMF in (5) exploits only partial nonlinearity of the data.
Formally, we make the following assumption.
Assumption 2. Suppose Y = X + E, where Y ∈ Rn1×n2 and X = U>QV . The columns of U ∈ Rm1×n1 and V ∈ Rm2×n2 are generated respectively by
u = f1(z) and v = f2(s),
where f1 : Rd1 → Rm1 and f2 : Rd2 → Rm2 are nonlinear mappings, d1 < min(m1, n1), d2 < min(m2, n2), and z ∈ Rd1 and s ∈ Rd2 denote latent variables. The entries of Q ∈ Rm1×m2 and E ∈ Rn1×n2 are drawn from N (0, σ2Q) and N (0, σ2E) respectively.
Suppose we have a matrix Y given by Assumption 2, we can factorize it as
Y ≈ Û>Q̂V̂ , ûi1 = fθ1(ẑi1), v̂i2 = fθ2(ŝi2), i1 = 1, . . . , n1, i2 = 1, . . . , n2. (8)
In (8), fθ1 : Rd1 → Rm1 and fθ2 : Rd2 → Rm2 are multi-layer neural networks with parameter sets θ1 = {W (1)1 , . . . ,W (1) L1 } and θ2 = {W (2)1 , . . . ,W (2) L2 }, where W (j)i denotes the weight matrix in the i-th layer of neural network j, and Lj denotes the number of layers of neural network j, j = 1, 2.
We call (8) two-mode or two-way nonlinear deep matrix factorization. Figure 2 shows an intuitive example of (8) with depth-three factorizations in both directions.
Suppose we observed a few entries of Y randomly (sampling without replacement). Based on Assumption 2, we propose the following matrix completion method
minimize Z,Q,S,θ1,θ2
1 2 ‖PΩ(Y − fθ1(Z)>Qfθ2(S))‖2F +R(Z,Q,S, θ1, θ2), (9)
where fθ1(Z) = g(W (1) L1 g(W (1) L1−1 · · · g(W (1) 1 Z) · · · )),
fθ2(S) = g(W (2) L2 g(W (2) L2−1 · · · g(W (2) 1 S) · · · )),
and Z ∈ Rd1×n1 , S ∈ Rd2×n2 , W (j)l ∈ R h (j) l ×h (j) l−1 , l = 1, . . . , Lj , h (j) Lj = mj , h (j) 0 = dj , h (j) −1 = nj , j = 1, 2, and Q ∈ Rm1×m2 . We have let all activation functions be the same without loss of generality. In (9), R(·) is a regularization operator to reduce overfitting. For instance, R(Z,Q,S, θ1, θ2) = λ1 2 (‖Z‖2F + ‖S‖2F ) + λ2 2 ( ‖Q‖2F + ∑ W∈θ1∪θ2 ‖W ‖ 2 F ) . Compared to (5), (9) exploits the full nonlinearity of the data. The following theorem provides the excess risk for (9).
Theorem 2. Suppose Z, S, {W (1)l } L1 l=1, {W (2) l } L2 l=1, and Q are given by (9). Let X̂ = fθ1(Z) >Qfθ2(S). Suppose ‖Q‖F ‖Z‖F ‖S‖F ∏2 j=1 ∏Lj l=1 ‖W (j) l ‖F ≤ κ, max(‖Y ‖∞, ‖X̂‖∞) ≤ ξ, and the Lipschitz constant of g is η. Let N = n1n2. Then with probability at least 1− 2N−1, there exists a numerical constant c such that
1√ N ‖X − X̂‖F ≤ 1√ |Ω| ‖PΩ(Y − X̂)‖F + 1√ N ‖E‖F
+ cξ
 ( m1m2 + ∑2 j=1 ∑Lj l=0 h (j) l h (j) l−1 ) log(ηL1+L2ξ−1κ)
|Ω| 1/4 . In Theorem 2, the last term of the RHS of the inequality can be written as
ψ3 = cξ ( |Ω|−1(n1d1 + n2d2 +m1m2 + 2∑ j=1 Lj∑ l=1 h (j) l h (j) l−1) log(η L−1ξ−1κ) )1/4 . (10)
When the hL−1 in (6) is much larger than the d1 in (10), we have ψ3 < ψ1. It means the two-mode nonlinear DMF (9) provides a tighter generalization bound than classical nonlinear DMF (5). It is worth mentioning that one may formulate DMF in a two-sided form, e.g., X ≈ g ( W
(2) L g(W (1) L ZW (1) R )W (2) R
) or X ≈ g ( W
(1) L g(ZLZR)W (1) R
) equivalently. However, this
model has two shortcomings compared to the one-sided DMF (4) and two-mode DMF (8). First, (4) and (8) are more compact (in terms of the number of parameters) than the two-sided DMF. For example, assume the sizes of X , ZL, and ZR are 500×500, 100×10, and 10×100 respectively, then the number of parameters in the two-sided DMF is larger than 100,000. In the one-sided DMF (4), let X = g(W2g(W1Z)) and let the sizes of W1, W2, and Z be 500×100, 100×10, and 10×500 respectively, then the number of parameters is 56,000, which is much smaller than that of the two-sided DMF. The other shortcoming is that the two-sided DMF is difficult to interpret in real applications. In contrast, (4) and (8) are much easier to interpret and more realistic.

heading: 4 MULTI-MODE NONLINEAR DEEP TENSOR FACTORIZATION
text: In this section, we extend the two-mode deep matrix factorization to multi-mode nonlinear deep tensor factorization. First, we make the following assumption. Assumption 3. Suppose Y = X +E , where Y ∈ Rn1×n2···×nk and X = C×1U1×2U2 · · ·×kUk. For j = 1, . . . , k, the rows of Uj ∈ Rnj×mj are generated by u(j) = fj(z(j)), where fj : Rdj → Rmj is a nonlinear mapping, dj < mj , and z(j) ∈ Rdj are latent variables. The entries of core tensor C ∈ Rm1×m2···×mk and noise tensor E ∈ Rn1×n2···×nk are drawn from N (0, σ2C) and N (0, σ2E) respectively.
Although {fj}kj=1, C, and {Z(j)}kj=1 are unknown, we can take advantage of neural networks and tensor decomposition to model the data, i.e.,
Y ≈ Ĉ ×1 Û1 ×2 Û2 · · · ×k Ûk, Ûj = fθj (Ẑ(j)), j = 1, . . . , k, (11)
where fθj : Rdj → Rdj×mj is a neural network with parameter set θj . Figure 3 shows an intuitive example of the proposed factorization method for a third-order tensor. Different from Tucker decomposition, we perform nonlinear DMF on the factor matrices and do not constrain the norms of the columns of the factor matrices. This factorization method can take advantage of the nonlinear structures of the tensor fully.
To recover X from a few entries of Y , we propose to solve the following problem
minimize {θj ,Z(j)}kj=1,C
1 2 ‖PΩ(Y−C×1fθ1(Z(1))×2fθ2(Z(2)) · · ·×k fθk(Z(k))‖2F +R({θj ,Z(j)}kj=1,C),
(12) where fθj (Z (j)) = g(W (j) Lj g(W (j) Lj−1 · · · g(W (j) 1 Z (j)) · · · )) and θj = {W (j)1 , . . . ,W (j) Lj }, j = 1, . . . , k. h denotes the regularization on the parameters, e.g.,
R({θj ,Z(j)}kj=1,C) = λ1 2 k∑ j=1 ‖Z(j)‖2F + λ2 2 ( ‖C‖2F + k∑ j=1 Lj∑ l=1 ‖W (j)l ‖ 2 F ) . (13)
The following theorem shows the excess risk bound for (12).
Theorem 3. Suppose C and {θj ,Z(j)}kj=1 are given by (12). Let X̂ = C ×1 fθ1(Z(1)) ×2 fθ2(Z (2)) · · · ×k fθk(Z(k)). Suppose ‖C‖F ≤ βC , ‖W (j) l ‖F ≤ β (j) l , l = 1, . . . , Lj , ‖Z(j)‖F ≤ β (j) 0 , j = 1, . . . , k, max(‖Y‖∞, ‖X̂‖∞) ≤ ξ, and the Lipschitz constant of g is η. Let N = ∏k j=1 nj . Then with probability at least 1 − 2N−1, there exists a numerical constant c such that 1√ N ‖X − X̂‖F ≤ 1√ |Ω| ‖PΩ(Y − X̂ )‖F + 1√ N ‖E‖F
+ cξ
 (∏k j=1mj + ∑k j=1 ∑Lj l=0 h (j) l h (j) l−1 ) log ( η ∑k j=1 Ljξ−1βC ∏k j=1 ∏Lj l=0 β (j) l ) |Ω| 1/4 . Remark 1. One can obtain slightly tighter bounds in Theorems 1, 2, and 3 when using spectral norms instead of Frobenius norms. But in (5), (9), and (12), the regularizations only control the Frobenius norms explicitly. In addition, (12) and Theorem 3 are also applicable to CP-like decomposition based tensor completion if C is super-diagonal.
We see that Theorem 2 is a special case of Theorem 3 when k = 2. Denote ψ4 the last term of the RHS of the inequality in Theorem 3. Since h(j)0 = dj and h (j) −1 = nj , we have
ψ4 = cξ |Ω|−1( k∏ j=1 mj + k∑ j=1 njdj + k∑ j=1 Lj∑ l=1 h (j) l h (j) l−1 ) log ( η ∑k j=1 Ljξ−1βC k∏ j=1 Lj∏ l=0 β (j) l )1/4 . (14) If g is linear and L = 0, (12) degenerates into the Tucker decomposition based LRTC, e.g. (3). Accordingly, ψ4 becomes
ψ5 = cξ |Ω|−1( k∏ j=1 mj + k∑ j=1 njmj ) log ( ξ−1βC k∏ j=1 β (j) 0 )1/4 . (15) Since mj dj , we have ψ5 > ψ4. It means the nonlinear method (12) has a tighter generalization error bound than linear tensor completion methods. Therefore, we expect that (12) outperforms linear LRTC methods when the data have nonlinear structures. For convenience, in this paper, we call (9) and (12) multi-mode deep matrix and tensor factorizations (M2DMTF).
Optimization The optimization for (9) and (12) are non-trivial because there are high-order tensors and neural networks with unknown inputs and incomplete outputs. We tried alternating minimization, nonlinear conjugate descent, L-BFGS (Liu & Nocedal, 1989), iRprop+ (Igel & Hüsken, 2000), and Adam (Kingma & Ba, 2014), where all variables are randomly initialized from Gaussian distribution. It is found that iRProp+ outperforms other algorithms especially on real datasets. The details of optimization and the analysis of time and space complexity are presented in Appendix A.

heading: 5 NUMERICAL RESULTS

heading: 5.1 SYNTHETIC DATA
text: We use the following nonlinear model to generate matrices:
U = W (1) 3 σ(W (1) 2 sin(W (1) 1 Z)), V = W (2) 3 exp(W (2) 2 σ(W (2) 1 S)), X = U >V , (16)
where Z ∈ R2×100, W (1)1 ∈ R10×2, W (1) 2 ∈ R20×10, W (1) 3 ∈ R30×20, S ∈ R2×100, W (2) 1 ∈ R10×2, W (2)2 ∈ R20×10, and W (2) 3 ∈ R30×20 are drawn from U(−1, 1). σ denotes the sigmoid fucntion. The exact rank of X is 20 but the smallest 10 nonzero singular values are close to zero. Such a synthetic matrix is realistic because the nonlinearity in real data may not be very strong and the matrices are often low-rank. Similarly, the following nonlinear model is used to generate third-order tensors:
U1 = σ(W (1) 2 sin(W (1) 1 Z (1))), U2 = exp(W (2) 2 σ(W (2) 1 Z (2))),
U3 = cos(W (3) 2 exp(W (3) 1 Z (3))), X = C ×1 U1 ×2 U2 ×3 U3, (17)
where Z(1) ∈ R2×50, W (1)1 ∈ R10×2, W (1) 2 ∈ R20×10, Z(2) ∈ R2×50, W (2) 1 ∈ R10×2, W (2) 2 ∈ R20×10, Z(3) ∈ R2×50, W (3)1 ∈ R10×2, and W (3) 2 ∈ R20×10 are drawn from U(−1, 1) and C ∈ R20×20×20 are drawn from N (0, 1). The exact n-rank of X is (20, 20, 20) but it can be well approximated by a tensor with n-rank (10, 10, 10) because here the nonlinearity is not very strong.
First, we compare the proposed method with linear matrix factorization (MF (1)), nuclear norm minimization (NNM) (Candès & Recht, 2009), TNNR (Hu et al., 2013), HM-IRLS (Kümmerle & Sigl, 2018), FGSR (Fan et al., 2019), (stacked) Autoencoders (SAE) (Sedhain et al., 2015), and DMF (Fan & Cheng, 2018). We use the tanh activation function in SAE, DMF, and M2DMTF. We determine the hyper parameters of all methods via cross-validation to ensure they perform as well as possible. Details about the optimizations and parameter settings are in Appendix B. We use the relative recovery error (Fan et al., 2019) ‖PΩ̄(X̂ −X)‖F /‖PΩ̄(X)‖F to evaluate the performance of matrix completion, where Ω̄ consists of the locations of the missing entries. Figure 4(a) shows the average performance of 20 repeated trials when the missing rate increases from 0.8 to 0.99. MF
and NNM have much higher recovery error than FGSR. This result is consistent with the results of (Fan et al., 2019; Arora et al., 2019). Although SAE is a nonlinear method, it still has high recovery error when the missing rate is less than 0.9. The reason is that filling the missing entries of the input of the encoder with zero introduces relatively large bias, which may not be eliminated effectively by the hidden layers. DMF and FGSR have quite similar performance, which is consistent with the analysis at the end of Section 2. Our M2DMTF outperformed other methods significantly especially when the missing rate is high. The reason is that M2DMTF is able to exploit the full nonlinearity of the data and has tighter generalization bound than other methods.
Figure 4(b) shows the performance (average of 20 trials) of FaLRTC (Liu et al., 2012), TenALS (Jain & Oh, 2014), TMac (Xu et al., 2013), KBR-TC (Xie et al., 2018), TRLRF (Yuan et al., 2019), CoSTCo (Liu et al., 2019), OITNN-O (Wang et al., 2020), and M2DMTF. Since the rank of the tensor is quite low, TMac, KBR-TC, and M2DMTF have high recovery accuracy when the missing rate is less than 0.96. Note that TMac and KBR-TC slightly outperformed M2DMTF when the missing rate is less than 0.96. The reason is that the optimization of M2DMTF is more difficult. In Figure 4(b), when the missing rate is higher than 0.97, the proposed M2DMTF outperformed all methods largely. The optimization curves of M2DMTF are shown in Figure 7 of Appendix B.

heading: 5.2 REAL MATRICES
text: We consider two benchmark datasets: MovieLens-100k and MovieLens-1M. They consist of the ratings (1 to 5) for movies given by users. The sizes (movies×users) of the rating matrices in the datasets are 1682 × 943 and 3706 × 6040 respectively. For each matrix (highly incomplete), we randomly split the known entries to a training set and a test set. The relative recovery errors (average of 10 repeated trials) on the two datasets are reported in Tables 1 and 2 respectively. FGSR, SAE, and M2DMTF outperformed MF and DMF in all cases. When the training-test ratio is 9/1 or 7/3, SAE has less recovery error than M2DMTF. This is quite different from the results on the synthetic data in Section 5.1. Currently, it is unknown why SAE has such a good performance on the MovieLens data. One possible reason is that the MovieLens data are quite noisy and SAE can successfully reduce the influence of the noise when the missing rate is relatively low. When the training-test ratio is smaller, i.e., 5/5, 3/7, or 1/9, M2DMTF outperformed SAE.

heading: 5.3 REAL TENSORS
text: We compare the proposed method M2DMTF with the baselines on the following datasets: Amino acid fluorescence (Bro, 1997) (5 × 201 × 61), Flow injection (Nørgaard & Ridder, 1994) (12 × 100× 89), and SW-NIR kinetic data (Bijlsma & Smilde, 2000) (301× 241× 8). More details about the experimental settings are in Appendix B.2. The relative recovery error (average of 10 trials) are reported in Table 3. The results of FaLRTC and CoSTCo are not reported because they have much higher recovery error than other methods in most cases and we want to fit the table to the width of the page. On Amino, our M2DMTF outperformed other methods significantly in all cases. On Flow, when the missing rate is 0.95 or 0.97, M2DMTF has much lower recovery error than other methods. On SW-NIR, M2DMTF has low recovery error consistently while most baselines failed totally.
Comparative studies were also conducted on three larger datasets. The results are reported in Appendix B.2.4. The experiments of tensor completion based recommendation system are in Section B.2.5. The visualization of the learned tensor and factors of M2DMTF are in Appendix B.2.6.

heading: 6 CONCLUSION
text: This paper has proposed a framework of multi-mode deep matrix and tensor factorizations. The proposed factorization methods can explore and exploit the full nonlinearity of the data effectively. We applied the factorization methods to matrix and tensor completion and provided the generalization error bounds, which verified the superiority of the proposed methods over conventional linear and nonlinear matrix factorization methods and linear tensor completion methods. The proposed framework reduced the gap between tensor decomposition and deep learning and showed that their combinations can provide better performance in tensor recovery. On some datasets, when the missing rates were low, our methods did not achieve the best performance. The reason is that the nonlinearity in the data is not very strong and the problem is not difficult for the baselines methods of which the optimizations are easier than our methods. Future work may focus on improving the optimization efficiency of M2DMTF.

heading: ACKNOWLEDGEMENTS
text: The work was supported by the research funding T00120210002 of Shenzhen Research Institute of Big Data and the Youth program 62106211 of National Natural Science Foundation of China.

heading: A OPTIMIZATION DETAILS AND COMPUTATIONAL COMPLEXITY
text: In our M2DMTF, the matrix model is a special case of the tensor model. Therefore, we only show the optimization detail and the computational complexity of the tensor case. The gradient-based optimization is shown in Algorithm 1. We suggest using iRprop+ as the optimizer in the algorithm because we found empirically that iRprop+ is more stable and efficient than other optimizers such as Adam in solving our M2DMTF.
Algorithm 1 Gradient-based optimization for M2DMTF (12) Input: Y ∈ Rn1×···×nk , Ω, λ1, λ2, g, network structures { mj , Lj , {h(j)l } Lj l=1 }k j=1 , T .
1: Initialization: C0, { Z (j) 0 , {W (j) l0 } Lj l=1 , θ (j) 0 }k j=1 , { U (j) 0 = fθ(j)0 (Z (j) 0 ) }k j=1
. 2: for t = 1, 2, . . . , T do 3: for j = 1, 2, . . . , k do{
∇θ(j) ,∇Z(j) } backprop.←−−−−− Lj , where Lj = 1
2 ∥∥∥PΩ(j) (Y(j) −U (j)t [Ct](j)V >t )∥∥∥2 F +
R({θ(j)t ,Z(j)}kj=1,Ct) and Vt = U (j−1) t ⊗ · · ·U (1) 1 ⊗U (k) t · · · ⊗U (j+1) t .
4: end for 5: ∇C ←− ∂Lk/∂C 6: { {θ(j)t ,Z (j) t }kj=1,Ct } ←− gradient-based optimizer ( {∇θ(j) ,∇Z(j)}kj=1,∇C ) 7: end for
Output: Ŷ = CT ×1 U (1)T ×2 U (2) T · · · ×k U (k) T .
Without loss of generality, we here assume that in (12), n1 = · · ·nk = n, m1 = · · ·mk = m, L1 = · · · = Lk = L, and h(1)l = · · · = h (k) l = hl, l = 0, 1, . . . , L. The time cost in each iteration of the gradient-based optimization is
O ( kmin{nmk + nkmk−1, nkm+ nk−1mk}+ kn
L∑ l=1 hl−1hl
) ,
in which the first term is from the Tucker-like decomposition while the second term is from the forward and backward propagations of the neural networks. Since the second term is much smaller than the first one and m < n, the time complexity per iteration is about O ( knkm+ knk−1mk ) .
As we need to store the observed elements of Y , core tensor C, latent factors Z(j), network parameter, their gradients, and the values of hidden units, the space complexity is O ( |Ω| + mk +
k ∑L l=1 hl−1hl + kn ∑L l=1 hl ) , which is further simplified to O ( |Ω|+mk + kn ∑L l=1 hl ) because n is larger than all hl. Now we see that the time (per iteration in the optimization) and space complexity of our M2DMTF are slightly higher than the Tucker decomposition based tensor completion methods. However, M2DMTF requires much more iterations than the baseline methods because it is composed of tensor decomposition and neural networks.

heading: B MORE ABOUT THE EXPERIMENTS

heading: B.1 MATRIX COMPLETION
text: In the proposed method (9), the regularization is
R(Z,Q,S, θ1, θ2) = λ1 2 (‖Z‖2F + ‖S‖2F ) + λ2 2
( ‖Q‖2F + 2∑ j=1 Lj∑ l=1 ‖W (j)l ‖ 2 F ) .
We let λ2 = λ′2 √ n1n2 and determine λ1 and λ′2 instead. In the optimization, S, Z, and {W (j) l } are initialized randomly, Q is initialized to an identity matrix.

heading: B.1.1 PARAMETER SETTING ON THE SYNTHETIC DATA
text: • In MF (problem (1) in the main paper), the factorization dimension d is 5 because it outperforms other choices. The λ is chosen from {0.01, 0.1, 1} and the optimizer is iRprop+. The maximum iteration is 2000.
• In FGSR (Fan et al., 2019), we use the variational form of Schatten-1/2 norm, set d = 20, and choose λ from {0.01, 0.1}. The maximum iteration is 2000.
• In SAE (Sedhain et al., 2015), the network structure is [100, 20, 10, 5, 10, 20, 100]. The weight decay parameter is chosen from {0.1, 0.5}.
• In DMF (Fan & Cheng, 2018), the network structure is [3, 10, 20, 100]. The weight decay parameters are chosen from {0.01, 0.1}.
• In M2DMTF, L = 2, d1 = d2 = 3, h (1) 1 = h (2) 1 = 10, m1 = m2 = 20, and λ1 = λ ′ 2 = 1.
In SAE, DMF, and M2DMTF, the activation function is the hyperbolic tangent function. The optimizer is iRprop+ and the maximum iteration is 3000.

heading: B.1.2 MORE RESULTS ON THE SYNTHETIC DATA
text: Figure 5 compares DMF with linear MF (FGSR) in the cases of different n2/n1. We see that when n2/n1 is too small, FGSR outperforms DMF. When n2/n1 increases, the superiority of DMF over FGSR becomes more significant. These empirical results are consistent with our theoretical analysis in Section 2 of the main paper.
Figure 6 shows the performance of M2DMTF with different number of hidden layers and different activation function. it can be found that tanh and sigmoid functions lead to lower recovery error than ReLU in all cases. The main reason is that the data were generated by smooth functions that can be
well approximated by tanh and sigmoid functions. When the tanh function or sigmoid function is used, increasing the depth of the neural networks can reduce the recovery error. But when L is too large, i.e. L = 4, the recovery error increased slightly, which might be caused by overfitting.
Figure 7 shows the values of the objective functions and relative recovery errors against the iterations in matrix completion and tensor completion when the missing rate is 0.95. The recovery error decreases slightly when the number of iterations is very large. In general, the optimization efficiency is not very high because M2DMTF is composed of multiple neural networks with unknown inputs and incomplete outputs.
We replace the sin, cos, exp, and σ in (16) and (17) with tanh such that the nonlinearity in the data generating models are consistent with the activation functions used in our matrix and tensor completion method M2DMTF. The average results of 10 trials are reported in Figure 8. We see that M2DMTF outperformed other methods significantly when the missing rate is high enough.

heading: B.1.3 PARAMETER SETTING ON THE MOVIELENS DATASETS
text: • In MF (problem (1) in the main paper), we set d = 10 and λ = 10. In stead of iRprop+, we use alternating minimization because it outperformed iRprop+ in this case. The maximum iteration is 2000.
• In FGSR Fan et al. (2019), we use the variational form of Schatten-1/2 norm, set d = 100, and set λ = 0.015 or 0.007 (for MovieLens-100k or MovieLens-1M). The maximum iteration is 2000.
• In SAE Sedhain et al. (2015), we use a single hidden layer of size 500. The weight decay parameter is 0.5.
• In DMF, Fan & Cheng (2018), the network structure is [10 or 20, 50, 100, 1682]. The weight decay parameters are chosen from {0.05, 0.1}.
• In M2DMTF, , we set L = 2, d = 10 or 20, h(1)1 = h (2) 1 = 50, m1 = m2 = 100, and
λ1 = λ ′ 2 = 1.

heading: B.1.4 MORE RESULTS ON THE MOVIELENS DATASETS
text: The RMSE (average of 10 trials) on the MovieLens-100k are reported in Table 4. The proposed M2DMTF outperformed MF, FGSR and DMF in all cases and outperformed SAE when the missing rate was high. Note that RMSE is a monotonic transformation of the relative recovery error used in the main paper.

heading: B.2 TENSOR COMPLETION
text: In the proposed method (12), the regularization is
R({θj ,Z(j)}kj=1,C) = λ1 2 k∑ j=1 ‖Z(j)‖2F + λ2 2 ( ‖C‖2F + k∑ j=1 Lj∑ l=1 ‖W (j)l ‖ 2 F ) . (18)
We let λ2 = λ′2 (∏k j=1 nj )1/k and determine λ1 and λ′2 instead. In the optimization, Uj is initialized by the singular value decomposition of the j-mode matricization of Y , C is initialized to zero, and {Z(j)}, {W (j)l } are initialized randomly.

heading: B.2.1 PARAMETER SETTING ON THE SYNTHETIC DATA
text: • In FaLRTC (Liu et al., 2012), we set α1 = α2 = α3 = 13 choose µ from {0.5, 1, 2, 5, 10}. • In TenALS (Jain & Oh, 2014), the initial rank is 5 or 10 and the max iteration is 200.
• In TMac (Xu et al., 2013), the rank is initialized to [10, 10, 10] or [5, 5, 5] and adjusted adaptively. The maximum iteration is 500.
• In KBR-TC (Xie et al., 2018), ρ = 1.01 and λ = 0.01. The maximum iteration is 300.
• TRLRF (Yuan et al., 2019), the rank is set to [5, 5, 5] or [10, 10, 10]. We choose λ from {0.01, 0.1, 1, 10, 50}.
• In CoSTCo (Liu et al., 2019), the learning rate is chosen from {0.01, 0.001, 0.001}, the rank is chosen from {2, 3, 5, 10, 20}, the epoch number is 50 or 200, and the batch size is 256× 2 or 256× 5.
• In OITNN-O (Wang et al., 2020), we choose α from {0.0001, 0.001, 0.01, 0.1}. The max iteration is 500.
• In M2DMTF, L = 2, d1 = d2 = d3 = 3, h (1) 1 = h (2) 1 = h (3) 1 = 10, m1 = m2 = m3 = 20,
and λ1 = λ′2 = 1. The optimizer is iRprop+ and the maximum iteration is 3000.

heading: B.2.2 TIME COST COMPARISON ON THE SYNTHETIC DATA
text: We report the average time costs (second) and the recovery errors of all methods in Table 5, which also includes the performance of M2DMTF with different iteration numbers. We see that our method with 3000 iterations is slower than other methods. However, when we use less iterations (e.g. 1000 or 200), the time cost of our method is comparable with or even lower than the baselines and its recovery error is still much lower than the baselines.

heading: B.2.3 EXPERIMENTAL SETTINGS ON REAL TENSORS
text: We normalize each tensor by X ← X/‖X‖∞. In M2DMTF, we set L = 1 and λ1 = λ′2 = 0.001. For i = 1, 2, 3, we set mi ≈ ni5 , di ≈ ni 10 if ni > 10 and set mi ≈ 0.8ni, di ≈ 0.7ni if ni ≤ 10. The parameters of other methods are carefully determined to provide their best performance.

heading: B.2.4 MORE RESULTS ON TENSOR DATASETS
text: In addition to the three datasets considered in the main paper, we test the proposed method on three more datasets: Porphyrin (Wold et al., 2006) 170 × 274 × 35, Sugar (Bro, 1999) 256 × 571 × 7, and Bonnie (Acar et al., 2008) 89× 97× 549. To reduce the computational costs of all methods on Bonnie, we only consider a subset of size 89 × 97 × 100. The results are reported in Table 6. We see that M2DMTF outperforms other methods when the missing rate is high enough.

heading: B.2.5 TENSOR-COMPLETION BASED RECOMMENDATION AS A DOWNSTREAM TASK
text: We split the MovieLens-100k dataset by time to form a tensor of size 943 × 1682 × 31, where 31 denotes the number of weeks as a dimension of time. Since the baselines KBR-TC, TRLRF, OITNN-O, and the proposed M2DMTF have high time cost on large tensors, we reduce the size of the tensor by removing the users giving 49 or less ratings and the movies rated by 49 or less users. Then we get a tensor of size 563 × 593 × 31. The relative recovery error (average of 5 trials) are reported in Table 7. Our M2DMTF outperformed other methods. In addition, comparing Table 7 with Table 1, we found that tensor completion methods are more effective than matrix completion methods in recommendation.
B.2.6 VISUALIZATION OF FACTORS LEARNED BY M2DMTF
Here we take the ‘Flow injection’ tensor (Nørgaard & Ridder, 1994) as an example. The data is from measuring samples with three different analytes on a flow injection analysis (FIA) system where a pH-gradient is imposed. The three analytes in twelve measured samples are 2-, 3-, and 4-hydroxy-benzaldehyde. They have different absorption spectra depending on whether they are in their acidic or basic form. The data set is a 12 (samples) × 100 (wavelengths) × 89 (times) array. Figure 9 presents the original data, incomplete data, and recovered data for the fifth horizontal slice (sample-5) of the ‘Flow’ tensor, where the missing rate is 0.97. We see that FaLRTC, TRLRF, and CoSTCo failed totally. TenALS, TMac, and KBR-TC have high recovery error visually. In contrast, the data recovered by our M2DMTF is very close to the original data. More importantly, the result of M2DMTF indicates that there are three analytes, while other methods failed to reveal the fact.
We analyze the characteristics of the low-dimensional factors learned by our M2DMTF. Figure 10 presents the values of U2 ∈ R100×20 and Z(2) ∈ R100×10 as well as the correlation coefficients
and mutual information. In the first column, we see that the values of U2 are more chaotic than the values Z(2). The reason is that the 20 variables in U2 are correlated mutually. In the second column of Figure 10, we see that the autocorrelation of the samples (rows) of U2 is similar to that of Z(2), because U2 relies on Z(2). The third column of the figure indicates that the correlation between the variables in U2 are higher than that in Z(2). The last column compares the mutual information that can quantify nonlinear relationship. Now comparing column 4 with column 3, we conclude that there exist significant nonlinear relationship between the variables in U2 and it is necessary to represent U2 by a lower-dimensional latent factor Z(2) via deep factorization. These results verified that the nonlinear latent factor assumption of our M2DMTF is practical and useful.

heading: C PROOF OF THEOREMS
text: First of all, we provide the following lemma.
Lemma 1. Let S be a set defined over tensors of size n1 × n2 × · · · × nk. Let |S| be the - covering number of S with respect to the Frobenius norm. Let N = ∏k i=1 ni. Suppose X ∈ S and max{‖Y‖∞, ‖X‖∞} ≤ ξ. Then with probability at least 1− 2N−1, there exists a constant c such
that
sup X∈S ∣∣∣∣∣ 1√N ‖Y −X‖F − 1√|Ω| ‖PΩ(Y −X )‖F ∣∣∣∣∣ ≤ 2 √|Ω| + ( 8ξ4 log (|S|N) |Ω| )1/4 .
The proof of the lemma is in Appendix D.1.

heading: C.1 PROOF FOR THEOREM 1
text: The following lemma (proved in Appendix 2) shows an upper bound for the covering number of the deep factorization matrix set. Lemma 2. Let S = {X ∈ Rm×n : X = WLg(WL−1 · · · g(W1Z)), ‖Z‖F ≤ β0, ‖Wl‖F ≤ βl, 1 ≤ l ≤ L}, where Z ∈ Rh0×n, Wl ∈ Rhl×hl−1 , l = 1, . . . , L, hL = m, h0 = d, and h−1 = n. Suppose the Lipschitz constant of g is η. Then the covering numbers of S with respect to the Frobenius norm satisfy
N (S, ‖ · ‖F , ) ≤
( 3(L+ 1)ηL−1 ∏L l=0 βl )∑L l=0 hlhl−1
.
Now using Lemma 1 and Lemma 2, we have 1√ N ‖X − X̂‖F = 1√ N ‖Y −E − X̂‖F
≤ 1√ N ‖Y − X̂‖F + 1√ N ‖E‖F ≤ 1√ |Ω| ‖PΩ(Y − X̂)‖F + 1√ N ‖E‖F + 2 √ |Ω| + ( 8ξ4 log (|S|N) |Ω| )1/4
≤ 1√ |Ω| ‖PΩ(Y − X̂)‖F + 1√ N ‖E‖F + 2 √ |Ω| + 8ξ4 ( logN + ∑L l=0 hlhl−1 log 3(L+ 1)ηL−1 ∏L l=0 βl ) |Ω|  1/4
≤ 1√ |Ω| ‖PΩ(Y − X̂)‖F + 1√ N ‖E‖F + 6ξ(L+ 1)√ |Ω| +
( 8ξ4 ( logN + ∑L l=0 hlhl−1 log(η L−1ξ−1 ∏L l=0 βl))
|Ω|
)1/4
≤ 1√ |Ω| ‖PΩ(Y − X̂)‖F + 1√ N ‖E‖F + cξ
(∑L l=0 hlhl−1 log(η L−1ξ−1 ∏L l=0 βl)
|Ω|
)1/4 .
In the fourth inequality, we have set = 3ξ(L + 1). The last inequality holds because N is much smaller than |S| and c can be a constant large enough. This finished the proof.

heading: C.2 PROOF FOR THEOREM 2
text: The following lemma (proved in Appendix D.3) provides an upper bound for the covering number of the two-mode deep factorization matrix set. Lemma 3. Let S = {X ∈ Rn1×n2 : X = fθ1(Z)>Qfθ2(S), ‖Q‖F ≤ βQ, ‖Z‖F ≤ β
(1) 0 , ‖S‖F ≤ β (2) 0 , ‖W (j) l ‖F ≤ β (j) l , 1 ≤ l ≤ Lj , j = 1, 2; fθ1(Z) =
g(W (1) L1 g(W (1) L1−1 · · · g(W (1) 1 Z))), fθ2(S) = g(W (2) L2 g(W (2) L2−1 · · · g(W (2) 1 S))), }, where Z ∈ Rd1×n1 , S ∈ Rd2×n2 , W (j)l ∈ R h (j) l ×h (j) l−1 , l = 1, . . . , Lj , h (j) Lj = mj , h (j) 0 = dj , h (j) −1 = nj , j = 1, 2, and Q ∈ Rm1×m2 . Suppose the Lipschitz constant of g is η. Then the covering numbers of S with respect to the Frobenius norm satisfy
N (S, ‖ · ‖F , ) ≤ ( 3(k + 1)(max(L1, L2) + 1)s )p ,
where p = m1m2 + ∑2 j=1 ∑Lj l=0 h (j) l h (j) l−1 and s = βQη L1+L2 (∏L1 l=0 β (1) l )(∏L2 l=0 β (2) l ) .
Similar to the proof of Theorem 1, using Lemma 1, we have
1√ N ‖X − X̂‖F
≤ 1√ |Ω| ‖PΩ(Y − X̂)‖F + 1√ N ‖E‖F + 2 √ |Ω| +
( 8ξ4 log (|S|N)
|Ω|
)1/4 ,
where N = n1n2. Let = 3ξ(k + 1)(max(L1, L2) + 1). According to Lemma 3 and we obtain
log(|S|) = ( m1m2 + 2∑ j=1 Lj∑ l=0 h (j) l h (j) l−1 ) log ( βQη L1+L2ξ−1 ( L1∏ l=0 β (1) l )( L2∏ l=0 β (2) l )) .
Now we have
1√ N ‖X − X̂‖F ≤ 1√ |Ω| ‖PΩ(Y − X̂)‖F + 1√ N ‖E‖F
+ cξ
 ( m1m2 + ∑2 j=1 ∑Lj l=0 h (j) l h (j) l−1 ) log ( βQη L1+L2ξ−1 (∏L1 l=0 β (1) l )(∏L2 l=0 β (2) l )) |Ω| 1/4 . This finished the proof.

heading: C.3 PROOF FOR THEOREM 3
text: The following lemma (proved in Appendix D.4) provides an upper bound for the covering number of the multi-mode deep factorization tensor set.
Lemma 4. Let S = {X ∈ Rn1×n2···×nk : X = C ×1 fθ1(Z(1)) ×2 fθ2(Z(2)) · · · ×k fθk(Z (k)); fθj (Z) = g(W (j) Lj g(W (j) Lj−1 · · · g(W (j) 1 Z
(j)))), j = 1, . . . , k; ‖C‖F ≤ βC , ‖Z(j)‖F ≤ β(j)0 , ‖W (j) l ‖F ≤ β (j) l , 1 ≤ l ≤ Lj , j = 1, . . . , k}, where Z(j) ∈ Rdj×nj , W (j) l ∈ R h (j) l ×h (j) l−1 , l = 1, . . . , Lj , h (j) Lj = mj , h (j) 0 = dj , h (j) −1 = nj , j = 1, . . . , k, and C ∈ Rm1×m2···×mk . Suppose the Lipschitz constant of g is η. Then the covering numbers of S with respect to the Frobenius norm satisfy
N (S, ‖ · ‖F , ) ≤ ( 3(k + 1)(Lmax + 1)s )p ,
where p = ∏k j=1mj + ∑k j=1 ∑Lj l=0 h (j) l h (j) l−1 and s = βCη ∑k j=1 Lj ∏k j=1 ∏Lj l=0 β (j) l .
The proof is similar to that for Theorem 2. Using Lemma 1 and Lemma 4 , we arrive at
1√ N ‖X − X̂‖F ≤ 1√ |Ω| ‖PΩ(Y − X̂ )‖F + 1√ N ‖E‖F
+ cξ
 (∏k j=1mj + ∑k j=1 ∑Lj l=0 h (j) l h (j) l−1 ) log ( βCη ∑k j=1 Ljξ−1 ∏k j=1 ∏Lj l=0 β (j) l ) |Ω| 1/4 , where c is a constant.

heading: D PROOF FOR LEMMAS

heading: D.1 PROOF FOR LEMMA 1
text: We give the following lemma.
Lemma 5 (Hoeffding inequality for sampling without replacement (Serfling, 1974)). Let X1, X2, . . . , Xs be a set of samples taken without replacement from a distribution {x1, x2, . . . , xN} of mean u and variance σ2. Denote a = mini xi and b = maxi xi. Then
P [∣∣∣1 s s∑ i=1 Xi − u ∣∣∣ ≥ t] ≤ 2 exp(− 2st2 (1− (s− 1)/N)(b− a)2 ) .
We define
L̂(X ) := 1 |Ω| ‖PΩ(Y −X )‖2F ,
L(X ) := 1 N ‖Y −X‖2F ,
where N = ∏k i=1 ni. Suppose max{‖Y‖∞, ‖X‖∞} ≤ ξ. According to Lemma 5, we have
P [ L̂(X )− L(X )| ≥ t ] ≤ 2 exp ( − 2|Ω|t 2
(1− (|Ω| − 1)/N)δ2
) ,
where δ = 4ξ2. Using union bound for all X̄ ∈ S yields
P [ sup X̄∈S |L̂(X̄ )− L(X̄ )| ≥ t ]
≤2|S| exp ( − 2|Ω|t 2
(1− (|Ω| − 1)/N)δ2
) .
Or equivalently, with probability at least 1− 2N−1,
sup X̄∈S |L̂(X̄ )− L(X̄ )| ≤
√ δ2 log (|S|N)
2
( 1
|Ω| − 1 N + 1 N |Ω|
)
≤
√ δ2 log (|S|N)
2|Ω| , ∆.
Since | √ u− √ v| ≤ √ |u− v| holds for any non-negative u and v, we have
sup X̄∈S ∣∣∣∣√L̂(X̄ )−√L(X̄ )∣∣∣∣ ≤ √∆. Recall that ≥ ‖X − X̄‖F ≥ ‖P(X − X̄ )‖F , we have∣∣∣∣√L(X )−√L(X̄ )∣∣∣∣
= 1√ N ∣∣∣‖Y −X‖F − ‖Y − X̄‖F ∣∣∣ ≤ √ N
and ∣∣∣∣√L̂(X )−√L̂(X̄ )∣∣∣∣ =
1√ |Ω| ∣∣∣‖PΩ(Y −X )‖F − ‖PΩ(Y − X̄ )‖F ∣∣∣ ≤ √|Ω| . It follows that
sup X∈S ∣∣∣∣√L̂(X )−√L(X )∣∣∣∣ ≤ sup
X∈S ∣∣∣∣√L̂(X )−√L̂(X̄ )∣∣∣∣+ ∣∣∣∣√L̂(X̄ )−√L(X̄ )∣∣∣∣+ ∣∣∣∣√L(X̄ )−√L(X )∣∣∣∣ ≤ √ |Ω| + √ ∆ + √ N
≤ 2 √ |Ω| +
( δ2 log (|S|N)
2|Ω|
)1/4 .
This finished the proof.

heading: D.2 PROOF FOR LEMMA 2
text: Let X = WLg(WL−1 · · · g(W1Z)), where Z ∈ Rd×n, Wl ∈ Rhl×hl−1 , l = 1, . . . , L, and h0 = d. Let Sab := {W ∈ Ra×b : ‖W ‖F ≤ β}. Then there exists an -net S̄ab obeying
N (Sab, ‖ · ‖F , ) ≤ ( 3β )ab such that ‖W − W̄ ‖F ≤ . Now replace with /γ and let ‖Wl − W̄l‖F ≤
γl , l = 1, . . . , L,
‖Z − Z̄‖F ≤ γ0 . Let γl = (L+ 1)ηL−1
∏ i 6=l βi, l = 0, . . . , L.
‖X − X̄‖F =‖WLg(WL−1 · · ·W2g(W1Z))− W̄Lg(W̄L−1 · · · W̄2g(W̄1Z̄))‖F =‖WLg(WL−1 · · ·W2g(W1Z))± W̄Lg(WL−1 · · ·W2g(W1Z)) ± W̄Lg(W̄L−1 · · ·W2g(W1Z)) · · · ± W̄Lg(W̄L−1 · · · W̄2g(W1Z)) ± W̄Lg(W̄L−1 · · · W̄2g(W̄1Z))− W̄Lg(W̄L−1 · · · g(W̄1Z̄))‖F ≤‖WLg(WL−1 · · ·W2g(W1Z))− W̄Lg(WL−1 · · ·W2g(W1Z))‖F
+ ‖W̄Lg(WL−1 · · ·W2g(W1Z))− W̄Lg(W̄L−1 · · ·W2g(W1Z))‖F + · · ·+ ‖W̄Lg(W̄L−1 · · · W̄2g(W1Z))− W̄Lg(W̄L−1 · · · W̄2g(W̄1Z))‖F + ‖W̄Lg(W̄L−1 · · · W̄2g(W̄1Z))− W̄Lg(W̄L−1 · · · W̄2g(W̄1Z̄))‖F
≤ηL−1‖WL − W̄L‖F ‖Z‖F L−1∏ l=1 ‖Wl‖F
+ ηL−1‖WL−1 − W̄L−1‖F ‖Z‖F ‖W̄L‖F L−2∏ l=1 ‖Wl‖F
+ · · ·+ ηL−1‖W1 − W̄1‖F ‖Z‖F L−1∏ l=2 ‖W̄l‖F
+ ηL−1‖Z − Z̄‖F L∏ l=1 ‖W̄l‖F
≤ηL−1 γL L−1∏ l=0 βl + η L−1 γL−1 ∏ l 6=L−1 βl + · · ·+ ηL−1 γ1 ∏ l 6=1 βl + η L−1 γ0 ∏ l 6=0 βl
= .
(19)
The second inequality utilized the Lipschitz continuity of g and the submultiplicativity of the Frobenius norm. Therefore, S̄ is an -cover of S. We have
N (S, ‖ · ‖F , )
≤ L∏ l=0 ( 3βlγl )hlhl−1
= L∏ l=0
( 3(L+ 1)ηL−1 ∏L l=0 βl )hlhl−1
=
( 3(L+ 1)ηL−1 ∏L l=0 βl )∑L l=0 hlhl−1
.

heading: D.3 PROOF FOR LEMMA 3
text: Lemma 3 is a special case of Lemma 4.

heading: D.4 PROOF FOR LEMMA 4
text: Let ρ0 = ‖C‖F . For j = 1, . . . , k, we have
‖Uj‖F ≤ ηLj‖Z(j)‖F Lj∏ l=1 ‖W (j)l ‖F = η Lj Lj∏ l=0 β (j) l , ρj . (20)
Denote τ = ∏k j=0 ρj . Similar values can be defined for C̄ and Ūj .
Let Sjl := {W (j)l ∈ Ra×b : ‖W (j) l ‖F ≤ β (j) l }. Then there exists an -net S̄jl obeying
N (Sjl, ‖ · ‖F , ) ≤
( 3β
(j) l )h(j)l h(j)l−1 such that ‖W (j)l − W̄ (j) l ‖F ≤ . Let Γ (j) l = (Lj + 1)η Lj ∏ i 6=l βi, l = 0, . . . , Lj , and Υj =
(k + 1) ∏ t 6=j ρt, j = 0, 1, . . . , k. Replace with
Γ (j) l Υj
, we have ‖W (j)l − W̄ (j) l ‖F ≤
Γ (j) l Υj
and
have
N (Sjl, ‖ · ‖F , Γ (j) l Υj ) ≤
( 3Γ
(j) l Υjβ (j) l )h(j)l h(j)l−1 .
It follows that
‖Uj − Ūj‖F =‖g(W (j)Lj g(W (j) Lj−1 · · ·W (j) 2 g(W (j) 1 Z (j))))− g(W̄ (j)Lj g(W̄ (j) Lj−1 · · · W̄ (j) 2 g(W̄ (j) 1 Z̄ (j))))‖F
=‖g(W (j)Lj g(W (j) Lj−1 · · ·W (j) 2 g(W (j) 1 Z (j))))± g(W̄ (j)Lj g(W (j) Lj−1 · · ·W (j) 2 g(W (j) 1 Z̄ (j))))
± g(W̄ (j)Lj g(W̄ (j) Lj−1 · · ·W (j) 2 g(W (j) 1 Z (j)))) · · · ± g(W̄ (j)Lj g(W̄ (j) Lj−1 · · · W̄ (j) 2 g(W (j) 1 Z (j))))
± g(W̄Lg(W̄L−1 · · · W̄2g(W̄1Z))− g(W̄ (j)Lj g(W̄ (j) Lj−1 · · · g(W̄ (j) 1 Z̄ (j))))‖F
≤‖g(WLjg(W (j) Lj−1 · · ·W (j) 2 g(W (j) 1 Z (j))))− g(W̄ (j)Lj g(W (j) Lj−1 · · ·W (j) 2 g(W (j) 1 Z (j))))‖F
+ ‖g(W̄ (j)LJ g(W (j) Lj−1 · · ·W (j) 2 g(W (j) 1 Z (j))))− g(W̄ (j)Lj g(W̄ (j) Lj−1 · · ·W (j) 2 g(W (j) 1 Zj)))‖F
+ · · ·+ ‖g(W̄ (j)Lj g(W̄ (j) Lj−1 · · · W̄ (j) 2 g(W (j) 1 Z (j))))− g(W̄ (j)Lj g(W̄ (j) Lj−1 · · · W̄ (j) 2 g(W̄ (j) 1 Zj)))‖F
+ ‖g(W̄ (j)Lj g(W̄ (j) Lj−1 · · · W̄ (j) 2 g(W̄ (j) 1 Z (j))))− g(W̄ (j)Lj g(W̄ (j) Lj−1 · · · W̄ (j) 2 g(W̄ (j) 1 Z̄ (j))))‖F
≤ηLj‖W (j)Lj − W̄ (j) Lj ‖F ‖Z(j)‖F Lj−1∏ l=1 ‖W (j)l ‖F
+ ηLj‖W (j)Lj−1 − W̄ (j) Lj−1‖F ‖Z (j)‖F ‖W̄ (j)Lj ‖F Lj−2∏ l=1 ‖W (j)l ‖F
+ · · ·+ ηLj‖W (j)1 − W̄ (j) 1 ‖F ‖Z(j)‖F L−1∏ l=2 ‖W̄ (j)l ‖F + ηLj‖Z(j) − Z̄(j)‖F Lj∏ l=1 ‖W̄ (j)l ‖F
≤ηLj Γ
(j) L Υj Lj−1∏ l=0 β (j) l + η Lj Γ (j) L−1Υj ∏ l 6=Lj−1 β (j) l + · · ·+ η Lj Γ (j) 1 Υj ∏ l 6=1 β (j) l + η Lj Γ (j) 0 Υj ∏ l 6=0 β (j) l
= Υj .
(21)
Now we have
‖X − X̄‖F =‖C ×1 U1 ×2 · · · ×k Uk − C̄ ×1 Ū1 ×2 · · · ×k Ūk‖F =‖C ×1 U1 ×2 · · · ×k Uk ± C ×1 U1 ×2 · · · ×k Ūk ± C ×1 U1 ×2 · · · ×k−1 Ūk−1 ×k Ūk ± · · · ± C ×1 Ū1 ×2 · · · ×k−1 Ūk−1 ×k Ūk − C̄ ×1 Ū1 ×2 · · · ×k Ūk‖F ≤‖C ×1 U1 ×2 · · · ×k (Uk − Ūk)‖F + ‖C ×1 U1 ×2 · · · ×k−1 (Uk−1 − Ūk−1)×k Ūk‖F
+ · · ·+ ‖C ×1 (U1 − Ū1)×2 Ū2 · · · ×k Ūk‖F + ‖(C − C̄)×1 Ū1 ×2 Ū2 · · · ×k Ūk‖F ≤‖C‖F ‖U1‖F · · · ‖Uk−1‖F ‖Uk − Ūk‖F + ‖C‖F ‖U1‖F · · · ‖Uk−1 − Ūk−1‖F ‖Ūk‖F
+ · · ·+ ‖C‖F ‖U1 − Ū1‖F ‖Ū2‖F · · · ‖Ūk‖F + ‖C − C̄‖F ‖Ū1‖F ‖Ū2‖F · · · ‖Ūk‖F ≤‖C‖F ‖U1‖F · · · ‖Uk−1‖F
Υk + ‖C‖F ‖U1‖F · · · Υk−1 ‖Ūk‖F
+ · · ·+ ‖C‖F Υ1 ‖Ū2‖F · · · ‖Ūk‖F + Υ0 ‖Ū1‖F ‖Ū2‖F · · · ‖Ūk‖F
≤  1 Υk ∏ j 6=k ρj + 1 Υk−1 ∏ j 6=k−1 ρj + · · ·+ 1 Υ1 ∏ j 6=1 ρj + 1 Υ0 ∏ j 6=0 ρj  =
Therefore S̄ is an -cover of S. Then the covering number of Sdnk can be bounded as
N (Sdnk, ‖ · ‖F , )
≤ ( 3Υ0ρ0 )∏k j=1mj k∏
j=1 Lj∏ l=0
( 3Γ
(j) l Υjβ (j) l )h(j)l h(j)l−1 .
=
( 3Υ0ρ0 )∏k j=1mj k∏
j=1 Lj∏ l=0 ( 3(Lj + 1)ρjΥj )h(j)l h(j)l−1
=
( 3(k + 1)τ )∏k j=1mj k∏
j=1 Lj∏ l=0 ( 3(k + 1)(Lj + 1)τ )h(j)l h(j)l−1 .
≤ ( 3(k + 1)τ )∏k j=1mj k∏
j=1
( 3(k + 1)(Lmax + 1)τ )∑Lj l=0 h (j) l h (j) l−1
=
( 3(k + 1)τ )∏k j=1mj ( 3(k + 1)(Lmax + 1)τ )∑k j=1 ∑Lj l=0 h (j) l h (j) l−1
≤ ( 3(k + 1)(Lmax + 1)τ )∏k j=1mj+ ∑k j=1 ∑Lj l=0 h (j) l h (j) l−1
=
3(k + 1)(Lmax + 1)ρ0∏kj=1 ( ηLj ∏Lj l=0 β (j) l )  ∏k j=1mj+ ∑k j=1 ∑Lj l=0 h (j) l h (j) l−1
.

Review summaries:
This paper summarizes the theoretical guarantee for the existing LRMC and LRTC algorithms, and provides the theoretical analysis for a new proposed Multi-Mode Nonlinear deep tensor factorization. The analytical results show that when n2 is larger than n1, the nonlinear DMF provides a tighter generalization bound than MF. Similar analysis has been extended to two-mode matrix factorization and multi-mode  tensor factorization. Experimental results in synthetic data and real data show better results of the proposed algorithm as compared to other algorithms in completion tasks. 

To bridge the gap between deep learning and tensor decomposition, this paper presents two novel approaches named as two mode non-linear deep matrix factorization and multi-mode nonlinear deep tensor factorization (extension of two mode model to multi-mode scenario). The main contribution of the methods lie in full exploration of non-linearity of data in matrix and tensor factorization. To better motivate the proposed models, the authors provide theoretical analysis for why and when nonlinear deep matrix factorization outperforms linear deep matrix factorization in matrix completion. The experimental evaluation demonstrates that in some datasets, the proposed models outperform the existing models on matrix and tensor completion tasks.

This paper studies the nonlinear low-rank completion of matrices and tensors. Specifically, it first presents the theoretical results showing why nonlinear deep matrix factorization is better than the ordinary matrix factorization model. Then, it proposes a model named two-mode nonlinear deep matrix factorization to make full use of the nonlinearity of the nearly square matrices. The authors also extend this method to tensor factorization by further factorizing the factor matrices in the Tucker decomposition using the deep factorization method. Impressive results are obtained using both synthetic and real-world datasets.

The paper provides a multi-mode framework for the deep learning based tensor decomposition, which could be useful for dealing with nonlinear high-dimensional data sets. In particular, it extends the deep matrix factorization (DMF) method and proposes a multi-mode deep matrix factorization method for matrix completion with convergence guarantee. Based on this, it also develops a multi-mode nonlinear deep tensor factorization method with convergence guarantee. The proposed models are solved by various optimization algorithms. Numerical experiments on synthetic and real data sets of the matrix/tensor form have shown that the proposed methods outperform other state of the arts.