Title: EXPLOITING INVARIANT STRUCTURES FOR COMPRES-SION IN NEURAL NETWORKS

Abstract: Modern neural networks often require deep compositions of high-dimensional nonlinear functions (wide architecture) to achieve high test accuracy, and thus can have overwhelming number of parameters. Repeated high cost in prediction at test-time makes neural networks ill-suited for devices with constrained memory or computational power. We introduce an efficient mechanism, reshaped tensor decomposition, to compress neural networks by exploiting three types of invariant structures: periodicity, modulation and low rank. Our reshaped tensor decomposition method exploits such invariance structures using a technique called tensorization (reshaping the layers into higher-order tensors) combined with higher order tensor decompositions on top of the tensorized layers. Our compression method improves low rank approximation methods and can be incorporated to (is complementary to) most of the existing compression methods for neural networks to achieve better compression. Experiments on LeNet-5 (MNIST), ResNet-32 (CI-FAR10) and ResNet-50 (ImageNet) demonstrate that our reshaped tensor decomposition outperforms (5% test accuracy improvement universally on CIFAR10) the state-of-the-art low-rank approximation techniques under same compression rate, besides achieving orders of magnitude faster convergence rates.

1. INTRODUCTION:

Modern neural networks achieve unprecedented accuracy over many difficult learning problems at the cost of deeper and wider architectures with overwhelming number of model parameters. The large number of model parameters causes repeated high cost in test-time as predictions require loading the network into the memory and repeatedly passing the unseen examples through the large network. Therefore, the model size becomes a practical bottleneck when neural networks are deployed on constrained devices, such as smartphones and IoT cameras.

Compressing a successful large network (i.e., reducing the number of parameters), while maintaining its performance, is non-trivial. Many approaches have been employed, including pruning, quantization, encoding and knowledge distillation (see appendix A for a detailed survey). A complementary compression technique, on top of which the aforementioned approaches can be used, is low rank approximation. For instance, singular value decomposition (SVD) can be performed on fully connected layers (weights matrices) and tensor decomposition on convolutional layers (convolutional kernels). Low rank approximation methods can work well and reduce the number of parameters by a factor polynomial in the dimension only when the weight matrices or convolutional kernels have low rank structures, which might not always hold in practice.

We propose to exploit additional invariant structures in the neural network for compression. A set of experiments on several benchmark datasets justified our conjecture (Section 4): large neural networks have some invariant structures, namely periodicity, modulation and low rank, which make part of the parameters redundant. Consider this toy example of a vector with periodic structure [1, 2, 3, 1, 2, 3, 1, 2, 3] or modulated structure [1, 1, 1, 2, 2, 2, 3, 3, 3] in Figure 1 . The number of parameters needed to represent this vector, naively, is 9. However if we map or reshape the vector into a higher order object, for instance, a matrix [1,1,1;2,2,2;3,3,3] where the columns of the matrix are repeated, then apparently this reshaped matrix can be decomposed into rank one without losing information. Therefore only 6 parameters are needed to represent the original length-9 vector. Although the invariant structures in large neural networks allow compression of redundant parameters, designing a sophisticated way of storing a minimal representation of the parameters (while maintaining the expressive power of the network) is nontrivial. To solve this problem, we proposed a new framework called reshaped tensor decomposition (RTD) which has three phases:

1. Tensorization. We reshape the neural network layers into higher-order tensors.

• For instance, consider a special square tensor convolutional kernel T ∈ R D×D×D×D , we reshape T into a higher m-order tensor

T ′ ∈ R m D ′ × • • • × D

′ , so called tensorization. Tensorization guarantees that the reshaped tensor has a smaller dimension, i.e., D ′ < D ( as D ′ m = D 4 and m > 4). Similar ideas apply to general convolutional kernels.

2. Higher-order tensor decomposition. We deploy tensor decomposition (a low rank approximation technique detailed in section 3) on the tensorized layers to exploit the periodic, modulated as well as low rank structures in the original layers.

• A rank-R tensor decomposition of the above 4-order tensor T will result in R number of components (each contains 4D parameters), and thus 4DR number of parameters in totalsmaller than the original D 4 number of parameters if R is small. • A rank-R tensor decomposition of the above reshaped m-order kernel tensor T ′ maps the layer into m + 1 narrower layers. The decomposition will result in R number of components with mD 4 m parameters and thus mD 4 m R in total -better than the 4DR number of parameters required by doing tensor decomposition on the original tensor T (D is usually large). Now the weights of the tensorized neural networks are the components of the tensor, i.e., result of the tensor decomposition. However, decomposing higher order tensors is challenging and known methods are not guaranteed to converge to the minimum error decomposition (Hillar & Lim, 2013) . Therefore fine tuning is needed to achieve high performance.

3. Data reconstruction-based sequential tuning. We fine-tune the parameters using a data reconstruction-based sequential tuning (Seq) method which minimizes the difference between training output of the uncompressed and compressed, layer by layer. Our Seq tuning is a novel approach inspired by a sequential training method proved to converge faster and achieve guaranteed accuracy using a boosting framework (Huang et al., 2017) . Unlike traditional end-to-end (E2E) backpropagation through the entire network, Seq tunes individual compressed "blocks" one at a time, reducing the memory and complexity required during compression.

Summary of Contributions:

• Novel compression schemes. We propose new reshaped tensor decomposition methods to exploit invariant structures for compressing the parameters in neural networks. By first tensorizing the kernel/weights into a higher-order tensor, our reshaped tensor decomposition discovers extra invariant structures and therefore outperform existing "low rank approximation methods". • Efficient computational framework. We introduce a system of tensor algebra that enables efficient training and inference for our compressed models. We show that a tensor decomposition on the parameters is equivalent to transforming one layer into multiple narrower sublayers in the compressed model. Therefore, other compression techniques (e.g. pruning) can be applied on top of our method by further compressing the sublayers returned by our method. • Sequential knowledge distillation. We introduce Seq tuning to transfer knowledge to a compressed network from its uncompressed counterpart by minimizing the data reconstruction error block by block. With our strategy, only one block of the network is loaded into the GPU "at each time", therefore allowing compression of large networks on moderate devices. Furthermore, we show empirically that our strategy converges much faster than normal end to end tuning.

• Comprehensive experiments. We perform extensive experiments to demonstrate that our reshaped tensor decomposition outperforms state-of-the-art low-rank approximation techniques (obtains 5% higher accuracy on CIFAR10 under same compression rates). Our experiments also show that our method scales to deep residual neural networks on large benchmark dataset, ImageNet.

Organization of the paper Section 2 introduces tensor operations, tensor decompositions and their representations in tensor diagrams. In Section 3, we introduce convolutional layer diagram, review existing low-rank approximation techniques, and propose three new schemes to exploit additional invariant structures. In Section 4 and Appendix B, we demonstrate by extensive experiments that our compression obtains higher accuracy than existing low-rank approximation techniques. Appendix A surveys compression techniques and discuss how our method is related or complementary to existing techniques. For simplicity, we will use tensor diagrams throughout the text. However we provide a detailed appendix where the tensor operations are mathematically defined. Notations An m-dimensional array T is defined as an m-order tensor T ∈ R I0×•••×Im-1 . Its (i 0 , • • • , i n-1 , i n+1 , • • • , i m-1 ) th mode-n fiber, a vector along the n th axis, is denoted as

2. TENSOR PRELIMINARIES:

T i0,••• ,in-1,:,in+1,••• ,im-1 .

Tensor Diagrams. Following the convention in quantum physics Cichocki et al. (2016) , Figure 2 introduces graphical representations for multi-dimensional objects. In tensor diagrams, an array (scalar/vector/matrix/tensor) is represented as a node in the graph, and its order is denoted by the number of edges extending from the node, where each edge corresponds to one mode (whose dimension is denoted by the number associated to the edge) of the multi-dimensional array.

= X Y T 1 I0 I1 I1 I2 I2 J0 J0 J1 J2 J2 I0 = J1 (a)

Mode-(0,1) tensor contraction X × 0 1 Y → T 1 T 1 i 1 ,i 2 ,j 0 ,j 2 = r Xr,i 1 ,i 2 Yj 0 ,r,j 2

= X M T 2 I0 I1 I1 I2 I2 J1 J1 J0 I0 = J0 (b)

Mode-0 tensor product/multiplication X ×0 M → T 2 : T 2 j 0 ,i 1 ,i 2 = r Xr,i 1 ,i 2 Mj 0 ,r

= X Y T 3 I0 I1 I1 I2 I2 J0 J0 J1 J2 J2 I ′ 0 (c)

Mode-(0,1) tensor convolution X * 0 1 Y → T 3 T 3 :,i 1 ,i 2 ,j 0 ,j 2 = X:,i 1 ,i 2 * Yj 0 ,:,j 2

= X Y T 4 I0 I1 I1 I2 I2 J0 J0 J1 J2 J2 I0 I0 = J1 (d)

Mode-(0,1) tensor partial outer product X ⊗ 0 1 Y → T 4 : T 4 r,i 1 ,i 2 ,j 0 ,j 2 = Xr,i 1 ,i 2 Yj 0 ,r,j 2 

T 1 ∈ R I1×I2×J0×J2 , T 2 ∈ R J0×I1×I2 , T 3 ∈ R I ′ 0 ×I1×I2×J0×J2

and T 4 ∈ R I0×I1×I2×J0×J2 are output tensors of corresponding operations. Similar definitions apply to general mode-(i, j) tensor operations.

Tensor Operations. In Figure 3 , we use some simple examples to introduce four types of tensor operations, which are higher-order generalization of their matrix/vector counterparts, on input tensors X and Y and input matrix M. In tensor diagram, an operation is represented by linking edges from the input tensors, where the type of operation is denoted by the shape of line that connects the nodes: solid line stands for tensor contraction / tensor multiplication, dashed line represents tensor convolution, and curved line is for tensor partial outer product. The rigorous definitions of high-order general tensor operations are defined in Appendix D.

Tensor Decompositions. We introduce generalized tensor decomposition as the reverse mapping of the general tensor operations (detailed in Appendix F): given a set of operations and a tensor, the generalized tensor decomposition recovers the factors/components such that the operations on these factors result in a tensor approximately equal to the original one. Several classical types of tensor decompositions (such as CANDECOMP/PARAFAC (CP), Tucker (TK) and Tensor-train (TT) decompositions) are introduced in Appendix F, and their applications on the convolutional kernel in Figure 4a (defined in Section 3) are illustrated as tensor diagrams in Figures 4b, 4c and 4d . 

K H W S T (a) Original kernel K C K S K T H W S T R (b) CP K S K C K T H W S T Rs Rt (c) Tucker (TK) K S K H K W K T H W S T Rs R Rt (d) Tensor-train (TT) K ′ H W S 0 S 1 S 2 T 0 T 1 T 2 (e) Tensorized kernel K C K 0 K 1 K 2 H W S0 S1 S2 T0 T1 T2 R (f) r-CP (m=3) C P 0 P 1 P 2 Q 0 Q 1 Q 2 H W S0 S0 S1 S2 T0 T1 T2 R s 0 R s 1 R s 2 R t 0 R t 1 R t 2 (g) r-TK (m=3) K 0 K 1 K 2 K C H W S0 S1 S1 S2 T0 T2 R0 R1 R2 (h) r-TT (m=3)

3. COMPRESSING CONVOLUTIONAL LAYER BY TENSOR DECOMPOSITIONS:

A standard convolutional layer in neural networks is parameterized by a 4-order kernel K ∈ R H×W ×S×T where H, W are height/width of the filters, and S, T are the numbers of input/output channels. The layer maps a 3-order input tensor U ∈ R X×Y ×S (with S number of feature maps of height X and width Y ) to another 3-order output tensor V ∈ R X ′ ×Y ′ ×T (with T number of feature maps of height X ′ and width Y ′ ) according to the following equation:

V x,y,t = S-1 s=0 i,j K i,j,s,t U i+dx,j+dy,s (3.1)

where d is the stride of the convolution. With HW ST parameters, it takes O(HW ST XY ) operations (FLOPs) to compute the output V. The diagram of the convolutional layer is in Figure 5a .

Plain Tensor Decomposition (PD) Traditional techniques compress a convolutional layer by directly factorizing the kernel K using tensor decompositions Jaderberg et al. (2014) ; Lebedev et al. (2014) ; Kim et al. (2015) , such as CANDECOMP/PARAFAC (CP), Tucker (TK) and Tensor-train (TT) decompositions. For example, consider a Tensor-train decomposition on K, the kernel can be factorized and stored as 4d . The decomposition is rigorously defined element-wisely as

K S ∈ R S×Rs , K H ∈ R Rs×H×R , K W ∈ R R×W ×Rt and K T ∈ R Rt×T , which only requires (SR s + HR s R + W R t R + T R t ) parameters as illustrated in Figure

K i,j,s,t = Rs-1 rs=0 R-1 r=0 Rt-1 rt=0 K S s,rs K H rs,i,r K W r,j,rt K T rt,t (3.2)

We defer the details of using CP and TK to Appendix G, although their tensor diagrams are illustrated in Figures 4b, 4c and 4d and their complexities are summarized in Tables 1 and 10 . 

EQUATION

Motivation of Tensorization -Invariant Structures:

Consider a matrix M ∈ R L 2 ×L 2 with M = a ⊗ b, a, b ∈ R L 2 , where a = [c; • • • ; c] is periodic repetition of the same vector c ∈ R L and b ⊤ = [d 0 1 ⊤ ; d 1 1 ⊤ ; • • • ; d L-1 1 ⊤ ] is a modulated version of another vector d ∈ R L .

Obviously, M as a rank-1 matrix can be represented by two length-L 2 vectors a and b, resulting in a total of 2L 2 parameters. However, if we reshape the matrix M into a 4-order tensor T ∈ R L×L×L×L , it can be factorized by CP decomposition as

T = 1 ⊗ c ⊗ d ⊗ 1 (1 ∈ R L )

, and represented by four length-L vectors, requiring only 4L parameters. We refer the process of reshaping an array into a higher-order tensor as tensorization, and the use of tensor decomposition following tensorization as reshaped tensor decomposition (RTD). Therefore, the example above demonstrates that RTD discovers additional invariant structures that baseline plain tensor decomposition (PD) fails to identify.

Reshaped Tensor Decomposition (RTD) Inspired by this intuition, we tensorize the convolutional kernel K into a higher-order tensor K ′ ∈ H×W ×S0× 

U ′ ∈ R X×Y ×S0×•••×Sm-1 and V ′ ∈ R X ′ ×Y ′ ×T0×•••×Tm-1 . V ′ x,y,t0,••• ,tm-1 = S0-1 s0=0 • • • Sm-1-1 sm-1=0 i,j K ′ i,j,s0••• ,sm-1,t0,••• ,tm-1 U ′ i+dx,j+dy,s0,••• ,sm-1 (3.3)

Now we can compress the convolutional layer by factorizing the tensorized kernel K ′ by tensor decompositions, and name the schemes using CP, Tucker and Tensor-train as reshaped CP (r-CP), reshaped Tucker (r-TK) and reshaped Tensor-train (r-TT) respectively. For example, consider a r-TT decomposition on K ′ , the tensorized kernel can now be stored in 4h ). The decomposition scheme is rigorously defined element-wisely as

