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.

summaries: 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.

summaries: 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.

summaries: 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.

abstractText: 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.

authors:
Evrim Acar
Rasmus Bro
Bonnie Schmidt

title: New exploratory clustering tool

venue: Journal of Chemometrics: A Journal of the Chemometrics Society,

authors:
Evrim Acar
Daniel M Dunlavy
Tamara G Kolda
Morten Mørup

title: Scalable tensor factorizations for incomplete data

venue: Chemometrics and Intelligent Laboratory Systems,

authors:
Sanjeev Arora
Nadav Cohen
Wei Hu
Yuping Luo

title: Implicit regularization in deep matrix factorization

authors:
Boaz Barak
Ankur Moitra

title: Noisy tensor completion via the sum-of-squares hierarchy

venue: In Conference on Learning Theory,

authors:
Sabina Bijlsma
Age K Smilde

title: Estimating reaction rate constants from a two-step reaction: a comparison between two-way and three-way methods

venue: Journal of Chemometrics,

title: Parafac. tutorial and applications

venue: Chemometrics and intelligent laboratory systems,

title: Exploratory study of sugar production using fluorescence spectroscopy and multiway analysis

authors:
Emmanuel J. Candès
Benjamin Recht

title: Exact matrix completion via convex optimization

venue: Foundations of Computational Mathematics,

authors:
Yudong Chen
Srinadh Bhojanapalli
Sujay Sanghavi
Rachel Ward

title: Coherent matrix completion

venue: In Proceedings of the 31st International Conference on Machine Learning

title: Approximation by superpositions of a sigmoidal function

venue: Mathematics of control, signals and systems,

authors:
Jicong Fan
Jieyu Cheng

title: Matrix completion by deep matrix factorization

venue: Neural Networks,

authors:
Jicong Fan
Tommy W.S. Chow

title: Matrix completion by least-square, low-rank, and sparse selfrepresentations

venue: Pattern Recognition,

authors:
Jicong Fan
Lijun Ding
Yudong Chen
Madeleine Udell

title: Factor group-sparse regularization for efficient low-rank matrix recovery

venue: In Advances in Neural Information Processing Systems

authors:
Jicong Fan
Mingbo Zhao
Tommy W.S. Chow

title: Matrix completion via sparse factorization solved by accelerated proximal alternating linearized minimization

venue: IEEE Transactions on Big Data,

authors:
Dylan J Foster
Andrej Risteski

title: Sum-of-squares meets square loss: Fast rates for agnostic tensor completion

authors:
Silvia Gandy
Benjamin Recht
Isao Yamada

title: Tensor completion and low-n-rank tensor recovery via convex optimization

venue: Inverse Problems,

authors:
Paris Giampouras
Rene Vidal
Athanasios Rontogiannis
Benjamin Haeffele

title: A novel variational form of the schatten-p quasi-norm

authors:
Suriya Gunasekar
Blake E Woodworth
Srinadh Bhojanapalli
Behnam Neyshabur
Nati Srebro

title: Implicit regularization in matrix factorization

authors:
Boris Hanin
Mark Sellke

title: Approximating continuous functions by relu nets of minimal width

venue: arXiv preprint arXiv:1710.11278,

title: Understanding alternating minimization for matrix completion

venue: IEEE 55th Annual Symposium on Foundations of Computer Science (FOCS),,

authors:
Christopher J Hillar
Lek-Heng Lim

title: Most tensor problems are np-hard

venue: Journal of the ACM (JACM),

authors:
Yao Hu
Debing Zhang
Jieping Ye
Xuelong Li
Xiaofei He

title: Fast and accurate matrix completion via truncated nuclear norm regularization

venue: IEEE Transactions on Pattern Analysis and Machine Intelligence,

authors:
Huyan Huang
Yipeng Liu
Jiani Liu
Ce Zhu

title: Provable tensor ring completion

venue: Signal Processing,

authors:
Christian Igel
Michael Hüsken

