On Model Compression for Neural Networks: Framework, Algorithm, and
  Convergence Guarantee
Model compression is a crucial part of deploying neural networks (NNs),
especially when the memory and storage of computing devices are limited in many
applications. This paper focuses on two model compression techniques: low-rank
approximation and weight pruning in neural networks, which are very popular
nowadays. However, training NN with low-rank approximation and weight pruning
always suffers significant accuracy loss and convergence issues. In this paper,
a holistic framework is proposed for model compression from a novel perspective
of nonconvex optimization by designing an appropriate objective function. Then,
we introduce NN-BCD, a block coordinate descent (BCD) algorithm to solve the
nonconvex optimization. One advantage of our algorithm is that an efficient
iteration scheme can be derived with closed-form, which is gradient-free.
Therefore, our algorithm will not suffer from vanishing/exploding gradient
problems. Furthermore, with the Kurdyka-{\L}ojasiewicz (K{\L}) property of our
objective function, we show that our algorithm globally converges to a critical
point at the rate of O(1/k), where k denotes the number of iterations. Lastly,
extensive experiments with tensor train decomposition and weight pruning
demonstrate the efficiency and superior performance of the proposed framework.
Our code implementation is available at https://github.com/ChenyangLi-97/NN-BCD
Chenyang Li, Jihoon Chung, Mengnan Du, Haimin Wang, Xianlian Zhou, Bo Shen
# Provable Convergence of Tensor Decomposition-Based Neural Network Training

# Provable Convergence of Tensor Decomposition-Based Neural Network Training

**Chenyang Li, Bo Shen* Department of Mechanical and Industrial Engineering, New Jersey Institute of Technology *Corresponding Author: bo.shen@njit.edu**

**Abstract**

Advanced tensor decomposition, such as tensor train (TT), has been widely studied for tensor decomposition-based neural network (NN) training, which is one of the most common model compression methods. However, training NN with tensor decomposition always suffers significant accuracy loss and convergence issues. In this paper, a holistic framework is proposed for tensor decomposition-based NN training by formulating TT decomposition-based NN training as a nonconvex optimization problem. This problem can be solved by the proposed tensor block coordinate descent (tenBCD) method, which is a _gradient-free_ algorithm. The global convergence of tenBCD to a critical point at a rate of \(\mathcal{O}(1/k)\) is established with the Kurdyka Lojasiewicz (KL) property, where \(k\) is the number of iterations. The theoretical results can be extended to the popular residual neural networks (ResNets). The effectiveness and efficiency of our proposed framework are verified through an image classification dataset, where our proposed method can converge efficiently in training and prevent overfitting.

**Keywords**

Model Compression, Tensor Train Decomposition, Global Convergence, Gradient-free Training.

## 1 Introduction

Neural network (NN) has revolutionized many facets of our modern society, such as image classification [1], object detection [2, 3], speech recognition [4], etc. These advances have become possible because of algorithmic advances, large amounts of available data, and modern hardware. Despite their widespread success and popularity, there still remains a significant challenge in executing NNs with many parameters on edge devices. For most embedded and Internet-of-Things (IoT) systems, the sizes of many state-of-the-art NN models are too large, thereby causing high storage and computational demands and severely hindering the practical deployment of NNs. For example, wearable robots [5, 6], such as exoskeletons, typically have limited processing power, memory, storage, and energy supply due to their small size and portability. In addition, these wearable devices rely on wireless communication with remote servers, as larger models would require more bandwidth and higher latency, leading to slower and less reliable performance.

To address this issue, numerous model compression techniques are proposed in the literature, which can be summarized into the following categories. (1) Pruning [7, 8, 9]: this technique involves removing unnecessary connections or neurons from a pre-trained model. This can result in a smaller network with similar performance. (2) Quantization [10, 11]: this involves reducing the number of bits required to represent the weights and activations in a neural network. For example, weights and activations may be represented using 8-bit integers instead of 32-bit floating-point numbers. (3) Structured sparsity [12]: this involves imposing a structured sparsity pattern on the weights of a model, such as by sparsifying entire rows or columns of weight matrices. (4) Knowledge
distillation [13]: this involves training a smaller model to mimic the behavior of a larger, more complex model, using the outputs of the larger model as labels. (5) Low-rank approximation [14]: this technique involves approximating the weight matrices/tensors of a deep learning model with low-rank matrices/tensors.

Among all model compression methods, low-rank approximation, especially tensor decomposition [15], is an extremely attractive NN model compression technique since it can reduce the number of parameters in a model while maintaining a high level of accuracy. Specifically, tensor decomposition is a mathematical tool that explores the low tensor rank characteristics of large-scale tensor data, which stands out by offering an ultra-high compression ratio. By utilizing advanced tensor decomposition techniques like tensor train (TT) [16], it is possible to achieve more than a 1,000\(\times\) reduction in parameters for the input-to-hidden layers of neural network models [17; 18]. Moreover, these compression methods can also enhance the classification accuracy in video recognition tasks significantly. Given such impressive compression performance, there has been a surge of interest in exploring the potential of tensor decomposition-based neural network models in prior research efforts [19]. Due to the benefits brought by the TT-based NN models, several TT-based NN hardware accelerators have been developed and implemented in different chip formats including digital CMOS ASIC [20], memristor ASIC [21] and IoT board [22].

Although tensor decomposition shows strong compression performance, the training of tensor decomposition-based NN is a quite challenging task [23; 19] because it involves tensor decomposition in NN training. In general, there are two ways to use tensor decomposition to obtain a compressed model: (1) Train from scratch in the decomposed format, and (2) Decompose a pre-trained uncompressed model and then retrain. In the first case, when the required tensor decomposition-based, e.g. TT-format model, is directly trained from scratch because the structure of the models is already pre-set to low tensor rank format before the training, the corresponding model capacity is typically limited as compared to the full-rank structure, thereby causing the training process being very sensitive to initialization and more challenging to achieve high accuracy. In the latter scenario, though the pre-trained uncompressed model provides a good initialization position, the straightforwardly decomposing full-rank uncompressed model into low tensor rank format causes inevitable and non-negligible approximation error, which is still very difficult to be recovered even after a long-time re-training period.

No matter which training strategy with tensor decomposition is adopted, the training of NN heavily relies on gradient-based methods, which make use of backpropagation [24] to compute gradients of network parameters. These gradient-based methods are based on the Stochastic Gradient Descent (SGD) method [25]. In recent years, a considerable amount of research has been dedicated to developing adaptive versions of the vanilla SGD algorithm. These adaptive variants include AdaGrad [26], RMSProp [27], Adam [28], and AMSGrad [29]. Despite the great success of these gradient-based methods, tensor decomposition always brings a linear increase in network depth, which implies training the tensor decomposition format NNs are typically more prone to the gradient vanishing problem [30] and hence being difficult to be trained well.

This paper aims to address the current limitations and fully unlock the potential of tensor decomposition-based NN training. To achieve this objective, a holistic framework for tensor decomposition-based NN training is proposed, which formulates tensor train decomposition-based NN training as a nonconvex optimization problem. This problem can be solved by the proposed tensor block coordinate descent (tenBCD) methods. BCD is a _gradient-free_ method that has been
recently adapted to NN training [31, 32]. The main reasons for the surge of attention of BCD algorithms are twofold. One reason is that they are _gradient-free_, and thus are able to deal with non-differentiable nonlinearities and potentially avoid the vanishing gradient issue. The other reason is that BCD can be easily implemented in a distributed and parallel manner, therefore in favor of distributed/federated scenarios. To summarize, the contributions of this paper are as follows:

* A holistic framework is proposed for tensor decomposition-based NN training, which involves a highly nonconvex optimization problem.
* An efficient tensor BCD (tenBCD) algorithm is implemented to solve the proposed optimization problem;
* Convergence of the iterative sequence generated by the tenBCD algorithm is analyzed, which is proved to be globally convergent to a critical point at a rate of \(\mathcal{O}(1/k)\).

## 2 Background and Preliminaries

In Section 2.1, the notation and basics of multi-linear/tensor algebra used in this paper are reviewed. Then, tensor train decomposition [16] is reviewed briefly in Section 2.2. Afterward, the tensor train fully-connected layer [33] is reviewed in Section 2.3.

### Notation and Tensor Basis

Throughout this paper, scalars are denoted by lowercase letters, e.g., \(x\); vectors are denoted by lowercase boldface letters, e.g., \(\mathbf{x}\); matrices are denoted by uppercase boldface, e.g., \(X\); and tensors are denoted by calligraphic letters, e.g., \(X\). The order of a tensor is the number of its modes or dimensions. A real-valued tensor of order-\(d\) is denoted by \(\mathbf{\mathcal{X}}\in\mathbb{R}^{n_{1}\times n_{2}\times\cdots\times n_{d}}\) and its entries by \(\mathbf{\mathcal{X}}(i_{1},\cdots,i_{d})\). The inner product of two same-sized tensors \(\mathbf{\mathcal{X}}\) and \(\mathbf{\mathcal{Y}}\) is the sum of the products of their entries, namely, \(\langle\mathbf{\mathcal{X}},\mathbf{\mathcal{Y}}\rangle=\sum_{i_{1}}\cdots\sum_{i_{d }}X\left(i_{1},\ldots,i_{d}\right)\cdot\mathbf{\mathcal{Y}}\left(i_{1},\ldots,i_{ d}\right)\). Following the definition of inner product, the Frobenius norm of a tensor \(\mathbf{\mathcal{X}}\) is defined as \(\left\|\mathbf{\mathcal{X}}\right\|_{F}=\sqrt{\langle\mathbf{\mathcal{X}},\mathbf{ \mathcal{X}}\rangle}\).

### Tensor Train (TT) Decomposition

Given a tensor \(\mathbf{\mathcal{A}}\in\mathbb{R}^{n_{1}\times n_{2}\times\cdots\times n_{d}}\), it can be decomposed to a sort of 3-order tensors via Tensor Train Decomposition (TTD) [16] as follows:

\[\begin{split}\mathbf{\mathcal{A}}(i_{1},i_{2},\cdots,i_{d})& =\mathbf{\mathcal{G}}_{1}(:,i_{1},:)\mathbf{\mathcal{G}}_{2}(:,i_{2},:) \cdots\mathbf{\mathcal{G}}_{d}(:,i_{d},:)\\ &=\sum_{\alpha_{0},\alpha_{1}\cdots\alpha_{d}}^{r_{0},r_{1}, \cdots r_{d}}\mathbf{\mathcal{G}}_{1}(\alpha_{0},i_{1},\alpha_{1})\mathbf{\mathcal{G} }_{2}(\alpha_{1},i_{2},\alpha_{2})\cdots\mathbf{\mathcal{G}}_{d}(\alpha_{d-1},i_{ d},\alpha_{d}),\end{split} \tag{1}\]

where \(\mathbf{\mathcal{G}}_{k}\in\mathbb{R}^{r_{k-1}\times n_{k}\times r_{k}}\) are called TT-cores for \(k=1,2,\cdots,d\), and \(\mathbf{r}=[r_{0},r_{1},\cdots,r_{d}],r_{0}=r_{d}=1\) are called TT-ranks, which determine the storage complexity of TT-format tensor. The representation of \(\mathbf{\mathcal{A}}\) via the explicit enumeration of all its entries requires storing \(\Pi_{k=1}^{d}n_{k}\) numbers compared with \(\sum_{k=1}^{d}n_{k}r_{k-1}r_{k}\) numbers if the tensor is stored in TT-format.

### Tensor Train Fully-Connected Layer

Consider a simple fully-connected layer with weight matrix \(\mathbf{W}\in\mathbb{R}^{M\times N}\) and input \(\mathbf{x}\in\mathbb{R}^{N}\), where \(M=\prod_{k=1}^{d}m_{k}\) and \(N=\prod_{k=1}^{d}n_{k}\), the output \(\mathbf{y}\in\mathbb{R}^{M}\) is obtained by \(\mathbf{y}=\mathbf{W}\mathbf{x}\). In order to transform this standard layer to TT fully-connected (TT-FC) layer, the weight matrix \(\mathbf{W}\) is first tensorized to a \(d\)-order weight tensor \(\mathbf{\mathcal{W}}\in\mathbb{R}^{(m_{1}\times n_{1})\times\cdots\times(m_{d} \times n_{d})}\) by reshaping and order transposing. Then \(\mathbf{\mathcal{W}}\) can be decomposed to TT-format:

\[\mathbf{\mathcal{W}}((i_{1},j_{1}),\cdots,(i_{d},j_{d}))=\mathbf{\mathcal{G}}_{1}(:,i_{ 1},j_{1},:)\cdots\mathbf{\mathcal{G}}_{d}(:,i_{d},j_{d},:) \tag{2}\]
Here, each TT-core \(\mathbf{\mathcal{G}}_{k}\in\mathbb{R}^{r_{k-1}\times m_{k}\times n_{k}\times r_{k}}\) is a 4-order tensor, which is one dimension more than the standard one (1) since the output and input dimensions of \(\mathbf{W}\) are divided separately. Hence, the forward propagation on the TT-FC layer can be expressed in the tensor format as follows (the bias term is ignored here):

\[\mathbf{\mathcal{Y}}(i_{1},\cdots,i_{d})=\sum_{j_{1},\cdots,j_{d}}\mathbf{\mathcal{G}}_{ 1}(:,i_{1},j_{1},:)\cdots\mathbf{\mathcal{G}}_{d}(:,i_{d},j_{d},:)\mathbf{\mathcal{X}}( j_{1},\cdots,j_{d})\]

where \(\mathbf{\mathcal{X}}\in\mathbb{R}^{m_{1}\times\cdots\times m_{d}}\) and \(\mathbf{\mathcal{Y}}\in\mathbb{R}^{n_{1}\times\cdots\times n_{d}}\) are the tensorized input and output corresponding to \(\mathbf{x}\) and \(\mathbf{y}\), respectively. The details about the TT-FC layer are introduced in [33]. As the TT-FC layer and the corresponding forward propagation schemes are formulated, standard stochastic gradient descent (SGD) algorithm can be used to update the TT-cores with the rank set \(\mathbf{r}\), which determines the target compression ratio. The initialization of the TT-cores can be either randomly set or obtained from directly TT-decomposing a pre-trained uncompressed model.

## 3 Proposed Methodology

Consider \(N\)-layer feedforward neural networks with \(N-1\) hidden layers of the neural networks. Particularly, let \(n_{i}\in\mathbb{N}\) be the number of hidden units in the \(i\)-th hidden layer for \(i=1,\ldots,N-1\). Let \(n_{0}\) and \(n_{N}\) be the number of units of input and output layers, respectively. Let \(\mathbf{W}_{i}\in\mathbb{R}^{n_{i}\times n_{i-1}}\) be the weight matrix between the \((i-1)\)-th layer and the \(i\)-th layer for any \(i=1,\ldots,N\). Let \(\mathcal{Z}:=\{(\mathbf{x}_{j},\mathbf{y}_{j})\}_{j=1}^{n}\subset\mathbb{R}^{n_{0}} \times\mathbb{R}^{n_{N}}\) be \(n\) samples, where \(\mathbf{y}_{j}\)'s are the one-hot vectors of labels. Denote \(\mathcal{W}:=\{\mathbf{W}_{i}\}_{i=1}^{N},\mathbf{X}:=(\mathbf{x}_{1},\mathbf{x}_{2},\ldots, \mathbf{x}_{n})\in\mathbb{R}^{n_{0}\times n}\) and \(\mathbf{Y}:=(\mathbf{y}_{1},\mathbf{y}_{2},\ldots,\mathbf{y}_{n})\in\mathbb{R}^{n_{N}\times n}\).

### Problem Formulation

As shown in Figure 1, the weight in \(i\)-th layer, namely, \(\mathbf{W}_{i}\), can be transformed into a tensor \(\mathbf{\mathcal{W}}_{i}\). The tensor can be further decomposed into TT-format. Therefore, tensor train decomposition-based NN training problem can be formulated as the following empirical risk (i.e., training loss) minimization:

\[\min_{\mathcal{W}}\mathcal{R}_{n}(\Phi(\mathbf{X};\mathcal{W}),\mathbf{Y}),\text{ subject to }\mathbf{\mathcal{W}}_{i}=\text{TTD}(\mathbf{r}_{i})\quad i=1,\ldots,N \tag{3}\]

Figure 1: The framework of tensor train decomposition-based NN training.
where \(\mathcal{R}_{n}(\Phi(\mathbf{X};\mathcal{W}),\mathbf{Y}):=\frac{1}{n}\sum_{j=1}^{n}\ell( \Phi(\mathbf{x}_{j};\mathcal{W}),\mathbf{y}_{j})\) with loss function \(\ell:\mathbb{R}^{n_{N}}\times\mathbb{R}^{n_{N}}\rightarrow\mathbb{R}_{+}\cup \{0\}\), \(\Phi(\mathbf{x}_{j};\mathcal{W})=\sigma_{N}(\mathbf{W}_{N}\sigma_{N-1}(\mathbf{W}_{N-1} \cdots\mathbf{W}_{2}\sigma_{1}(\mathbf{W}_{1}\mathbf{x}_{j})))\) is the neural network model with \(N\) layers. TTD\((\mathbf{r}_{i})\) is the tensor train decomposition with rank \(\mathbf{r}_{i}\) in (2) for weight tensor \(\mathbf{\mathcal{W}}_{i}\) and \(\sigma_{i}\) is the activation function of the \(i\)-th layer (generally, \(\sigma_{N}\equiv\text{Id}\) is the identity function).

Note that the NN training model (3) is highly nonconvex as the variables are coupled via the NN architecture, which brings many challenges for the design of efficient training algorithms and also its theoretical analysis. To make Problem (3) more computationally tractable, variable splitting is one of the most commonly used ways [31, 32]. The main idea of variable splitting is to transform a complicated problem (where the variables are coupled nonlinearly) into a relatively simpler one (where the variables are coupled much looser) by introducing some additional variables.

Considering general NN architectures, the regularized NN training model is applied here, which can reduce the original NN training model (3). Specifically, the variable splitting model is:

\[\min_{\mathcal{W},\mathcal{V}}\mathcal{L}_{0}(\mathcal{W},\mathcal{ V}) :=\mathcal{R}_{n}(\mathbf{V}_{N};\mathbf{Y})+\sum_{i=1}^{N}\tau_{i}(\mathbf{W} _{i})+\sum_{i=1}^{N}s_{i}(\mathbf{V}_{i})\] (4) subject to \[\mathbf{U}_{i} =\mathbf{W}_{i}\mathbf{V}_{i-1},\mathbf{V}_{i}=\sigma_{i}(\mathbf{U}_{i}),\mathbf{ \mathcal{W}}_{i}=\text{TTD}(\mathbf{r}_{i})\quad i=1,\ldots,N,\]

where \(\mathcal{R}_{n}(\mathbf{V}_{N};\mathbf{Y}):=\frac{1}{n}\sum_{j=1}^{n}\ell((\mathbf{V}_{N})_ {:j},\mathbf{y}_{j})\) denotes the empirical risk, \(\mathcal{V}:=\{\mathbf{V}_{i}\}_{i=1}^{N},(\mathbf{V}_{N})_{:j}\) is the \(j\)-th column of \(\mathbf{V}_{N}\). In addition, \(\tau_{i}\) and \(s_{i}\) are extended-real-valued, nonnegative functions revealing the priors of the weight variable \(\mathbf{W}_{i}\) and the state variable \(\mathbf{V}_{i}\) (or the constraints on \(\mathbf{W}_{i}\) and \(\mathbf{V}_{i}\) ) for each \(i=1,\ldots N\), and define \(\mathbf{V}_{0}:=\mathbf{X}\). To solve the formulation in (4), the following alternative minimization problem was considered:

\[\min_{\mathcal{W},\mathcal{V},\mathcal{U},\mathcal{G}} \mathcal{L}(\mathcal{W},\mathcal{V},\mathcal{U},\mathcal{G}):= \mathcal{L}_{0}(\mathcal{W},\mathcal{V})+\frac{\gamma}{2}\sum_{i=1}^{N}\|\mathbf{V }_{i}-\sigma_{i}(\mathbf{U}_{i})\|_{F}^{2} \tag{5}\] \[+\frac{\rho}{2}\sum_{i=1}^{N}\|\mathbf{U}_{i}-\mathbf{W}_{i}\mathbf{V}_{i-1} \|_{F}^{2}+\frac{\tau}{2}\sum_{i=1}^{N}\|\mathbf{\mathcal{W}}_{i}-\text{TTD}(\mathbf{r} _{i})\|_{F}^{2}\]

where \(\gamma,\rho,\tau>0\) are hyperparameters for different regularization terms, \(\mathcal{U}:=\{\mathbf{U}_{i}\}_{i=1}^{N}\), and \(\mathcal{G}:=\{\mathcal{G}_{i}\}_{i=1}^{N}\) is the set of TT-cores \(\mathcal{G}_{i}\) from \(i\)-th layer. The NN training model (5) can be very general, where: (a) \(\ell\) can be the squared, logistic, hinge, cross-entropy or other commonly used loss functions; (b) \(\sigma_{i}\) can be ReLU, leaky ReLU, sigmoid, linear, polynomial, softplus or other commonly used activation functions; (c) \(\tau_{i}\) can be the squared \(\ell_{2}\) norm, the \(\ell_{1}\) norm, the elastic net, the indicator function of some nonempty closed convex set (such as the nonnegative closed half-space or a closed interval \([0,1]\)); (d) \(s_{i}\) can be the \(\ell_{1}\) norm, the indicator function of some convex set with simple projection. Particularly, if there is no regularizer or constraint on \(\mathbf{W}_{i}\) (or \(\mathbf{V}_{i}\)), then \(\tau_{i}\) (or \(s_{i}\)) can be zero. The network architectures considered in this paper exhibit generality to various types of NNs, including but not limited to the fully (or sparse) connected MLPs, convolutional neural networks (CNN) and residual neural networks (ResNets) [34].

**Remark 1** (Advantage of Formulation (5)).: _As mentioned before, an existing TT-format NN is either 1) trained from randomly initialized tensor cores; or 2) trained from a direct decomposition of a pre-trained model. For the first strategy, it does not utilize any information related to the high-accuracy uncompressed model; while other model compression methods, e.g. pruning and knowledge distillation, have shown that proper utilization of the pre-trained models is very critical for NN compression. For
the second strategy, though the knowledge of the pre-trained model is indeed utilized, because the pre-trained model generally lacks low TT-rank property, after direct low-rank tensor decomposition the approximation error is too significant to be properly recovered even using long-time re-training. Such inherent limitations of the existing training strategies, consequently, cause significant accuracy loss for the compressed TT-format NN models. To overcome these limitations, it is to maximally retain the knowledge contained in the uncompressed model, or in other words, minimize the approximation error after tensor decomposition with given target tensor ranks. In our formulation (5), \(\mathcal{L}_{0}(\mathcal{W},\mathcal{V})\) is the loss function of the uncompressed model while the regularization term \(\|\boldsymbol{\mathcal{W}}_{i}-\mathrm{TTD}(\boldsymbol{r}_{i})\|_{F}^{2}\) can encourage the uncompressed DNN models to gradually exhibit low tensor rank properties._

### Tensor BCD Algorithms

Note that (5) is a nonconvex optimization problem with multi-block variables. BCD is a Gauss-Seidel type method for a minimization problem with multi-block variables to update all the variables cyclically while fixing the remaining blocks at their last updated values [35]. A tensor BCD (tenBCD) algorithm is developed for solving (5). In this paper, proximal terms are added to some sub-problems arising from the tenBCD algorithm for two major reasons: (1) To practically stabilize the training process; (2) To yield the desired "sufficient descrease" property for theoretical justification. At each iteration \(k\), the tenBCD method with the backward order is considered for the updates of variables, i.e., the variables are updated from the output layer (layer \(N\)) to the input layer (layer 1). For each layer, the variables \(\{\boldsymbol{V}_{i},\boldsymbol{U}_{i},\boldsymbol{W}_{i},\mathcal{G}_{i}\}\) are updated cyclically for Problem (5). Since \(\sigma_{N}\equiv\mathrm{Id}\), the output layer is paid special attention. The tenBCD algorithms for (5) can be summarized in Algorithm 1.

#### 3.2.1 Optimization over \(\boldsymbol{V}_{i}\)

At iteration \(k\), \(\boldsymbol{V}_{N}\) can be updated through the following optimization problem

\[\boldsymbol{V}_{N}^{k}=\operatorname{argmin}_{\boldsymbol{V}_{N}}\left\{s_{N} (\boldsymbol{V}_{N})+\mathcal{R}_{n}(\boldsymbol{V}_{N};\boldsymbol{Y})+ \frac{\gamma}{2}\|\boldsymbol{V}_{N}-\boldsymbol{U}_{N}^{k-1}\|_{F}^{2}+\frac {\alpha}{2}\|\boldsymbol{V}_{N}-\boldsymbol{V}_{N}^{k-1}\|_{F}^{2}\right\}, \tag{6}\]