m+1 factors {K 0 , • • • , K m-1 , K C } (a special example of m = 3 is illustrated in Figure

K ′ i,j,s0,••• ,sm-1,t0,••• ,tm-1 = R0-1 r0=1 • • • Rm-1-1 rm-1=0 K 0 s0,r0,t0 K 1 r0,s1,t1,r1 • • • K C i,j,rm-1 (3.4)

where

K l ∈ R R l-1 ×S l ×T l ×R l , ∀l ∈ [m] and K C ∈ R Rm-1×H×W . Assuming S l = S 1 m , T l = T 1 m , R l = R, ∀l ∈ [m], we require O(m(ST ) 1 m R 2 + HW R) parameters.

We defer the detailed descriptions of r-CP and r-TK to Appendix I, but we illustrate their tensor diagrams in Figures 4f and 4g and summarize their complexities in Tables 1 and 12 .

Sequential Tuning Tensor decompositions provide weight estimates in the tensorized convolutional layers. However, decomposing higher order tensors is challenging and known methods are not guaranteed to converge to the minimum error decompositions (Hillar & Lim, 2013) . Therefore fine tuning is needed to restore high performance. Analogous to Huang et al. (2017) , our strategy of data reconstruction-based sequential tuning (Seq) sequentially fine-tunes the parameters layer by layer, using backpropagation to minimize the difference between the outputs from uncompressed layer V ′ and the tensorized compressed layer V ′ .

Computational Complexity As shown in Figure 5c , reshaped tensor decomposition maps a layer into multiple (and thus deeper) narrower layers, each of which has width R l that is usually smaller than the original width T . We design efficient algorithms for forward/backward propagations for prediction/fine-tuning on these modified layers using tensor algebra. (1) Forward propagation: computing the output V ′ , given the input U ′ and kernel factors {K l } m-1 l=0 .

(2) Backward propagation: computing the derivative of loss function L with respect to (w.r.t.) the input ∂L/∂U ′ and kernel factors {∂L/∂K l } m-1 l=0 , given the derivative w.r.t. the output ∂L/∂V ′ . A naive forward and backward propagation mechanism is to explicitly reconstruct the original kernel K ′ using the factors {K l } m-1 l=0 , which however makes propagations highly inefficient as shown in Appendix F, G and I. Alternatively, we propose a framework where both propagations are evaluated efficiently without explicitly forming or computing the original kernel. The key idea is to interact the input U ′ with each of the factors K l individually. Taking r-TT as an example, we plug the decomposition 3.4 into 3.3, then the computation of V ′ is reduced into m + 1 steps:

U l x,y,••• ,r l = R l-1 -1 r l-1 =0 S l -1 s l =0 K l-1 r l-1 ,s l ,t l ,r l U l-1 x,y,s l ,••• ,sm-1,t0,••• ,t l-1 ,r l-1 , ∀l = 1, • • • , m (3.5) V x,y,t0,••• ,tm-1 = Rm-1-1 rm-1=0 i,j K m rm-1,i,j U m i+dx,j+dy,t0,••• ,tm-1,rm-1 (3.6)

where U l is the intermediate result after interacting with K l-1 , and U 0 = U. Each step in 3.5 takes O(max(S, T ) 1+ 1 m RXY ) operations, while the last step in 3.6 requires O(HW T RXY ) 1 . Therefore, the time complexity for the forward pass is O((m max(S, T ) 1+ 1 m R + HW T )RXY ), more efficient than uncompressed O(HW ST XY ) as R ≪ S, T . Backpropagation is derived and analyzed in Appendix I, and the analyses of other decomposition are in Appendix G and I, but we summarize their computational complexities in Table 1 , 10 and 12.

Parallel Computational Complexity Tensor algebra allows us to implement the propagations 3.5 and 3.6 in parallel given enough computational resources, further speeding up prediction. The parallel tine complexities of prediction2 with our RTD implementation is displayed in Table 2 . The prediction time complexity of RTD outperforms the baseline PD, whereas the PD outperforms the original convolutional layers as R ≪ N and m ≥ 3.

Decomp.:

O(# of parameters) O(# of forward ops.) 

O(# of backward ops.) original k 2 N 2 k 2 N 2 D 2 N 2 D 4 CP (k 2 + 2N )R (k 2 + 2N )RD 2 (D 2 + 2N )RD 2 TK (k 2 R + 2N )R (k 2 R + 2N )RD 2 (D 2 R + 2N )RD 2 TT 2(kR + N )R 2(kR + N )RD 2 2(DR + N )RD 2 r-CP (k 2 + mN 2 m )R (mN 1+ 1 m + k 2 N )RD 2 (mN 1+ 1 m + N D 2 )RD 2 r-TK (k 2 R 2m-1 + 2mN )R (k 2 R 2m-1 + 2mN )RD 2 (R 2m-1 D 2 + 2mN )RD 2 r-TT (mN 2 m R + k 2 )R (mN 1+ 1 m R + k 2 N )RD 2 (mN 1+ 1 m R + N D 2 )RD 2

X = Y = X ′ = Y ′ = D, S = T = N , H = W = k and D ≫ k).

General settings are summarized in Tables 10 and 12 . 

Decomp. O(parallel complexity) Decomp. O(parallel complexity) CP N + k 2 + R r-CP mN 1 m + k 2 R TK N + k 2 R + R r-TK 2mN 1 m + k 2 R m TT N + 2kR + R r-TT mN 1 m R + k 2 R

4. EXPERIMENTS:

We evaluate our reshaped tensor decomposition method on the state-of-art networks for a set of benchmark datasets: we evaluate convolutional layer compression on ResNet-32 He et al. (2016) for CIFAR-10; we evaluate fully-connected layer compression on MNIST; and we evaluate the scalability of our compression method on ResNet-50 He et al. (2016) for ImageNet (2012) dataset. The baseline we compare against is the state-of-the-art low-rank approximation methods called plain tensor decomposition (PD), as other compression methods are complementary and can be used on top of our reshaped tensor decomposition (RTD) method. All types of tensor decomposition (CP, TK, and TT) in baseline PD will be evaluated and compared with corresponding types of tensor decomposition (r-CP, r-TK and r-TT) in our RTD method.

Our primary contribution is to introduce a new framework, reshaped tensor decomposition, that picks out additional invariance structure such as periodicity and modulation, which the low rank approximation baseline, plain tensor decomposition, fails to find. Now we demonstrate that our RTD maintains high accuracy even when the networks are highly compressed on CIFAR-10. We refer to traditional backpropogation-based tuning of the network as end-to-end (E2E) tuning, and to our proposed approach that trains each block individually as data reconstruction-based sequential (Seq) tunning.

Compression Our algorithm achieves 5% higher accuracy than baseline on ResNet-34 CIFAR10. As in Table 3, using baseline CP decomposition with end-to-end tuning, ResNet-34 is compressed to 10% of its original size, reducing the accuracy from 93.2% to 86.93%. Our reshaped tensor decomposition using r-CP, paired with Seq tuning, increases the accuracy to 91.28% with the same 10% compression rate -a performance loss of 2% with only 10% of the number of parameters. It achieves further aggressive compression -a performance loss of 6% with only 2% of the number of parameters. We observe similar trends (higher compression and higher accuracy) for Tensor-train decomposition. The structure of the Tucker decomposition (see section I) makes it less effective with very high compression, since the "internal structure" of the network reduces to very low rank, which may lose necessary information. Increasing the network size to 20% of the original provides reasonable performance on CIFAR-10 for Tucker as well. 

5. CONCLUSION AND PERSPECTIVES:

We describe an efficient mechanism for compressing neural networks by tensorizing network layers. We implement tensorized decompositions to find approximations of the tensorized kernel, potentially preserving invariance structures missed by implementing decompositions on the original kernels. We extend vector/matrix operations to their higher order tensor counterparts, providing systematic notations and libraries for tensorization of neural networks and higher order tensor decompositions. As a future step, we will explore optimizing the parallel implementations of the tensor algebra.

Appendix: Exploiting Invariant Structures for Compression in Neural Networks:

A RELATED WORKS A recent survey Cheng et al. (2017) reviews state-of-the-art techniques for compressing neural networks, in which they group the methods into four categories: (1) low-rank factorization;

(2) design of compact filters; (3) knowledge distillation; and 4) parameters pruning, quantization and encoding. Generally, our decomposition schemes fall into the category of low-rank factorization, but these schemes also naturally lead to novel designs of compact filters in the sublayers of the compressed network. On the other hand, our strategy of sequential tuning is an advanced scheme of knowledge distillation that transfers information from pre-trained teacher network to compressed student network block by block. Furthermore, our method is complementary to the techniques of parameters pruning, quantization and encoding, which can be applied on top of our method by further compressing the parameters in the sublayers returned by tensor decomposition.

• Low-rank Factorization. Low-rank approximation techniques have been used for a long time to reduce the number of parameters in both fully connected and convolutional layers. Pioneering papers propose to flatten/unfold the parameters in convolutional layer into matrices (a.k.a matricization), followed by (sparse) dictionary learning or matrix decomposition (Jaderberg et al., 2014; Denton et al., 2014; Zhang et al., 2015) . Subsequently in Lebedev et al. (2014) ; Kim et al. (2015) , the authors show that it is possible to compress the tensor of parameters directly by standard tensor decomposition (in particular CP or Tucker decomposition). The groundbreaking work (Novikov et al., 2015) demonstrates that the parameters in fully connected layer can be efficiently compressed by tensor decomposition by first reshaping the matrix of parameters into higher-order tensor, and the idea is later extended to compress LSTM and GRU layers in recurrent neural networks (Yang et al., 2017) . Concurrent to Wang et al. (2018) , our paper extends this basic idea to convolutional layer by exploiting the invariant structures among the filters. Different from Wang et al. (2018) that only focuses Tensor-ring decomposition, we investigates, analyzes and implements a boarder range of decomposition schemes, besides other benefits discussed below. • Design of Compact Filters. These techniques reduce the number of parameters by imposing additional constraints on linear layers (fully connected or convolutional). For example, the matrix of parameters in fully connected layer is restricted to circular (Cheng et al., 2015) , Toeplitz/Vandermonde/Cauchy (Sindhwani et al., 2015) , or multiplication of special matrices (Yang et al., 2015) . Historically, convolutional layer is proposed as a compact design of fully connected layer, where spatial connections are local (thus sparse) with repeated weights.

Recent research further suggests to use more compact convolutional layers, such as 1 × 1 convolutional layer (Szegedy et al., 2017; Wu et al., 2017) (where each filter is simply a scalar), and depthwise convolutional layer (Chollet, 2016) (where connections between the feature maps are also sparse). In our paper, we show that the sublayers returned by our decomposition schemes are in fact 1 × 1 depthwise convolutional layers, combing advantages from both designs above. • Knowledge Distillation. The algorithms of knowledge distillation aim to transfer information from a pre-trained teacher network to a smaller student network. In Ba & Caruana (2014) ; Hinton et al. (2015) , the authors propose to train the student network supervised by the logits (the vector before softmax layer) of the teacher network. Romero et al. (2014) extends the idea to matching the outputs from both networks at each layer, up to an affine transformation. Our Seq tuning strategy is therefore similar to Romero et al. (2014) , but we use identical mapping instead of affine transformation, and train the compressed network block by block. • Pruning, Quantization and Encoding. Han et al. (2015) proposes a three-step pipeline to compress a pre-trained network, by (1) pruning uninformative connections, (2) quantizing the remaining weights and (3) encoding the discretized parameters. Since our decomposition schemes effectively transform one layer in the original network into the multiple sublayers, this pipeline can be applied by further compressing all sublayers. Therefore, our method is complementary (and can be used independently) to the techniques in this pipeline.

B SUPPLEMENTARY EXPERIMENTS:

Convergence Rate Compared to end-to-end, an ancillary benefit of Seq tuning is much faster and leads to more stable convergence. Figure 6 plots compression error over number of gradient updates for various methods. (This experiment is for PD with 10% compression rate.) There are three salient points: first, Seq tuning has very high error in the beginning while the "early" blocks of the network are being tuned (and the rest of the network is left unchanged to tensor decomposition values). However, as the final block is tuned (around 2 × 10 11 gradient updates) in the figure, the errors drop to nearly minimum immediately. In comparison, end-to-end tuning requires 50-100% more gradient updates to achieve stable performance. Finally, the result also shows that for each block, Seq tuning achieves convergence very quickly (and nearly monotonically), which results in the stair-step pattern since extra tuning of a block does not improve (or appreciably reduce) performance.

Performance on Fully-Connected Layers An extra advantage of reshaped tensor decomposition compression is that it can apply flexibly to fully-connected as well as convolutional layers of a neural network. Table 7 shows the results of applying reshaped tensor decomposition compression to various tensor decompositions on a variant of LeNet-5 network LeCun et al. (1998) . The convolutional layers of the LeNet-5 network were not compressed, trained or updated in these experiments. The uncompressed network achieves 99.31% accuracy. Table 7 shows the fully-connected layers can be compressed to 0.2% losing only about 2% accuracy. In fact, compressing the dense layers to 1% of their original size reduce accuracy by less then 1%, demonstrating the extreme efficacy of reshaped tensor decomposition compression when applied to fully-connected neural network layers. 

C NOTATIONS:

Symbols: Lower case letters (e.g. v) are used to denote column vectors, while upper case letters (e.g. M) are used for matrices, and curled letters (e.g. T ) for multi-dimensional arrays (tensors). For a tensor T ∈ R I0×•••×Im-1 , we will refer to the number of indices as order, each individual index as mode and the length at one mode as dimension. Therefore, we will say that T ∈ R I0×•••×Im-1 is an m-order tensor which has dimension I k at mode-k. Tensor operations are extensively used in this paper: The tensor (partial) outer product is denoted as ⊗, tensor convolution as * , and finally × denotes either tensor contraction or tensor multiplication. Each of these operators will be equipped with subscript and superscript when used in practice, for example × m n denotes mode-(m, n) tensor contraction (defined in Appendix D). Furthermore, the symbol • is used to construct compound operations. For example, ( * • ⊗) is a compound operator simultaneously performing tensor convolution and tensor partial outer product between two tensors.

Indexing: In this paragraph, we explain the usages of subscripts/superscripts for both multidimensional arrays and operators, and further introduce several functions that are used to alter the layout of multi-dimensional arrays.

• Nature indices start from 0, but reversed indices are used occasionally, which start from -1.

Therefore the first entry of a vector v is v 0 , while the last one is v -1 . • For multi-dimensional arrays, the subscript is used to denote an entry or a subarray within an object, while superscript is to index among a sequence of arrays. For example, M i,j denotes the entry at i th row and j th column of a matrix M, and M (k) is the k th matrix in a set of N matrices

{M (0) , M (1) , • • • M (N -1) }.

For operators, as we have seen, both subscript and superscript are used to denote the modes involved in the operation.

• The symbol colon ':' is used to slice a multi-dimensional array. For example, M :,k denotes the k th column of M, and T :,:,k denotes the k th frontal slice of a 3-order tensor T .

• Big-endian notation is adopted in conversion between multi-dimensional array and vectors. Specifically, the function vec(•) flattens (a.k.a. vectorize)

a tensor T ∈ R I0×•••×Im-1 into a vector v ∈ R m-1 l=0 I l such that T i0,••• ,im-1 = v im-1+im-2Im-1+•••+i0I1•••Im-1 .

• The function swapaxes(•) is used to permute ordering of the modes of a tensor as needed. For example, given two tensors U ∈ R I×J×K and V ∈ R K×J×I , the operation V = swapaxes(U) convert the tensor U into V such that V k,j,i = U i,j,k .

• The function flipaxis(•, •) flips a tensor along a given mode. For example, given a tensor U ∈ R I×J×K and V = flipaxis(U, 0), the entries in V is defined as

V i,j,k = U I-1-i(mod I),j,k . D TENSOR OPERATIONS Operator Notation Definition mode-(k, l) Tensor Contraction T (0) = X × k l Y T (0) i0,••• ,i k-1 ,i k+1 ,••• ,im-1,j0,••• ,j l-1 ,j l+1 ,••• ,jn-1 = X i0,••• ,i k-1 ,:,i k+1 ,••• ,im-1 , Y j0,••• ,j l-1 ,:,j l+1 ,••• ,jn-1 inner product of mode-k fiber of X and mode-l fiber of Y mode-k Tensor Multiplication T (1) = X × k M T (1) i0,••• ,i k-1 ,r,i k+1 ,••• ,im-1 = X i0,••• ,i k-1 ,:,i k+1 ,••• ,im-1 , M :,r inner product of mode-k fiber of X and r th column of M mode-(k, l) Tensor Convolution T (2) = X * k l Y T (2) i0,••• ,i k-1 ,:,i k+1 ,••• ,im-1,j0,••• ,j l-1 ,j l+1 ,••• ,jn-1 = X i0,••• ,i k-1 ,:,i k+1 ,••• ,im-1 * Y j0,••• ,j l-1 ,:,j l+1 ,••• ,jn-1 convolution of mode-k fiber of X and mode-l fiber of Y mode-(k, l) Partial- Outer Product T (3) = X ⊗ k l Y T (3) i0,••• ,i k-1 ,r,i k+1 ,••• ,im-1,j0,••• ,jn-1 = X i0,••• ,i k-1 ,r,i k+1 ,••• ,im-1 Y j0,••• ,j l-1 ,r,j l+1 ,••• ,jn-1

Hadamard product of mode-k fiber of X and mode-l fiber of Y 

∈ R I0×•••×Im-1 , Y ∈ R J0×•••×Jn-1 and matrix M ∈ R I k ×J . Mode-(k, l)

tensor contraction and mode-(k, l) tensor partial-outer product are legal only if 1) is an m-order tensor, T (2) is an (m + n -1)-order tensor and T (3) is an (m + n -1)-order tensor.

EQUATION

In this section, we introduce a number of tensor operations that serve as building blocks of tensorial neural networks. To begin with, we describe several basic tensor operations that are natural generalization to their vector/matrix counterparts. Despite their simplicity, these basic operations can be combined among themselves to construct complicated compound operators that are actually used in all designs. We will analyze their theoretical sequential time complexities in details, and point out the implementational concerns along the way. Although all these operations can in principle be implemented by parallel programs, the degree of parallelism depends on their particular software and hardware realizations. Therefore, we will use the sequential time complexity as a rough estimate of the computational expense in this paper.

Tensor contraction Given a m-order tensor T (0) ∈ R I0×•••×Im-1 and another n-order tensor

T (1) ∈ R J0×•••×Jn-1

, which share the same dimension at mode-k of T (0) and mode-l of T (1) ( i.e.

I k = J l ), the mode-(k, l) contraction of T (0) and T (1) , denoted as T T (0) × k l T (1) , returns a (m + n -2)-order tensor T ∈ R I0×•••×I k-1 ×I k+1 ×•••×Im-1×J0×•••×J l-1 ×J l+1 ×•••×Jn-1 , whose entries are computed as T i0,••• ,i k-1 ,i k+1 ,••• ,im-1,j0,••• ,j l-1 ,j l+1 ,••• ,jn-1 = I k -1 r=0 T (0) i0,••• ,i k-1 ,r,i k+1 ,••• ,im-1 T (1) j0,••• ,j l-1 ,r,j l+1 ,••• ,jn-1 (D.1a) = T (0) i0,••• ,i k-1 ,:,i k+1 ,••• ,im-1 , T (1) j0,••• ,j l-1 ,:,j l+1 ,••• ,jn-1 (D.1b)

Notice that tensor contraction is a direct generalization of matrix multiplication to higher-order tensor, and it reduces to matrix multiplication if both tensors are 2-order (and therefore matrices). As each entry in T can be computed as inner product of two vectors, which requires I k = J l multiplications, the total number of operations to evaluate a tensor contraction is therefore O((

m-1 u=0 I u )( n-1 v=0,v =l J v )), taking additions into account.

Notice that the analysis of time complexity only serves as a rough estimate for actual execution time, because we do not consider the factors of parallel computing and computer systems. In practice, (1) the modes that are not contracted over can be computed in parallel, and summations can be computed in logarithmic instead of linear time; (2) The spatial locality of the memory layout plays a key role in speeding up the computation of tensor operations. These arguments equally apply to all tensor operations in this paper, but we will not repeat them in the analyses for simplicity.

Tensor multiplication (Tensor product) Tensor multiplication (a.k.a. tensor product) is a special case of tensor contraction where the second operant is a matrix. Given a m-order tensor U ∈ R I0×•••×Im-1 and a matrix M ∈ R I k ×J , where the dimension of U at mode-k agrees with the number of the rows in M, the mode-k tensor multiplication of U and M, denoted as

V U × k M, yields another m-order tensor V ∈ R I0×•••×I k-1 ×J×I k+1 ×•••Im-1 , whose entries are computed as V i0,••• ,i k-1 ,j,i k+1 ,••• ,im-1 = I k -1 r=0 U i0,••• ,i k-1 ,r,i k+1 ,••• ,im-1 M r,j (D.2a) = U i0,••• ,i k-1 ,:,i k+1 ,••• ,im-1 , M :,j (D.2b)

