A vector, a matrix, and a tensor are indicated by a bold lowercase letter, a ∈ ℝ^I, a bold uppercase letter, B ∈ ℝ^{I × J}, and a bold calligraphic letter, C ∈ ℝ^{J₁ × J₂ × ⋯ × J_N}, respectively. An Nth-order tensor, X ∈ ℝ^{I₁ × I₂ × ⋯ × I_N}, can be transformed into a vector, which is denoted by using the same character, x=vec(X)∈ℝ∏n=1NIn. An (i₁, i₂, ..., i_N)-element of X is denoted by x_{i1, i2, ..., iN} or [X]_{i1, i2, ..., iN}. The operators ⊡ represent the Hadamard product. ·^† represents the Moore-Penrose pseudoinverse. sign(X) ∈ ℝ^{I₁ × I₂ × ⋯ × I_N} and abs(X) ∈ ℝ^{I₁ × I₂ × ⋯ × I_N} are, respectively, operations that return the sign and absolute value for each entry of X ∈ ℝ^{I₁ × I₂ × ⋯ × I_N}.

Tensor decomposition (TD) is a mathematical model to represent a tensor as a product of tensors/matrices. There are many TD models such as canonical polyadic decomposition (CPD) [41–43, 55], Tucker decomposition (TKD) [44–46], tensor-train decomposition (TTD) [49, 50], tensor-ring decomposition (TRD) [51], block-term decomposition [56–58], coupled tensor decomposition [59–61], hierarchical Tucker decomposition [62], tensor wheel decomposition [63], fully connected tensor network [64], t-SVD [65], Hankel tensor decomposition [66–69], and convolutional tensor decomposition [21].

Figure 1 shows typical models of tensor decompositions in graphical notation [10]. In graphical notation, each node represents a core tensor and the edges connecting the nodes represent tensor products. As can be seen in the figure, all TD models are expressed by tensor products of core tensors.

In this paper, we will write any TD model as

Note that Equation 1 is a somewhat abstract expression which only represents that the entire tensor X is reconstructed by multiplying the L core tensors G₁, G₂, ..., G_L−1 and G_L. The least squares (LS) based tensor decomposition problem of a given tensor V ∈ ℝ^{I₁ × I₂ × ⋯ × I_N} is given by

This squared error is a non-convex function with respect to overall optimization parameters (G₁, G₂, ..., G_L), but it is a convex quadratic function with respect to only one core tensor G_l for any l ∈ {1, 2, ..., L}. If we focus on optimizing a certain core tensor G_l and temporarily consider all other core tensors as constant variables and can be summarized as a single tensor, denoted as G_−l, then the sub-optimization problem is reduced to an LS-based matrix optimization problem:

where we used abstract notations 〈〈G_l, G_−l〉〉 = 〈〈G₁, G₂, ..., G_L〉〉, and mat_lm, m ∈ {1, 2, 3}, are the corresponding appropriate matricization operators for the certain l-th core tensor. Solving sub-optimization shown in Equation 3 has a closed-form solution and updating G_l by

always reduces (non-increases) the squared error. In other words, putting θ = (G₁, G₂, ..., G_L) and f(θ)=∥V-〈〈G1,G2,...,GL〉〉∥F2, the following algorithm

provides the result f(θ_k) ≥ f(θ_k+1) for any l ∈ {1, 2, ..., L}. Therefore, the LS-based tensor decomposition problem in Equation 2 can be solved by repeating the steps shown in Equation 6 for every l ∈ {1, 2, ..., L}, and this is called the alternating least squares (ALS) algorithm [42, 43, 47, 50, 51]. ALS is a workhorse algorithm in TDs because it has no hyper-parameter (e.g., step-size of gradient descent), and the objective function can be stably reduced (monotonically non-increasing).

Tensor completion is a task to estimate missing values in an observed incomplete tensor by using the low-rank structure of a tensor. In case of low-rank matrix completion, theories, algorithms, and applications are well studied [30, 35, 70–73]. Unlike matrices, which have a unique rank definition, tensors have various ranks (e.g., CP rank, Tucker rank, and TT rank), and the meaning of low rank is different for decomposition models [26, 27, 64, 74–77]. Since the appropriate TD model depends on the application, it is desirable to have an environment in which various TD models can be freely selected and tested.

Here, we introduce a basic formulation of LS-based tensor completion as follows:

where T ∈ ℝ^{I₁ × I₂ × ⋯ × I_N} is an observed incomplete tensor and W is a mask tensor of the same size as T. The entries of W are given by

To solve the problem in Equation 14, many algorithms have been studied [21, 26, 28, 30, 78–80]. The gradient descent-based optimization algorithm [26] is proposed for the CP model. For the Tucker model, which imposes orthogonality on factor matrices, algorithms based on optimization on a manifold are proposed [28, 81]. However, these algorithms are designed specifically for specific TD models and it is difficult to generalize them to various TD models.

Expectation-maximization alternating least squares (EM-ALS) [82] is a versatile tensor completion algorithm that is less dependent on differences in TD models. In fact, EM-ALS is incorporated into various tensor completion algorithms such as TMac [29], TMac-TT [76], MTRD [77], MDT-Tucker [66, 67], and SPC [34]. The algorithm of EM-ALS can be derived from majorization-minimization (MM) [40, 83] which iteratively minimizes an auxiliary function g(X, X′) that serves as an upper bound on the objective function f(X). In more detail, an auxiliary function holds the following conditions

and the update rule of MM is given by

This ensures that the objective function is monotonically decreasing as follows:

Specifically, the objective function and its auxiliary functions are given by

where W¯=1-W is a tensor that flips the 0 and 1 of W. Looking at the additional terms in Equation 20, we can see that they clearly satisfy the condition shown in Equation 16. Furthermore, update rule can be transformed as

where Vk=W⊡T+W¯⊡Xk. Since the structure of Equation 21 is the same as Equation 8, it can be solved by ALS. Finally, EM-ALS can be given by the following algorithm:

where X₀ ∈ 𝕊_t is required for the initialization. The steps shown in Equations 22, 23 are called the E-step and the M-step, respectively. Since M-step is an iterative algorithm, EM-ALS becomes a double iterative algorithm and is inefficient. Since the auxiliary function can also be decreased by one step of ALS Xk+1=UVk,𝕊tALS(Xk), the M-step can often be replaced with this. In EM-ALS, the operation that depends on the TD model is only the part in Equation 23. This means that any update rule of various LS-based TDs that decreases the objective function can be used as is.

Remark 2. Plug-and-play of TD algorithms is possible for tensor completion in EM-ALS.

Robust tensor decomposition (RTD) is a task to reconstruct a low-rank tensor from an observed tensor with outliers. Typically, the outlier components are assumed to be sparse (or follow a Laplace distribution) and the problem is formulated as a tensor decomposition based on ℓ₁ loss as follows:

where B ∈ ℝ^{I₁ × I₂ × ⋯ × I_N} is an observed tensor that includes outlier components. Various TD models and algorithms for RTD have been proposed, such as CPD [84], TKD [85], TRD [86], t-SVD [87], and more sophisticated models for specific domains [37].

Zhang and Ding [85] have proposed Tucker-based RTD using the alternating direction method of multipliers (ADMM). Zhang and Ding's ADMM formulation is also versatile, like EM-ALS, and we review it briefly here. First, we introduce a new variable E ∈ ℝ^{I₁ × I₂ × ⋯ × I_N} and add a constraint E = B−X to the RTD problem shown in Equation 24. For the constrained optimization problem, the augmented Lagrangian is given by

where Λ ∈ ℝ^{I₁ × I₂ × ⋯ × I_N} is a Lagrange multiplier, and β > 0 is a hyperparameter. In each step of ADMM, X is updated as a minimizer of LRTD(X,E,Λ) with respect to X, E is updated as a minimizer of LRTD(X,E,Λ) with respect to E, and Λ is updated by the method of multipliers. ADMM algorithm is given by

where initializations X₀, E₀, and Λ₀ are required. The sub-optimization problems shown in Equations 26, 27 can be solved by LS-based TD algorithm and soft-thresholding, respectively. In practice, Equations 26, 27 can be replaced as

where soft_ρ(·) with ρ > 0 is a soft-thresholding operator:

In a similar way to EM-ALS, Proj𝕊tALS(Vk,Xk) often be replaced with UVk,𝕊tALS(Xk) in Equation 30, and the any update rule of various LS-based TD can be used as is.

Remark 3. Plug-and-play of TD algorithms is possible for robust tensor decomposition in ADMM.

In this paper, we consider solving a more versatile tensor reconstruction problem than tensor completion and RTD. First, we assume the following linear observation model