where \(s_{N}(\boldsymbol{V}_{N})+\mathcal{R}_{n}(\boldsymbol{V}_{N};\boldsymbol{Y})\) is regarded as a new proximal function \(\tilde{s}_{N}(\boldsymbol{V}_{N})\). When \(i<N\), \(\boldsymbol{V}_{i}\) can be updated through the following optimization problem

\[\boldsymbol{V}_{i}^{k}=\operatorname{argmin}_{\boldsymbol{V}_{i}}\left\{s_{i} (\boldsymbol{V}_{i})+\frac{\gamma}{2}\|\boldsymbol{V}_{i}-\sigma_{i}( \boldsymbol{U}_{i}^{k-1})\|_{F}^{2}+\frac{\rho}{2}\|\boldsymbol{U}_{i+1}^{k}- \boldsymbol{W}_{i+1}^{k}\boldsymbol{V}_{i}\|_{F}^{2}\right\}. \tag{7}\]

For subproblem (6), \(\frac{\alpha}{2}\|\boldsymbol{V}_{N}-\boldsymbol{V}_{N}^{k-1}\|_{F}^{2}\) is the proximal term, where \(\alpha>0\) is the positive coefficient.

The above two problems (6) and (7) are simple proximal updates [35, 36] (or just least squares problems), which usually have closed-form solutions to many commonly used NNs. For \(\boldsymbol{V}_{N}^{k}\)-update, \(s_{N}(\boldsymbol{V}_{N})+\mathcal{R}_{n}(\boldsymbol{V}_{N};\boldsymbol{Y})\) is regarded as a new proximal function \(\tilde{s}_{N}(\boldsymbol{V}_{N})\). Some typical examples leading to the closed-form solutions include: (a) \(s_{i}\) are \(0\) (i.e., no regularization), or the squared \(\ell_{2}\) norm, or the indicator function of a nonempty closed convex set with a simple projection like the nonnegative closed half-space and the closed interval \([0,1]\); (b) the loss function \(\ell\) is the squared loss or hinge loss.1

Footnote 1: The \(\boldsymbol{V}_{N}\)-update with hinge loss and other smooth losses is provided in Appendix A.1.

#### 3.2.2 Optimization over \(\boldsymbol{U}_{i}\)

At iteration \(k\), \(\boldsymbol{U}_{N}\) can be updated through the following optimization problem

\[\boldsymbol{U}_{N}^{k}=\operatorname{argmin}_{\boldsymbol{U}_{N}}\left\{\frac {\gamma}{2}\|\boldsymbol{V}_{N}^{k}-\boldsymbol{U}_{N}\|_{F}^{2}+\frac{\rho}{2} \|\boldsymbol{U}_{N}-\boldsymbol{W}_{N}^{k-1}\boldsymbol{V}_{N-1}^{k-1}\|_{F}^ {2}\right\} \tag{8}\]
\(\mathbf{U}_{i},i<N\) can be updated through the following optimization problem

\[\mathbf{U}_{i}^{k}=\operatorname*{argmin}_{\mathbf{U}_{i}}\left\{\frac{\gamma}{2}\|\mathbf{V} _{i}^{k}-\sigma_{i}(\mathbf{U}_{i})\|_{F}^{2}+\frac{\rho}{2}\|\mathbf{U}_{i}-\mathbf{W}_{i} ^{k-1}\mathbf{V}_{i-1}^{k-1}\|_{F}^{2}+\frac{\alpha}{2}\|\mathbf{U}_{i}-\mathbf{U}_{i}^{k- 1}\|_{F}^{2}\right\} \tag{9}\]

For subproblem (9), \(\frac{\alpha}{2}\|\mathbf{U}_{i}-\mathbf{U}_{i}^{k-1}\|_{F}^{2}\) is the proximal term. The subproblem (8) is a least-square optimization where the closed-form solution can be derived. The subproblem (9) is a nonlinear and nonsmooth optimization where \(\sigma_{i}\) is ReLU or leaky ReLU. Accordingly, the closed-form solution to solve the subproblem (9) is provided in Appendix A.2.

#### 3.2.3 Optimization over \(\mathbf{W}_{i}\)

At iteration \(k\), \(\mathbf{W}_{i},i=1,\ldots,N\) can be updated through the following optimization problem

\[\mathbf{W}_{i}^{k}=\operatorname*{argmin}_{\mathbf{W}_{i}}\left\{\tau_{i}(\mathbf{W}_{i}) +\frac{\rho}{2}\|\mathbf{U}_{i}^{k}-\mathbf{W}_{i}\mathbf{V}_{i-1}^{k-1}\|_{F}^{2}+\frac{ \tau}{2}\|\mathbf{\mathcal{W}}_{i}-\operatorname*{TTD}(\mathbf{r}_{i})\|_{F}^{2} \right\}, \tag{10}\]

The closed-form solution to solve the above optimization problem can be obtained when \(\tau_{i}\) is \(0\) (i.e., no regularization), or the squared \(\ell_{2}\) norm (i.e., weight decay), or the indicator function of a nonempty closed convex set with a simple projection like the nonnegative closed half-space and the closed interval \([0,1]\).

#### 3.2.4 Optimization over \(\mathcal{G}_{i}\)

At iteration \(k\), \(\mathcal{G}_{i},i=1,\ldots,N\) can be updated through the following optimization problem

\[\mathcal{G}_{i}^{k}=\operatorname*{argmin}_{\mathcal{G}_{i}}\left\{\frac{\tau }{2}\|\mathbf{\mathcal{W}}_{i}^{k}-\operatorname*{TTD}(\mathbf{r}_{i})\|_{F}^{2}+ \frac{\alpha}{2}\|\mathcal{G}_{i}-\mathcal{G}_{i}^{k-1}\|_{F}^{2}\right\} \tag{11}\]

where \(\frac{\alpha}{2}\|\mathcal{G}_{i}-\mathcal{G}_{i}^{k-1}\|_{F}^{2}\) is the proximal terms. This subproblem is implemented in TensorLy package [37].

```
0: Sample \(\mathbf{X}\in\mathbb{R}^{n_{0}\times n}\) and \(\mathbf{Y}\in\mathbb{R}^{n_{N}\times n}\), \(\gamma,\rho,\tau,\alpha>0\) Initialization:\(\{\mathbf{W}_{i}^{0},\mathbf{V}_{i}^{0},\mathbf{U}_{i}^{0},\mathcal{G}_{i}^{0}\}_{i=1}^{N}\), \(\mathbf{V}_{0}^{k}\equiv\mathbf{V}_{0}:=\mathbf{X}\)
1:for\(k=1,\ldots\)do
2: Get \(\mathbf{V}_{N}^{k}\) by solving (6)
3: Get \(\mathbf{U}_{N}^{k}\) by solving (8)
4: Get \(\mathbf{W}_{N}^{k}\) by solving (10)
5: Get \(\mathcal{G}_{N}^{k}\) by solving (11)
6:for\(i=N-1,\ldots,1\)do
7: Get \(\mathbf{V}_{i}^{k}\) by solving (7)
8: Get \(\mathbf{U}_{i}^{k}\) by solving (9)
9: Get \(\mathbf{W}_{i}^{k}\) by solving (10)
10: Get \(\mathcal{G}_{i}^{k}\) by solving (11)
11:endfor
12:endfor
13:\(\{\mathcal{G}_{i}\}_{i=1}^{N}\)
```

**Algorithm 1** tenBCD Algorithm

### Global Convergence Analysis of tenBCD

In this section, the global convergence of Algorithm 1 for Problem (5) is established. Firstly, let \(h:\mathbb{R}^{p}\rightarrow\mathbb{R}\cup\{+\infty\}\) be an extended-real-valued function, its graph is defined by \(Graph(h):=\{(\mathbf{x},y)\in\mathbb{R}^{p}\times\mathbb{R}:y=h(\mathbf{x})\}\), and its domain by \(\operatorname{dom}(h):=\{\mathbf{x}\in\mathbb{R}^{p}:h(\mathbf{x})<+\infty\}\). The subdifferential of a function is defined as follows.
**Definition 1** (Subdifferentials [38, 39]).: _Assume that \(f:\ \mathbb{R}^{p}\to(-\infty,+\infty)\) is a proper and lower semicontinuous function._

1. _The domain of_ \(f\) _is defined and denoted by_ \(\mathrm{dom}f\coloneqq\{\mathbf{x}\in\mathbb{R}^{p}:f(\mathbf{x})<+\infty\}\)__
2. _For a given_ \(\mathbf{x}\in\mathrm{dom}f\)_, the Frechet subdifferential of_ \(f\) _at_ \(\mathbf{x}\)_, written_ \(\hat{\partial}f(\mathbf{x})\)_, is the set of all vectors_ \(\mathbf{u}\in\mathbb{R}^{p}\) _that satisfy_ \[\lim_{\mathbf{y}\neq\mathbf{x}}\inf_{\mathbf{y}\to\mathbf{x}}\frac{f(\mathbf{y})-f(\mathbf{x})-\langle \mathbf{u},\mathbf{y}-\mathbf{x}\rangle}{\|\mathbf{y}-\mathbf{x}\|}\geq 0.\]
3. _The limiting-subdifferential, or simply the subdifferential, of_ \(f\) _at_ \(\mathbf{x}\)_, written_ \(\partial f(\mathbf{x})\) _is defined through the following closure process_ \[\partial f(\mathbf{x}):=\{\mathbf{u}\in\mathbb{R}^{p}:\exists\mathbf{x}^{k}\to\mathbf{x},f( \mathbf{x}^{k})\to f(\mathbf{x})\ \text{and}\ \mathbf{u}^{k}\in\hat{\partial}f(\mathbf{x}^{k})\to\mathbf{u}\ \text{as}\ k\to\infty\}.\]

Now, our first main lemma about the sufficient decrease property of the iterative sequence \(\{\mathcal{P}^{k}:=(\{\mathbf{W}_{i}^{k}\}_{i=1}^{N},\{\mathbf{V}_{i}^{k}\}_{i=1}^{N}, \{\mathbf{U}_{i}^{k}\}_{i=1}^{N}),\{\mathcal{G}_{i}^{k}\}_{i=1}^{N}\}_{k\in \mathbb{N}}\) from Algorithm 1 is ready to be introduced.

**Lemma 2** (Sufficient Decrease Property).: _Given that \(\alpha,\gamma,\rho,\tau>0\), \(\left\{\mathcal{P}^{k}\right\}_{k\in\mathbb{N}}\) is the sequence generated by the tenBCD algorithm 1, then the sequence satisfies_

\[\mathcal{L}(\mathcal{P}^{k})\leq\mathcal{L}(\mathcal{P}^{k-1})-\lambda\| \mathcal{P}^{k}-\mathcal{P}^{k-1}\|_{F}^{2}. \tag{12}\]

_For the case that \(\mathbf{V}_{N}\) is updated via the proximal strategy, \(\lambda:=\min\left\{\frac{\alpha}{2},\frac{\gamma+\rho}{2},\frac{\tau}{2}\right\}\). For the case that \(\mathbf{V}_{N}\) is update via the prox-linear strategy, \(\lambda:=\min\left\{\frac{\alpha}{2},\frac{\gamma+\rho}{2},\frac{\tau}{2}, \alpha+\frac{\gamma-LR}{2}\right\}\), where \(\nabla\mathcal{R}_{n}\) is Lipschitz continuous with a Lipschitz constant \(L_{R}\) and \(\alpha>\max\{0,\frac{L_{R}-\gamma}{2}\}\)._

Proof.: The inequality (12) can be developed by considering the descent quantity along the update of each block variable, i.e., \(\{\mathbf{V}_{i}\}_{i=1}^{N}\), \(\{\mathbf{U}_{i}\}_{i=1}^{N}\), \(\{\mathbf{W}_{i}\}_{i=1}^{N}\), and \(\{\mathcal{G}_{i}\}_{i=1}^{N}\). To begin with, the following notations are introduced. Specifically, \(\mathbf{W}_{<i}:=\left(\mathbf{W}_{1},\mathbf{W}_{2},\ldots,\mathbf{W}_{i-1}\right),\mathbf{W}_{>i }:=\left(\mathbf{W}_{i+1},\mathbf{W}_{i+1},\ldots,\mathbf{W}_{N}\right)\), and \(\mathbf{V}_{<i},\mathbf{V}_{>i},\mathbf{U}_{<i},\mathbf{U}_{>i},\mathcal{G}_{<i},\mathcal{G}_ {>i}\) are defined similarly. We will consider each case separately.

**Optimization over \(\mathbf{V}_{i}\)**

\(\mathbf{V}_{N}^{k}\)**-block:** at iteration \(k\), there are two ways to update the variable: (1) proximal update with closed-form solution: the following inequality can be derived