Following the convention of multi-linear algebra, the mode for J now substitutes the location originally for I k (which is different from the definition of tensor contraction). Regardlessly, the number of operations for tensor multiplication follows tensor contraction exactly, that is O((

m-1 u=0 I u )J).

Tensor convolution Given a m-order tensor T (0) ∈ R I0×I1×•••×Im-1 and another n-order tensor

T (1) ∈ R J0×J1×•••×Jn-1 . The mode-(k, l) convolution of T (0) and T (1) , denoted as T T (0) * k l T (1) , returns a (m + n -1)-order tensor T ∈ R I0×•••×I ′ k ×•••×Im-1×J0×•••×J l-1 ×J l+1 ×•••×Jn-1 .

The entries of T can be computed using any convolution operation * that is defined for two vectors.

T i0,••• ,i k-1 ,:,i k+1 ,••• ,im-1,j0,••• ,j l-1 ,j l+1 ,••• ,jn-1 = T (0) i0,••• ,i k-1 ,:,i k+1 ,••• ,im-1 * T (1) j0,••• ,j l-1 ,:,j l+1 ,••• ,jn-1 (D.3a) = T (1) j0,••• ,j l-1 ,:,j l+1 ,••• ,jn-1 * T (0) i0,••• ,i k-1 ,:,i k+1 ,••• ,im-1 (D.3b)

Here we deliberately omit the exact definition of vector convolution * , as it can be defined in multiple forms depending on the user case (Interestingly, the "convolution" in convolutional layer indeed computes correlation instead of convolution). Correspondingly, the resulted dimension I ′ k at modek is determined by the chosen type of convolution. For example, the "convolution" in convolutional layer typically yields I ′ k = I k (with same padding) or I ′ k = I k -J l + 1 (with valid padding). Notice that vector convolution itself is generally asymmetric , i.e. u * v = v * u (except for the case of circular convolution). For convenience, we can define its conjugate as * such that u * v = v * u. With this notations, Equation D.3a can also be written as D.3b.

Generally speaking, Fast Fourier Transform (FFT) plays a critical role to lower the computational complexities for all types of convolution. In the case of tensor convolution, the number of required operations without FFT is O((

m-1 u=0 I u )( n-1 v=0 J v )), while FFT is able to reduce it to O(( m-1 u=0,u =k I u )( n-1 v=0,v =l J v ) max(I k , J l ) log(max(I k , J l ))

). That being said, FFT is not always necessary: if min(I k , J l ) < log (max(I k , J l )) (which is typical in convolutional layers, where I k is the height/width of the feature maps and J l is the side length of the square filters), computing the convolution without FFT is actually faster. Furthermore, FFT can be difficult to implement (thus not supported by popular software libraries) if convolution is fancily defined in neural networks (e.g. dilated, atrous). Therefore, we will assume that tensor convolutions are computed without FFT in subsequent sections unless otherwise noted.

Tensor outer product Given a m-order tensor T (0) ∈ R I0×I1×•••×Im-1 and another n-order tensor T (1) ∈ R J0×J1×•••×Jn-1 , the outer product of T (0) and T (1) , denoted T T (0) ⊗ T (1) , concatenates all the indices of T (0) and T (1) , and returns a

(m + n)-order tensor T ∈ R I0×•••×Im-1×J0×•••×Jn-1 whose entries are computed as T i0,••• ,im-1,j0,••• ,jn-1 = T (0) i0,••• ,im-1 T (1) j0,••• ,jn-1 (D.4)

It is not difficult to see that tensor outer product is a direct generalization for outer product for two vectors M = u⊗ v = u v ⊤ . Obviously, the number of operations to compute a tensor outer product explicitly is O((

m-1 u=0 I u )( n-1 v=0 J v )).

Tensor outer product is rarely calculated alone in practice because it requires significant amounts of computational and memory resources.

Tensor partial outer product Tensor partial outer product is a variant of tensor outer product defined above, which is widely used in conjunction with other operations. Given a m-order tensor T (0) ∈ R I0×I1×•••×Im-1 and another n-order tensor T (1) ∈ R J0×J1×•••×Jn-1 , which share the same dimension at mode-k of T (0) and mode-l of T (1) (i.e. I k = J l ), the mode-(k, l) partial outer product of T (0) and T (1) , denoted as T T

(0) ⊗ k l T (1) , returns a (m + n -1)-order tensor T ∈ R I0×•••×Im-1×J0×•••×J l-1 ×J l+1 ×•••×Jn-1 , whose entries are computed as T i0,••• ,i k-1 ,r,i k+1 ,••• ,im-1,j0,••• ,j l-1 ,j l+1 ,••• ,jn-1 = T (0) i0,••• ,i k-1 ,r,i k+1 ,••• ,im-1 T (1) j0,••• ,j l-1 ,r,j l+1 ,••• ,jn-1 (D.5a) = T (0) ••• ,r,••• ⊗ T (1) ••• ,r,••• (D.5b)

The operation bears the name "partial outer product" because it reduces to outer product once we fix the indices at mode-k of T (0) and mode-l of T (1) . Referring to the computational complexity of tensor outer product, the number of operations for each fixed index is

O(( m-1 u=0,u =k I u )( n-1 v=0,v =l J v ))

, therefore the total time complexity for the tensor partial outer product is O((

m-1 u=0 I u )( n-1 v=0,v =l J v ))

, the same as tensor contraction.

Precautions of usage:

• Similar to matrix multiplication, the operants in tensor operations are not commutative in general.

For example, neither 0) holds even if the dimensions at the specified modes happen to match.