title: Improving the rprop learning algorithm. In Proceedings of the second international ICSC symposium on neural computation (NC

authors:
Prateek Jain
Sewoong Oh

title: Provable tensor factorization with missing data

venue: In Advances in Neural Information Processing Systems,

authors:
Chi Jin
Sham M Kakade
Praneeth Netrapalli

title: Provable efficient online matrix completion via non-convex stochastic gradient descent

venue: In Proceedings of the 30th International Conference on Neural Information Processing Systems,

authors:
Misha E Kilmer
Karen Braman
Ning Hao
Randy C Hoover

title: Third-order tensors as operators on matrices: A theoretical and computational framework with applications in imaging

venue: SIAM Journal on Matrix Analysis and Applications,

authors:
Diederik P Kingma
Jimmy Ba

title: Adam: A method for stochastic optimization

venue: arXiv preprint arXiv:1412.6980,

authors:
Tamara G Kolda
Brett W Bader

title: Tensor decompositions and applications

venue: SIAM review,

authors:
Daniel Kressner
Michael Steinlechner
Bart Vandereycken

title: Low-rank tensor completion by riemannian optimization

venue: BIT Numerical Mathematics,

authors:
Christian Kümmerle
Juliane Sigl

title: Harmonic mean iteratively reweighted least squares for lowrank matrix recovery

venue: The Journal of Machine Learning Research,

authors:
Dong C. Liu
Jorge Nocedal

title: On the limited memory bfgs method for large scale optimization

venue: Mathematical Programming,

authors:
Hanpeng Liu
Yaguang Li
Michael Tsang
Yan Liu

title: Costco: A neural tensor completion model for sparse tensors

venue: In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining,

authors:
Ji Liu
Przemyslaw Musialski
Peter Wonka
Jieping Ye

title: Tensor completion for estimating missing values in visual data

venue: IEEE transactions on pattern analysis and machine intelligence,

authors:
Zhou Lu
Hongming Pu
Feicheng Wang
Zhiqiang Hu
Liwei Wang

title: The expressive power of neural networks: a view from the width

venue: In Proceedings of the 31st International Conference on Neural Information Processing Systems,

authors:
Feiping Nie
Heng Huang
Chris Ding

title: Low-rank matrix recovery via efficient Schatten p-norm minimization

venue: In Proceedings of the Twenty-Sixth AAAI Conference on Artificial Intelligence,

authors:
Lars Nørgaard
Carsten Ridder

title: Rank annihilation factor analysis applied to flow injection analysis with photodiode-array detection

title: Approximation theory of the mlp model

venue: Acta Numerica 1999: Volume

authors:
Noam Razin
Asaf Maman
Nadav Cohen

title: Implicit regularization in tensor factorization

title: A simpler approach to matrix completion

venue: J. Mach. Learn. Res.,

authors:
Jasson DM Rennie
Nathan Srebro

title: Fast maximum margin matrix factorization for collaborative prediction

venue: In Proceedings of the 22nd international conference on Machine learning,

authors:
Suvash Sedhain
Aditya Krishna Menon
Scott Sanner
Lexing Xie

title: Autorec: Autoencoders meet collaborative filtering

venue: In Proceedings of the 24th international conference on World Wide Web,

title: Probability inequalities for the sum in sampling without replacement

venue: The Annals of Statistics, pp

authors:
Ohad Shamir
Shai Shalev-Shwartz

title: Matrix completion with the trace norm: learning, bounding, and transducing

authors:
Fanhua Shang
Yuanyuan Liu
James Cheng

title: Tractable and scalable Schatten quasi-norm approximations for rank minimization

venue: In Artificial Intelligence and Statistics,

authors:
Nathan Srebro
Ruslan R Salakhutdinov

title: Collaborative filtering in a non-uniform world: Learning with the weighted trace norm

venue: Advances in Neural Information Processing Systems

authors:
Nathan Srebro
Adi Shraibman

title: Rank, trace-norm and max-norm

venue: In International Conference on Computational Learning Theory,

authors:
Nathan Srebro
Jason Rennie
Tommi S. Jaakkola

title: Maximum-margin matrix factorization

venue: In Advances in neural information processing systems,

authors:
Ruoyu Sun
Zhi-Quan Luo

title: Guaranteed matrix completion via non-convex factorization

venue: IEEE Transactions on Information Theory,

authors:
George Trigeorgis
Konstantinos Bousmalis
Stefanos Zafeiriou
Björn W Schuller

title: A deep matrix factorization method for learning attribute representations

authors:
Andong Wang
Chao Li
Zhong Jin
Qibin Zhao

title: Robust tensor decomposition via orientation invariant tubal nuclear norms

venue: In AAAI,

authors:
Qi Wang
Mengying Sun
Liang Zhan
Paul Thompson
Shuiwang Ji
Jiayu Zhou

title: Multimodality disease modeling via collective deep matrix factorization

venue: In Proceedings of the 23rd ACM SIGKDD international conference on knowledge discovery and data mining,

authors:
Zaiwen Wen
Wotao Yin
Yin Zhang

title: Solving a low-rank factorization model for matrix completion by a nonlinear successive over-relaxation algorithm

venue: Mathematical Programming Computation,

authors:
JP Wold
Rasmus Bro
A Veberg
F Lundby
AN Nilsen
J Moan

title: Active photosensitizers in butter detected by fluorescence spectroscopy and multivariate curve resolution

venue: Journal of agricultural and food chemistry,

authors:
Q. Xie
Q. Zhao
D. Meng
Z. Xu

title: Kronecker-basis-representation based tensor sparsity and its applications to tensor recovery

authors:
Yangyang Xu
Ruru Hao
Wotao Yin
Zhixun Su

title: Parallel matrix factorization for low-rank tensor completion

venue: arXiv preprint arXiv:1312.1254,

authors:
Hong-Jian Xue
Xinyu Dai
Jianbing Zhang
Shujian Huang
Jiajun Chen

title: Deep matrix factorization models for recommender systems

venue: In IJCAI,

authors:
Longhao Yuan
Chao Li
Danilo Mandic
Jianting Cao
Qibin Zhao

title: Tensor ring decomposition with rank minimization on latent space: An efficient approach for tensor completion

venue: In Proceedings of the AAAI Conference on Artificial Intelligence,

authors:
Ming Yuan
Cun-Hui Zhang

title: On tensor completion via nuclear norm minimization

authors:
Zemin Zhang
Shuchin Aeron

title: Exact tensor completion using t-svd

venue: IEEE Transactions on Signal Processing,

authors:
Handong Zhao
Zhengming Ding
Yun Fu

title: Multi-view clustering via deep matrix factorization

authors:
Qibin Zhao
Guoxu Zhou
Shengli Xie
Liqing Zhang
Andrzej Cichocki

title: Tensor ring decomposition

venue: arXiv preprint arXiv:1606.05535,

authors: DMF Fan Cheng

title: 2018), the network structure

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
.