\[\mathcal{L}\left(\{\mathbf{W}_{i}^{k-1}\}_{i=1}^{N},\mathbf{V}_{i<N}^{k-1 },\mathbf{V}_{N}^{k},\{\mathbf{U}_{i}^{k-1}\}_{i=1}^{N},\{\mathcal{G}_{i}^{k-1}\}_{i=1 }^{N}\right) \tag{13}\] \[\leq \mathcal{L}\left(\{\mathbf{W}_{i}^{k-1}\}_{i=1}^{N},\mathbf{V}_{i<N}^{k-1 },\mathbf{V}_{N}^{k-1},\{\mathbf{U}_{i}^{k-1}\}_{i=1}^{N},\{\mathcal{G}_{i}^{k-1}\}_{i= 1}^{N}\right)-\frac{\alpha}{2}\|\mathbf{V}_{N}^{k}-\mathbf{V}_{N}^{k-1}\|_{F}^{2}.\]

The above inequality (13) is due to the fact that \(\mathbf{V}_{N}^{k}\) is the optimal solution for subproblem (6). (2) proximal-linear case: let \(h^{k}(\mathbf{V}_{N}):=s_{N}(\mathbf{V}_{N})+\mathcal{R}_{n}(\mathbf{V}_{N};\mathbf{Y})+\frac{ \gamma}{2}\|\mathbf{V}_{N}-\mathbf{U}_{N}^{k-1}\|_{F}^{2}\) and \(\bar{h}^{k}(\mathbf{V}_{N}):=s_{N}(\mathbf{V}_{N})+\mathcal{R}_{n}(\mathbf{V}_{N}^{k-1}; \mathbf{Y})+\langle\nabla\mathcal{R}_{n}(\mathbf{V}_{N}^{k-1};\mathbf{Y}),\mathbf{V}_{N}-\mathbf{V }_{N}^{k-1}\rangle+\frac{\alpha}{2}\|\mathbf{V}_{N}-\mathbf{V}_{N}^{k-1}\|_{F}^{2}+\frac {2}{2}\|\mathbf{V}_{N}-\mathbf{U}_{N}^{k-1}\|_{F}^{2}\). By the optimality of \(\mathbf{V}_{N}^{k}\) and the strong convexity2 of \(\bar{h}^{k}(\mathbf{V}_{N})\) with modulus at least \(\alpha+\gamma\), the following holds

Footnote 2: The function \(h\) is called a strongly convex function with parameter \(\gamma>0\) if \(h(u)\geq h(v)+\langle\nabla h(v),u-v\rangle+\frac{\gamma}{2}\|u-v\|^{2}\).

\[\bar{h}^{k}(\mathbf{V}_{N}^{k})\leq\bar{h}^{k}(\mathbf{V}_{N}^{k-1})-\frac{\alpha+\gamma} {2}\|\mathbf{V}_{N}^{k}-\mathbf{V}_{N}^{k-1}\|_{F}^{2}, \tag{14}\]

which implies

\[h^{k}(\mathbf{V}_{N}^{k})\leq h^{k}(\mathbf{V}_{N}^{k-1})+\mathcal{R}_{n}(\mathbf{V}_{N}^{k}; \mathbf{Y})-\mathcal{R}_{n}(\mathbf{V}_{N}^{k-1};\mathbf{Y})-\langle\nabla\mathcal{R}_{n}(\mathbf{V }_{N}^{k-1};\mathbf{Y}),\mathbf{V}_{N}^{k}-\mathbf{V}_{N}^{k-1}\rangle\]
\[\begin{split}&-(\alpha+\frac{\gamma}{2})\|\mathbf{V}_{N}^{k}-\mathbf{V}_{N}^{k- 1}\|_{F}^{2}\\ &\leq h^{k}(\mathbf{V}_{N}^{k-1})-(\alpha+\frac{\gamma-L_{R}}{2})\|\bm {V}_{N}^{k}-\mathbf{V}_{N}^{k-1}\|_{F}^{2},\end{split}\] (15a) where inequality ( 15a ) is due to the inequality ( 14 ), the relationship between \[h^{k}(\mathbf{V}_{N}^{k-1})\] and \[\bar{h}^{k}(\mathbf{V}_{N}^{k-1})\], and the relationship between \[h^{k}(\mathbf{V}_{N}^{k})\] and \[\bar{h}^{k}(\mathbf{V}_{N}^{k})\]. The inequality ( 15b ) holds for the \[L_{R}\] -Lipschitz continuity of \[\nabla\mathcal{R}_{n}\], i.e., the following inequality by [40] \[\mathcal{R}_{n}(\mathbf{V}_{N}^{k};\mathbf{Y})\leq\mathcal{R}_{n}(\mathbf{V}_{N}^{k-1}; \mathbf{Y})+\langle\nabla\mathcal{R}_{n}(\mathbf{V}_{N}^{k-1};\mathbf{Y}),\mathbf{V}_{N}^{k}- \mathbf{V}_{N}^{k-1}\rangle+\frac{L_{R}}{2}\|\mathbf{V}_{N}^{k}-\mathbf{V}_{N}^{k-1}\|_{F }^{2}.\] According to the relationship between \[h^{k}(\mathbf{V}_{N})\] and \[\mathcal{L}\left(\{\mathbf{W}_{i}^{k-1}\}_{i=1}^{N},\mathbf{V}_{i<N}^{k-1},\mathbf{V}_{N}, \{\mathbf{U}_{i}^{k-1}\}_{i=1}^{N},\{\mathcal{G}_{i}^{k-1}\}_{i=1}^{N}\right)\], and the inequality ( 15 ), \[\begin{split}&\mathcal{L}\left(\{\mathbf{W}_{i}^{k-1}\}_{i=1}^{N}, \mathbf{V}_{i<N}^{k-1},\mathbf{V}_{N}^{k},\{\mathbf{U}_{i}^{k-1}\}_{i=1}^{N},\{\mathcal{G} _{i}^{k-1}\}_{i=1}^{N}\right)\\ \leq&\mathcal{L}\left(\{\mathbf{W}_{i}^{k-1}\}_{i=1}^{N}, \mathbf{V}_{i<N}^{k-1},\mathbf{V}_{N}^{k-1},\{\mathbf{U}_{i}^{k-1}\}_{i=1}^{N},\{\mathcal{ G}_{i}^{k-1}\}_{i=1}^{N}\right)-(\alpha+\frac{\gamma-L_{R}}{2})\|\mathbf{V}_{N}^{k}-\mathbf{V}_{ N}^{k-1}\|_{F}^{2}.\end{split} \tag{16}\]

\(\mathbf{V}_{i}^{k}\)**-block (\(i<N\)):**\(\mathbf{V}_{i}^{k}\) is updated according to the following

\[\mathbf{V}_{i}^{k}\leftarrow\operatorname*{argmin}_{\mathbf{V}_{i}}\left\{s_{i}(\mathbf{V} _{i})+\frac{\gamma}{2}\|\mathbf{V}_{i}-\sigma_{i}(\mathbf{U}_{i}^{k-1})\|_{F}^{2}+\frac {\rho}{2}\|\mathbf{U}_{i+1}^{k}-\mathbf{W}_{i+1}^{k}\mathbf{V}_{i}\|_{F}^{2}\right\}.\]

Let \(h^{k}(\mathbf{V}_{i})=s_{i}(\mathbf{V}_{i})+\frac{\gamma}{2}\|\mathbf{V}_{i}-\sigma_{i}(\bm {U}_{i}^{k-1})\|_{F}^{2}+\frac{\rho}{2}\|\mathbf{U}_{i+1}^{k}-\mathbf{W}_{i+1}^{k}\mathbf{V }_{i}\|_{F}^{2}\). By the convexity of \(s_{i}\), the function \(h^{k}(\mathbf{V}_{i})\) is a strongly convex function with modulus no less than \(\gamma\). By the optimality of \(\mathbf{V}_{i}^{k}\), the following holds

\[h^{k}(\mathbf{V}_{i}^{k})\leq h^{k}(\mathbf{V}_{i}^{k-1})-\frac{\gamma}{2}\|\mathbf{V}_{i} ^{k}-\mathbf{V}_{i}^{k-1}\|_{F}^{2}. \tag{17}\]

Based on the inequality (17), it yields for

\[\begin{split}&\mathcal{L}(\mathbf{W}_{\leq i}^{k-1},\mathbf{W}_{>i}^{k},\mathbf{V}_{< i}^{k-1},\mathbf{V}_{i}^{k},\mathbf{V}_{>i}^{k},\mathbf{U}_{\leq i}^{k-1},\mathbf{U}_{>i}^{k}, \mathcal{G}_{\leq i}^{k-1},\mathcal{G}_{>i}^{k})\\ \leq&\mathcal{L}(\mathbf{W}_{\leq i}^{k-1},\mathbf{W}_{>i}^{k}, \mathbf{V}_{<i}^{k-1},\mathbf{V}_{i}^{k-1},\mathbf{V}_{>i}^{k},\mathbf{U}_{\leq i}^{k-1}, \mathcal{U}_{>i}^{k},\mathcal{G}_{\leq i}^{k-1},\mathcal{G}_{>i}^{k})-\frac{ \gamma}{2}\|\mathbf{V}_{i}^{k}-\mathbf{V}_{i}^{k-1}\|_{F}^{2}\end{split} \tag{18}\]

for \(i=1,\ldots,N-1\), where

\[\begin{split} h^{k}(\mathbf{V}_{i}^{k})-h^{k}(\mathbf{V}_{i}^{k-1})& =\mathcal{L}(\mathbf{W}_{\leq i}^{k-1},\mathbf{W}_{>i}^{k},\mathbf{V}_{<i}^{ k-1},\mathbf{V}_{i}^{k},\mathbf{V}_{>i}^{k},\mathbf{U}_{\leq i}^{k-1},\mathbf{U}_{>i}^{k}, \mathcal{G}_{\leq i}^{k-1},\mathcal{G}_{>i}^{k})\\ &-\mathcal{L}(\mathbf{W}_{\leq i}^{k-1},\mathbf{W}_{>i}^{k},\mathbf{V}_{<i}^{k-1 },\mathbf{V}_{i}^{k-1},\mathbf{V}_{>i}^{k},\mathbf{U}_{\leq i}^{k-1},\mathbf{U}_{>i}^{k}, \mathcal{G}_{\leq i}^{k-1},\mathcal{G}_{>i}^{k}).\end{split}\]

**Optimization over \(\mathbf{U}_{i}\)**

\(\mathbf{U}_{N}^{k}\)**-block**: similar to the inequality (18), the descent quantity is established as follows

\[\begin{split}&\mathcal{L}(\mathbf{W}_{\leq N}^{k-1},\mathbf{V}_{<N}^{k-1}, \mathbf{V}_{N}^{k},\mathbf{U}_{<N}^{k-1},\mathbf{U}_{N}^{k},\mathcal{G}_{\leq N}^{k-1})\\ \leq&\mathcal{L}(\mathbf{W}_{\leq N}^{k-1},\mathbf{V}_{<N}^{k-1},\mathbf{V}_{N}^{k}, \mathbf{U}_{<N}^{k-1},\mathbf{U}_{N}^{k-1},\mathcal{G}_{\leq N}^{k-1})-\frac{\gamma+ \rho}{2}\|\mathbf{U}_{N}^{k}-\mathbf{U}_{N}^{k-1}\|_{F}^{2},\end{split} \tag{19}\]

where the above inequality is because the objective function in subproblem (8) is a strongly convex function with modulus at least \(\gamma+\rho\).
\(\mathbf{U}_{i}^{k}\)**-block (\(i<N\))**: the following can be obtained

\[\begin{split}&\mathcal{L}(\mathbf{W}_{\leq i}^{k-1},\mathbf{W}_{>i}^{k},\mathbf{V} _{<i}^{k-1},\mathbf{V}_{\geq i}^{k},\mathbf{U}_{<i}^{k-1},\mathbf{U}_{i}^{k},\mathbf{U}_{>i}^{k},\mathcal{G}_{\leq i}^{k-1},\mathcal{G}_{>i}^{k})\\ \leq&\mathcal{L}(\mathbf{W}_{\leq i}^{k-1},\mathbf{W}_{>i}^{k },\mathbf{V}_{<i}^{k-1},\mathbf{V}_{\geq i}^{k},\mathbf{U}_{<i}^{k-1},\mathbf{U}_{i}^{k-1},\bm {U}_{>i}^{k},\mathcal{G}_{\leq i}^{k-1},\mathcal{G}_{>i}^{k})-\frac{\alpha}{2} \|\mathbf{U}_{i}^{k}-\mathbf{U}_{i}^{k-1}\|_{F}^{2}\end{split} \tag{20}\]