T (0) × k l T (1) = T (1) × k l T (0) nor T (0) × k l T (1) = T (1) × l k T (

• Different from matrix multiplication, the law of associative also fails in general. For example,

(T (0) × k l T (1) ) × p q T (2) = T (0) × k l (T (1) × p q T (2)

), mainly because tensor operations can change the locations of modes in a tensor.

• However, both problems are not fundamental, and can be fixed by adjusting the superscripts and subscripts of the operators carefully (and further permute ordering of the modes in the result accordingly). For example, T (0) × k l T (1) = swapaxes(T (1) × l k T (0) ) holds if swapaxes(•) is properly performed. Due to space limits, we can not develop general rules in this paper, and will derive such identities as needed. In general, the take away message is a simple statement: Given an expression that contains multiple tensor operations, these operations need to be evaluated from left to right unless a bracket is explicitly supplied.

Compound operations: As building blocks, the basic tensor operations defined above can further combined to construct compound operations that perform multiple operations on multiple tensors simultaneously. For simplicity, we illustrate their usage using two representative examples in this section. More examples will arise naturally when we discuss the derivatives and backpropagation rules for compound operations in Appendix E.

• Simultaneous multi-operations between two tensors. For example, given two 3-order tensors T (0) ∈ R R×X×S and T (1) ∈ R R×H×S , we can define a compound operation ⊗ 0 0 • * 1 1 • × 2 2 between T (0) and T (1) , where mode-(0, 0) partial outer product, mode-(1, 1) convolution and mode-(2, 2) contraction are performed simultaneously, which results in a 2-order tensor T of R R×X ′ (it is indeed a matrix, though denoted as a tensor). The entries of T T

(0) ⊗ 0 0 • * 1 1 • × 2 2 T (1) are computed as T r,: = S-1 s=0 T (0) r,:,s * T (1)

r,:,s

(D.6)

For commonly used vector convolution, it is not difficult to show that number of operations required to compute the result T is O (R max(X, H) log(max(X, H))S) with FFT and O(RXHS) without FFT, as each of the R vectors in T is computed with a sum of S vector convolutions.

• Simultaneous operations between a tensor and a set of multiple tensors. For example, given a 3-order tensor U ∈ R R×X×S and a set of three tensors

T (0) ∈ R R×P , T (1) ∈ R K×Q and T (2) ∈ R S×T , we can define a compound operation on U as V U(⊗ 0 0 T (0) * 1 0 T (1) × 2 0 T (2)

), which performs mode-(0, 0) partial outer product with T (0) , mode-(1, 0) convolution with T (1) and mode-(2, 0) contraction with T (2) simultaneously. In this case, a 5-order tensor V ∈ R R×X ′ ×P ×Q×T is returned, with entries calculated as

V r,:,p,q,t = S-1 s=0 T (0) r,p U r,:,s * T (1) :,q T (2) s,t (D.7)

The analysis of time complexity of a compound operation with multiple tensors turns out to be a non-trivial problem. To see this, let us first follow the naive way to evaluate the output according to the expression above: (1) each vector in the result V r,:,p,q,t can be computed with a sum of S vector convolutions, which requires O(XHS) operations; (2) and with RP QT such vectors in the result V, the time complexity for the whole compound operation is therefore O(RXHP RST ). However, it is obviously not the best strategy to perform these operations. In fact, the equations can be equivalently rewritten as

V r,:,p,q,t = S-1 s=0 U r,:,s T (2) s,t * T (1) :,q T (0) r,p (D.8)

If we follows the supplied brackets and break the evaluation into three steps, it is not difficult to verify that these steps take O(RXST ), O(RXHP T ) and O(RX ′ HP T ) operations respectively, and result in a total time complexity of O(RXST + RXHP T + RX ′ P QT ) for the compound operation, which is far lower than the one with the naive way. Unfortunately, it is an NP-hard problem to determine the best order (with minimal number of operations) to evaluate a compound operation over multiple tensors, therefore in practice the order is either determined by exhaustive search (if there are only a few tensors) or follows a heuristic strategy (if the number of tensors is large).

The examples provided above are by no mean comprehensive, and in fact more complicated compound operations simultaneously perform multiple operations on multiple tensors can be defined, and we will see examples of them in the next section when we derive the backpropagation equations for the compound operations above. Generally, compound operations over multiple tensors are difficult to flatten into mathematical expressions without introducing tedious notations. Therefore, these operations are usually described by graphical representations, which are usually called tensor network in the physics literature (not to confuse with tensorial network in this paper). Interested readers are referred to the monograph Cichocki et al. (2016) , which serves a comprehensive introduction to the application of tensor network in the field of machine learning.

E DERIVATIVES AND BACKPROPAGATION OF TENSOR OPERATIONS:

All operations introduced in the last section, both basic and compound, are linear in their operants. Therefore, the derivatives of the result with respect to its inputs are in principle easy to calculate. In this section, we will explicitly derive the derivatives for all operations we have seen in Appendix D.

These derivatives can be further combined with classic chain rule to obtain the corresponding backpropagation equations (i.e. how gradient of the loss function propagates backward through tensor operations), which are the cornerstones of modern feed-forward neural networks. In the section, we show that these backpropagation equations can also be characterized by (compound) tensor operations, therefore their computational complexities can be analyzed similarly as in Appendix D.

Interestingly, the backpropagation equations associated with a tensor operation, though typically appear to be more involved, share the same asymptotic complexities as in the forward pass (with tensor convolution as an exception). This observation is extremely useful in the analyses of tensorial neural networks in Appendix G, H and I, which allows us to reuse the same number in the forward pass in the analysis of backward propagation.

In this section, we will assume for simplicity the loss function L is differentiable. However, all derivatives and backpropagation equations equally apply when L is only sub-differentiable (piecewise smooth). Also, we will focus on one step of backpropagation, therefore we assume the gradient of the loss function is known to us in prior.

Tensor contraction Recall the definition of tensor contraction in Equation D.1a, the partial derivatives of the result T with respect to its operants T (0) , T (1) can be computed at the entries level:

∂T i0,••• ,i k-1 ,i k+1 ,••• ,im-1,j0,••• ,j l-1 ,j l+1 ,••• ,jn-1 ∂T (0) i0,••• ,i k-1 ,r,i k+1 ,••• ,im-1 = T (1) j0,••• ,j l-1 ,r,j l+1 ,••• ,jn-1 (E.1a) ∂T i0,••• ,i k-1 ,i k+1 ,••• ,im-1,j0,••• ,j l-1 ,j l+1 ,••• ,jn-1 ∂T (1) j0,••• ,j l-1 ,r,j l+1 ,••• ,jn-1 = T (0) i0,••• ,i k-1 ,r,i k+1 ,••• ,im-1 (E.1b)

With classic chain rule, the derivatives of L with respect to T (0) and T (1) can be obtained through the derivative of Lwith respect to T .

∂L ∂T:

(0) i0,••• ,i k-1 ,r,i k+1 ,••• ,im-1 = J0-1 j0=0 • • • J l-1 -1 j l-1 =0 J l+1 -1 j l+1 =0 • • • Jn-1-1 jn-1=0 ∂L ∂T i0,••• ,i k-1 ,i k+1 ,••• ,im-1,j0,••• ,j l-1 ,j l+1 ,••• ,jn-1 T (1) j0,••• ,j l-1 ,r,j l+1 ,••• ,jn-1 (E.2a) ∂L ∂T (1) j0,••• ,j l-1 ,r,j l+1 ,••• ,jn-1 = I0-1 i0=0 • • • I k-1 -1 i k-1 =0 I k+1 -1 i k+1 =0 • • • Im-1-1 im-1=0 ∂L ∂T i0,••• ,i k-1 ,i k+1 ,••• ,im-1,j0,••• ,j l-1 ,j l+1 ,••• ,jn-1 T (0) i0,••• ,i k-1 ,r,i k+1 ,••• ,im-1 (E.2b)

Though tedious at entries level, it can be simplified with tensor notations in this paper.

(0) = swapaxes ∂L ∂T × m-1 0 • • • • • × m+l-2 l-1 • × m+l-1 l+1 • • • • • × m+n-3 n-1 T (1) (E.3a) ∂L ∂T (1) = swapaxes ∂L ∂T × 0 0 • • • • • × k-1 k-1 • × k k+1 • • • • • × m-2 m-1 T (0) (E.3b)

where swapaxes(•) is used to align the modes of outputs. Notice that the backpropagation equations are compound operations, even if the original operation is a basic one. It is not difficult to show that the number of operations required for both backpropagation equations are O((

m-1 u=0 I u )( n-1 v=0,v =l J v ))

, which are exactly the same as in the forward pass in Equation D.1a. The result should not surprise us however, since the tensor contraction is a direct generalization to matrix multiplication (where backward propagation has exactly the same time complexity as the matrix multiplication itself).

Tensor multiplication (Tensor product) As a special case of tensor contraction, the derivatives and backpropagation equations for tensor multiplication can be obtained in the same manner. To begin with, the derivatives of V with respect to U and M can be computed from the definition in Equation D.2a.

∂V i0,••• ,i k-1 ,j,i k+1 ,••• ,im-1 ∂U i0,••• ,i k-1 ,r,i k+1 ,••• ,im-1 = M r,j (E.4a) ∂V i0,••• ,i k-1 ,j,i k+1 ,••• ,im-1 ∂M r,j = U i0,••• ,i k-1 ,r,i k+1 ,••• ,im-1 (E.4b)

Subsequently, the derivatives of L with respect to U and M can be computed as with chain rule,

∂L ∂V i0,••• ,i k-1 ,r,i k+1 ,••• ,im-1 = J-1 j=0 ∂L ∂U i0,••• ,i k-1 ,j,i k+1 ,••• ,im-1 M r,j (E.5a) ∂L ∂M r,j = I0-1 i0=0 • • • I k-1 -1 i k-1 =0 I k+1 -1 i k+1 =0 • • • Im-1-1 im-1=0 ∂L ∂V i0,••• ,i k-1 ,j,i k+1 ,••• ,im-1 U i0,••• ,i k-1 ,r,i k+1 ,••• ,im-1 (E.5b)

Again, the backpropagation equations above can be succinctly written in tensor notations.

∂L ∂V = ∂L ∂U × k M ⊤ (E.6a) ∂L ∂M = U × 0 0 • • • • • × k-1 k-1 • × k+1 k+1 • • • • • × m-1 m-1 ∂L ∂V (E.6b)

where the time complexities for both equations are O(( m-1 u=0 I u )J), which is identical to the forward pass in Equation D.2a (obviously since tensor multiplication is a special of tensor contraction).

Tensor convolution Recall in the definition of tensor convolution in Equation D.3a, we deliberately omit the exact definition of vector convolution for generality. For simplicity, we temporarily limit ourselves to the special case of circular convolution. In this case, tensor convolution can be concretely defined by either equation below:

T i0,••• ,i k-1 ,:,i k+1 ,••• ,im-1,j0,••• ,j l-1 ,j l+1 ,••• ,jn-1 = Cir T (0) i0,••• ,i k-1 ,:,i k+1 ,••• ,im-1 T (1) j0,••• ,j l-1 ,:,j l+1 ,••• ,jn-1 (E.7a) T i0,••• ,i k-1 ,:,i k+1 ,••• ,im-1,j0,••• ,j l-1 ,j l+1 ,••• ,jn-1 = Cir T (1) j0,••• ,j l-1 ,:,j l+1 ,••• ,jn-1 T (0) i0,••• ,i k-1 ,:,i k+1 ,••• ,im-1 (E.7b)

where Cir(•) returns a circular matrix of the input vector. Concretely, given a vector v ∈ R I , the circular matrix Cir(v) is defined as Cir(v) i,j = v i-j(modI) . Now, the derivatives of the result tensor T with respect to T (0) and T (1) can be obtained by matrix calculus.

∂T i0,••• ,i k-1 ,:,i k+1 ,••• ,im-1,j0,••• ,j l-1 ,j l+1 ,••• ,jn-1 ∂T (0) i0,••• ,i k-1 ,:,i k+1 ,••• ,im-1 = Cir T (1) j0,••• ,j l-1 ,:,j l+1 ,••• ,jn-1 ⊤ (E.8a) ∂T i0,••• ,i k-1 ,:,i k+1 ,••• ,im-1,j0,••• ,j l-1 ,j l+1 ,••• ,jn-1 ∂T (1) j0,••• ,j l-1 ,:,j l+1 ,••• ,jn-1 = Cir T (0) i0,••• ,i k-1 ,:,i k+1 ,••• ,im-1 ⊤ (E.8b)

Applying chain rule to the equations above, we arrive at two lengthy equations:

∂L ∂T (0) i0,••• ,i k-1 ,:,i k+1 ,••• ,im-1 = J0-1 j0=0 • • • J l-1 -1 j l-1 =0 J l+1 -1 j l+1 =0 • • • Jn-1-1 jn-1=0 Cir T (1) j0,••• ,j l-1 ,:,j l+1 ,••• ,jn-1 ⊤ ∂L ∂T i0,••• ,i k-1 ,:,i k+1 ,••• ,im-1,j0,••• ,j l-1 ,j l+1 ,••• ,jn-1 (E.9a) ∂L ∂T (1) j0,••• ,j l-1 ,r,j l+1 ,••• ,jn-1 = I0-1 i0=0 • • • I k-1 -1 i k-1 =0 I k+1 -1 i k+1 =0 • • • Im-1-1 im-1=0 Cir T (0) i0,••• ,i k-1 ,:,i k+1 ,••• ,im-1 ⊤ ∂L ∂T i0,••• ,i k-1 ,:,i k+1 ,••• ,im-1,j0,••• ,j l-1 ,j l+1 ,••• ,jn-1 (E.9b)

With notations of tensor operations, they can be greatly simplified as

∂L ∂T (0) = swapaxes ∂L ∂T × m 0 • • • × m+l-1 l-1 • * k l • × m+l-2 l+1 • • • • • × m+n-2 n-1 flipaxis(T (1) , l) (E.10a) ∂L ∂T (1) = swapaxes ∂L ∂T × 0 0 • • • • • × k-1 k-1 • * k k • × k+1 k+1 • • • • • × m-1 m-1 flipaxis(T (0) , k) (E.10b)

Although these backpropagation equations are derived for the special case of circular convolution, they hold for general convolution if we replace * k l by its corresponding adjoint operator ( * k l ) ⊤ (and use the unflipped versions of T (0) and T (1) ). With adjoint of convolution, Equations E.10a and E.10b can be rewritten as

∂L ∂T (0) = swapaxes ∂L ∂T × m 0 • • • × m+l-1 l-1 •( * k l ) ⊤ • × m+l-2 l+1 • • • × m+n-2 n-1 T (1) (E.11a) ∂L ∂T (1) = swapaxes ∂L ∂T × 0 0 • • • • • × k-1 k-1 • ( * k k ) ⊤ • × k+1 k+1 • • • • • × m-1 m-1 T (0) (E.11b)

where the exact form of the adjoint operator ( * k l ) ⊤ depends on the original definition of vector convolution. Generally, the trick to start with circular convolution and generalize to general cases is very useful to derive backpropagation equations for operations that convolution plays a part.

Despite varieties in the definitions of tensor convolution, the analyses of their time complexities of backpropagation equations are identical, since the numbers of operations only differ by a constant for different definitions (therefore asymptotically the same). With FFT, the number of operations for these two backpropagation equations are O((

m-1 u=0,u =k I u )( n-1 v=0,v =l J v ) max(I ′ k , J l ) log(max(I ′ k , J l ))) and O(( m-1 u=0,u =k I u )( n-1 v=0,v =l J v ) max(I ′ k , I k ) log(max(I ′ k , I k ))), and with- out FFT O(( m-1 u=0,u =k I u )( n-1 v=0,v =l J v )I ′ k J l ) and O(( m-1 u=0,u =k I u )( n-1 v=0,v =l J v )I ′ k I k ) respec- tively.

Different from other operations, the time complexities for forward and backward passes are different (with circular convolution as an exception). This asymmetry can be utilized in neural networks (whereI ′ k ≈ I k ≫ J l ) to accelerate backpropagation with FFT.

Tensor outer product The derivatives of T with respect to T (0) and T (1) can be directly observed from the definition of tensor outer product in Equation D.4.

∂T i0,••• ,im-1,j0,••• ,jn-1 ∂T (0) i0,••• ,im-1 = T (1) j0,••• ,jn-1 (E.12a) ∂T i0,••• ,im-1,j0,••• ,jn-1 ∂T (1) j0,••• ,jn-1 = T (0) i0,••• ,im-1 (E.12b)

Similar to all previously discussed operations, we first derive the backpropagation equations for T (0) and T (1) at the entries level using standard chain rule.

(0) i0,••• ,im-1 = J0-1 j0=0 • • • Jn-1-1 jn-1=0 ∂L ∂T i0,••• ,im-1,j0,••• ,jn-1 T (1) j0,••• ,jn-1 (E.13a) ∂L ∂T (1) j0,••• ,jn-1 = I0-1 i0=0 • • • Im-1-1 im-1=0 ∂L ∂T i0,••• ,im-1,j0,••• ,jn-1 T (0) i0,••• ,im-1 (E.13b)

which can then be converted to tensor notations in this paper:

∂L ∂T (0) = ∂L ∂T × m 0 • • • • • × m+n-1 n-1 T (1) (E.14a) ∂L ∂T (1) = ∂L ∂T × 0 0 • • • • • × m-1 m-1 T (0) (E.14b)

The number of operations required for both equations are O((

m-1 u=0,u =k I u )( n-1 v=0 J v ))

, which are again identical to one in the forward pass in Equation D.4.

Tensor partial outer product Finally, the derivatives of T with respect to T (0) and T (1) can be obtained from the definition of tensor partial outer product in Equation D.5a.

∂T i0,••• ,i k-1 ,r,i k+1 ,••• ,im-1,j0,••• ,j l-1 ,j l+1 ,••• ,jn-1 ∂T (0) i0,••• ,i k-1 ,r,i k+1 ,••• ,im-1 = T (1) j0,••• ,j l-1 ,r,j l+1 ,••• ,jn-1 (E.15a) ∂T i0,••• ,i k-1 ,r,i k+1 ,••• ,im-1,j0,••• ,j l-1 ,j l+1 ,••• ,jn-1 ∂T (1) j0,••• ,j l-1 ,r,j l+1 ,••• ,jn-1 = T (0) i0,••• ,i k-1 ,r,i k+1 ,••• ,im-1 (E.15b)

Again with chain rule, the backpropagation equations for T (0) and T (1) at the entries level are

∂L ∂T (0) i0,••• ,i k-1 ,r,i k+1 ,••• ,im-1 = J0-1 j0=0 • • • J l-1 -1 j l-1 =0 J l+1 -1 j l+1 =0 • • • Jn-1-1 jn-1=0 ∂L ∂T i0,••• ,i k-1 ,r,i k+1 ,••• ,im-1,j0,••• ,j l-1 ,j l+1 ,••• ,jn-1 T (1) j0,••• ,j l-1 ,r,j l+1 ,••• ,jn-1 (E.16a) ∂L ∂T (1) j0,••• ,j l-1 ,r,j l+1 ,••• ,jn-1 = I0-1 i0=0 • • • I k-1 -1 i k-1 =0 I k+1 -1 i k+1 =0 • • • Im-1-1 im-1=0 ∂L ∂T i0,••• ,i k-1 ,r,i k+1 ,••• ,im-1,j0,••• ,j l-1 ,j l+1 ,••• ,jn-1 T (0) i0,••• ,i k-1 ,r,i k+1 ,••• ,im-1 (E.16b)

Though the backpropagation equations above appear very similar to the ones for tensor contraction in Equations E.1a and E.1b, written in tensor notations, they are almost the same as the ones for tensor convolution in Equations E.10a and E.10b, except that ( * k l ) ⊤ 's are now replaced by

⊗ k l . ∂L ∂T (0) = swapaxes ∂L ∂T × m 0 • • • • • × m+l-1 l-1 • ⊗ k l • × m+l-2 l+1 • • • • • × m+n-2 n-1 T (1) (E.17a) ∂L ∂T (1) = swapaxes ∂L ∂T × 0 0 • • • • • × k-1 k-1 • ⊗ k k • × k+1 k+1 • • • • • × m-1 m-1 T (0) (E.17b)

It is not difficult to recognize the time complexity for the two equations above are O((

m-1 u=0 I u ) ( n-1 v=0,v =l J v ))

, which are identical to the ones in forward pass in Equation D.5a.

Compound operations Up to this point, we have developed the derivatives and backpropagation equations for all basic operations. In this part, we will continue to show similar techniques above equally apply to compound operations, though slightly more involved, and derive the backpropagation equations for the examples we used in Appendix D. Though these equations are not immediately useful in later sections, the techniques to derive them are useful for all other compound operations. Furthermore, these induced equations, which are more complicated than their original definitions, serve as complementary examples of compound operations to the ones in the last section.

• Simultaneous multi-operations between two tensors. In Appendix D, we introduced a compound operation (⊗ 0 0 • * 1 1 • × 2 2 ) on two tensors T (0) ∈ R R×X×S and T (1) ∈ R R×H×S , which returns a tensor T ∈ R R×X ′ . Here, we recap its definitions as follows:

T T (0) ⊗ 0 0 • * 1 1 • × 2 2 T (1) (E.18a) T r,: S-1 s=0 T (0) r,:,s * T (1)

r,:,s

(E.18b)

To ease the derivation, we use the trick to start with circular convolution: directly apply the chain rule, the backpropagation equations at entries level are obtained as follows:

∂L ∂T (0) r,:,s = Cir T (1)

r,:,s

⊤ ∂L T r,: (E.19a) ∂L ∂T (1) r,:,s = Cir T (0)

r,:,s

⊤ ∂L T r,: (E.19b)

Now we convert the equations above to tensor notations, and replace the circular convolutions with their adjoints to obtain general backpropagation rules:

∂L ∂T (0) = ∂L ∂T ⊗ 0 0 • * 1 1 flipaxis(T (1) ) = ∂L ∂T ⊗ 0 0 • ( * 1 1 ) ⊤ T (1) (E.20a) ∂L ∂T (1) = ∂L ∂T ⊗ 0 0 • * 1 1 flipaxis(T (0) ) = ∂L ∂T ⊗ 0 0 • ( * 1 1 ) ⊤ T (0) (E.20b)

For simplicity, we assume FFT is not used to accelerate the backpropagation equations. In this case, the derivatives with respect to T (0) and T (1) can be computed in O(RHX ′ ST ) and O(RXX ′ ST ) operations respectively. Again, the time complexities of forward and backward passes are not the same when a (compound) tensor operation contains convolution.

• Simultaneous operations between a tensor and a set of multiple tensors. Another compound operation presented in Appendix D is defined between a tensor U ∈ R R×X×S and a set of tensors T (0) ∈ R R×P , T (1) ∈ R H×Q and T (2) ∈ R S×T , which returns a tensor V ∈ R R×X ′ ×P ×Q×T . Again, we recap its definitions in the following:

V U ⊗ 0 0 T (0) * 1 0 T (1) × 2 0 T (2) (E.21a) V r,:,p,q,t S-1 s=0 T (0) r,p U r,:,s * T (1) :,q T (2) s,t (E.21b)

In order to derive the backpropagation rule for the core tensor U, we follow the standard procedure to (1) first obtain its entries level representation, and (2) explicitly convert it to tensor notations subsequently. Concretely, the backpropagation equation in both formats are displayed as follows:

∂L ∂U r,:,s = P -1 p=0 Q-1 q=0 T -1 t=0 T (0) r,p Cir T (1) :,q ⊤ ∂L ∂V r,:,p,q,t T (2) s,t (E.22a) ∂L ∂U = ∂L ∂V ⊗ 0 0 • × 1 2 T (0) ( * 0 1 ) ⊤ • × 1 3 T (1) × 0 2 • × 1 4 T (2) (E.22b)

Notice that the equation above is indeed simultaneous multi-operations between a tensor and a set of multiple tensors, which combines two types of "basic" compound operations introduced in Appendix D. In principle, we can obtain backpropagation equations for {T (0) , T (1) , T (2) } in the same manner. However, there is a simpler way to derive them by rewriting the definition as:

V = swapaxes U * 1 0 T (1) × 2 0 T (2) ⊗ 0 0 T (0) = swapaxes U (0) ⊗ 0 0 T (0) (E.23a) = swapaxes U ⊗ 0 0 T (0) × 2 0 T (2) * 1 0 T (1) = swapaxes U (1) * 1 1 T (1) (E.23b) V = swapaxes U ⊗ 0 0 T (0) * 1 0 T (1) × 2 0 T (2) = swapaxes U (2) × 2 2 T (2) (E.23c)

where U (0) , U (1) and U (2) are short-hand notations for U( * 1 0 T (1) × 2 0 T (2) ), U(⊗ 0 0 T (0) × 2 0 T (2) ) and U(⊗ 0 0 T (0) * 1 0 T (1) ). With these notations, we are able to reduce these complex expressions to basic ones, by which we can reuse the backpropagation rules derived in this section:

∂L ∂T (0) = ∂L ∂V ⊗ 0 0 • × 1 1 • × 3 2 • × 4 3 U (0) = ∂L ∂V ⊗ 0 0 • × 1 1 • × 3 2 • × 4 3 U * 1 0 T (1) × 2 0 T (2) (E.24a) ∂L ∂T (1) = ∂L ∂V × 0 0 • ( * 1 1 ) ⊤ • × 2 2 • × 4 3 U (1) = ∂L ∂V × 0 0 • ( * 1 1 ) ⊤ • × 2 2 • × 4 3 U ⊗ 0 0 T (0) × 2 0 T (2) (E.24b) ∂L ∂T (2) ⊤ = ∂L ∂V × 0 0 • × 1 1 • × 2 3 • × 3 4 U (2) = ∂L ∂V × 0 0 • × 1 1 • × 2 3 • × 3 4 U ⊗ 0 0 T (0) * 1 0 T (1) (E.24c)

The complicity of tensor operation culminates at this point: the equations above are examples of simultaneous multi-operations on multiple tensors, which we omitted in the discussion in Appendix D due to their complexity. Although the expressions themselves suggest particular orderings to evaluate the compound operations, they are merely the traces of the techniques used in deriving them. It is completely reasonable to reorganize the equations such that they can be computed with more efficient strategies: for instance, one can verify that the following set of equations is actually equivalent to the one above:

∂L ∂T (0) = ∂L ∂V ( * 1 0 ) ⊤ • × 3 1 T (1) × 4 1 T (2) ⊗ 0 0 • × 1 1 • × 3 2 U (E.25a) ∂L ∂T (1) = ∂L ∂V ⊗ 0 0 • × 2 1 T (0) × 4 1 T (2) × 0 0 • ( * 1 1 ) ⊤ • × 3 2 U (E.25b) ∂L ∂T (2) ⊤ = ∂L ∂V ⊗ 0 0 • × 2 1 T (0) ( * 1 0 ) ⊤ • × 3 1 T (1) × 0 0 • × 1 1 U (E.25c)

As discussed in Appendix D, the problem to find the optimal order to evaluate a compound operation over multiple tensors is NP-hard in general and usually we need to resort to heuristics to obtain a reasonably efficient algorithm. Indeed, one can verify that the second set of equations is more efficient than the first one. For this example, interested readers are encouraged to find the most efficient way by combinatoric search.

F TENSOR DECOMPOSITIONS:

Tensor decompositions are natural extensions of matrix factorizations for multi-dimensional arrays. In this section, we will review three commonly used tensor decompositions, namely CANDE-COMP/PARAFAC (CP) decomposition, Tucker (TK) decomposition and Tensor-train (TT) decomposition. CP and Tucker decompositions are classic methods to factorize a tensor and are thoroughly reviewed in an early survey Kolda & Bader (2009) , while Tensor-train decomposition is recently proposed in Oseledets (2011), in which the authors show that it has superior numerical stability and scalability than the classic ones. More advanced schemes of tensor decompositions (known as tensor networks) are reviewed in Cichocki et al. (2016) .

For each of these decompositions, we will present their forms both at the entries level and in tensor notations introduced in Appendix D. When tensor decompositions are used in neural networks, a natural question to ask is how the backpropagation algorithm adapts to the decomposition schemes, i.e. how the gradient of the original tensor backpropagates to its factors. In this section, we will follow the standard procedure in Appendix E to derive the corresponding backpropagation equation for each tensor decomposition. Different from previous works (Novikov et al., 2015; Kossaifi et al., 2017b ) that use matrix calculus following matricization, we present the backpropagation equations directly in tensor notations, which makes our presentation concise and easy to analyze. As we will see in the analyses, backpropagation equations through the original tensor to its factors are computationally expensive for all decomposition schemes, therefore it is preferable to avoid explicit computation of these equations in practice.

CP decomposition CP decomposition is a direct generalization of singular value decomposition (SVD) which decomposes a tensor into additions of rank-1 tensors (outer product of multiple vectors). Specifically, given an m-order tensor T ∈ R I0×I1×•••×Im-1 , CP decomposition factorizes it into m factor matrices {M (0) } m-1 l=0 , where

M (l) ∈ R R×I l , ∀l ∈ [m],

where R is called the canonical rank of the CP decomposition, which is allowed to be larger than the I l 's.

T i0,••• ,im-1 R-1 r=0 M (0) r,i0 • • • M (m-1) r,im-1 (F.1a) T R-1 r=0 M (0) r,: ⊗ • • • ⊗ M (m-1)

r,:

= 1 × 0 0 M (0) ⊗ 0 0 • • • ⊗ 0 0 M (m-1) (F.1b)

where 1 ∈ R R is an all-ones vector of length R. With CP decomposition, T can be represented with only ( m-1 l=0 I l )R entries instead of ( m-1 l=0 I l ) as in the original tensor. Now we proceed to derive the backpropagation rules for CP decomposition, i.e. the equations relating ∂L/∂T to {∂L/∂M (l) } m-1 l=0 . In order to avoid deriving these equations from the entries level, we first isolate the factor of interest and rewrite the definition of CP decomposition as:

T = swapaxes M (0) • • • ⊗ 0 0 M (l-1) ⊗ 0 0 M (l+1) • • • ⊗ 0 0 M (m-1) × 0 0 M (l) (F.2) = swapaxes A (l) × 0 0 M (l)

where we treat the first term as a constant tensor A (l) . Once we reduce the compound operation to a basic one, we can simply refer to the rule derived in Appendix E, which gives us

∂L ∂M (l) = ∂L ∂T × 0 1 • • • • × l-1 l • × l+1 l+1 • • • × m-1 m-1 A (l) = ∂L ∂T × 0 1 • • • • × l-1 l • × l+1 l+1 • • • × m-1 m-1 M (0) ⊗ 0 0 • • • ⊗ 0 0 M (l-1) ⊗ 0 0 M (l+1) ⊗ 0 0 • • • ⊗ 0 0 M (m-1) (F.3)

The number of operations to compute one such equation both is O((

m-1 l=0 I l )R), therefore a total number of O(m( m-1 l=0 I l )R

) is required for all m equations. Therefore, evaluating these equations are computationally expensive (which takes O(mR) order as many operations as the size ( m-1 l=0 I l ) of the original tensor T ), and should be avoided whenever possible.

Tucker decomposition Tucker decomposition provides more general factorization than CP decomposition. Given an m-order tensor T ∈ R I0×I1×•••×Im-1 , Tucker decomposition factors it into m factor matrices {M (l) } m-1 l=0 , where

M (l) ∈ R R l ×I l , ∀l ∈ [m]

and an additional m-order core tensor C ∈ R R0×R1×•••×Rm-1 , where the Tucker ranks R l 's are required to be smaller or equal than the dimensions at their corresponding modes, i.e. 

R l ≤ I l , ∀l ∈ [m]. T i0,••• ,im-1 R0-1 r0=0 • • • Rm-1-1 rm-1=0 C r0,...,rm-1 M (0) r0,i0 • • • M (m-1) rm-1,im-1 (F.4a) T C × 0 M (0) × 1 M (1) • • • × m-1 M (m-1) (F.4b) Notice that when R 0 = • • • = R m-1 = R

( m-1 l=0 R l + m-1 l=0 I l R l ) entries.

The backpropagation equations relating ∂L/∂T to ∂L/∂C and {∂L/M (l) } m-1 l=0 can be derived similarly as in CP decomposition. First, we derive the equation for C at the entries level:

∂L ∂C r0,••• ,rm-1 = I0-1 i0=0 • • • Im-1-1 im-1=0 ∂L ∂T i0,••• ,im-1 M (0) r0,i0 • • • M (m-1) rm-1,im-1 (F.5)

The equation above, written in tensor notations, reveals an expression in "reversed" Tucker form:

∂L ∂C = ∂L ∂T × 0 (M (0) ) ⊤ • • • × m-1 (M (m-1) ) ⊤ (F.6)

Although the number of operations to evaluate the equation depends on the particular order of tensor multiplications between ∂L/∂T and {M (l) ) ⊤ } m-1 l=0 , the time complexity can be bounded by