summaries: The paper works on a novel problem of interpreting and explaining structured output models. The paper utilizes an energy based model to account for correlations between structured outputs and learns an interpretability block which given as input an image learns to mask it such that the energy based model would assign a similar score to the ground truth output and input as well as the ground truth output and perturbed input. In essence, the energy based model is a proxy for the actual deep neural network performing the structured prediction task. Results on a couple of datasets demonstrate that the work does better than baselines like LIME which do not utilize the correlations in the outputs modeled by the energy based model.

summaries: The paper proposes a technique for identifying what input variables are most relevant for determining the value of a single, given output variable in structured-output (MAP) inference.  The idea is to learn an energy model that predicts which input variables are relevant to a particular structured prediction (x, y).  The authors propose to implement the energy model using a neural network followed by a Gumbel-softmax activation, and to train it by maximizing a structured hinge loss.  The proposed approach is evaluated on three datasets and compared to standard attribution techniques (LIME, SHAP, L2X).

summaries: The authors propose an energy-based training method for achieving model interpretability, which performs instance-wise feature selection. The proposed model adopts a similar approach of one of the pre-existing interpretable methods by calculating feature-level importance score with regard to each instance. The authors validate their method on synthetic and public datasets.

summaries: This paper propose a method for interpreting structured output model. The key idea is to find an "interpretation" which explains an "target" output random variable based on subset of rest of the output variables. The training objective is on finding a small subset which keeps the target output random variable invariant. The proposed methodology is applied to explain a synthetic energy function and structured prediction energy networks.

abstractText: We propose a novel interpretation technique to explain the behavior of structured output models, which learn mappings between an input vector to a set of output variables simultaneously. Because of the complex relationship between the computational path of output variables in structured models, a feature can affect the value of output through other ones. We focus on one of the outputs as the target and try to find the most important features utilized by the structured model to decide on the target in each locality of the input space. In this paper, we assume an arbitrary structured output model is available as a black box and argue how considering the correlations between output variables can improve the explanation performance. The goal is to train a function as an interpreter for the target output variable over the input space. We introduce an energy-based training process for the interpreter function, which effectively considers the structural information incorporated into the model to be explained. The effectiveness of the proposed method is confirmed using a variety of simulated and real data sets.

authors:
S. Bach
A. Binder
G. Montavon
F. Klauschen
K.R. Mller
W. Samek

title: On pixel-wise explanations for non-linear classifier decisions by layer-wise relevance propagation

venue: PloS one,

authors:
D. Belanger
A. McCallum

title: Structured prediction energy networks

venue: In Proceedings of the 33th International Conference on Machine Learning (ICML 2016),

authors:
D. Belanger
B. Yang
A. McCallum

title: End-to-end learning for structured prediction energy networks

venue: arXiv preprint arXiv:1703.05667,

authors:
J. Chen
L. Song
M.J. Wainwright
M.I. Jordan

title: Learning to explain: An information-theoretic perspective on model interpretation

venue: In Proceedings of the 35th International Conference on Machine Learning,