for \(i=1,\ldots,N-1\) since \(\mathbf{U}_{i}^{k}\) is the optimal solution for subproblem (9).

**Optimization over \(\mathbf{W}_{i}\)**

\(\mathbf{W}_{i}^{k}\)**-block (\(i\leq N\))**: \(\mathbf{W}_{i}^{k}\) is updated according to the following

\[\mathbf{W}_{i}^{k}\rightarrow\operatorname*{argmin}_{\mathbf{W}_{i}}\left\{r_{i}(\bm {W}_{i})+\frac{\rho}{2}\|\mathbf{U}_{i}^{k}-\mathbf{W}_{i}\mathbf{V}_{i-1}^{k-1}\|_{F}^{2} +\frac{\tau}{2}\|\mathbf{\mathcal{W}}_{i}-\operatorname{TTD}(\mathbf{r}_{i})\|_{F}^{ 2}\right\},\]

where \(h^{k}(\mathbf{W}_{i})=r_{i}(\mathbf{W}_{i})+\frac{\rho}{2}\|\mathbf{U}_{i}^{k}-\mathbf{W}_{i} \mathbf{V}_{i-1}^{k-1}\|_{F}^{2}+\frac{\tau}{2}\|\mathbf{\mathcal{W}}_{i}-\operatorname {TTD}(\mathbf{r}_{i})\|_{F}^{2}\) is a strongly convex function with modulus at least \(\tau\). Accordingly, the following holds

\[\begin{split}&\mathcal{L}(\mathbf{W}_{<i}^{k-1},\mathbf{W}_{i}^{k},\mathbf{W}_{> i}^{k},\mathbf{V}_{<i}^{k-1},\mathbf{V}_{\geq i}^{k},\mathbf{U}_{<i}^{k-1},\mathbf{U}_{\geq i }^{k},\mathcal{G}_{\leq i}^{k-1},\mathcal{G}_{>i}^{k})\\ \leq&\mathcal{L}(\mathbf{W}_{<i}^{k-1},\mathbf{W}_{i}^{k-1}, \mathbf{W}_{>i}^{k},\mathbf{V}_{<i}^{k-1},\mathbf{V}_{\geq i}^{k},\mathbf{U}_{<i}^{k-1},\mathbf{U} _{\geq i}^{k},\mathcal{G}_{\leq i}^{k-1},\mathcal{G}_{>i}^{k})-\frac{\tau}{2} \|\mathbf{W}_{i}^{k}-\mathbf{W}_{i}^{k-1}\|_{F}^{2},\end{split} \tag{21}\]

which is due to the relationship between \(h^{k}(\mathbf{W}_{i})\) and \(\mathcal{L}(\mathbf{W}_{<i}^{k-1},\mathbf{W}_{i},\mathbf{W}_{>i}^{k},\mathbf{V}_{<i}^{k-1},\bm {V}_{\geq i}^{k},\mathbf{U}_{<i}^{k-1},\mathbf{U}_{\geq i}^{k},\mathcal{G}_{\leq i}^{k -1},\mathcal{G}_{>i}^{k})\).

**Optimization over \(\mathcal{G}_{i}\)**

\(\mathcal{G}_{i}\)**-block (\(i\leq N\))**: the descent quantity for \(\mathcal{G}_{i}\) can be derived as follows

\[\begin{split}&\mathcal{L}(\mathbf{W}_{<i}^{k-1},\mathbf{W}_{\geq i}^{k}, \mathbf{V}_{<i}^{k-1},\mathbf{V}_{\geq i}^{k},\mathbf{U}_{<i}^{k-1},\mathbf{U}_{\geq i}^{k}, \mathcal{G}_{\leq i}^{k-1},\mathcal{G}_{i}^{k},\mathcal{G}_{>i}^{k})\\ \leq&\mathcal{L}(\mathbf{W}_{<i}^{k-1},\mathbf{W}_{\geq i}^{k },\mathbf{V}_{<i}^{k-1},\mathbf{V}_{\geq i}^{k},\mathbf{U}_{<i}^{k-1},\mathbf{U}_{\geq i}^{k}, \mathcal{G}_{<i}^{k-1},\mathcal{G}_{i}^{k-1},\mathcal{G}_{>i}^{k})-\frac{ \alpha}{2}\|\mathcal{G}_{i}^{k}-\mathcal{G}_{i}^{k-1}\|_{F}^{2},\end{split} \tag{22}\]

where the above inequality (22) is due to the fact that \(\mathcal{G}_{i}^{k}\) is the optimal solution for subproblem (11).

By summing up inequalities (13) (or (16)), (19), (20), (21), and (22), it yields the

\[\mathcal{L}(\mathcal{P}^{k})\leq\mathcal{L}(\mathcal{P}^{k-1})-\lambda\| \mathcal{P}^{k}-\mathcal{P}^{k-1}\|_{F}^{2},\]

where \(\lambda:=\min\left\{\frac{\alpha}{2},\frac{\gamma+\rho}{2},\frac{\tau}{2}\right\}\) (or \(\lambda:=\min\left\{\frac{\alpha}{2},\frac{\gamma+\rho}{2},\frac{\tau}{2}, \alpha+\frac{\gamma-L_{R}}{2}\right\}\)). 

From Lemma 2, the Lagrangian sequence \(\left\{\mathcal{L}(\mathcal{P}^{k})\right\}_{k\in\mathbb{N}}\) is monotonically decreasing, and the descent quantity of each iterate can be lower bounded by the discrepancy between the current iterate and its previous iterate. This lemma is crucial for the global convergence of a nonconvex algorithm. It tells at least the following four important items: (i) \(\{\mathcal{L}(\mathcal{P}^{k})\}_{k\in\mathbb{N}}\) is convergent if \(\mathcal{L}\) is lower bounded; (ii) \(\{\mathcal{P}^{k}\}_{k\in\mathbb{N}}\) itself is bounded if \(\mathcal{L}\) is coercive and \(\mathcal{P}^{0}\) is finite; (iii) \(\{\mathcal{P}^{k}\}_{k\in\mathbb{N}}\) is square summable, i.e., \(\sum_{k=1}^{\infty}\|\mathcal{P}^{k}-\mathcal{P}^{k-1}\|_{F}^{2}<\infty\), implying its asymptotic regularity, i.e., \(\|\mathcal{P}^{k}-\mathcal{P}^{k-1}\|_{F}\to 0\) as \(k\rightarrow\infty\); and (iv) \(\frac{1}{K}\sum_{k=1}^{K}\|\mathcal{P}^{k}-\mathcal{P}^{k-1}\|_{F}^{2}\to 0\) at a rate of \(\mathcal{O}(1/K)\). Leveraging Lemma 2, we can establish the global convergence (i.e., the whole sequence convergence) of tenBCD algorithm 1 in NN training settings. In contrast, [41] only establish the subsequence convergence of SGD in NN training settings. Such a gap between the subsequence convergence of SGD in [41] and the whole sequence convergence of tenBCD algorithm 1 in this paper exists mainly because SGD can only achieve the descent property but not the sufficient descent property.
It can be noted from Lemma 2 that neither multiconvexity and differentiability nor Lipschitz differentiability assumptions are imposed on the NN training models to yield this lemma, as required in the literature [42, 35, 43, 36]. Instead, we mainly exploit the proximal strategy for all nonstrongly convex subproblems in Algorithm 1 to establish this lemma.

Our second main lemma is about the subgradient lower bound.

**Lemma 3** (Subgradient Lower Bound).: _Under the same conditions of Lemma 2, let \(\mathcal{B}\) be an upper bound of \(\mathcal{P}^{k-1}\) and \(\mathcal{P}^{k}\) for any positive integer \(k,L_{\mathcal{B}}\) be a uniform Lipschitz constant of \(\sigma_{i}\) on the bounded set \(\{\mathcal{P}:\|\mathcal{P}\|_{F}\leq\mathcal{B}\}\), and_

\[\delta:=\max\{\gamma,\alpha+\rho\mathcal{B},\alpha+\gamma L_{\mathcal{B}},2 \rho\mathcal{B}+2\rho\mathcal{B}^{2},\alpha+\sqrt{N}\tau\mathcal{B}^{N-1}\}\]

_(or, for the prox-linear case, \(\delta:=\max\{\gamma,L_{R}+\alpha+\rho\mathcal{B},\alpha+\gamma L_{\mathcal{B} },2\rho\mathcal{B}+2\rho\mathcal{B}^{2},\alpha+\sqrt{N}\tau\mathcal{B}^{N-1}\}\)), then for any positive integer \(k\), there holds,_

\[\begin{split}\operatorname{dist}(\mathbf{0},\partial\mathcal{L} (\mathcal{P}^{k}))&\leq\delta\sum_{i=1}^{N}\left[\|\mathbf{W}_{i}^{k -}-\mathbf{W}_{i}^{k-1}\|_{F}+\|\mathbf{V}_{i}^{k}-\mathbf{V}_{i}^{k-1}\|_{F}+\|\mathbf{U}_{i}^ {k}-\mathbf{U}_{i}^{k-1}\|_{F}+\|\mathcal{G}_{i}^{k}-\mathcal{G}_{i}^{k-1}\|_{F} \right]\\ &\leq\bar{\delta}\|\mathcal{P}^{k}-\mathcal{P}^{k-1}\|_{F}\end{split} \tag{23}\]

_where \(\bar{\delta}:=\delta\sqrt{4N}\), \(\operatorname{dist}(\mathbf{0},\mathcal{S}):=\inf_{\mathbf{s}\in\mathcal{S}}\|\bm {s}\|_{F}\) for a set \(\mathcal{S}\), and_

\[\partial\mathcal{L}(\mathcal{P}^{k}):=(\{\partial_{\mathbf{W}_{i}}\mathcal{L}\}_ {i=1}^{N},\{\partial_{\mathbf{V}_{i}}\mathcal{L}\}_{i=1}^{N},\{\partial_{\mathbf{U}_{i }}\mathcal{L}\}_{i=1}^{N},\{\partial_{\mathcal{G}_{i}}\mathcal{L}\}_{i=1}^{N} )(\mathcal{P}^{k}).\]

Proof.: The inequality (23) is established via bounding each term of \(\partial\mathcal{L}(\mathcal{P}^{k})\). Specifically, the following holds

\[\mathbf{0} \in\partial s_{N}(\mathbf{V}_{N}^{k})+\partial\mathcal{R}_{n}(\mathbf{V}_ {N}^{k};\mathbf{Y})+\gamma(\mathbf{V}_{N}^{k}-\mathbf{U}_{N}^{k-1})+\alpha(\mathbf{V}_{N}^{k}- \mathbf{V}_{N}^{k-1}), \tag{24a}\] \[\mathbf{0} \in\partial s_{N}(\mathbf{V}_{N}^{k})+\nabla\mathcal{R}_{n}(\mathbf{V}_ {N}^{k-1};\mathbf{Y})+\gamma(\mathbf{V}_{N}^{k}-\mathbf{U}_{N}^{k-1})+\alpha(\mathbf{V}_{N}^{k }-\mathbf{V}_{N}^{k-1}),\ \text{(proximal-linear)}\] (24b) \[\mathbf{0} =\gamma(\mathbf{U}_{N}^{k}-\mathbf{V}_{N}^{k})+\rho(\mathbf{U}_{N}^{k}-\mathbf{W} _{N}^{k-1}\mathbf{V}_{N-1}^{k-1}),\] (24c) \[\mathbf{0} \in\partial\tau_{N}(\mathbf{W}_{N}^{k})+\rho(\mathbf{W}_{N}^{k}\mathbf{V}_{N -1}^{k-1}-\mathbf{U}_{N}^{k})\mathbf{V}_{N-1}^{k-1}+\tau\left(\mathbf{W}_{N}^{k}-\mathrm{ TTD}^{k-1}(\mathbf{r}_{N})\right),\] (24d) \[\mathbf{0} \in\partial\left(\frac{\tau}{2}\|\mathbf{\mathcal{W}}_{N}^{k}- \mathrm{TTD}^{k}(\mathbf{r}_{N})\|_{F}^{2}\right)+\alpha(\mathcal{G}_{N}^{k}- \mathcal{G}_{N}^{k-1}), \tag{24e}\]

where (24a), (24b), (24c), (24d), and (24e) are due to the optimality conditions of all updates in (6), (33), (8), (10), and (11), respectively.

For \(i=N-1,\ldots,1\), the following holds