O(( m-1 l=0 R l )( m-1 l=0 I l ))

, where the worst case is achieved when

R l = I l , ∀l ∈ [m].

Regarding the backpropagation equations for the factors {M (l) } m-1 l=0 , we again isolate the factor of interest and rewrite the definition of Tucker decomposition as:

T = C × 0 M (0) • • • × l-1 M (l-1) × l+1 M (l+1) • • • × m-1 M (m-1) × l M (l) (F.7) = A (l) × l M (l)

where the first term is abbreviated as a tensor A (l) . Subsequently, we apply the standard backpropagation rule of tensor multiplication in Appendix E and obtain the following equation:

∂L ∂M (l) ⊤ = ∂L ∂T × 0 0 • • • • • × l-1 l-1 • × l+1 l+1 • • • • • × m-1 m-1 A (l) = ∂L ∂T × 0 (M (0) ) ⊤ • • • × l-1 (M (l-1) ) ⊤ × l+1 (M (l+1) ) ⊤ • • • × m-1 (M (m-1) ) ⊤ × 0 0 • • • • • × l-1 l-1 • × l+1 l+1 • • • • • × m-1 m-1 C (F.8)

where the second expression is equivalent to the first one, but requires fewer operations. Though the exact number of operations depends on the order of tensor multiplications, it can be (again) bounded by O((

m-1 l=0 R l )( m-1 l=0 I l )).

Therefore, the total time complexity for all (m + 1) equations (m for {M (l) } m-1 l=0 and one for C) is bounded by O((m + 1)(

m-1 l=0 R l )( m-1 l=0 I l ))

, which is also highly inefficient and should be avoided in practice.

Tensor-train decomposition Tensor-train decomposition factorizes a m-order tensor into m interconnected low-order core tensors {T (l) } m-1 l=0 , where

T (l) ∈ R R l ×I l ×R l+1 , l = 1, • • • , m -2 with T (0) ∈ R I0×R0 , and T (m-1) ∈ R Rm-1×Im-1 such that T i0,...,im-1 R0-1 r0=1 • • • Rm-2-1 rm-2=1 T (0) i0,r0 T (1) r0,i1,r1 • • • T (m-1) rm-2,im-1 (F.9a) T T (0) × -1 0 T (1) × -1 0 • • • × -1 0 T (m-1) (F.9b)

where the R l 's are known as Tensor-train ranks, which controls the tradeoff between the number of parameters and accuracy of representation. With Tensor-train decomposition, a tensor is represented by

(R 0 I 0 + m-2 l=1 R l I l R l+1 + R m-1 I m-1

) entries. The backpropagation equations are derived following the paper Novikov et al. (2015) , although we reformat them in tensor notations. To begin with, we introduce two sets of auxiliary tensors {P (l) } m-1 l=0 and {Q (l) } m-1 l=0 as follows:

P (l) = T (0) × -1 0 • • • × -1 0 T (l) = P (l-1) × -1 0 T (l) , ∀l = 1, • • • , m -1 (F.10a) Q (l) = T (l) × -1 0 • • • × -1 0 T (m-1) = T (l) × -1 0 Q (l+1) , ∀l = 0, • • • , m -2 (F.10b)

with corner cases as P (0) = T (0) and Q (m-1) = T (m-1) . Note that both {P (l) } m-1 l=0 and {Q (l) } m-1 l=0 can be computed using dynamic programming (DP) using the recursive definitions above. With these auxiliary tensors, the definition of Tensor-train decomposition in Equation F.9b can be rewritten as:

T = P (l-1) × -1 0 T (l) × -1 0 Q (l+1) , ∀l = 0, m -1 (F.11) = T (0) × -1 0 Q (1) = P (m-2) × -1 0 T (m-1)

Applying the backpropagation rule for tensor contraction twice, the backpropagation equations can be obtained in tensor notations as:

∂L ∂T (l) = P (l-1) × 0 0 • • • • • × l-1 l-1 ∂L ∂T × 2 1 • • • • • × m-l m-l-1 Q (l+1) , ∀l = 0, m -1 (F.12) ∂L ∂T (0) = ∂L ∂T × 1 1 • • • • • × m-1 m-1 Q (1) , ∂L ∂T (m-1) = P (m-2) × 0 0 • • • • • × m-2 m-2

Neglecting the expense of precomputing the auxiliary tensors {P (l) } m-1 l=0 and {Q (l) } m-1 l=0 , the time complexity for the

l th equation is O(R l R l+1 ( m-1 l=0 I l )), therefore the total number of operations for all m equations is O((R 0 + m-2 l=0 R l R l+1 + R m-1 )( m-1 l=0 I l ))

, which is again prohibitively large for practical use.

Variants of standard decompositions:

In this paper, tensor decompositions are usually used in flexible ways, i.e. we will not stick to the standard formats defined in the previous paragraphs. Indeed, we consider tensor decomposition as a reverse mapping of tensor operations: given a tensor T and a set of operations, the corresponding tensor decomposition aims to recover the input factors {T (l) } m-1 l=0 such that the operations on these factors return a tensor approximately equal to the given one. In the following, we demonstrate some possibilities using examples:

• The ordering of the modes can be arbitrary. Therefore, CP decomposition of 3-order tensor T ∈ R I0×I1×I2 can be factorized as

T i0,i1,i2 = R-1 r=0 M (0) i0,r M (1) i1,r M (2) i2,r or T i0,i1,i2 = R-1 r=0 M (0) r,i0 M (1) r,i1 M (2) r,i2 , or even T i0,i1,i2 = R-1 r=0 M (0) i0,r M (1)

r,i1 M

(2) i2,r . It is easy to observe these decompositions are equivalent to each other if factor matrices are properly transposed.

• A tensor may be partially factorized over a subset of modes. For example, we can define a partial Tucker decomposition which factors only the last two modes of a 4-order tensor T ∈ R I0×I1×I2×I3 into a core tensor C ∈ R I0×I1×R2×R3 and two factor matrices 3) if written in our tensor notations.

M (2) ∈ R R2×I2 , M (3) ∈ R R3×I3 such that T i0,i1,i2,i3 = R2-1 r2=0 R3-1 r3=0 C i0,i1,r2,r3 M (2) r2,i2 M (3) r3,i3 or alternatively T = C × 2 M (2) × 3 M (

• Multiple modes can be grouped into supermode and decomposed like a single mode. For example, given a 6-order tensor T ∈ R I0×I1×I2×J0×J1×J2 can be factorized into three factors

T (0) ∈ R I0×J0×R0 , T (1) ∈ R R0×I1×J1×R1 and T (2) ∈ R R1×I2×J2 such that T i0,j0,i1,j1,i2,j2 = R2-1 r2=0 R3-1 r3=0 T (0) i0,j0,r0 T (1) r0,i1,j1,r1 T (2)

r1,i2,j2 , or more succinctly as 2) . where I 0 and J 0 are grouped into a supermode (I 0 , J 0 ) and similarly for (I 1 , J 1 ) and (I 2 , J 2 ).