authors:
Liang-Chieh Chen
Alexander Schwing
Alan Yuille
Raquel Urtasun

title: Learning deep structured models

venue: In International Conference on Machine Learning,

authors:
Colin Graber
Ofer Meshi
Alexander Schwing

title: Deep structured prediction with nonlinear output transformations

venue: arXiv preprint arXiv:1811.00539,

authors:
M. Gygli
M. Norouzi
A. Angelova

title: Deep value networks learn to evaluate and iteratively refine structured outputs

venue: arXiv preprint arXiv:1703.04363,

authors:
E. Jang
S. Gu
B. Poole

title: Categorical reparameterization with gumbel-softmax

venue: In stat,

authors:
S.M. Lundberg
S.I. Lee

title: A unified approach to interpreting model predictions

authors:
W.J. Murdoch
C. Singh
K. Kumbier
R. Abbasi-Asl
B. Yu

title: Interpretable machine learning: definitions, methods, and applications

authors:
Jian Peng
Liefeng Bo
Jinbo Xu

title: Conditional neural fields

venue: Advances in neural information processing systems,

authors:
M.T. Ribeiro
S. Singh
C. Guestrin

title: Why should i trust you?: Explaining the predictions of any classifier

venue: In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining,

authors:
Alexander G Schwing
Raquel Urtasun

title: Fully connected deep structured networks

venue: arXiv preprint arXiv:1503.02351,

authors:
Avanti Shrikumar
Peyton Greenside
Anshul Kundaje

title: Learning important features through propagating activation differences

venue: In Proceedings of the 34th International Conference on Machine Learning,

authors:
K. Simonyan
V. Andrea
Z. Andrew

title: Deep inside convolutional networks: Visualising image classification models and saliency maps

venue: arXiv preprint arXiv:1312.6034,

authors:
M.D. Zeiler
R. Fergus

title: Visualizing and understanding convolutional networks

venue: arXiv preprint arXiv:1311.2901,

authors:
J. Zhou
O.G. Troyanskaya

title: Predicting effects of noncoding variants with deep learning-based sequence model

venue: Nature Methods, pp

authors:
M.L. Zintgraf
S.T. Cohen
T. Adel

title: Visualizing deep neural network decisions: Prediction difference analysis