\[\mathbf{0} \in\partial s_{i}(\mathbf{V}_{i}^{k})+\gamma(\mathbf{V}_{i}^{k}-\sigma_{i }(\mathbf{U}_{i}^{k-1}))+\rho{\mathbf{W}_{i+1}^{k}}^{\top}(\mathbf{W}_{i+1}^{k}\mathbf{V}_{i}^ {k}-\mathbf{U}_{i+1}^{k}), \tag{25a}\] \[\mathbf{0} \in\gamma\left[(\sigma_{i}(\mathbf{U}_{i}^{k})-\mathbf{V}_{i}^{k})\odot \partial\sigma_{i}(\mathbf{U}_{i}^{k})\right]+\rho(\mathbf{U}_{i}^{k}-\mathbf{W}_{i}^{k-1} \mathbf{V}_{i-1}^{k-1})+\alpha(\mathbf{U}_{i}^{k}-\mathbf{U}_{i}^{k-1}),\] (25b) \[\mathbf{0} \in\partial\tau_{i}(\mathbf{W}_{i}^{k})+\rho(\mathbf{W}_{i}^{k}\mathbf{V}_ {i-1}^{k-1}-\mathbf{U}_{i}^{k})\mathbf{V}_{i-1}^{k-1}\top+\tau\left(\mathbf{W}_{i}^{k}- \mathrm{TTD}^{k-1}(\mathbf{r}_{i})\right),\] (25c) \[\mathbf{0} \in\partial\left(\frac{\tau}{2}\|\mathbf{\mathcal{W}}_{i}^{k}- \mathrm{TTD}^{k}(\mathbf{r}_{i})\|_{F}^{2}\right)+\alpha(\mathcal{G}_{i}^{k}- \mathcal{G}_{i}^{k-1}), \tag{25d}\]

where (25a), (25b), (25c), and (25d) are due to the optimality conditions of all updates in (7), (9), (10), and (11), respectively. \(\mathbf{V}_{0}^{k}\equiv\mathbf{V}_{0}=\mathbf{X}\) for all \(k\), and \(\odot\) is the Hadamard product. Through the
above relationship (24), we have

\[-\alpha(\mathbf{V}_{N}^{k}-\mathbf{V}_{N}^{k-1})-\gamma(\mathbf{U}_{N}^{k}-\mathbf{U} _{N}^{k-1})\in\partial s_{N}(\mathbf{V}_{N}^{k})+\partial\mathcal{R}_{n}(\mathbf{V}_{N}^ {k};\mathbf{Y})+\gamma(\mathbf{V}_{N}^{k}-\mathbf{U}_{N}^{k})=\partial_{\mathbf{V}_{N}}\mathcal{ L}(\mathcal{P}^{k}),\] \[\Big{(}\nabla\mathcal{R}_{n}(\mathbf{V}_{N}^{k};\mathbf{Y})-\nabla \mathcal{R}_{n}(\mathbf{V}_{N}^{k-1};\mathbf{Y})\Big{)}-\alpha(\mathbf{V}_{N}^{k}-\mathbf{V}_ {N}^{k-1})-\gamma(\mathbf{U}_{N}^{k}-\mathbf{U}_{N}^{k-1})\in\partial_{\mathbf{V}_{N}} \mathcal{L}(\mathcal{P}^{k}),\ \text{(proximal-linear)}\] \[-\rho(\mathbf{W}_{N}^{k}-\mathbf{W}_{N}^{k-1})\mathbf{V}_{N-1}^{k}-\rho\mathbf{W }_{N}^{k-1}(\mathbf{V}_{N-1}^{k}-\mathbf{V}_{N-1}^{k-1})=\gamma(\mathbf{U}_{N}^{k}-\mathbf{V}_ {N}^{k})+\rho(\mathbf{U}_{N}^{k}-\mathbf{W}_{N}^{k}\mathbf{V}_{N-1}^{k})=\partial_{\mathbf{U}_ {N}}\mathcal{L}(\mathcal{P}^{k}),\] \[\rho\mathbf{W}_{N}^{k}[\mathbf{V}_{N-1}^{k}(\mathbf{V}_{N-1}^{k}-\mathbf{V}_{N-1}^ {k-1})^{\top}+(\mathbf{V}_{N-1}^{k}-\mathbf{V}_{N-1}^{k-1})\mathbf{V}_{N-1}^{k-1}\top]- \rho\mathbf{U}_{N}^{k}(\mathbf{V}_{N}^{k}-\mathbf{V}_{N}^{k-1})^{\top}+\tau\left(\mathrm{ TTD}^{k}(\mathbf{r}_{N})-\mathrm{TTD}^{k-1}(\mathbf{r}_{N})\right)\] \[\in\partial r_{N}(\mathbf{W}_{N}^{k})+\rho(\mathbf{W}_{N}^{k}\mathbf{V}_{N-1}^ {k}-\mathbf{U}_{N}^{k})\mathbf{V}_{N-1}^{k}\top+\tau\left(\mathbf{W}_{N}^{k}-\mathrm{TTD}^ {k}(\mathbf{r}_{i})\right)=\partial_{\mathbf{W}_{N}}\mathcal{L}(\mathcal{P}^{k}),\] \[-\alpha(\mathcal{G}_{N}^{k}-\mathcal{G}_{N}^{k-1})\in\partial_{ \mathcal{G}_{N}}\mathcal{L}(\mathcal{P}^{k}). \tag{26}\]

For \(i=N-1,\ldots,1\), the relationship (25) implies

\[-\gamma(\sigma_{i}(\mathbf{U}_{i}^{k})-\sigma_{i}(\mathbf{U}_{i}^{k-1})) \in\partial s_{i}(\mathbf{V}_{i}^{k})+\rho(\mathbf{V}_{i}^{k}-\sigma_{i}(\mathbf{U}_{i}^ {k}))+\gamma\mathbf{W}_{i+1}^{k}\top(\mathbf{W}_{i+1}^{k}\mathbf{V}_{i}^{k}-\mathbf{U}_{i+1}^{ k})=\partial_{\mathbf{V}_{i}}\mathcal{L}(\mathcal{P}^{k}),\] \[-\rho\mathbf{W}_{i}^{k-1}(\mathbf{V}_{i-1}^{k}-\mathbf{V}_{i-1}^{k-1})-\rho( \mathbf{W}_{i}^{k}-\mathbf{W}_{i}^{k-1})\mathbf{V}_{i-1}^{k}-\alpha(\mathbf{U}_{i}^{k}-\mathbf{U}_ {i}^{k-1})\] \[\in\gamma\left[(\sigma_{i}(\mathbf{U}_{i}^{k})-\mathbf{V}_{i}^{k})\odot \partial\sigma_{i}(\mathbf{U}_{i}^{k})\right]+\rho(\mathbf{U}_{i}^{k}-\mathbf{W}_{i}^{k}\bm {V}_{i-1}^{k})=\partial_{\mathbf{U}_{i}}\mathcal{L}(\mathcal{P}^{k}),\] \[\rho\mathbf{W}_{i}^{k}[\mathbf{V}_{i-1}^{k}(\mathbf{V}_{i-1}^{k}-\mathbf{V}_{i-1}^ {k-1})^{\top}+(\mathbf{V}_{i-1}^{k}-\mathbf{V}_{i-1}^{k-1})\mathbf{V}_{i-1}^{k-1}]-\rho \mathbf{U}_{i}^{k}(\mathbf{V}_{i-1}^{k}-\mathbf{V}_{i-1}^{k-1})^{\top}+\tau\left(\mathrm{ TTD}^{k}(\mathbf{r}_{i})-\mathrm{TTD}^{k-1}(\mathbf{r}_{i})\right)\] \[\in\partial r_{i}(\mathbf{W}_{i}^{k})+\rho(\mathbf{W}_{i}^{k}\mathbf{V}_{i-1}^ {k}-\mathbf{U}_{i}^{k})\mathbf{V}_{i-1}^{k}\top=\partial_{\mathbf{W}_{i}}\mathcal{L}( \mathcal{P}^{k}),\] \[-\alpha(\mathcal{G}_{i}^{k}-\mathcal{G}_{i}^{k-1})\in\partial_{ \mathcal{G}_{i}}\mathcal{L}(\mathcal{P}^{k}). \tag{27}\]

Based on the above relationships, and by the Lipschitz continuity of the activation function on the bounded set \(\{\mathcal{P}:\|\mathcal{P}\|_{F}\leq\mathcal{B}\}\) and the bounded assumption of both \(\mathcal{P}^{k-1}\) and \(\mathcal{P}^{k}\), we have

\[\|\xi_{\mathbf{V}_{N}}^{k}\|_{F}\leq\alpha\|\mathbf{V}_{N}^{k}-\mathbf{V}_{N}^ {k-1}\|_{F}+\gamma\|\mathbf{U}_{N}^{k}-\mathbf{U}_{N}^{k-1}\|_{F}, \xi_{\mathbf{V}_{N}}^{k}\in\partial_{\mathbf{V}_{N}}\mathcal{L}(\mathcal{P}^{k}), \tag{28}\] \[(\omega\|\xi_{\mathbf{V}_{N}}^{k}\|_{F}\leq(L_{R}+\alpha)\|\mathbf{V}_{N}^ {k}-\mathbf{V}_{N}^{k-1}\|_{F}+\gamma\|\mathbf{U}_{N}^{k}-\mathbf{U}_{N}^{k-1}\|_{F})\ \text{proximal-linear}\] \[\|\xi_{\mathbf{V}_{N}}^{k}\|_{F}\leq\rho\mathcal{B}\|\mathbf{W}_{N}^{k}- \mathbf{W}_{N}^{k-1}\|_{F}+\rho\mathcal{B}\|\mathbf{V}_{N}^{k}-\mathbf{V}_{N}^{k-1}\|_{F}, \xi_{\mathbf{U}_{N}}^{k}\in\partial_{\mathbf{U}_{N}}\mathcal{L}(\mathcal{P}^{k}),\] \[\|\xi_{\mathbf{W}_{N}}^{k}\|_{F}\leq 2\rho\mathcal{B}^{2}\|\mathbf{V}_{N-1}^ {k}-\mathbf{V}_{N-1}^{k-1}\|_{F}+\rho\mathcal{B}\|\mathbf{V}_{N}^{k}-\mathbf{V}_{N}^{k-1}\|_{F}\] \[+\tau\|\mathrm{TTD}^{k}(\mathbf{r}_{N})-\mathrm{TTD}^{k-1}(\mathbf{r}_{N}) \|_{F}, \xi_{\mathbf{W}_{N}}^{k}\in\partial_{\mathbf{W}_{N}}\mathcal{L}(\mathcal{P}^{k}),\] \[\|\xi_{\mathcal{G}_{N}}^{k}\|_{F}\leq\alpha\|\mathcal{G}_{N}^{k}- \mathcal{G}_{N}^{k-1}\|_{F}, \xi_{\mathcal{G}_{i}}^{k}\in\partial_{\mathcal{G}_{i}}\mathcal{L}(\mathcal{P}^{k}).\]

and for \(i=N-1,\ldots,1\),