T = T (0) × -1 0 T (1) × -1 0 T (

Notations O(# of params.) O(# of backprop ops.) original

T I 0 CP 1 × 0 0 (M (0) ⊗ 0 0 • • • ⊗ 0 0 M (m-1) ) mI 1 m R mIR TK C(× 0 M (0) • • • × m-1 M (m-1) ) R m + mI 1 m R m 2 IR TT T (0) × -1 0 • • • × -1 0 T (m-1) mI 1 m R 2 mIR 2

Table 9: Summary of tensor decompositions. In this table, we summarize three types of tensor decompositions in tensor notations, and list their numbers of parameters and time complexities to backpropagate the gradient of a tensor T ∈ R I0×I1×•••Im-1 to its m factors (and an additional core tensor C for Tucker decomposition). For simplicity, we assume all dimensions I l 's of T are equal, and denote the size of T as the product of all dimensions I = m-1 l=0 I l . Furthermore, we assume all ranks R l 's (in Tucker and Tensor-train decompositions) share the same number R.

G PLAIN TENSOR DECOMPOSITION S ON CONVOLUTIONAL LAYER:

In this section, we will show how tensor decomposition is able to compress (and accelerate) the standard convolutional layer in neural networks. In order to achieve this, we first represent the operation of a standard convolutional layer in tensor notations. By factorizing the tensor of parameters (a.k.a. kernel) into multiple smaller factors, compression is achieved immediately.

As we discussed in Appendix F, learning the factors through the gradient of the original tensor of parameters is highly inefficient. In this section, we provide an alternative strategy that interacts the input with the factors individually, in which explicit reference to the original kernel is avoided. Therefore, our strategy also reduces the computational complexity along with compression.

For simplicity, we assume FFT is not used in computing convolutions, although we show in Appendix E that FFT can possibly speed up the backward pass Mathieu et al. (2013) .

Standard convolutional layer:

In modern convolutional neural network (CNN), a standard convolutional layer is parameterized by a 4-order kernel K ∈ R H×W ×S×T , where H and W are height and width of the filters (which are typically equal), S and T are the number of input and output channels respectively. A convolutional layer maps a 3-order tensor U ∈ R X×Y ×S to another 3order tensor V ∈ R X ′ ×Y ′ ×T , where X and Y are the height and width for the input feature map, while X ′ and Y ′ are the ones for the output feature map, with the following equation:

V x,y,t = S-1 s=0 i,j K i,j,s,t U i+dx,j+dy,s (G.1)

where d is the stride of the convolution layer and the scopes of summations over i and j are determined by the boundary conditions. Notice that the number of parameters in a standard convolutional layer is HW ST and the number of operations needed to evaluate the output V is O(HW ST XY ).

With tensor notations in this paper, a standard convolutional layer can be defined abstractly as

V = U * 0 0 • * 1 1 • × 2 2 K (G.2)

which states that the standard convolutional layer in fact performs a compound operation of two tensor convolutions and one tensor contraction simultaneously between the input tensor U and the kernel of parameters K. Following the standard procedure in Appendix E, we obtain both backpropagation equations in tensor notations as follows:

∂L ∂U = ∂L ∂V ( * 0 0 ) ⊤ • ( * 1 1 ) ⊤ • × 2 3 K, ∂L ∂K = U ( * 0 0 ) ⊤ • ( * 1 1 ) ⊤ ∂L ∂V (G.3)

It is not difficult to verify that the numbers of operations to compute these two backpropagation equations are O(HW ST X ′ Y ′ ) and O(XY ST X ′ Y ′ ) respectively.

In the next few paragraphs, we will apply various decompositions in Appendix F as well as singular value decomposition (SVD) on the kernel K, and derive the steps to evaluate Equation G.1 that interact the input with the factors individually. Interestingly, these steps are themselves (non-standard) convolutional layers, therefore tensor decomposition on the parameters is equivalent to decoupling a layer in the original model into several sublayers in the compressed network, which can be implemented efficiently using modern deep learning libraries. For simplicity, we assume in the analyses of these decomposition schemes that the output feature maps have approximately the same size as the input ones, i.e. X ≈ X ′ , Y ≈ Y ′ .

SVD-convolutional layer:

Many researchers propose to compress a convolutional layer using singular value decomposition, under the name of dictionary learning (Jaderberg et al., 2014; Denton et al., 2014; Zhang et al., 2015) . These methods differ in their matricization of the tensor of parameters K, i.e. how to group the four modes into two and flatten the kernel K into a matrix. By simple combinatorics, it is not difficult to show there are seven different types of matricization in total. Here, we only pick to present the one by Jaderberg et al. (2014) , which groups filter height and input channels as a supermode (H, S) and filter width and output channels (W, T ) as another.

K i,j,s,t = R-1 r=0 K (0) i,s,r K (1) j,r,t (G.4a) K = swapaxes K (0) × 2 1 K (1) (G.4b)

where K (0) ∈ R H×S×R and K (1) ∈ R W ×R×T are the two factor tensors. It is easy to see an SVDconvolutional layer has (HS + W T )R parameters in total (HSR in K (0) and W T R in K (1) ). Now we plug the Equation G.4a into G.1, and break the evaluation of V into two steps such that only one factor is involved at each step.

U (0) x,j+dy,r = S-1 s=0 i K (0) i,s,r U i+dx,j+dy,s (G.5a) V x,y,t = R-1 r=0 j K (1) j,r,t U (0) x,j+dy,r (G.5b)

where U (0) ∈ R X ′ ×Y ×R is the intermediate result after the first step. This two-steps procedure requires only requires O((HSXY +W T X ′ Y )R) operations in the forward pass (O(HSRXY ) at the first step and O(W T RX ′ Y ) at the second). This number is smaller than O(HW ST XY ) as in the standard convolutional layer, since

(HS + W T )R ≤ HW ST implies (HSXY + W T X ′ Y )R ≤ (HS + W T )RXY ≤ HW ST XY .

Therefore, SVD-convolutional layer also outperforms the standard one in efficiency. (Notice that if K is reconstructed explicitly, the forward pass requires at least O(HW ST XY ) operations.) In tensor notations, these steps can be written concisely as follows:

U (0) = U * 0 0 • × 2 1 K (0) (G.6a) V = U (0) * 1 1 • × 2 1 K (1) (G.6b)

After decomposition, each operation is still a compound operation of tensor convolution and tensor contraction, and therefore itself a convolutional layer whose filters have size either H × 1 and 1 × W . Effectively, SVD-convolutional layer is in fact a concatenation of two convolutional layers without nonlinearity in between. Now we proceed to derive the corresponding backpropagation equations for these two steps following the procedure in Appendix E, which are presented in the following:

∂L ∂U (0) = ∂L ∂V ( * 1 0 ) ⊤ • × 2 2 K (1) , ∂L ∂K (1) = U (0) × 0 0 • ( * 1 1 ) ⊤ ∂L ∂V (G.7a) ∂L ∂U = ∂L ∂U (0) ( * 0 0 ) ⊤ • × 2 2 K (0) , ∂L ∂K (0) = U ( * 0 0 ) ⊤ • × 1 1 ∂L ∂U (0) (G.7b)

It is not hard to show the number of operations required to obtain the derivatives with respect to U and

U (0) is O((HSX ′ Y + W T X ′ Y ′ )R

), and the one for the factors is

O((XSX ′ Y + Y T X ′ Y ′ )R).

CP-convolutional layer Both Lebedev et al. (2014); Denton et al. (2014) propose to decompose the kernel K using CP decomposition, differing at whether the height H and width W of the filters are grouped into a supermode. For simplicity, we follow the scheme in Denton et al. ( 2014):

K i,j,s,t = R-1 r=0 K (0) s,r K (1) i,j,r K (2) r,t (G.8a) K = 1 × 0 2 K (1) ⊗ 2 1 K (0) ⊗ 2 0 K (2) (G.8b)

where K (0) ∈ R S×R , K (1) ∈ R H×W ×R and K (2) ∈ R R×T are three factor tensors, which contain (HW + S + T )R parameters in total. Again, plugging Equation G.8a into G.1 yields a three-steps

Under review as a conference paper at ICLR 2019 procedure to evaluate V:

U (0) i+dx,j+dy,r = S-1 s=0 K (0) s,r U i+dx,j+dy,s (G.9a) U (1) x,y,r = i,j K (1) i,j,r U (0) i+dx,j+dy,r (G.9b) V x,y,t = R-1 r=0 K (2) r,t U (1)

x,y,r

(G.9c)

where U (0) ∈ R X×Y ×R and U (1) ∈ R X ′ ×Y ′ ×R are two intermediate tensors. Written in tensor notations, these equations are represented as:

U (0) = U × 2 K (0) (G.10a) U (1) = U (0) * 0 0 • * 1 1 • ⊗ 2 2 K (1) (G.10b) V = U (1) × 2 K (2)

(G.10c) After CP decomposition, the first and third steps are basic tensor multiplications on the input/intermediate tensor, which are usually named (weirdly) as 1 × 1 convolutional layers despite that no convolution is involved at all, while the second step is a compound operation of two tensor convolutions and one partial outer product, which is known as depth-wise convolutional layer (Chollet, 2016) . The number of operations for these three steps are O(SRXY ), O(HW RXY ) and O(T RX ′ Y ′ ) respectively, resulting in a time complexity of O((SXY + HW XY + T X ′ Y ′ )R) for the forward pass, which is faster than the standard convolutional layer, since (HW + S + T )R ≤ HW ST implies (SXY + HW XY + T X ′ Y ′ )R ≤ HW ST XY . Now we proceed to obtain their backpropagation equations following the procedure in Appendix E:

∂L ∂U (1) = ∂L ∂V × 2 (K (2) ) ⊤ , ∂L ∂K (2) = U (1) × 0 0 • × 1 1 ∂L ∂V (G.11a) ∂L ∂U (0) = ∂L ∂U (1) ( * 0 0 ) ⊤ • ( * 1 1 ) ⊤ • ⊗ 2 2 K (1) , ∂L ∂K (1) = ∂L ∂U (1) ( * 0 0 ) ⊤ • ( * 1 1 ) ⊤ • ⊗ 2 2 U (1) ∂L ∂U = ∂L ∂U (0) × 2 (K (0) ) ⊤ , ∂L ∂K (0) = U × 0 0 • × 1 1 ∂L ∂U (0) (G.11b)

The number of operations in all three steps to calculate the derivatives with respect to input/intermediate tensors can be counted as O((SXY + HW X ′ Y ′ + T X ′ Y ′ )R), while the one for the factors as O((SXY

+ XY X ′ Y ′ + T X ′ Y ′ )R).

Tucker-convolutional layer The use of Tucker decomposition to compress and accelerate convolutional layers is proposed in Kim et al. (2015) . Despite the name of Tucker decomposition, they in fact suggest a partial Tucker decomposition, which only factorizes the modes over the numbers of input/output filters and keeps the other two modes for filter height/width untouched.

EQUATION

i,j,rs,rt K

(2) rt,t (G.12a)

K = K (1) × 2 (K (0) ) ⊤ × 3 K (2) (G.12b)

where K (0) ∈ R S×Rs , K (1) ∈ R H×W ×Rs×Rt and K (2) ∈ R Rt×T are three factor tensors, with a total of (SR s + HW R s R t + R t T ) parameters. All that follow are identical to the ones for SVD and CP layers. A three-steps forward pass procedure is obtained by plugging Equation G.12a in G.1.

EQUATION

i,j,rs,rt U (0) i+dx,j+dy,rs

(G.13b) V x,y,t = Rt-1 rt=0 K (2) rt,t U (1)

x,y,rt

(G.13c)

where U (0) ∈ R X×Y ×Rs and U (1) ∈ R X ′ ×Y ′ ×Rt are two intermediate tensors. These three steps, with number of operations of

O(R s SXY ), O(HW R s R t XY ) and O(R t T X ′ Y ′ ) respectively, leads to a total time complexity of O(SR s XY + HW R s R t XY + R t T X ′ Y ′ )

for the forward pass. Like CP and SVD convolutional layers, Tucker-convolutional layer is faster than the standard convolutional layer, since

SR s + HW R s R t + R t T ≤ HW ST implies SR s XY + HW R s R t XY + R t T X ′ Y ′ ≤ HW ST XY .

These equations, again, can be concisely written in tensor notations:

U (0) = U × 2 K (0) (G.14a) U (1) = U (0) * 0 0 • * 1 1 • × 2 2 K (1) (G.14b) V = U (1) × 2 K (2) (G.14c)

where the first and the third steps are two 1 × 1 convolutional layers, and the second step is itself a standard convolutional layer, which only differs from CP-convolutional layer at the second step. For completeness, we summarize all backpropagation equations in the following:

∂L ∂U (1) = ∂L ∂V × 2 (K (2) ) ⊤ , ∂L ∂K (2) = U (1) × 0 0 • × 1 1 ∂L ∂V (G.15a) ∂L ∂U (0) = ∂L ∂U (1) ( * 0 0 ) ⊤ • ( * 1 1 ) ⊤ • × 2 3 K (1) , ∂L ∂K (1) = ∂L ∂U (1) ( * 0 0 ) ⊤ • ( * 1 1 ) ⊤ U (1) (G.15b) ∂L ∂U = ∂L ∂U (0) × 2 (K (0) ) ⊤ , ∂L ∂K (0) = U × 0 0 • × 1 1 ∂L ∂U (0) (G.15c)

Referring to the CP-convolutional layer, the time complexity for the backward pass is obtained with slight modification: the number of operations for the input/intermediate tensors is

O(SR s XY + HW R s R t X ′ Y ′ + R t T X ′ Y ′ ), and the one for factors is O(SR s XY + XY HW X ′ Y ′ + R t T X ′ Y ′ ).

Tensor-train-convolutional layer Lastly, we propose to apply Tensor-train decomposition to compress a convolutional layer. However, naive Tensor-train decomposition on the kernel K may give inferior results (Garipov et al., 2016) , and careful reordering of the modes is necessary. In this paper, we propose to reorder the modes as (input channels S, filter height H, filter width W , output channels T ), and decompose the kernel as

EQUATION

rs,i,r K

(2) r,j,rt K

(3) rt,t (G.16a)

K = swapaxes K (0) × -1 0 K (1) × -1 0 K (2) × -1 0 K (3) (G.16b)

where 

K (0) ∈ R S×Rs , K (1) ∈ R Rs×H×R , K (2) ∈ R R×W ×Rt and K (3) ∈ R

RXY + W RR T X ′ Y + R t T X ′ Y ′ ) in the forward pass. U (0) i+dx,j+dy,rs = S-1 s=0 K (0) s,rs U i+dx,j+dy,s (G.17a) U (1)

x,j+dy,r =

Rs-1

rs=0 i K (1)

rs,i,r U (0) i+dx,j+dy,rs (G.17b)

U (2) x,y,rt = R-1 r=0 j K (2)

r,j,rt U

(1)

x,j+dy,r (G.17c)

V x,y,t = Rt-1 rt=0 K (3) rt,t U (2)

x,y,rt (G.17d) where U (0) ∈ R X×Y ×Rs , U (1) ∈ R X ′ ×Y ×R and U (2) ∈ R X ′ ×Y ′ ×Rt are three intermediate tensors.

In tensor notations, these equations can be rewritten as as follows:

U (0) = U × 2 K (0) (G.18a) U (1) = U (0) ( * 0 1 • × 2 0 ) K (1) (G.18b) U (2) = U (1) ( * 1 1 • × 2 0 ) K (2) (G.18c) V = U (2) × 2 K (3) (G.18d)

Tensor-train-convolutional layer is concatenation of four sub-layers, where the first and the last ones are 1 × 1 convolutional layers, while the other two in between are convolutional layers with rectangular kernels. In fact, Tensor-train-convolutional layer can either be interpreted as (1) a Tucker-convolutional layer where the second sublayer is further compressed by a SVD, or (2) a SVD-convolutional layer where both factors are further decomposed again by SVD. Referring to the previous results, the corresponding backpropagation equations are easily derived as

∂L ∂U (2) = ∂L ∂V × 2 (K (3) ) ⊤ , ∂L ∂K (3) = U (2) × 0 0 • × 1 1 ∂L ∂V (G.19a) ∂L ∂U (1) = ∂L ∂U (2) ( * 1 0 ) ⊤ • × 2 2 K (2) , ∂L ∂K (2) = U (1) × 0 0 • ( * 1 1 ) ⊤ ∂L ∂U (2) (G.19b) ∂L ∂U (0) = ∂L ∂U (1) ( * 0 0 ) ⊤ • × 2 2 K (1) , ∂L ∂K (2) = U (0) ( * 0 0 ) ⊤ • × 1 1 ∂L ∂U (1) (G.19c) ∂L ∂U = ∂L ∂U (0) × 2 (K (0) ) ⊤ , ∂L ∂K (0) = U × 0 0 • × 1 1 ∂L ∂U (0) (G.19d)

Similar to all previous layers, the time complexities for input/intermediate tensors and factors can be calculated as 

O(SR s XY +HR s RX ′ Y +W R t T X ′ Y ′ +R t T X ′ Y ′ ) and O(SR s XY +XR s RX ′ Y + Y R t RX ′ Y ′ + R t T X ′ Y ′ ) respectively.

′ Y ′ XY ST X ′ Y ′ SVD (HS + W T )R HSXY + W T X ′ Y )R (HSX ′ Y + W T X ′ Y ′ )R (XSX ′ Y + Y T X ′ Y ′ )R CP (HW + S + T )R (SXY + HW XY + T X ′ Y ′ )R (SXY + HW X ′ Y ′ + T X ′ Y ′ )R (SXY + XY X ′ Y ′ + T X ′ Y ′ )R TK (HW R s R t + SR s + R t T ) (HW R s R t XY + SR s XY + R t T X ′ Y ′ ) (HW R s R t X ′ Y ′ + SR s XY + R t T X ′ Y ′ ) (XY R s R t X ′ Y ′ + SR s XY + R t T X ′ Y ′ ) TT (SR s + HR s R+ W R t R + R t T ) (SR s XY + HR s RXY + W R t RX ′ Y + R t T X ′ Y ′ ) (SR s XY + HR s RX ′ Y + W R t T X ′ Y ′ + R t T X ′ Y ′ ) (SR s XY + XR s RXY + Y R t RX ′ Y ′ + R t T X ′ Y ′ )

Table 10: Summary of plain tensor decomposition on convolutional layer. We list the number of parameters and the number of operations required by forward/backward passes for various plain tensor decomposition on convolutional layer. For reference, a standard convolutional layer maps a set of S feature maps with height X and width Y , to another set of T feature maps with height X ′ and width Y ′ . All filters in the convolutional layer share the same height H and width W .

H RESHAPED TENSOR DECOMPOSITION S ON DENSE LAYER:

The operation of dense layer (a.k.a. fully connected layer) in neural network can be simply characterized by a matrix-vector multiplication, which maps a vector u ∈ R S to another vector v ∈ R T , where S and T are the number of units for the input and output respectively.

v t = S-1 s=0 K s,t u s (H.1a) v = K u (H.1b)

It is easy to see that a dense layer is parameterized by a matrix K with ST parameters, and evaluating the output v requires O(ST ) operations. With a matrix at hand, the simplest compression is via singular value decomposition (SVD), which decomposes K into multiplication of two matrices K = P Q, where P ∈ R S×R , Q ∈ R R×T with R ≤ min(S, T ). With SVD decomposition, the number of parameters is reduced from ST to ((S +T )R) and time complexity from O(ST ) to O((S +T )R).

Inspired by the intuition that invariant structures can be exploited by tensor decompositions in Section 3, we tensorize the matrix K into a tensor

K ∈ R S0×•••×Sm-1×T0×•••×Tm-1 such that S = m-1 l=0 S i , T = m-1 l=0 T l and vec(K) = vec(K). Correspondingly, we reshape the in- put/output u, v into U ∈ R S0×•••×Sm-1 , V ∈ R T0×•••×Tm-1 such that vec(U) = u, vec(V) = v

and present an (uncompressed) tensorized dense layer as follows:

V t0,••• ,tm-1 = S0-1 s0=0 • • • Sm-1-1 sm-1=0 K s0,••• ,sm-1,t0,••• ,tm-1 U s0,••• ,sm-1 (H.2a) V = U × 0 0 • • • • • × m-1 m-1 K (H.2b)

Therefore, a tensorized dense layer, parameterized by a 2m-order tensor K, maps an m-order tensor U to another m-order tensor V. It is straightforward to observe that the tensorized dense layer is mathematically equivalent is to the dense layer in Equation H.1a. Correspondingly, its backpropagation equations can be obtained by simply reshaping the ones for standard dense layer:

∂L ∂u = K ⊤ ∂L ∂v , ∂L ∂K = u ∂L ∂v ⊤ = u ⊗ ∂L ∂v (H.3) ∂L ∂U = K × m 0 • • • × 2m-1 m-1 ∂L ∂V , ∂L ∂K = U ⊗ ∂L ∂V (H.4)

In the section, we will compress the tensorized dense layer by decomposing the kernel K into multiple smaller factors. As we will see, the schemes of tensor decompositions used in this section are not as straightforward as in Appendix F and G, and our principles in their designs are analyzed at the end of this section. Again, learning the the factors through the gradient of the original tensor is extremely costly, therefore a multi-steps procedure to compute the output by interacting the input with the factors individually is desirable. For simplicity of analyses, we will assume for the rest of the paper that S and T are factored evenly, that is

S l ≈ S 1 m , T l ≈ T 1 m , ∀l ∈ [m]

and all ranks are equal to a single number R. r-CP-dense layer Obviously, the simplest way to factorize K is to perform naive CP decomposition over all 2m modes without grouping any supermode. However, such naive decomposition leads to significant loss of information, as we discuss at the end of this section. In this paper, we instead propose to factor the kernel K by grouping (S l , T l )'s as supermodes. Concretely, the tensor of parameters K is decomposed as:

K s0,••• ,sm-1,t0,••• ,tm-1 = R-1 r=0 K (0) r,s0,t0 • • • K (m-1) r,sm-1,tm-1 (H.5a) K = swapaxes 1 × 0 0 (K (0) ⊗ 0 0 • • • ⊗ 0 0 K (m-1) (H.5b)

where K (l) ∈ R R×S l ×T l , ∀l ∈ [m] are m factors and R is the canonical rank that controls the tradeoff between the number of parameters and the fidelity of representation. Therefore, the total number of parameters of a r-CP-dense layer is approximately

m-1 l=0 S l T l R ≈ m(ST ) 1

m R, which is significantly smaller than ST given that R is reasonably small. The next step to derive the sequential procedure mirrors the ones in all schemes in Appendix G, by plugging the Equation H.5a into H.2a, we arrive at a multi-steps procedure for the forward pass.

U (0) r,s0,••• ,sm-1 = U s0,••• ,sm-1 (H.6a) U (l+1) r,s l+1 ,••• ,sm-1,t0,••• ,t l = S l -1 s l =0 K (l) r,s l ,t l U (l) r,s l ,••• ,sm-1,t0,••• ,t l-1 (H.6b) V t0,t1,••• ,tm-1 = R-1 r=0 U (m) r,t0,••• ,tm-1 (H.6c)

where the m intermediate results

U (l) ∈ R R×S l ×•••Sm-1×T0×•••×T l-1 , ∀l ∈ [m]

are (m+1)-order tensors, and the (l + 1) th step above (i.e. the equation interacting U (l) with K (l) ) requires O(( m R). These steps in tensor notations are presented as follows:

l k=0 S k ) ( m-1 k=l T k )R)

U (0) = 1 ⊗ U (H.7a) U (l+1) = U (l) ⊗ 0 0 • × 1 1 K (l) (H.7b) V = 1 × 0 0 U (m)

(H.7c) Following the procedure in Appendix E, their backpropagation equations are obtained as:

∂L ∂U (l) = K (l) ⊗ 0 0 • × 2 -1 ∂L ∂U (l+1) (H.8a) ∂L ∂K (l) = U (l) ⊗ 0 0 • × 2 1 • • • • • × m m-1 ∂L ∂U (l+1) (H.8b) ∂L ∂U (m) = 1 ⊗ ∂L ∂V , ∂L ∂U = 1 × 0 0 ∂L ∂U (0)

As we discussed in Appendix E, if a compound operation does not contain convolution, the time complexity of backpropagation is identical to the forward pass. Therefore, we claim the number of operations required for backward pass is also bounded by O(m max(S, T ) 1+ 1 m R).

r-Tucker-dense layer The application of TK decomposition is rather straightforward, which factors the tensor of parameters K exactly the same as in Appendix F.

K s0,••• ,sm-1,t0,••• ,tm-1 = R s 0 -1 r s 0 =0 • • • R s m-1 -1 r s m-1 =0 R t 0 -1 r t 0 =0 • • • R t m-1 -1 r t m-1 =0 P (0) s0,r s 0 • • • P (m-1) sm-1,r s m-1 C r s 0 ,••• ,r s m-1 ,r t 0 ,••• ,r t m-1 Q (0) r t 0 ,t0 • • • Q (m-1) r t m-1 ,tm-1 (H.9a) K = C × 0 (P (0) ) ⊤ • • • × m-1 (P (m-1) ) ⊤ × m Q (0) • • • × 2m-1 Q (m-1) (H.9b)

where

P (l) ∈ R S l ×R s l , ∀l ∈ [m] are named as input factors, C ∈ R R s 0 ×•••×R s m-1 ×R t 0 ×•••×R t

m-1 as core factor, and lastly

Q (l) ∈ R R t l ×T l , ∀l ∈ [m]

as output factors. Similar to Tucker-convolutional layer, the procedure to evaluate the result V can be broken into three steps, by plugging Equation H.9a in H.2a.

U (0) r s 0 ,••• ,r s m-1 = S0-1 s0=0 • • • Sm-1-1 sm-1=0 P (0) s0,r s 0 • • • P (m-1) sm-1,r s m-1 U s0,••• ,sm-1 (H.10a) U (1) r t 0 ,••• ,r t m-1 = R s 0 -1 r s 0 =0 • • • R s m-1 -1 r s m-1 =0 C r s 0 ,••• ,r s m-1 ,r t 0 ,•••r t m-1 U (0) r s 0 ,••• ,r s m-1 (H.10b) V t0,••• ,tm-1 = R t 0 -1 r t 0 =0 • • • R t m-1 -1 r t m-1 =0 Q (0) r t 0 ,t0 • • • Q (m-1) r t m-1 ,tm-1 U (1) r t 0 ,••• ,r t m-1 (H.10c)