for j ∈ {1, 2, ..., J}, where an observation b_j is obtained from innerproduct between a given tensor A_j ∈ ℝ^{I₁ × I₂ × ⋯ × I_N} and an unknown low-rank tensor X ∈ ℝ^{I₁ × I₂ × ⋯ × I_N} with noise component n_j. This observation model can be transformed into vector-matrix form as

where b=[b1,b2,...,bJ]⊤∈ℝJ is an observed signal, n=[n1,n2,...,nJ]⊤∈ℝJ is noise component, x = vec(X) is a vector form of low-rank tensor, and

is a design matrix. Introducing a low-rank TD constraint x ∈ 𝕊 instead of X ∈ 𝕊_t and loss function based on the noise assumption of n, the optimization problem can be given as

where 𝕊: = {vec(〈〈G₁, G₂, ..., G_L〉〉) | G₁ ∈ D₁, G₂ ∈ D₂, ..., G_L ∈ D_L} is a space of vectorization of low-rank tensors, D(·, ·) stands for loss function such as ℓ₂ loss, ℓ₁ loss and generalized Kullback–Leibler (KL) divergence. This study aims to solve the problem shown in Equation 36. Note that x ∈ ℝ^I is a I-dimensional vector in the form, and we consider that x represents a tensor X∈ℝI1×I2×⋯×IN by x = vec(X), where I=∏i=1NIi.

The problem in Equation 36 includes low-rank tensor completion and RTD. When the loss function is ℓ₂, b = vec(W⊡T), and A = diag(vec(W)), then the problem in Equation 36 reduces to the tensor completion problem in Equation 14. When the loss function is ℓ₁, b = vec(B), and A = I, then the problem in Equation 36 reduces to the RTD problem in Equation 24.

Figure 2 shows the concept of the tensor reconstruction problem considered in this study. In our problem formulation in Equation 36 includes various other patterns of tensor reconstruction. Tensor completion and RTD are just a few of them, and there are still many missing patterns. We aim to solve the problem in Equation 36 for any design matrix A ∈ ℝ^{J × I}, and other loss functions such as the generalized Kullback-Leibler (KL) divergence. This generalization is important for various applications. For the design matrix A, the Toeplitz matrix is used in image deblurring, the downsampling matrix is used in the super-resolution task, the Radon transform matrix is used in computed tomography and the random projection matrix is used in compressed sensing [1]. For loss functions, ℓ₂ loss is used in Gaussian noise setting, ℓ₁ loss is used in Laplace noise (or sparse noise) setting, and generalized KL divergence is used in Poisson noise setting.

Furthermore, we consider a penalized version as follow:

where α ≥ 0 and p(θ) is a penalty function for core tensors θ = (G₁, G₂, ..., G_L). This is a generalization of the problem in Equation 36, since Equation 37 with α = 0 is equivalent to Equation 36. This plays an important role for the introduction of Tikhonov regularization, sparse regularization, low-rank regularization (e.g., nuclear norm), and smoothing. In particular, Tikhonov regularization of the core tensors reduces the problem of non-uniqueness of scale and improves the convergence of the TD algorithm [88].

The key idea to solve a diverse set of problems formulated in Equation 36 is to employ the plug-and-play (PnP) approach of TD algorithms such as EM-ALS and ADMM. In this study, we first solve the case of ℓ₂ loss using LS-based TD with the MM framework, which is a generalization of EM-ALS. Furthermore, we use it to solve the cases of ℓ₁ loss and KL divergence with the ADMM framework. Thus, we call the proposed algorithm ADMM-MM. By replacing the LS-based TD module in a PnP manner, various types of TD can be easily generalized and applied to various applications.

In this section, we explain one-by-one how to optimize the problem in Equation 37 for three loss functions: ℓ₂ loss, ℓ₁ loss, and generalized KL divergence. Then, we explain how to combine these three cases as the ADMM-MM algorithm.

Finally, the proposed ADMM-MM algorithm that supports three typical loss functions can be summarized in Algorithm 1. Our ADMM-MM allows for the immediate application of sophisticated TD models to other tasks such as completion, robust reconstruction, and compressive sensing. By simply selecting the update rule in the 7th line, we can accommodate the three loss functions. The majorization step on the 9th line is important to allow for accommodating a variety of design matrices A. Many LS-based TD with penalty can be applied into the ADMM-MM algorithm in a plug-and-play manner at the 10th line. The basic expectation of plug-and-play module, Proj_𝕊γ(v, x) or Uv,𝕊γ(x), is that it has a monotonically decreasing property with respect to the penalized squared error shown in Equation 39. Note that our algorithm can also involve tensor nuclear norm regularization [65] although not a direct tensor decomposition. That is, a proximal mapping given in closed form of nuclear norm regularization can be replaced directly with Proj_𝕊γ(v, x).