\[\|\xi_{\mathbf{V}_{i}}^{k}\|_{F}\leq\gamma L_{\mathcal{B}}\|\mathbf{U}_{i}^ {k}-\mathbf{U}_{i}^{k-1}\|_{F}, \xi_{\mathbf{V}_{i}}^{k}\in\partial_{\mathbf{V}_{i}}\mathcal{L}(\mathcal{P}^{k}),\] (29) \[\|\xi_{\mathbf{U}_{i}}^{k}\|_{F}\leq\rho\mathcal{B}\|\mathbf{V}_{i-1}^{k}- \mathbf{V}_{i-1}^{k-1}\|_{F}+\rho\mathcal{B}\|\mathbf{W}_{i}^{k}-\mathbf{W}_{i}^{k-1}\|_{F}+ \alpha\|\mathbf{U}_{i}^{k}-\mathbf{U}_{i}^{k-1}\|_{F}, \xi_{\mathbf{U}_{i}}^{k}\in\partial_{\mathbf{U}_{i}}\mathcal{L}(\mathcal{P}^{k}),\] \[\|\xi_{\mathbf{W}_{i}}^{k}\|_{F}\leq(2\rho\mathcal{B}^{2}+\rho \mathcal{B})\|\mathbf{V}_{i-1}^{k}-\mathbf{V}_{i-1}^{k-1}\|_{F}+\tau\|\mathrm{TTD}^{k}( \mathbf{r}_{i})-\mathrm{TTD}^{k-1}(\
where

\[\delta:=\max\{\gamma,\alpha+\rho\mathcal{B},\alpha+\gamma L_{\mathcal{B}},2\rho \mathcal{B}+2\rho\mathcal{B}^{2},\alpha+\sqrt{N}\tau\mathcal{B}^{N-1}\},\]

(or, for the prox-linear case, \(\delta:=\max\{\gamma,L_{R}+\alpha+\rho\mathcal{B},\alpha+\gamma L_{\mathcal{B}},2 \rho\mathcal{B}+2\rho\mathcal{B}^{2},\alpha+\sqrt{N}\tau\mathcal{B}^{N-1}\})\). 

**Definition 4** (Critical point [38, 39]).: _A necessary condition for \(\mathbf{x}\) to be a minimizer of a proper and lower semicontinuous (PLSC) function \(f\) is that_

\[\mathbf{0}\in\partial f(\mathbf{x}). \tag{31}\]

_A point that satisfies (31) is called limiting-critical or simply critical._

**Definition 5** (Global convergence [44, 45]).: _Any iterative algorithm for solving an optimization problem over a set \(X\), is said to be **globally convergent** if for any starting point \(\mathbf{x}_{0}\in X\), the sequence generated by the algorithm always has an accumulation critical point._

To build the global convergence of our iterative sequence \(\{\mathcal{P}^{k}\}_{k\in\mathbb{N}}\) from Algorithm 1, the function \(\mathcal{L}(\mathcal{W},\mathcal{V},\mathcal{U},\mathcal{G})\) needs to have the Kurdyka Lojasiewicz (KL) property as follows

**Definition 6** (KL property [35, 36]).: _A real function \(f:\mathbb{R}^{p}\to(-\infty,+\infty]\) has the Kurdyka Lojasiewicz (KL) property, namely, for any point \(\bar{\mathbf{u}}\in\mathbb{R}^{p}\), in a neighborhood \(N(\bar{\mathbf{u}},\sigma)\), there exists a desingularizing function \(\phi(s)=cs^{1-\theta}\) for some \(c>0\) and \(\theta\in[0,1)\) such that_

\[\phi^{\prime}(|f(\mathbf{u})-f(\bar{\mathbf{u}})|)\mathrm{d}(0,\partial f(\mathbf{u}))\geq 1 \tag{32}\]

_for any \(\mathbf{u}\in N(\bar{\mathbf{u}},\sigma)\) and \(f(\mathbf{u})\neq f(\bar{\mathbf{u}})\)._

The real analytic and semi-algebraic functions, which are related to KL property, are introduced below.

**Definition 7** (Real analytic [46]).: _A function \(h\) with domain an open set \(U\subset\mathbb{R}\) and range the set of either all real or complex numbers, is said to be real analytic at \(u\) if the function \(h\) may be represented by a convergent power series on some interval of positive radius centered at \(u\), i.e., \(h(x)=\sum_{j=0}^{\infty}\alpha_{j}(x-u)^{j}\), for some \(\{\alpha_{j}\}\subset\mathbb{R}\). The function is said to be real analytic on \(V\subset U\) if it is real analytic at each \(u\in V\)[46, Definition 1.1.5]. The real analytic function \(f\) over \(\mathbb{R}^{p}\) for some positive integer \(p>1\) can be defined similarly._

**Definition 8** (Semi-algebraic [36]).: _A subset \(S\) of \(\mathbb{R}^{p}\) is a real **semi-algebraic set** if there exists a finite number of real polynomial functions \(g_{ij},h_{ij}\): \(\mathbb{R}^{p}\to\mathbb{R}\) such that \(S=\cup_{j=1}^{q}\cap_{i=1}^{m}\left\{\mathbf{u}\in\mathbb{R}^{p}:g_{ij}(\mathbf{u})=0\right.\) and \(h_{ij}(\mathbf{u})<0\}.\) In addition, a function \(h:\mathbb{R}^{p+1}\to\mathbb{R}\cup+\infty\) is called **semi-algebraic** if its graph \(\{(\mathbf{u},t)\in\mathbb{R}^{p+1}:h(\mathbf{u})=t\}\) is a real semi-algebraic set._

Based on the above definitions, the following lemma can be obtained.

**Lemma 9**.: _Most of the commonly used NN training models (5) can be verified to satisfy the following_

* _the loss function_ \(\ell\) _is a proper lower semicontinuous and nonnegative function. For example, the squared, logistic, hinge, or cross-entropy losses._
* _the activation functions_ \(\sigma_{i}(i=1\ldots,N-1)\) _are Lipschitz continuous on any bounded set. For example, ReLU, leaky ReLU, sigmoid, hyperbolic tangent, linear, polynomial, or softplus activations._
_._
3. _the regularizers_ \(\tau_{i}\) _and_ \(s_{i}(i=1,\ldots,N)\) _are nonegative lower semicontinuous convex functions._ \(\tau_{i}\) _and_ \(s_{i}\) _are the squared_ \(\ell_{2}\) _norm, the_ \(\ell_{1}\) _norm, the elastic net, the indicator function of some nonempty closed convex set (such as the nonnegative closed half-space, box set or a closed interval_ \([0,1]\)_), or 0 if no regularization._
4. _all these functions_ \(\ell,\sigma_{i},\tau_{i}\) _and_ \(s_{i}(i=1,\ldots,N)\) _are either real analytic or semialgebraic, and continuous on their domains._

_Accordingly, the objective function \(\mathcal{L}(\mathcal{W},\mathcal{V},\mathcal{U},\mathcal{G})\) in (5) has **Kurdyka Lojasiewicz (KL)** property._

Proof.: **On the loss function \(\ell\):** Since these losses are all nonnegative and continuous on their domains, they are proper lower semicontinuous and lower bounded by 0. In the following, we only verify that they are either real analytic or semialgebraic.

1. If \(\ell(t)\) is the squared \((t^{2})\) or exponential \((e^{t})\) loss, then according to [46], they are real analytic.
2. If \(\ell(t)\) is the logistic loss (\(\log(1+\mathrm{e}^{-t})\)), since it is a composition of logarithm and exponential functions which both are real analytic, thus according to [46], the logistic loss is real analytic.
3. If \(\ell(\boldsymbol{u};\boldsymbol{y})\) is the cross-entropy loss, i.e., given \(\boldsymbol{y}\in\mathbb{R}^{d_{N}},\ell(\boldsymbol{u};\boldsymbol{y})=- \frac{1}{d_{N}}[\langle\boldsymbol{y},\log\widehat{\boldsymbol{y}}(\boldsymbol {u})\rangle+\langle\boldsymbol{1}-\boldsymbol{y},\log(\boldsymbol{1}- \widehat{\boldsymbol{y}}(\boldsymbol{u}))\rangle]\), where \(\log\) is performed elementwise and \((\widehat{\boldsymbol{y}}(\boldsymbol{u})_{i})_{1\leq i\leq d_{N}}:=((1+ \mathrm{e}^{-u_{i}})^{-1})_{1\leq i\leq d_{N}}\) for any \(\boldsymbol{u}\in\mathbb{R}^{d_{N}}\), which can be viewed as a linear combination of logistic functions, then by (a2) and [46], it is also analytic.
4. If \(\ell\) is the hinge loss, i.e., given \(\boldsymbol{y}\in\mathbb{R}^{d_{N}},\ell(\boldsymbol{u};\boldsymbol{y}):= \max\{0,1-\langle\boldsymbol{u},\boldsymbol{y}\rangle\}\) for any \(\boldsymbol{u}\in\mathbb{R}^{d_{N}}\), by [47], it is semialgebraic, because its graph is \(\mathrm{cl}(\mathcal{D})\), the closure of the set \(\mathcal{D}\), where \(\mathcal{D}=\{(\boldsymbol{u},z):1-\langle\boldsymbol{u},\boldsymbol{y} \rangle-z=0,\boldsymbol{1}-\boldsymbol{u}\succ 0\}\cup\{(\boldsymbol{u},z):z=0, \langle\boldsymbol{u},\boldsymbol{y}\rangle-1>0\}\).

**On the activation function \(\sigma_{i}\):** Since all the considered specific activations are continuous on their domains, they are Lipschitz continuous on any bounded set. In the following, we only need to check that they are either real analytic or semialgebraic.

1. If \(\sigma_{i}\) is a linear or polynomial function, then according to [46] is real analytic.
2. If \(\sigma_{i}(t)\) is sigmoid, \((1+\mathrm{e}^{-t})^{-1}\), or hyperbolic tangent, \(\tanh(t):=\frac{\mathrm{e}^{t}-\mathrm{e}^{-t}}{\mathrm{e}^{t}+\mathrm{e}^{-t}}\), then the sigmoid function is a composition \(g\circ h\) of these two functions where \(g(u)=\frac{1}{1+u},u>0\) and \(h(t)=\mathrm{e}^{-t}\) (resp. \(g(u)=1-\frac{2}{u+1},u>0\) and \(h(t)=\mathrm{e}^{2t}\) in the hyperbolic tangent case). According to [46], \(g\) and \(h\) in both cases are real analytic. Thus, sigmoid and hyperbolic tangent functions are real analytic.
3. If \(\sigma_{i}\) is ReLU, i.e., \(\sigma_{i}(u):=\max\{0,u\}\), then we can show that ReLU is semialgebraic since its graph is \(\mathrm{cl}(\,.\mathcal{D})\), the closure of the set \(\mathcal{D}\), where \(\mathcal{D}=\{(u,z):u-z=0,u>0\}\cup\{(u,z):z=0,-u>0\}\).
4. Similar to the ReLU case, if \(\sigma_{i}\) is leaky ReLU, i.e., \(\sigma_{i}(u)=u\) if \(u>0\), otherwise \(\sigma_{i}(u)=au\) for some \(a>0\), then we can similarly show that leaky ReLU is semialgebraic since its graph is \(\mathrm{cl}(\mathcal{D})\), the closure of the set \(\mathcal{D}\), where \(\mathcal{D}=\{(u,z):u-z=0,u>0\}\cup\{(u,z):au-z=0,-u>0\}\).
5. If \(\sigma_{i}\) is polynomial, then according to [46], it is real analytic.
6. If \(\sigma_{i}\) is softplus, i.e., \(\sigma_{i}(u)=\frac{1}{t}\log(1+\mathrm{e}^{tu})\) for some \(t>0\), since it is a composition of two analytic functions \(\frac{1}{t}\log(1+u)\) and \(\mathrm{e}^{tu}\), then according to [46], it is real analytic.

**On \(\tau_{i}(\boldsymbol{W}_{i}),s_{i}(\boldsymbol{V}_{i})\):** By the specific forms of these regularizers, they are nonnegative, lower semicontinuous and continuous on their domain. In the following, we only need to verify they are either real analytic and semialgebraic.
* the squared \(\ell_{2}\) norm \(\|\cdot\|_{2}^{2}\): According to [47], the \(\ell_{2}\) norm is semialgebraic, so is its square where \(g(t)=t^{2}\) and \(h(\boldsymbol{W})=\|\boldsymbol{W}\|_{2}\).
* the squared Frobenius norm \(\|\cdot\|_{F}^{2}\): The squared Frobenius norm is semialgebraic since it is a finite sum of several univariate squared functions.
* the elementwise \(1\)-norm \(\|\cdot\|_{1,1}\): Note that \(\|\boldsymbol{W}\|_{1,1}=\sum_{i,j}|\boldsymbol{W}_{ij}|\) is the finite sum of absolute functions \(h(t)=|t|\). According to [47], the absolute value function is semialgebraic since its graph is the closure of the following semialgebraic set \(\mathcal{D}=\{(t,s):t+s=0,-t>0\}\cup\{(t,s):t-s=0,t>0\}\). Thus, the elementwise \(1\)-norm is semialgebraic.
* the elastic net: Note that the elastic net is the sum of the elementwise \(1\)-norm and the squared Frobenius norm. Thus, by (c2), (c3), and [47], the elastic net is semialgebraic.
* If \(\tau_{i}\) or \(s_{i}\) is the indicator function of nonnegative closed half-space or a closed interval (box constraints), by [47], any polyhedral set is semialgebraic such as the nonnegative orthant \(\mathbb{R}_{+}^{p\times q}=\{\boldsymbol{W}\in\mathbb{R}^{p\times q}, \boldsymbol{W}_{ij}\geq 0,\forall i,j\}\), and the closed interval. Thus, \(\tau_{i}\) or \(s_{i}\) is semialgebraic in this case.

We first verify the KL property of \(\mathcal{L}\). From (5), we have

\[\mathcal{L}(\mathcal{W},\mathcal{V},\mathcal{U},\mathcal{G}):= \mathcal{R}_{n}\left(\boldsymbol{V}_{N};\boldsymbol{Y}\right)+\sum_{i=1}^{N}r_ {i}\left(\boldsymbol{W}_{i}\right)+\sum_{i=1}^{N}s_{i}\left(\boldsymbol{V}_{i}\right)\] \[+\frac{\gamma}{2}\sum_{i=1}^{N}\|\boldsymbol{V}_{i}-\sigma_{i}( \boldsymbol{U}_{i})\|_{F}^{2}+\frac{\rho}{2}\sum_{i=1}^{N}\|\boldsymbol{U}_{i }-\boldsymbol{W}_{i}\boldsymbol{V}_{i-1}\|_{F}^{2}+\frac{\tau}{2}\sum_{i=1}^{N }\|\boldsymbol{\mathcal{W}}_{i}-\mathrm{TTD}(\boldsymbol{r}_{i})\|_{F}^{2},\]

which mainly includes the following types of functions, i.e.,

\[\mathcal{R}_{n}\left(\boldsymbol{V}_{N};\boldsymbol{Y}\right),\tau_{i}\left( \boldsymbol{W}_{i}\right),s_{i}\left(\boldsymbol{V}_{i}\right),\left\| \boldsymbol{V}_{i}-\sigma_{i}\left(\boldsymbol{U}_{i}\right)\right\|_{F}^{2}, \left\|\boldsymbol{U}_{i}-\boldsymbol{W}_{i}\boldsymbol{V}_{i-1}\right\|_{F}^ {2},\sum_{i=1}^{N}\|\boldsymbol{\mathcal{W}}_{i}-\mathrm{TTD}(\boldsymbol{r}_ {i})\|_{F}^{2}.\]

To verify the KL property of the function \(\mathcal{L}\), we consider the above functions one.

On \(\mathcal{R}_{n}(\boldsymbol{V}_{N};\boldsymbol{Y})\): Note that given the output data \(\boldsymbol{Y},\mathcal{R}_{n}(\boldsymbol{V}_{N};\boldsymbol{Y}):=\frac{1}{n }\sum_{j=1}^{n}\ell((\boldsymbol{V}_{N})_{:j},\boldsymbol{y}_{j})\), where \(\ell:\mathbb{R}^{d_{N}}\times\mathbb{R}^{d_{N}}\rightarrow\mathbb{R}_{+}\cup \{0\}\) is some loss function. If \(\ell\) is real analytic (resp. semialgebraic), then \(\mathcal{R}_{n}(\boldsymbol{V}_{N};\boldsymbol{Y})\) is real-analytic (resp. semialgebraic).

On \(\|\boldsymbol{V}_{i}-\sigma_{i}(\boldsymbol{U}_{i})\|_{F}^{2}\) : Note that \(\|\boldsymbol{V}_{i}-\sigma_{i}(\boldsymbol{U}_{i})\|_{F}^{2}\) is a finite sum of simple functions of the form, \(|v-\sigma_{i}(u)|^{2}\) for any \(u,v\in\mathbb{R}\). If \(\sigma_{i}\) is real analytic (resp. semialgebraic), then \(v-\sigma_{i}(u)\) is real analytic (resp. semialgebraic), and further \(|v-\sigma_{i}(u)|^{2}\) is also real analytic (resp. semialgebraic) since \(|v-\sigma_{i}(u)|^{2}\) can be viewed as the composition \(g\circ h\) of these two functions where \(g(t)=t^{2}\) and \(h(u,v)=v-\sigma_{i}(u)\).

On \(\|\boldsymbol{U}_{i}-\boldsymbol{W}_{i}\boldsymbol{V}_{i-1}\|_{F}^{2}\): Note that the function \(\|\boldsymbol{U}_{i}-\boldsymbol{W}_{i}\boldsymbol{V}_{i-1}\|_{F}^{2}\) is a polynomial function with the variables \(\boldsymbol{U}_{i},\boldsymbol{W}_{i}\) and \(\boldsymbol{V}_{i-1}\), and thus according to [46] and [47], it is both real analytic and semialgebraic.

On \(\tau_{i}(\boldsymbol{W}_{i}),s_{i}(\boldsymbol{V}_{i}):\) All \(\tau_{i}\)'s and \(s_{i}\)'s are real analytic or semialgebraic.

On \(\|\boldsymbol{\mathcal{W}}_{i}-\mathrm{TTD}(\boldsymbol{r}_{i})\|_{F}^{2}:\) Note that the function \(\|\boldsymbol{\mathcal{W}}_{i}-\mathrm{TTD}(\boldsymbol{r}_{i})\|_{F}^{2}\) is a polynomial function with the variables \(\boldsymbol{W}_{i},\mathcal{G}_{i}\).
Since each part of the function \(\mathcal{L}\) is either real analytic or semialgebraic, \(\mathcal{L}\) is a subanalytic function [48, p.43]. Furthermore, by the continuity, \(\mathcal{L}\) is continuous in its domain. Therefore, \(\mathcal{L}\) is a KL function according to [49, Theorem 3.1].3 

Footnote 3: Let \(h:\mathbb{R}^{p}\to\mathbb{R}\cup\{+\infty\}\) be a subanalytic function with closed domain, and assume that \(h\) is continuous on its domain, then \(h\) is a KL function.

Based on Lemmas 2, 3, and 9 and conclusions in [36, Section 3.2], the following main theorem can be obtained.

**Theorem 10** (Global Convergence).: _Let \(\{\mathcal{P}^{k}:=(\{\mathbf{W}_{i}^{k}\}_{i=1}^{N},\{\mathbf{V}_{i}^{k}\}_{i=1}^{N}, \{\mathbf{U}_{i}^{k}\}_{i=1}^{N}),\{\mathcal{G}_{i}^{k}\}_{i=1}^{N}\}_{k\in\mathbb{N}}\) be the sequences generated from Algorithm 1. Suppose that \(\tau_{i}\) and \(\mathcal{L}\) are coercive for any \(i=1,\ldots,N\). Then for any \(\alpha,\gamma,\rho,\tau>0\) and any finite initialization \(\mathcal{P}^{0}\), the following hold_

1. \(\{\mathcal{L}(\mathcal{P}^{k})\}_{k\in\mathbb{N}}\) _converges to_ \(\mathcal{L}^{*}\)_._
2. \(\{\mathcal{P}^{k}\}_{k\in\mathbb{N}}\) _converges to a critical point of_ \(\mathcal{L}\)_._
3. _If further the initialization_ \(\mathcal{P}^{0}\) _is sufficiently close to some global minimum_ \(\mathcal{P}^{*}\) _of_ \(\mathcal{L}\)_, then_ \(\mathcal{P}^{k}\) _converges to_ \(\mathcal{P}^{*}\)_._
4. _Let_ \(\theta\) _be the KL exponent of_ \(\mathcal{L}\) _at_ \(\mathcal{P}^{*}\)_. There hold: (a) if_ \(\theta=0\)_, then_ \(\{\mathcal{P}^{k}\}_{k\in\mathbb{N}}\) _converges in a finite number of steps; (b) if_ \(\theta\in(0,\frac{1}{2}]\)_, then_ \(\|\mathcal{P}^{k}-\mathcal{P}^{*}\|_{F}\leq C\eta^{k}\) _for all_ \(k\geq k_{0}\)_, for certain_ \(k_{0}>0,C>0,\eta\in(0,1)\)_; and (c) if_ \(\theta\in(\frac{1}{2},1)\)_, then_ \(\|\mathcal{P}^{k}-\mathcal{P}^{*}\|_{F}\leq Ck^{-\frac{1-\theta}{2\theta-1}}\) _for_ \(k\geq k_{0}\)_, for certain_ \(k_{0}>0,C>0\)_._
5. \(\frac{1}{K}\sum_{k=1}^{K}\|\mathbf{g}^{k}\|_{F}^{2}\to 0\) _at the rate_ \(\mathcal{O}(1/K)\) _where_ \(\mathbf{g}^{k}\in\partial\mathcal{L}(\mathcal{P}^{k})\)_._

Lispchitz differentiable property is a required for nonconvex optimizations with multi-block variables to build the convergence in the existing literature [50]. However, the NN training problem (5) in this paper generally does not satisfy such a condition. For example, when ReLU activation is used. Theorem 10 establishes the global convergence under a very mild condition that most NN models satisfy.

**Extension to ResNets [34]:** the theoretical results in Theorem 10 can be extended to ResNets by considering the following optimization problem

\[\min_{\mathcal{W},\mathcal{V}}\mathcal{L}_{0}(\mathcal{W},\mathcal{V})\text{ subject to }\mathbf{U}_{i}=\mathbf{W}_{i}\mathbf{V}_{i-1},\mathbf{V}_{i}-\mathbf{V}_{i-1}=\sigma_{i}(\mathbf{U}_{i} ),\mathbf{\mathcal{W}}_{i}=\text{TTD}(\mathbf{r}_{i})\quad i=1,\ldots,N,\]

where the residual term \(\mathbf{V}_{i}-\mathbf{V}_{i-1}\) is considered instead of \(\mathbf{V}_{i}\). The corresponding algorithm can be easily modified from Algorithm 1.

## 4 Case Study

In this experiment, to evaluate the effectiveness and efficiency of our proposed method, NN model (5) training with different compression ratios (determined by TT-rank \(\mathbf{r}_{i}\)) is conducted on the image classification task. In terms of the NN model, ReLU activation, the squared loss, and the network architecture being an MLP with ten hidden layers are considered here. The number of hidden units in each layer is \(2^{9}=512\). The neural network is trained on the MNIST dataset, which is a handwritten digits dataset. The size of each input image is \(d_{0}=28\times 28=784\) and the output dimension is \(d_{11}=10\). The numbers of training and test samples are 60,000 and 10,000, respectively. For comparison, SGD is also considered as a benchmark method, where the learning rate is 0.001.

For each experiment, the same mini-batch sizes (512) and initializations for all algorithms. All the experiments are repeated ten times to obtain the average performance. Specifically, all the
weights \(\{\mathbf{W}_{i}\}_{i=1}^{N}\) are initialized from a Gaussian distribution with a standard deviation of 0.01. The auxiliary variables \(\{\mathbf{U}_{i}\}_{i=1}^{N}\), state variables \(\{\mathbf{V}_{i}\}_{i=1}^{N}\), TT-cores \(\{\mathcal{G}_{i}\}_{i=1}^{N}\) are initialized by a single forward pass [50]. Under these settings, the training loss, training accuracy, and test accuracy are shown in Table 1. With a smaller CR (\(\frac{\#}{\#\theta\text{parameters after compression}}\)), a higher training loss is observed. Our proposed method with CR<1 can outperform the uncompressed method and SGD. In addition, the curves of the training loss and test accuracy are plotted in Figure 2.

Figure 2(a) shows that the proposed method converges with different compression rates. The training loss of our proposed method also shows the monotone decreasing trend, which verified the statements in Theorem 10. Figure 2(b) shows that, for different CR (<1), the test accuracy of our proposed method keeps increasing as the number of iterations increases. When CR=1 (the model without compression), the test accuracy increases first and then decreases. This result demonstrates that model compression can prevent overfitting. In addition, our proposed method with CR<1 can outperform SGD significantly in terms of test accuracy.

\begin{table}
\begin{tabular}{l r r r r r r} \hline \hline  & \multicolumn{2}{c}{Total Loss} & \multicolumn{2}{c}{Training Accuracy} & \multicolumn{2}{c}{Test Accuracy} \\ \hline CR & Mean & Std & Mean & Std & Mean & Std \\ \hline
1.00 & 0.6726 & 0.0170 & 0.8921 & 0.0017 & 0.8794 & 0.0034 \\
0.77 & 0.9804 & 0.0331 & 0.9768 & 0.0003 & 0.9673 & 0.0014 \\
0.34 & 1.3564 & 0.0440 & 0.9677 & 0.0006 & 0.9623 & 0.0009 \\
0.09 & 3.4612 & 0.0264 & 0.9439 & 0.0009 & 0.9443 & 0.0014 \\
0.02 & 4.5591 & 0.0331 & 0.9046 & 0.0027 & 0.9090 & 0.0021 \\ SGD & NA & NA & 0.8940 & 0.0295 & 0.9002 & 0.0274 \\ \hline \hline \end{tabular}
\end{table}
Table 1: Results of tenBCD algorithm with different compression ratios.

Figure 2: The convergence analysis of tenBCD algorithm with different compression ratios: (a) training loss; (b) test accuracy.
## 5 Conclusion

In this paper, a holistic framework is proposed for tensor decomposition-based NN model compression by formulating TT decomposition-based NN training as a nonconvex optimization problem. The framework can be extended to other formats of tensor decomposition such as Tucker decomposition, and CP decomposition. For the first time in the literature on tensor decomposition-based NN model compression, global convergence is guaranteed for the proposed tensor BCD (tenBCD) algorithm. Specifically, tenBCD converges to a critical point at a rate of \(\mathcal{O}(1/k)\), where \(k\) is the number of iterations. The empirical experiment shows that the proposed method can converge and run efficiently in practice. Compared with SGD, the proposed method can maintain a high compression rate and high accuracy simultaneously.