where the first and last steps are compound operations between a tensor and a set of multiple tensors, while the middle step is a multi-operations between two tensors. Under the assumptions that

S l ≈ S 1 m , T l ≈ T 1 m , R s l = R t l = R, ∀l ∈ [m]

, the order to contract the factors in the first and last steps makes no difference, therefore we assume the order follows the indices without loss of generality. With this strategy, the contraction with the l th factor takes O(( by the definition of Tucker decomposition: therefore, the time complexity to contract all m factors is at most O(mSR). Likewise, the number of operations for the last step can also be bounded by O(mT R). Lastly, it is easy to see the middle step needs O(R 2m ) operations, therefore leads to a total time complexity of O(m(S + T )R + R 2m ) for the three-step procedure. In tensor notations, these equations can be concisely written as

l k=0 R k )( m-1 k=l S k )) operations,

U (0) = U × 0 P (0) • • • × m-1 P (m-1) (H.11a) U (1) = U (0) × 0 0 • • • • • × m-1 m-1 C (H.11b) V = U (1) × 0 Q (0) • • • × m-1 Q (m-1) (H.11c)

Though compound in nature, the procedure to derive their backpropagation rules are pretty straightforward: notice the equations for the first and last steps have the exactly the same form as standard Tucker decomposition in Appendix F. Therefore, we can simply modify the variable names therein to obtain the backpropagation equations for these two steps.

∂L ∂U (1) = ∂L ∂V × 0 (Q (0) ) ⊤ • • • × m-1 (Q (m-1) ) ⊤ (H.12a) ∂L ∂Q (l) ⊤ = ∂L ∂V × 0 0 • • • • • × l-1 l-1 • × l+1 l+1 • • • • • × m-1 m-1 U (1) × 0 Q (0) • • • × l-1 Q (l-1) × l+1 Q (l+1) • • • × m-1 Q (m-1) (H.12b) ∂L ∂U = ∂L ∂U (0) × 0 (P (0) ) ⊤ • • • × m-1 (P (m-1) ) ⊤ (H.12c) ∂L ∂P (l) ⊤ = ∂L ∂U (0) × 0 0 • • • • • × l-1 l-1 • × l+1 l+1 • • • • • × m-1 m-1 U × 0 P (0) • • • × l-1 P (l-1) × l+1 P (l+1) • • • × m-1 P (m-1) (H.12d)

The step in the middle is itself a tensorized layer defined in Equation H.2b, therefore its backpropagation rules can be obtained by renaming the variable in Equations H.4.

∂L ∂U:

(0) = C × m 0 • • • • • × 2m-1 m-1 ∂L ∂U (1) (H.13a) ∂L ∂C = U (0) ⊗ ∂L ∂U (1) (H.13b)

Despite their technical complicity, we can resort to conclusion that the complexities for forward and backward passes are the same for operations without convolution, and claim that the number of operations required for the backpropagation equations above is bounded by O(m(S + T )R + R 2m ).

r-Tensor-train-dense layer The layer presented in this part follows closely the pioneering work in compressing network using tensor decompositions (Novikov et al., 2015) , except that we replace the backpropagation algorithm in the original paper (as discussed in Appendix F) with a multistep procedure similar to all other layers in this paper. With the replacement, the efficiency of the backward is greatly improved compared to original design. Similar to r-CP-dense layer, we will group (S l , T l )'s as supermodes and decomposed the kernel K by Tensor-train decomposition following the order of their indices:

K s0,••• ,sm-1,t0,••• ,tm-1 = R0-1 r0=0 • • • Rm-2-1 rm-2=0 K (0) s0,t0,r0 K (1) r0,s1,t1,r1 • • • K (m-1) rm-2,sm-1,tm-1 (H.14a) K = swapaxes K (0) × -1 0 K (1) × -1 0 • • • × -1 0 K (m-1) (H.14b)

where the factor tensors are K

(l) ∈ R R l-1 ×S l ×T l ×R l , ∀l = 1, • • • , m -2, with two corner cases K (0) ∈ R S0×T0×R0 and K (m-1) ∈ R Rm-2×Sm-1×Tm-1 . The number of parameters in its l th factor K (l) is approximately (ST ) 1 m R 2 , therefore the layer has O(m(ST ) 1 m R 2 ) parameters in total. For aestheticism, we add singleton mode R -1 = 1 to U, K (0) such that U ∈ R S0×•••×Sm-1×R-1 , K (0) ∈ R R-1×S0×T0×R0

, and rename U as U (0) . Now insert the Equation H.14a into H.2a and expand accordingly, we obtain an (m + 1)-steps procedure to evaluate the output V:

U (l+1) s l+1 ,••• ,sm-1,t0,••• ,t l ,r l = R l-1 -1 r l-1 =0 S l -1 s l =0 K (l) r l-1 ,s l ,t l ,r l U (l) s l ,••• ,sm-1,t0,••• ,t l-1 ,r l-1 (H.15)

where the last result U (m) is equal to the final output V, and the other m tensors U (l) 's are intermediate results. Similar to r-CP-dense layer, the number of operations at the

l th step is O(( l k=0 S k ) ( m-1 k=l T l )R l-1 R l )

and can be bounded by O(max(S, T ) 1+ 1 m R 2 ), therefore the time complexity for all m-steps is bounded by O(m max(S, T ) 1+ 1 m R 2 ). These steps very simple in tensor notations.

U (l+1) = U (l) × 0 1 • × -1 0 K (l) (H.16)

Observe that the operations in the forward pass are entirely tensor contractions, therefore their backpropagation equations are easily derived following the procedure in Appendix E.

(l) = K (l) × 1 -2 • × 2 -1 ∂L ∂U (l+1) (H.17a) ∂L ∂K (l) = swapaxes U (l) × 1 0 • • • • • × m-1 m-2 ∂L ∂U (l+1) (H.17b)

In the analysis of the backward pass, we can again take advantage of the argument that the forward and backward passes share the same number of operations. Therefore, we claim the time complexity for backpropagation is bounded by O(m max(S, T ) 1+ 1 m R 2 ).

Relation to tensor contraction layer: In Kossaifi et al. (2017a) , the authors propose a novel tensor contraction layer, which takes a tensor of arbitrary order as input and return a tensor of the same order. Formally, a tensor contraction layer, parameterized by a set of m matrices {M (l) }) m-1 l=0 , with

M (l) ∈ R S l ×T l , ∀l ∈ [m], maps a m-order tensor U ∈ R S0×•••×Sm-1 to another m-order tensor V ∈ R T0×•••Tm-1 such that V t0,••• ,tm-1 = S0-1 s0=0 • • • Sm-1-1 sm-1=0 U s0,••• ,sm-1 M (0) s0,t0 • • • M (m-1) sm-1,tm-1 (H.18a) V = U × 0 M (0) • • • × m-1 M (m-1) (H.18b)

It is not difficult to observe that the tensor contraction layer is in fact special case of r-CP-dense layer where the kernel is restricted to rank-1, that is 1) in our tensor notations.