The proposed algorithm can be regarded as a unified and extended algorithm of EM-ALS and Zhang and Ding's ADMM. When A is the diagonal matrix of vectorization of the mask tensor W (e.g., A = diag(vec(W))) and ℓ₂ loss is assumed, the ADMM-MM algorithm results in EM-ALS. When A = I and ℓ₁ loss is assumed, the ADMM-MM algorithm results in Zhang and Ding's ADMM.

From the perspective of an optimization algorithm that uses ADMM with MM, the proposed algorithm can be regarded as a special case of linearized ADMM [92, 93]. However, the linearized ADMM work [92] considers only convex optimizations, while the other work [93] considers only non-convex matrix norms which have effective proximal mappings. The major difference in the proposed method from the above methods is that the update rule of the iterative algorithm is applied to penalized TD instead of proximal mapping. Our proposal to plug-and-play various TDs is a new and meaningful attempt, which will greatly improve the applicability of TDs.

The computational complexity of a single iteration in the ADMM-MM algorithm is often dominated by the LS-based TD at the 10th line. The complexity of lines 6 through 9 is O(Ω), where Ω is the number of nonzero elements in the design matrix A ∈ ℝ^{J × I}. We often assume that the design matrix is sparse and J ≤ Ω ≪ IJ. On the other hand, the complexity of LS-based TD for N-th order tensor V ∈ ℝ^{I₁ × I₂ × ⋯ × I_N} depends on the model and algorithm. CPD/NNCPD can be solved by ALS and HALS [2, 8]. The complexity of CP-ALS is O(NIRCP)+O(NRCP3), while the complexity of HALS is O(NIRCP), assuming that the CP rank is denoted as R_CP and I=∏i=1NIi. ALS is computationally more expensive than HALS because it requires matrix inversion. TKD is usually solved by the higher-order orthogonal iteration (HOOI) algorithm [48], which solves N eigenvalue problems alternately. Since the symmetric matrix of size (I_n, I_n) by tensor-matrix multiplication is usually more expensive than its eigenvalue problem, then the complexity of HOOI is O(NIRTK), assuming that the Tucker rank is (R_TK, R_TK, ..., R_TK). TTD is usually solved by TT-ALS with QR orthogonalization [50]. The advantage of this algorithm is that QR orthogonalization allows each least-squares solution to be found without matrix inversion. The complexity of TT-ALS is O(NIRTT), assuming that the TT-rank is (R_TT, R_TT, ..., R_TT). TRD is usually solved by TR-ALS [51]. The complexity of TR-ALS is O(NIRTR2)+O(NRTR6), assuming that the TR-rank is (R_TR, R_TR, ..., R_TR). From a complexity perspective, TKD and TTD are highly efficient. TRD has a higher complexity with respect to TR rank. CPD often requires a larger CP rank, making it more expensive than other models.

Unfortunately, this study does not include new theoretical results on the convergence of the ADMM-MM algorithm using Ls-based TD algorithms for sub-optimization. Instead, this section introduces existing theoretical results on convergence related to TD, EM/MM, and ADMM, and discusses their relevance to this study.

Tensor reconstruction tasks in this experiment include tensor denoising, completion, de-blurring, and super-resolution. We used an RGB image named “facade” represented as a third-order tensor of size 256 × 256 × 3 for ground truth in denoising, completion, de-blurring, and super-resolution tasks.

For tensor denoising task, we set A = I. Gaussian noise, salt-and-pepper noise, and Poisson noise are added for individual tasks. For tensor completion tasks, we set A = diag(vec(W)) with a randomly generated mask tensor W ∈ {0, 1}^{256 × 256 × 3}, where 90% of the entries are missing. For de-blurring tasks, we used a motion blur window of size 21 × 21, and the constructed block Toeplitz matrix is used as A. For super-resolution tasks, we used a Lanczos2 kernel for downsampling to 1/4 size, and the constructed downsampling matrix is used as A.

In this experiment, we demonstrate that the proposed algorithm connects various TD models with various image processing tasks.