venue: In Proceedings of the 5th International Conference on Learning Representations (ICLR

heading: 1 INTRODUCTION

text:
The impressive prediction performance of novel machine learning methods has motivated researchers of different fields to apply these models in challenging problems. However, their complex and non-linear inherence limit the ability to explain what they have learned. Interpretation gets more attention when we want to discover the reasons behind the model’s decision and be sure about the trustworthiness and fairness of a trained machine learning model in areas such as medicine, finance, and judgment. Additionally, interpreting a model with a satisfying prediction accuracy in a scientific problem, which results in understanding relationships behind the data, leads to new knowledge about the problem domain.Murdoch et al. (2019)
In many real-world applications, the goal is to map an input variable to a high-dimensional structured output, e.g., image segmentation and sequence labeling. In such problems, the output space includes a set of statistically related random variables. As considering these dependencies can increase the prediction accuracy, many structured output models have been introduced. Many of these methods use graphical models, including random fields, to capture the structural relations between variables. Most define an energy function over these random fields, with a global minimum at the ground truth. Therefore, an inference is needed to find the best configuration of output variables for input by minimizing the energy function in the prediction step. Early efforts to utilize deep neural networks in structured output problems adopt deep networks to extract high-level features from the input vector to incorporate them in calculating the energy function Peng et al. (2009); Chen et al. (2015); Schwing & Urtasun (2015). The computational complexity of the inference step in models that use random fields limits their ability to incorporate complex structures and interactions between output variables. Recent works in Belanger & McCallum (2016); Gygli et al. (2017); Belanger et al. (2017); Graber et al. (2018) propose to adopt deep networks instead of random fields to model the structure of the output space. Nevertheless, complex interactions between problem variables in such models make their interpretation too challenging, specifically when we focus on the model behavior in predicting a single output variable.
This paper attempts to interpret a structured output model by focusing on each output variable separately. Our approach to model interpretation is based on instance-wise feature selection. Its goal is to find the relative importance of each input feature in predicting a single output variable. The subset of important features can vary across the input space. The complicated interactions between computational paths of output variables in structured output models cause critical challenges for finding a subset of important features associated with each output variable. A feature may not be used directly in the computational path of output but affects its value through relations with other outputs. To compute the importance of a feature for a target output variable, we should aggregate its effect on all output variables that are correlated to this target.
Existing approaches of model interpretation can be divided into two groups, model-based and post hoc analysis Murdoch et al. (2019). The model-based interpretation approach encourages machine learning methods that readily provide insight into what the model learned. However, it leads to simple models that are not sufficiently effective for complex structured output problems. Here we follow the post hoc analysis and try to explain the behavior of a trained, structured output model provided as a black box. Many interpretation techniques to find the importance of features as a post hoc analysis have been introduced. Works in Zhou & Troyanskaya (2015); Zeiler & Fergus (2013); Zintgraf et al. (2017) make perturbations to some features and observe their impact on the final prediction. These techniques are computationally inefficient in situations where we search for the most valuable features. Since we should perform a forward propagation for all possible perturbations, in another trend, works in Simonyan et al. (2013); Bach et al. (2015) back-propagate an importance signal from the target output through the network to calculate the critical signal of features by calculating the gradient of the target w.r.t the input features. These models are computationally more efficient than perturbation-based techniques because they need only one pass of propagating. However, they need the structure of the network to be known. As this approach may cause a saturation problem, DeepLIFT Shrikumar et al. (2017) proposes that instead of propagating a gradient signal, the difference of the output from a reference value in terms of the difference of features from a reference value to be considered. In addition to these approaches, other ideas have also been introduced in model interpretation. Authors in Ribeiro et al. (2016) introduce LIME which trains a local interpretable surrogate model to simulate the behavior of a black box model in the vicinity of a sample. It randomly selects a set of instances of the input space around that sample and obtains the black box prediction for them, and trains the surrogate model by this new dataset. Therefore this interpretable model is a good approximation of the black box around the locality of the selected sample. Shapley value, a concept from the game theory, explains how to distribute an obtained payout between coalition players fairly. The work in Lundberg & Lee (2017) proposes the kernel SHAP for approximating the shapely value for each feature as its importance for a prediction. As an information-theoretic perspective on interpretation, the work in Chen et al. (2018) proposes to find a subset of features with the highest mutual information with the output. This subset is expected to involve the most important features for the output.
Existing interpretation techniques can be applied to explain the behavior of a structured model, w.r.t. a single output, by ignoring other output variables. However, none of these approaches consider possible correlations between output variables and only analyze the marginal behavior of the black box on the target. In this paper, we attempt to incorporate the structural information between output variables during training the interpreter. As our goal is to present a local interpreter, which is trained globally as Chen et al. (2018), we train a function over the input space which returns the index of most important features for decision making about the target output. Since the value of other output variables affects the value of the target, incorporating them into the training procedure of an interpreter function may lead to higher performance and decrease our uncertainty about the black box behavior. To the best of our knowledge, this is the first time an interpreter is designed mainly for structured output models, and dependencies between output variables are considered during training the interpreter.

heading: 2 PRELIMINARIES AND MOTIVATION

text:
Structured output prediction models map an arbitrary n-dimensional feature vector x ∈ X to the output y ∈ Y where y = [y1, y2, . . . , yd] includes a set of correlated variables with known and unknown complex relationships and Y shows a set of valid configurations.
Now we explain our intuition about an interpreter, which explains the behavior of a structured output model in predicting a single output variable. We assume a structured model is available as a black box about which we do not know. Our goal is to find indices of k important features of x which affect the black box prediction about the target output yt. As for different localities of the input space these indices may vary, the proposed interpreter is a function IN t(x;α) : X → {0, 1}n over the input space with a set of parameters α which returns an n-dimensional k-hot vector. In this vector, the value of 1 shows the indices of selected k important features for target output yt.
We define Θsb as the set of all parameters and hidden variables inside the structured black box. The probabilistic graphical model of Fig. 1 describes dependencies between problem variables. In this figure, x shows the input variable, ysbt and y sb −t = {ysbi |i 6= t} show black box predictions and αIN t is the set of parameters of IN t. The bidirectional edge between ysbt and ysb−t emphasizes the correlation between the outputs of a structured model. In fact αIN t is determined based on Θsb and the black box architecture, and the final prediction of the ysbt does not directly affect its value. However, here, Θsb is a latent variable which makes active paths between αIN t and output values ysbt and y sb −t. Therefore αIN t and y sb −t are dependent random variables and we have:
H(αIN t |x,ysbt ) > H(αIN t |x,ysbt ,ysb−t) (1)
where H(.|.) shows the conditional entropy. We use the strict inequality because αIN t and ysb−t are dependent random variables. The left term measures our uncertainty when we train the interpreter only by observing the target output ysbt . This inequality confirms that the uncertainty is decreased when we consider observed ysb−t during estimating αIN t . Motivated by this fact we propose a training procedure for an interpreter IN t which incorporates the structural information of the output space by observing the black box prediction on all output variables.
We call our method SOInter as we propose it to train an Interpreter specifically for Structured Output models.

heading: 3 PROPOSED METHOD

text:
We consider psb(y|x) as the distribution by which the structured black box predicts the output as follows,
ysb = arg max y psb(y|x). (2)
Our goal is to train the interpreter IN t(x;α) which explores a subset of most important features that affects the value of black box prediction on the target output yt in each locality of the input space. As the desired interpreter detects the subset of most important features, we expect perturbating other ones does not change the black box prediction of the target yt. Motivated by this statement, we are encouraged to compare the black box prediction for the target output when a sample and its perturbated version are passed through the black box. The interpreter IN t(x;α) returns a k-hot vector in which the value of 1 shows the index of a selected feature. We define x̃, a perturbated
version of x, as follows,
x̃ = x IN t(x;α) (3) in which only features selected by the interpreter are included, and other ones are filled with zeros. Therefore by passing x and x̃ through the black box, we expect the value of the tth element of predictions to be the same, and we can define a penalty over the value of the target in these two situations. However, since the structure of the black box is unknown and the arg max operation in equation 2 is non-differentiable, a loss function that directly compares these two output values can not be used to find the optimal interpreter. Therefore, in the following subsection, we try to achieve a penalty according to the difference between these values for the target, which can transfer the gradient through an interpreter block.

heading: 3.1 OBTAINING A TRACTABLE LOSS FUNCTION

text:
We consider ỹ as the black box prediction when masked inputs x̃ are given to the black box i.e.,
ỹ = arg max y
psb(y|x̃). (4)
We define distribution qα to estimate the probability density function of ỹ. We assume qα is from the exponential family and can be parameterized as follows,
qα(y|x) = exp(−Eq(x̃, y))∫
y′ exp(−Eq(x̃, y′)) (5)
As we define qα the distribution over the output variables when the masked input is given to the black box, Eq can be considered as an energy function describing the structural information incorporated into the black box. We train a deep neural network to estimate this energy block, using a set of samples from the input space and their associated outputs predicted by the black box. For an input feature x, the obtained energy function Eq(x, y) has the minimum value when y is equal to the black box prediction for x. The detailed architecture of this block depends on the data inherence. Different techniques to train an energy network have been introduced recently Belanger & McCallum (2016); Belanger et al. (2017); Gygli et al. (2017) which can be used to train Eq in a pre-training phase. Here we use the work in Gygli et al. (2017). Since after the pre-training phase, the only unknown parameters of qα are parameters of IN t, we show it by the index of α. Now to find the optimum value for α, we define the following penalty over the log-likelihood of qα,
SL(x, ysbt , ỹ; α) = max{0, (log qα([ỹt, ỹ−t] |x)− log qα( [ ysbt , ỹ−t ] |x) + L(ỹt,ysbt ))}.
(6)
When the target values for an input x and its perturbated version x̃ are the same, i.e. ỹt = ysbt , the penalty value is zero. Otherwise, as ỹ maximizes qα, SL(x, ysbt , ỹ; α) has a positive value which is proportional to the difference between the desired target value ysbt and the most offending answer ỹt. Using the log-likelihood of the joint distribution qα(.|x) in calculating the penalty, which ỹ−t is the most probable value for other outputs according to the qα, leads to incorporate the structural information between output variables.
Now we attempt to achieve a more simpler form of equation 6. According to equation 5 the loglikelihood of qα equals to,
log qα(y|x) = −Eq(x̃, y)− log ∫
y′ exp(−Eq(x̃, y′)). (7)
Because the second term of equation 7 for both of log qα( [ ysbt , ỹ−t ] |x) and log qα([ỹt, ỹ−t] |x) is the same, we can write the equation 6 as follows,
SL(x, ysbt , ỹ; α) = max{0, (Eq(x̃, [ ysbt , ỹ−t ] )− Eq(x̃, [ỹt, ỹ−t]) + L(ỹt,ysbt ))}. (8)
We propose to iteratively calculate the most offending answer ỹ for the current qα and then minimize the penalty function of equation 8 according to α. It is worth mentioning that we consider a constraint only over the ysbt as we intend to find the best interpreter for the target. As Eq is a deep neural network and x̃ is an element-wise multiplication of x and IN t(x;α), the gradient over the penalty of equation 8 can be back-propagated through the interpreter block.
The final proposed iterative optimization procedure for an interpreter block can be expressed as follows,
α(k) = arg min α Ep(x)[SL(x, ysbt , ỹ(k−1); α)]
ỹ(k) = arg max y psb(y|x IN t(x;α(k))) (9)
where ysb = arg maxy psb(y|x). In the first step, the expected penalty of equation 8 is minimized to find the optimum value for the parameter α. The penalty function SL is minimized when we find a value for α which makes [ ysbt ,y−t ] the strict maximum of qα(.|x) with the margin proportional to L. In the second step, the black box prediction for the input which is masked by the output of the current interpreter is calculated.
The initial value α(0) is selected randomly and its associated ỹ(0) is obtained using the second step of equation 9. The algorithm is continued until the value of the penalty SL does not considerably change which is usually obtained in less than 100 iterations. Here we consider L = I(ysbt , ỹt) which is an indicator function

heading: 3.2 THE INTERPRETER BLOCK

text:
The interpreter IN t includes a deep neural network, with a set of parameters α, followed by a Gumbel-Softmax Jang et al. (2017) unit. The detailed architecture of the deep network depends on the inherence of x. The dimension of the interpreter output is the same as the feature vector x. Fig. 2 describes the architecture used for training the interpreter IN t(x, α). The output of the deep networkWα shows the importance of the elements of the feature vector x. To encourage the interpreter to find top k important features associated with the target output yt, we use the GumbelSoftmax trick as proposed in Chen et al. (2018). To obtain top k important features, we consider the output ofWα(x) as parameters of the categorical distribution. Then we can independently draw a sample for k times. Each sample is a one-hot vector in which the element with the value of 1 shows the selected feature. To have a k-hot vector, we can simply get the element-wise maximum of these one-hot vectors. However this sampling process is not differentiable and we use its continuous approximation introduced by the Gumbel-Softmax trick. Considering following random variables,
gi = − log(− log(ui)) (10)
where ui ∼ Uniform(0, 1), we can use the re-parameterization trick instead of direct sampling fromWα(x) as follows:
ci = exp{log(Wα(x)i + gi)/τ}
Σj exp{log(Wα(x)j + gj)/τ} (11)
The vector c is the continuous approximation of the sampled one-hot vector. To have k selected features, we draw k vectors cj , j = 1 . . . k and obtain their element-wise maximum as follows Chen et al. (2018),
IN t(x, α)i = maxj {c j i , j = 1 . . . k} (12)
Finally, Algorithm 1 of Appendix A summarizes what we propose to train the interpreter IN t.

heading: 4 EXPERIMENTS

text: We evaluate the performance of our proposed interpreter on both synthetic and real datasets. In section 4.1, we define two arbitrary energy functions to synthesize structured data. We compare the performance of SOInter with two well-known interpretation techniques, Lime and Kernel-Shap, which are frequently used to evaluate the performance of interpretation methods, and L2X Chen et al. (2018) which proposes an information-Theoretic method for interpretation. None of these techniques are specifically designed for structured models. Indeed, they only consider the target output and ignoring other ones. In section 4.2 and 4.3, the efficiency of SOInter is shown with two real text and image datasets.

heading: 4.1 SYNTHETIC DATA

text:
Here we define two arbitrary energy functions on input vector x and output variables y, E1 and E2 in equation 13, which are linear and non-linear functions respectively according to the input features.
E1 = (x1y1 + x4)(1− y2) + (x2(1− y1) + x3)y2 E2 = (sin(x1)y1y3 + |x4|) (1− y2)y4 + ( exp(
x2 10 − 1)(1− y1)(1− y3) + x3
) y2(1− y4)
(13)
Input features are randomly generated using the standard normal distribution. Output variables are considered as binary discrete variables. For each input vector x, we found the corresponding output by the following optimization,
y∗ = arg min y E(x, y) (14)
where E shows the energy function from which we attempt to generate data. E1 describes the energy value over a structured output of size 2 and E2 describes an output of size 4. For each scenarios, we simulate input vectors with the dimension of 5, 10, 15 and 20.
For each generated dataset, we train a structured prediction energy network introduced in Gygli et al. (2017). As it has the sufficient ability to learn energy functions in equation 13, we can assume it has successfully captured the important features with a negligible error rate. We adopt each interpretation techniques to explain trained energy networks. According to equation 13 first four features affect the value of outputs. Fig. 3 compares the accuracy of results obtained by each method during the interpretation. Diagrams of Fig. 3 show results for target outputs y1 and y2 in E1 and two arbitrary outputs y3 and y4 in E2. As there may be randomness in interpretation methods, we run each interpreter five times for each dataset. Each line in the diagrams shows the average value, and the highlighted area shows the standard deviation. SOInter has an overall better performance compared to others.
As the accuracy measures the exact match of the subset of important features with ground truths, we consider an order for features and report the median rank of important ones in Fig. 4 as proposed in Chen et al. (2018). As the first four features are the solution, the desired median rank is 2.5 in all situations. As shown, SOInter has the nearest median rank to 2.5 nearly in all cases.
As the number of input features is increased, the performance of methods is generally degraded. This is because the ratio of important features compared to the size of the input vector is decreased, which can lead to confusion of the interpreter.
However, the obtained results confirm the robustness of SOInter for the more significant number of input features. Thus the proposed method is more reliable when the size of the input vector is large.

heading: 4.2 MULTI-LABEL CLASSIFICATION ON BIBTEX DATASET

text: Bibtex is a standard dataset for the multi-label classification of texts. Each sample in Bibtex involves an input feature vector corresponding to 1836 words mapped to a 159-dimensional output vector. Elements of the output vector are associated with a set of tags that describes the sample subject. We train a structured prediction energy network (SPEN) as a multi-label classifier on Bibtex with a desirable accuracy as shown in Gygli et al. (2017). A SPEN as a structured black box is a challenging benchmark for an interpreter because of its ability to capture more complicated relations between output variables. During interpreting this classifier with SOInter, we select an output variable, i.e., a tag, as a target of explanation and find the top 30 features related to this tag for each sample. According to SPEN decisions, we aggregate those top features over all samples for each tag and find the top 30 features which are expected to be correlated to this tag. Table 1 shows the general top 30 features for different 3 tags. More results are provided in Table 3 of Appendix A. As shown, word sets are meaningfully correlated to their corresponding tags. In addition, we highlight bold words that are correlated to each tag confirmed by human experts.

heading: 4.3 IMAGE SEGMENTATION ON WEIZMANN-HORSE DATASET

text:
Image segmentation is another structured output learning task in which the image is partitioned into semantic regions. Here we again train a SPEN for segmentation of 24×24 Weizmann-horse dataset. Each image of this dataset is partitioned into two regions which determines the borders of the horse. During the interpretation, we consider a pixel as the target and find pixels of the image that affect the target’s output.
In Fig. 5 the pixel in [10, 10] is considered as target and the interpretation results for arbitrary images are shown. We do experiments for different numbers of important features of 5, 10, 50, and 100. The red pixel shows the target, and green ones are obtained important input features as expected green pixels are placed in the locality of the target.
Table 2 compares the accuracy and post-accuracy of the SPEN in finding the label of the target during the segmentation process. Post-accuracy measures the model performance when only the selected features are passed through the SPEN and un-selected ones are masked with zero. These measures are shown for different choices for the number of important features. The negligible difference between the accuracy and post-accuracy confirms that the interpreter selects pixels which are the most effective ones in labeling the target. In addition, the increased value of the post-accuracy in some situations shows that omitted features confuse the segmentation model, i.e. SPEN, as a noise signal. These observations show that SOInter successfully selects important features which affect the SPEN decision for the target output.

heading: 5 CONCLUSION

text: We have presented SOInter, an interpreter for explaining structured output models. We focused on a single output variable of a structured model, available as a black box, as the target. Then we train a function over the input space, which returns a subset of important features for the black box to decide on the target. This is the first time an interpreter has been designed explicitly for structured output models to the best of our knowledge. These models learn complex relations between output variables which ignoring them while interpreting a single output can decline the explanation performance. We used an energy model to learn the structural information of the black box and utilize it during the interpreter’s training. The effectiveness of SO-Inter is confirmed using synthetic and real structured datasets.

heading: A APPENDIX

text:
Algorithm 1 The training procedure of the interpreter IN t(x;α) 1: function INTERPRETER-BLOCK(x, α) 2: return Gumbel-Softmax (Wα(x)) 3: end function 4: 5: function SL-CALCULATOR(x, ỹ, α) 6: x̃← x Interpreter-Block(x, α) 7: ysb ← arg maxy psb(y|x) 8: L ← |ỹt − ysbt | 9: SL ← max{0, Eq(x̃, [ ysbt , ỹ−t ] )− Eq(x̃, [ỹt, ỹ−t]) + L} 10: return SL 11: end function 12: 13: function TRAINING-INTERPRETER(t) 14: α(0) ← Sample Standard Normal distribution 15: ỹi
(0) ← arg maxy psb(y|x Interpreter-Block(xi, α(0)), for i = 1...N 16: k = 1 17: while the penalty does not change considerably do 18: α(k) = Adam-minimizerα ΣNi SL-calculator(xi, ỹi
(k−1), α) 19: ỹi
(k) = Adam-minimizery − psb(y|xi Interpreter-Block(xi, α(k))), for i = 1...N 20: k = k + 1 21: end while 22: end function

title: SOINTER: A NOVEL DEEP ENERGY-BASED INTERPRETATION METHOD FOR EXPLAINING STRUCTURED OUTPUT MODELS