K s0,••• ,sm-1,t0,••• ,tm-1 = M (0) s0,t0 • • • M (m-1) sm-1,tm-1 , or equivalently K = M (0) ⊗ • • • ⊗ M (m-

Relation to tensor regression layer: Along with the tensor contraction layer, tensor regression layer is also proposed in Kossaifi et al. (2017b) , which takes a tensor of arbitrary order as input and maps it to a scalar. Formally, given an m-order tensor U ∈ R S0×•••×Sm-1 , it is reduced to a scalar v by contracting of all the modes with another tensor of the same size K

∈ R S0×•••×Sm-1 . v = S0-1 s0=0 • • • Sm-1-1 sm-1=0 U s0,••• ,sm-1 K s0,••• ,sm-1 (H.19a) v = U × 0 0 • × 1 1 • • • • • × m-1 m-1 K (H.19b)

where the tensor K is stored in Tucker-format as in Equation F.4a. Therefore, the tensor regression layer is effectively parameterized by a set of matrices

{M (l) } m-1 l=0 , M (l) ∈ R S l ×R l , ∀l ∈ [m], with an additional core tensor C ∈ R R0×•••×Rm-1

. Therefore the definition of tensor regression layer in Equation H.19a can also be rephrased as

v = S0-1 s0=0 • • • Sm-1-1 sm-1=0 U s0,••• ,sm-1 M (0) s0,r0 • • • M (m-1) sm-1,rm-1 C r0,••• ,rm-1 (H.20a) = U × 0 M (0) • • • × m-1 M (m-1) × 0 0 • • • • • × m-1 m-1 C (H.20b)

Now we are able to observe that the tensor regression layer is indeed a special case of r-Tuckerdense layer where the input factors P (l) = M (l) , ∀l ∈ [m], while the output factors Q (l) 's are simply scalar 1's with

T l = 1, ∀l ∈ [m].

Comments on the designs: As we can observe, the design of r-Tucker-dense is different from the other two layers using CP and Tensor-train decompositions, in which (S l , T l )'s are first grouped into supermodes before factorization. Indeed, it is a major drawback in the design of r-Tucker-dense layer: notice that the first intermediate result U (0) obtained in the forward pass is a tensor of size

R s 0 ×R s 1 • • •×R s m-1

, which becomes very tiny if the kernel K is aggressively compressed. Therefore, the size of the intermediate tensor poses an "information bottleneck", causing significant loss during the forward pass, which is verified by our experimental results in Section app:experiments. Therefore, the use of r-Tucker-dense layer is not recommended when we expect excessive compression rate. On the other hand, by grouping (S l , T l )'s as supermodes in r-CP-dense layer and r-Tensor-traindense layer, all intermediate tensors U (l) 's have similar size as the input, therefore the bottleneck in r-Tucker-dense layer is completely avoided.

But why do we not group (S l , T l )'s together in the design of Tucker-dense layer at the beginning? In theory, we are for sure able to factorize the kernel as

K = R0-1 r0=1 • • • Rm-1-1 rm-1=0 C r0,••• ,rm-1 K (0) r0,s0,t0 • • • K (m-1) rm-1,sm-1,tm-1 (H.21)

However, the contractions among the input and the factors become problematic: (1) interacting the input with the factors K (l) 's yields an intermediate tensor of size

T 0 × • • • × T m-1 × R 0 × • • • × R m-1

, which is too large to fit into memory; (2) while reconstructing the kernel K from K (l) 's and subsequently invoking the Equation H.2a will make the time complexity for backward pass intractable as we discussed in Appendix F. Therefore, we have to abandon this attempting design, in order to maintain a reasonable time complexity. As a compensation for the possible loss, the current design of Tucker-dense layer actually has one benefit over the other two layers: the numbers of operations for the backward pass remains at the same order as the number of parameters, while the number of operations required by r-CP-dense and r-Tensor-train-dense layers are orders higher. As a result, r-Tucker-dense layer is much faster than r-CP-dense and r-Tensor-train-dense layers at the same compression rate. Therefore, r-Tucker-dense layer is more desirable if we value speed over accuracy and the compression rate is not too high.

O(# of params.)

O(# of forward/backprop ops.) original ST ST SVD (S + T )R (S + T )R r-CP m(ST ) 1 m R m max(S, T ) 1+ 1 m R r-TK m(S 1 m + T 1 m )R + R 2m m(S + T )R + R 2m r-TT m(ST ) 1 m R 2 m max(S, T ) 1+ 1 m R 2

Table 11: Summary of reshaped tensor decomposition on dense layer. In this table, we list the numbers of parameters and time complexities of forward/backward passes required by various reshaped tensor decompositions on dense layer. For simplicity, we assume that the number of input units S and outputs units T are factorized evenly, i.e. S l = S pass are easily obtained through the results of r-CP-dense layer: compared to r-CP-dense layer, each of the existing m steps will be scaled by a factor of XY , while the additional last step requires O(HW T RXY ) operations. Therefore, the total number of operations for r-CP-convolutional layer is O(m max(S, T ) 1+ 1 m RXY + HW T RXY ). In tensor notations, these steps are rephrased as: U (0) = 1 ⊗ U (I.5a)

U (l+1) = U (l) ⊗ 0 0 • × 3 1 K (l) (I.5b) V = U (m) × 0 0 • * 1 1 • * 2 2 K (m) (I.5c)

The backpropagation equations are also very similar to their r-CP-dense layer counterparts. For completeness, we list all of them in the following: l+1) (I.6d)

∂L ∂U (m) = K (m) ( * 1 0 ) ⊤ • ( * 2 1 ) ⊤ ∂L ∂V (I.6a) ∂L ∂K (m) = U (m) ( * 1 0 ) ⊤ • ( * 2 1 ) ⊤ × 3 2 • • • • • × m+2 m+1 ∂L ∂V (I.6b) ∂L ∂U (l) = swapaxes K (l) ⊗ 0 0 • × 2 -1 ∂L ∂U (l+1) (I.6c) ∂L ∂K (l) = U (l) ⊗ 0 0 • ⊗ 1 1 • ⊗ 2 2 • × 4 3 • • • • • × m+3 m+2 ∂L ∂U (

∂L ∂U = 1 × 0 0 ∂L ∂U (0) (I.6e)

In the analyses of these backpropagation equations, it is again convenient to make connections to their r-CP-dense layer counterparts: the time complexities for the first m steps are scaled by XY , while the last step of 3D-convolutional layer requires O(HW T RX ′ Y ′ ) operations for the derivative with respect to U (m) , and O(XY T RX ′ Y ′ ) for the gradient of K (m) .

r-Tucker-convolutional layer Incorporating the features from both Tucker-convolutional layer in Appendix G and r-Tucker-dense layer in Appendix H, we propose to apply partial Tucker decomposition on the tensorized kernel K over all modes except the filter height H and width W . Concretely, the tensorized kernel K is factorized as:

K i,j,s0,••• ,sm-1,t0,••• ,tm-1 = R s 0 -1 r s 0 =0 • • • R s m-1 -1 r s m-1 =0 R t 0 -1 r t 0 =0 • • • R t m-1 -1 r t m-1 =0 P (0) s0,r s 0 • • • P (m-1) sm-1,r s m-1 C i,j,r s 0 ,••• ,r s m-1 ,r t 0 ,••• ,r t m-1 Q (0) r t 0 ,t0 • • • Q (m-1) r t m-1 ,tm-1 (I.7a) K = C × 2 (P (0) ) ⊤ • • • × m+1 (P (m-1) ) ⊤ × m+2 Q (0) • • • × 2m+1 Q (m-1) (I.7b)

where [m] are again named as input factors, core factor and output factors respectively. The reason that (S l , T l )'s are not grouped into supermodes follows exactly the same arguments as in Appendix H. Compared to r-Tensor-train-dense layer, the only difference is that the core tensor now has two extra modes for filter height and width, and therefore the number of parameters is magnified by a factor of HW . Similar to Tucker-convolutional and r-Tucker-dense layers, the procedure to evaluate V can be sequentialized into three steps: (I.8b)

P (l) ∈ R S l ×R s l , ∀l ∈ [m], C ∈ R H×W ×R s 0 ×•••×R s m-1 ×R t 0 ×•••×R t m-1 and Q (l) ∈ R R t l ×T l , ∀l ∈

V x,y,t0,••• ,tm-1 = R t 0 -1 r t 0 =0 • • • R t m-1 -1 r t m-1 =0 Q (0) r t 0 ,t0 • • • Q (m-1) r t m-1 ,tm-1 U (1)

x,y,r t 0 ,••• ,r t m-1 (I.8c)

where U (0) ∈ R X×Y ×R s 0 ×•••×R s m-1 and U (1) ∈ R X ′ ×Y ′ ×R t 0 ×•••×R t m-1 are two intermediate tensors. Referring to r-Tucker-dense layer, the number of operations of the first step is scaled up by XY , the number for the middle step by HW XY and the last step by X ′ Y ′ . Therefore, the time complexity for the three-steps process is O(mS 1 m XY + HW R 2m XY + mT 1 m X ′ Y ′ ). Subsequently, we can rewrite these steps concisely in tensor notations. U (0) = U × 2 P (0) • • • × m+1 P (m-1) (I.9a)

U (1) = U (0) * 0 0 • * 1 1 × 2 2 • • • • • × m-1 m-1 C (I.9b) V = U (1) × 2 Q (0) • • • × m+1 Q (m-1) (I.9c)

In principle, the backpropagation rules for the sequential steps can be derived almost identically as in the r-Tucker-dense layer. For reference, we list all equations for the first and last steps as follows: The analyses for the backpropagation equations mimic the ones in the forward pass, again by comparison against the ones in r-Tucker-dense layer: the time complexity to obtain the derivatives in the first step is magnified by XY , while the ones for the middle step and the last step are scaled by HW X ′ Y ′ and X ′ Y ′ respectively. Therefore, the total number of operations for the derivatives with respect to input/intermediate results is O(mS

∂L ∂U (1) = ∂L ∂V × 2 (Q (0) ) ⊤ • • • × m+1 (Q (m-1) ) ⊤ (I.10a) ∂L ∂Q (l) ⊤ = ∂L ∂V × 0 0 • • • • • × l+1 l+1 • × l+3 l+3 • • • • • × m+1 m+1 U (1) × 2 Q (0) • • • × l+1 Q (l-1) × l+3 Q (l+1) • • • × m+1 Q (m-1) (I.10b) ∂L ∂U = ∂L ∂U (0) × 2 (P (0) ) ⊤ • • • × m+1 (P (m-1) ) ⊤ (I.10c) ∂L ∂P (l) ⊤ = ∂L ∂U (0) × 0 0 • • • • • × l+1 l+1 • × l+3 l+3 • • • • • × m+1

1 m RXY + HW R 2m X ′ Y ′ + mT 1 m RXY )

, and the number for the factors is O(mS

1 m RXY + XY R 2m X ′ Y ′ + mT 1 m RX ′ Y ′ ).

r-Tensor-train-convolutional layer Following the r-CP-convolutional layer, we propose to apply Tensor-train decomposition to the tensorized kernel K by grouping (S l , T l )'s and filter height/width (H, W ) as supermodes. In Tensor-train decomposition, these supermodes are ordered by their indices, with the extra supermode (H, W ) appended to the end. Concretely, the tensorized kernel K is decomposed as:

K i,j,s0,••• ,sm-1,t0,••• ,tm-1 = R0-1 r0=0 • • • Rm-1-1 rm-1=0 K (0) s0,t0,r0 K (1) r0,s1,t1,r1 • • • K (m) rm-1,i,j (I.12a) K = swapaxes K (0) × -1 0 K (1) × -1 0 • • • × -1 0 K (m) (I.12b)

where K (0) ∈ R S0×T0×R0 , K (l) ∈ R R l-1 ×S l ×T l ×R l and K (m) ∈ R Rm-1×H×W are (m + 1) factor tensors. Compared to r-Tensor-train-dense layer, the r-Tensor-train-convolutional layer has an additional factor K (m) that contains RHW parameters, which leads to a total number of O((m(ST ) 1 m + HW )R). For conciseness, we follow the preprocessing steps to add singleton mode R -1 = 1 to U and K (0) such that U ∈ R X×Y ×S0×•••×Sm-1×R-1 and K (0) ∈ R R-1×S0×T0×R0 and rename U as U (0) . As we shall expect, the multi-steps procedure to evaluate V now has (m + 1) steps, with the last step as a 3D-convolutional layer:

U (l) i+dx,j+dy,••• ,r l = R l-1 -1 r l-1 =0 S l -1 s l =0 K (l)

r l-1 ,s l ,t l ,r l U 

U (l) = U (l) × 2 1 • × -1 0 K (l) (I.14a) V = U (m) * 0 1 • * 1 2 • × -1 0 K (m) (I.

14b) The backpropagation rules are easily obtained by modifying the ones in r-Tensor-train-dense layer. For completeness, we list all backpropagation equations in the following: 

∂L ∂U (m) = ∂L ∂V ( * 0 1 ) ⊤ • ( * 2 1 ) ⊤ K (m) (I.15a) ∂L ∂K (m) = ∂L ∂V ( * 0 0 ) ⊤ • ( * 1 1 ) ⊤ • × 2 2 • • • × m+1 m+1 U (m) (I.15b) ∂L ∂U (l) = swapaxes K (l) × 2 -2 • × 3 -1 ∂L ∂U (l+1) (I.

:

The optimal complexity of tensor algebra is NP complete in generalLam et al. (1997), therefore the complexity presented in this paper is the complexity of our implementation.

2 Assuming adding n terms takes n time in parallel for memory efficiency, although it could be O(log n).

† The reshaped tensor ring (r-TR) results are cited fromWang et al. (2018), and the accuracy of 80.8% is achieved by 6.67% compression rate.

Note that Tucker remains an exception for aggressive compression due to the low rank internal structure that we previously discussed.

Figures and Tables:

Figure 1: A toy example of invariant structures. The periodic and modulated structures are picked out by exploiting the low rank structure in the reshaped matrix.

Figure 2: Diagram of a ∈ R, v ∈ R I , M ∈ R I×J and T ∈ R I×J×K .

Figure 3: Tensor operation illustration. Examples of tensor operations in which M ∈ R J0×J1 , X ∈ R I0×I1×I2 and Y ∈ R J0×J1×J2 are input matrix/tensors, andT 1 ∈ R I1×I2×J0×J2 , T 2 ∈ R J0×I1×I2 , T 3 ∈ R I ′ 0 ×I1×I2×J0×J2and T 4 ∈ R I0×I1×I2×J0×J2 are output tensors of corresponding operations. Similar definitions apply to general mode-(i, j) tensor operations.

Figure 4: Diagrams of kernel decompositions.Figure (b) (c) and (d) are three types of plain tensor decomposition for a traditional convolutional kernel K in (a).Figure (f) (g) and (h) are three types of reshaped tensor decomposition for our tensorized kernel K ′ in (e) where the reshaping order m ∈ Z is chosen to be 3 for illustrative simplicity.

Figure 4: Diagrams of kernel decompositions.Figure (b) (c) and (d) are three types of plain tensor decomposition for a traditional convolutional kernel K in (a).Figure (f) (g) and (h) are three types of reshaped tensor decomposition for our tensorized kernel K ′ in (e) where the reshaping order m ∈ Z is chosen to be 3 for illustrative simplicity.

Figure 5: Convolutional layer diagram. Input U is passed through the layer kernel K. The forward propogation operation of an uncompressed layer, a plain tensor decomposition compressed layer and our reshaped tensor decomposition compressed layer are illustrated in (a), (b) and (c) respectively.

Figure 6: Convergence rate for Seq vs. E2E tuning on CIFAR10.

of params.) O(# of forward ops.) O(# of backward ops. for inputs) O(# of backward ops. for params.

operations. The number can be bounded by O(max(S, T ) 1+ 1 m R), and equality is achieved when S l = T l = S 1 m = T 1 m . Since the first and last steps require O(SR) and O(T R) that are negligible compared to the other m steps, therefore the total number of operations is bounded by O(m max(S, T ) 1+ 1

which is roughly O(R l+1 S m-l m ) and can be further bounded by O(SR) since R ≤ S 1 m

, T l = T 1 m , ∀l ∈ [m] and all ranks (in r-TK and r-TT) share the same number R, i.e. R l = R, ∀l ∈ [m].

=0 i,j C i,j,r s 0 ,••• ,r s m-1 ,r t 0 ,•••r t m-1 U (0) i+dx,j+dy,r s 0 ,••• ,r s m-1

(0) • • • × l+1 P (l-1) × l+3 P (l+1) • • • × m+1 P (m-1) (I.10d)Notice that the middle step itself is a tensorized convolutional layer defined in Equation I.1b, therefore its backpropagation equations are exactly the same as the ones in Equations I.2a and I.2b.

,s l ,••• ,sm-1,t0,••• ,t l-1 ,r l-1 l) ∈ R X×Y ×S l ×•••×Sm-1×T0×•••×T l-1 ×R l-1 , ∀l ∈ [m]are the intermediate results. The number of operations for the first m steps is O(m max(S, T ) 1+ 1m RXY ) by comparing against the corresponding steps in r-Tensor-train-dense layer, and the last step of 3D-convolutional layer requires O(HW RT XY ) operations in the forward pass. In tensor notations, these steps are nicely represented as:

the backward steps again follow by building connections with r-Tensor-train-dense layer: each of the first m steps is scaled by XY , while the last step requires O(HW T RX ′ Y ′ ) for U(m) and O(XY T RX ′ Y ′ ) for K(m) .Decomp.O(# of params) O(# of forward ops.) O(# of back ops. for inputs) O(# of back ops. for params.)original HW ST HW ST XY HW ST X ′ Y ′ XY ST X ′ Y ′ r-CP #Params CP + HW R #Ops CP XY + HW T RXY #Ops CP XY + HW T RX ′ Y ′ #Ops CP XY + XY T RX ′ Y ′ r-TK (HW #Params TKC + #Params TKI + #Params TKO (HW #Ops TKC XY + #Ops TKI XY + #Ops TKO X ′ Y ′ ) (HW #Ops TKC X ′ Y ′ + #Ops TKI XY + #Ops TKO X ′ Y ′ ) (XY #Ops TKC X ′ Y ′ + #Ops TKI XY + #Ops TKO X ′ Y ′ ) r-TT #Params TT + HW R #Ops TT XY + HW T RXY #Ops TT XY + HW T RX ′ Y ′ #Ops TT XY + XY T RX ′ Y ′ †In order to compare against reshaped tensor decompositions on dense layer in Table11, we denote the numbers for the dense layers as:I = mSR #OpsTK C = R 2m #OpsTK O = mT R

, T l = T 1 m , ∀l ∈ [m], and all ranks are equal to R, i.e R l = R, ∀l ∈ [m].

Number of parameters and operations required by a compressed convolutional layer by various types of tensor decompositions (for the special case of

Table content: <table/>

Parallel time complexity of forward pass using various types of tensor decompositions on convolutional layers. The uncompressed parallel complexity of forward pass is O(k 2 N ).

Table content: <table/>

Percentage test accuracy of baseline PD with E2E tuning vs. our RTD with Seq tuning on CIFAR10. The uncompressed ResNet-32 achieves 93.2% accuracy with 0.46M parameters.

Table content: <table><tr><td/><td>5%</td><td>10%</td><td>rate 20%</td><td>40%</td><td>Decomp.</td><td>2%</td><td colspan="2">Compression rate 5% 10%</td><td>20%</td></tr><tr><td>SVD</td><td colspan="4">83.09 87.27 89.58 90.85</td><td>r-TR  †</td><td>-</td><td>80.80  †</td><td>-</td><td>90.60</td></tr><tr><td>CP</td><td colspan="4">84.02 86.93 88.75 88.75</td><td>r-CP</td><td>85.7</td><td colspan="2">89.86 91.28</td><td>-</td></tr><tr><td>TK</td><td colspan="4">83.57 86.00 88.03 89.35</td><td>r-TK</td><td colspan="3">61.06 71.34 81.59 87.11</td></tr><tr><td>TT</td><td colspan="4">77.44 82.92 84.13 86.64</td><td>r-TT</td><td colspan="3">78.95 84.26 87.89</td><td>-</td></tr></table>

These results confirm the existence of extra invariant structure in the parameter space of deep neural networks. Such invariant structure is picked up by our proposed aproach, tensorization combined with low rank approximation (i.e., our RTD), but not by low rank approximation itself (i.e., baseline PD). Therefore, our results show that RTD and Seq tuning are symbiotic, and both are necessary to simultaneously obtain a high accuracy and a high compression rate.

Table content: <table><tr><td>shows that RTD with Seq tuning outperforms PD</td></tr><tr><td>with end-to-end tuning. Now we address the following</td></tr><tr><td>question: is one factor (Seq tuning or RTD) primarily re-</td></tr><tr><td>sponsible for increased performance, or is the benefit due</td></tr><tr><td>to synergy between the two?</td></tr><tr><td>Seq tuning, reshaped tensor decomposition, or both?</td></tr><tr><td>(1) We present the effect of different tuning methods on</td></tr><tr><td>accuracy in Table 4. Other than at very high compres-</td></tr><tr><td>sion rate (5% column in Table 4), Seq tuning (Seq) consis-</td></tr><tr><td>tently outperforms end-to-end (E2E) tuning. In addition,</td></tr><tr><td>Seq tuning is also much faster and leads to more stable</td></tr><tr><td>convergence compared to end-to-end tuning. Figure 6 plots the compression error over the number</td></tr><tr><td>of gradient updates for various tuning methods. (2) We present the effect of different compression</td></tr><tr><td>methods on accuracy in Table 5. Interestingly, if our RTD is used, the test accuracy is restored for</td></tr></table>

Percentage accuracy of our Seq vs. baseline E2E tuning using PD on CIFAR10.

Table content: <table><tr><td>Decomp.</td><td colspan="3">Compression rate Decomp. 5% 10%</td><td colspan="2">Compression rate 5% 10%</td></tr><tr><td>CP</td><td>83.19</td><td>88.50</td><td>r-CP</td><td>89.86</td><td>91.28</td></tr><tr><td>TK</td><td>80.11</td><td>86.73</td><td>r-TK</td><td>71.34</td><td>81.59</td></tr><tr><td>TT</td><td>80.77</td><td>87.08</td><td>r-TT</td><td>84.26</td><td>87.89</td></tr></table>

Percentage accuracy of our RTD vs. baseline PD using Seq tuning on CIFAR10.

Table content: <table><tr><td colspan="5">Scalability Finally, we show that our methods scale to state-of-the-art large networks, by evaluat-</td></tr><tr><td colspan="5">ing performance on the ImageNet 2012 dataset on a 50-layer ResNet (uncompressed with 76.05%</td></tr><tr><td colspan="5">accuracy). Table 6 shows the accuracy of RTD (TT decomposition) with Seq tuning compared to</td></tr><tr><td colspan="5">plain tensor decomposition with E2E tuning and the uncompressed network, on ResNet-50 with</td></tr><tr><td colspan="5">10% compression rate. Table 6 shows that Seq tuning of RTD is faster than the alternative. This is</td></tr><tr><td colspan="5">an important result because it empirically validates our hypotheses that (1) our RTD compression</td></tr><tr><td colspan="5">captures the invariance structure of the ResNet (with few redundancies) better and faster than the</td></tr><tr><td colspan="5">baseline PD compression, (2) data reconstruction Seq tuning is effective even on the largest networks</td></tr><tr><td colspan="5">and datasets, and (3) our proposed efficient RTD compression methods scale to the state-of-the-art</td></tr><tr><td>neural networks.</td><td/><td/><td/><td/></tr><tr><td/><td/><td>Uncompressed</td><td>TT (E2E)</td><td>r-TT (Seq)</td></tr><tr><td colspan="5"># epochs # samples # params. = 25M # params. = 2.5M # params. = 2.5M</td></tr><tr><td>0.2</td><td>0.24M</td><td>4.22</td><td>2.78</td><td>44.35</td></tr><tr><td>0.3</td><td>0.36M</td><td>6.23</td><td>3.99</td><td>46.98</td></tr><tr><td>0.5</td><td>0.60M</td><td>9.01</td><td>7.48</td><td>49.92</td></tr><tr><td>1.0</td><td>1.20M</td><td>17.3</td><td>12.80</td><td>52.59</td></tr><tr><td>2.0</td><td>2.40M</td><td>30.8</td><td>18.17</td><td>54.00</td></tr></table>

Convergence of percentage accuracies of uncompressed vs. PD (TT decomposition) vs.

Table content: <table><tr><td>RTD (r-TT decomposition) achieving 10% compression rate for ResNet-50 ImageNet.</td></tr></table>

Reshaped

Table content: <table><tr><td/><td colspan="2">Compression rate</td></tr><tr><td colspan="2">Method 0.2% 0.5%</td><td>1%</td></tr><tr><td>r-CP</td><td colspan="2">97.21 97.92 98.65</td></tr><tr><td>r-TK</td><td colspan="2">97.71 98.56 98.52</td></tr><tr><td>r-TT</td><td colspan="2">97.69 98.43 98.63</td></tr></table>



Table content: <table/>

Rt×T are factors, which require (SR s +HR s R+W RR t +R t T ). Once we plug the decomposition scheme in Equation G.16a into G.1, the evaluation of V is decoupled into four steps , with number of operations asO(R s SXY ), O(HR s RXY ), O(W RR t X ′ Y ) and O(R t T X ′ Y ′ ) respectively,resulting in a total time complexity as O(R s SXY + HR s

Table content: <table/>

Summary of reshaped tensor decomposition on convolutional layer. In this table, we list the number of parameters and time complexities required by various reshaped tensor decomposition s on convolutional layer. Recall that a convolutional layer, composed with ST filters of size H × W , maps a set of S feature maps of size X × Y to another set of T feature maps of size X ′ × Y ′ . For simplicity, we assume the numbers of input/output feature maps S, T are factorized evenly, i.e. S

Table content: <table/>

Peer review (iclr2019):
 The submission suggests a methodology compressing the parameters of neural network (NN) architectures, by first converting them to higher-order tensors and then factorizing them with various tensor decompositions. The main assumption of this work is that there may be structure in lower-order tensors (e.g., matrices or low-order tensors) contained as parameters of a NN such as: periodicity (e.g., the vector [1,2,3,1,2,3,1,2,3]) or modularity (e.g., the vector [1,1,1,2,2,2,3,3,3]). The main argument of this work is that such invariant structures in lower-order objects often go unexploited in previous low-rank approximation techniques tackling the problem of compressing the parameters of a deep neural network. The authors argue that such structures may be revealed by reshaping lower-order tensors into higher-order ones and factorizing the resulting tensors. Despite the fact that the results seem promising, the paper suffers from the following issues: 1. Tensorizing to a higher-order tensor with lower mode sizes and factorizing through tensor network approaches is not a new approach, as suggested in the main contributions. The authors should more clearly position this part of the work with respect to existing literature. I provide more details on this point below. 2. There is no experiment or discussion to back up the assumption that: 1) the invariant structures suggested by the authors do appear within the parameters of a neural network architecture; 2) such invariant structures are captured by the higher-order tensor decomposition proposed. As such, it seems that there is no solid intuition presented to justify why the approach works. Even if the authors stress out in the experiments that: "These results confirm the existence of extra invariant structure in the parameter space of deep neural networks", this is really not clear/well-justified. 3. The experiment setup is missing crucial discussions on the hyper-parameter search. For example, as the tensor order increases, there are many more potential choices for setting the low rank parameters of a tensor train decomposition, as well as the corresponding mode sizes of the input tensor to be decomposed. How does this increased space of hyper-parameters to search for affects performance seems a crucial performance factor, but it is not discussed at all. 4. The description of the sequential fine-tuning procedure seems really hand-waving. Despite the lengthy appendix, there is not even a pointer from the sequential fine-tuning section so that the reader can go over this process in more detail with precise notation. The lack of any detailed description of this step raises doubts regarding this step*s novelty/methodological significance. Additional comments w.r.t. first issue: The idea of reshaping a lower-order object into a higher-order tensor with low mode sizes and then factorize it does exist in the literature. Especially, in combination with tensor train decomposition, the framework has been called as the Quantized-Tensor Train (QTT). For example, the authors could check Section 2.9 of the survey: Grasedyck, Lars, Daniel Kressner, and Christine Tobler. "A literature survey of low‐rank tensor approximation techniques." GAMM‐Mitteilungen 36.1 (2013): 53-78. and references therein. I think the authors should mention the connections existing with this line of work and argue about how their approach differs. Admittedly, it seems that this is the first effort of applying those ideas in the context of compressing neural network (NN) parameters, as well as deriving the update rules w.r.t. a NN; still, appropriate credit should be given to the QTT line of work. Another issue w.r.t. comparing with related work is that the authors should definitely try to include the comparison and discussions with other low-rank approximation-based NN compression in the main text, rather than the appendix. Minor comments: Tensor partial outer product should be more clearly defined in the main text. Currently, it is really not clear to the reader what this operation signifies. Also, it is not clear on how to interpret CP tensor network diagrams as they are presented by the authors. It seems that the current presentation does not take into account the diagonal core tensor that is present in CP.