A Novel Log-Based Tensor Completion Algorithm

2.1 Notations and definitions

Throughout this paper, each scalar is denoted by a normal letter, e.g. a, b, c, ...; each vector is denoted by a boldfaced lower-case letter, e.g. a, b, c, ...; each matrix is denoted by an upper-case letter, e.g. A, B, C, ...; each tensor is written as an italic letter, e.g. an N-order tensor is denoted as $\mathcal{X} \in R^{I_{1} \times \cdots \times I_{N}}$, where I_k is the dimensional size corresponding to mode-k, k=1, ..., N; each element of tensor $\mathcal{X}$ is denoted as $x_{i_{1} \cdots i_{N}} \in \mathcal{X}$.

The Frobenius norm of $\mathcal{X}$ is defined as the square root of the inner product of twofold tensors:

$\|\mathcal{X}\|_{F}=\sqrt{\langle \chi, \chi\rangle}=\sqrt{\sum_{i_{1}=1}^{I_{1}} \cdots \sum_{i_{N}=1}^{I_{N}} x_{i_{1} \cdots i_{N}}^{2}}$   (1)

The CP decomposition factorizes the target tensor into a sum of component one-rank tensors. For example, a third-order tensor $\mathcal{X} \in \mathbb{R}^{I \times J \times K}$ can be written as:

$\chi \approx \sum_{r=1}^{R} a_{r} \circ b_{r} \circ c_{r}$  (2)

where, R is a positive integer; $a_{r} \in \mathbb{R}^{I}$, $b_{r} \in \mathbb{R}^{J}$ and $c_{r} \in \mathbb{R}^{K}$ for r=1, ..., R. Elementwise, formula (2) can be written as

$x_{i j k} \approx \sum_{r=1}^{R} a_{i r} b_{j r} c_{k r}$

for $i=1, \ldots, I, j=1, \ldots, J, k=1, \ldots, K$  (3)

The Tucker decomposition, a higher-order principal component analysis (PCA), decomposes the target tensor into a core tensor multiplied (or transformed) by a matrix along each mode. In the three-way case where $\mathcal{X} \in \mathbb{R}^{I \times J \times K}$, the Tucker decomposition can be expressed as:

$\begin{aligned} \chi & \approx \mathcal{G} \times_{1} A \times_{2} B \times_{3} C \\=& \sum_{p=1}^{P} \sum_{q=1}^{Q} \sum_{r=1}^{R} g_{p q r} a_{r} \circ b_{r} \circ c_{r} \end{aligned}$  (4)

where, $A \in \mathbb{R}^{I \times P}$, $B \in \mathbb{R}^{J \times Q}$, and $C \in \mathbb{R}^{K \times R}$ are the factor matrices (usually orthogonal), i.e. the principal components in each mode. The tensor $\mathcal{G} \in \mathbb{R}^{P \times Q \times R}$ is called the core tensor, whose entries reflect the level of interaction between different components. Elementwise, formula (4) can be written as:

$x_{i j k} \approx \sum_{r=1}^{R} \sum_{q=1}^{Q} \sum_{r=1}^{R} g_{p q r} a_{i p} b_{j q} c_{k r}$

for $i=1, \ldots, I, j=1, \ldots, J, k=1, \ldots, K$  (5)

where, P, Q, and R are the number of components (columns) in factor matrices A, B and C, respectively.

The Tensor-Train decomposition represents a d-way tensor $\chi$ as d3-way tensors $\chi^{(1)}, \dots, \chi^{(d)}$, which are called the TT-cores. The k-th TT-core has dimensions r_k-1, n_k, r_k. where r_k-1, r_k are called the TT-ranks. Specifically, each entry of $\mathcal{X} \in \mathbb{R}^{I_{1} \times \cdots \times I_{d}}$ can be determined by:

$a_{i_{1} \cdots i_{d}}=\chi_{i_{1}}^{(1)} \cdots \mathcal{X}_{i_{d}}^{(d)}$  (6)

Definition 2.1 ([32]) Let $\chi \in R^{I_{1} \times \cdots \times I_{N}}$ be an N-order tensor. The k-th mode matrix unfolding is defined as the matrix $X_{(k)} \in R^{I_{k} \times J}$, $\mathrm{J}=\prod_{i=1, i \neq k}^{N} I_{i}$. The tensor element (i₁, ..., i_N) is mapped to the matrix element (i_k, j), whereNext, several important concepts were presented with the following definitions:

$\mathrm{j}=1+\sum_{m=1, m \neq k}^{N}\left(i_{m}-1\right) J_{i}$

with $J_{i}=\prod_{m=1, m \neq i}^{i-1} I_{m}$  (7)

Definition 2.2 ([32]) Let $X \in R^{I} k^{\times J}$ be a matrix with $J=\prod_{i=1, i \neq k}^{N} I_{i}$. The k-mode (I₁, ..., I_N) tensor folding $\chi$ is defined as the tensor $\chi \in R^{I_{1} \times \cdots \times I_{N}}$. The matrix element (i_k, j) is mapped to the tensor element (i₁, ..., i_N), where j is defined as in formula (7).

Definition 2.3 ([33]) The higher-order singular value decomposition (SVD) of a tensor $\chi \in R^{I_{1} \times \cdots \times I_{N}}$ can be defined as

$\chi=\mathcal{S} \times_{1} U_{1} \times_{2} U_{2} \times_{3} \cdots \times_{N} U_{N}$  (8)

or elementwise as

$\chi\left(i_{1}, \cdots, i_{N}\right)=\sum_{j_{1}=1}^{I_{1}} \cdots \sum_{j_{N}=1}^{I_{N}} \mathcal{S}\left(j_{1}, \cdots, j_{N}\right)$

$U_{1}\left(i_{1}, j_{1}\right) \cdots U_{N}\left(i_{N}, j_{N}\right)$  (9)

where, for 1≤k≤n, the U_k are unitary n_k×n_k matrices, the core tensor $\mathcal{S} \in R^{I_{1} \times \cdots \times I_{N}}$ is such that the (N-1)-order subtensor $\mathcal{S}^{i_{k}=p}$ defined by

$\mathcal{S}_{i_{1} \cdots i_{k-1} i_{k+1} \cdots i_{N}}^{i_{k}=p}:=s_{i_{1} \cdots i_{k-1} p i_{k+1} \cdots i_{N}}$  (10)

satisfies:

(1) all-orthogonality:

$\left\langle\mathcal{S}^{i_{k}=p}, \mathcal{S}^{i_{k}=q}\right\rangle=0, \forall p \neq q, 1 \leq k \leq N$   (11)

(2) the subtensors of the core tensor $\mathcal{S}$ are ordered by their Frobenius norm, i.e., for any 1≤k≤n,

$\left\|\mathcal{S}^{i_{k}=1}\right\|_{F} \geq \cdots \geq\left\|\mathcal{S}^{i_{k}=I_{N}}\right\|_{F} \geq 0$  (12)

N-rank of a tensor is the straightforward generalization of the column (row) rank for matrices. Let r_k denote the k-th rank of an N-order tensor $\chi \in R^{I_{1} \times \cdots \times I_{N}}$. Then, it is the rank of the mode-k unfolding matrix X_(k):

$r_{k}=\operatorname{rank}_{k}(\mathcal{X})=\operatorname{rank}\left(X_{(k)}\right)$  (13)

A tensor of which the n-rank are equal to r_n is called a rank-(r₁, ..., r_N) tensor. For any tensor $\mathcal{X} \in R^{I_{1} \times \cdots \times I_{N}}$, vec($\mathcal{X}$) denotes the vectorization of $\mathcal{X}$.

Next is to introduce the SVD of a matrix. Let there be an n₁×n₂ matrix of rank r, and $X=U(\operatorname{Diag} \sigma(X)) V^{T}$ be the SVD of X, where U and V are respectively n₁×n and n₂×n matrices with orthonormal columns, and the function $\sigma: C^{n_{1} \times n_{2}} \rightarrow R^{n}\left(\mathrm{n}=\min \left\{n_{1}, n_{2}\right\}\right)$ has components $\sigma_{1}(X) \geq \sigma_{2}(X) \geq \cdots \geq \sigma_{n}(X) \geq 0$, i.e. the singular values of the matrix X.

The necessary results in Lewis’s research [34] will be briefly described for subsequent formula derivation.

Definition 2.4 ([34]) For any vector γ in Rⁿ, vector $\tilde{\gamma}$ can be written with components |γ_i| arranged in nonincreasing order. A function f: Rⁿ→[-∞,+∞] is absolutely symmetric, if f(γ)= $f(\tilde{\gamma})$ for any vector γ in Rⁿ.

Theorem 2.1 ([34, Corollary 2.5]) Suppose function f: Rⁿ→(-∞,+∞] is absolutely symmetric, and n₁×n₂ matrix X has σ(X) in domf. Then, the n₁×n₂ matrix Y lies in $\partial(f \circ \sigma)(X)$, if and only if σ(Y) falls in $\partial f(\sigma(\mathrm{X}))$ and there exists a simultaneous SVD of the form U(Diagσ(X))V^T and Y=U(Diagσ(Y))V^T, where U and V are n₁×n and n₂×n matrices with orthonormal columns, respectively. In fact,

$\partial(f \circ \sigma)(X)$$=\left\{U(\operatorname{Diag} \mu) V^{T} | \mu \in \partial f(\sigma(x)), \mathrm{X}=U(\operatorname{Diag} \sigma(X)) V^{T}\right\}$  (14)

2.2 Low n-rank tensor completion problem

This subsection briefly illustrates the low n-rank tensor completion problem. Suppose there exists a partially observed tensor $\mathcal{T}$ under the low n-rank assumption. Then, the Low n-rank tensor completion problem can be expressed as:

$\min _{x} \sum_{i=1}^{N} \operatorname{rank}\left(X_{(i)}\right)$ s.t. $\mathcal{X}_{\Omega}=\mathcal{T}_{\Omega}$  (15)

where, $\chi \in R^{I_{1} \times \cdots \times I_{N}}$ is the decision variable whose size is identical with $\mathcal{T}$; Ω is a subset of index, indicating that the entries of $\mathcal{T}$ in the set Ω are given while the remaining entries are missing. This problem is the special case of the following tensor n-rank minimization problem:

$\min _{x} \sum_{i=1}^{N} \operatorname{rank}\left(X_{(i)}\right)$ s.t. $\mathbb{A}(\mathcal{X})=\boldsymbol{b}$  (16)

where, the linear map $\mathbb{A}: R^{I_{1} \times \cdots \times I_{N}} \rightarrow R^{p}$ with $\mathrm{p} \leq \prod_{i=1}^{N} n_{i}$, and vector $\mathbf{b} \in R^{p}$ is given. The rank minimization problem is NP-hard and computationally intractable ([10]). For low-rank matrix minimization problem, the rank function is often replaced with the nuclear norm of a matrix. As an elegant extension of the matrix case, many scholars have attempted to identify a possible convex relaxation as an alternative to the rank minimization problem. For example, Liu et al. proposed the sum of metricized nuclear norm of a tensor:

$\min _{x} \sum_{i=1}^{N} \alpha_{i}\left\|X_{(i)}\right\|_{*}$ s.t. $\mathbb{A}(\mathcal{X})=\boldsymbol{b}$  (17)

where, α_i≥0 (i=1, ..., N) are constants satisfying $\sum_{i=1}^{N} \alpha_{i}=1$; $\|Y\|_{*}=\sum_{i=1}^{n} \sigma_{i}(Y)$ is the nuclear norm of matrix Y; σ_i(Y) is the i-th largest singular value of matrix Y.

In previous work, the model (17) has successfully imputed missing tensor data. Recent studies suggest that the results can be significantly improved, using certain nonconvex functions [35-37].

2.3 DC programming and DC algorithm (DCA)

This subsection briefly outlies DC programming and DCA, facilitating subsequent model solving by DC programming.

Definition 2.4 [38] Let C be a convex subset of Rⁿ. A real-valued function f: C→R is called DC on C, if there exist two convex functions g,h: C→R such that f can be expressed as:

$f(x)=g(x)-h(x)$  (18)

Definition 2.5 [38] For an arbitrary function g:Rⁿ$\cup$R{+∞}, the function g^*:Rⁿ→R∪{+∞}, defined by

$g^{*}(y):=\sup \left\{\langle x, y\rangle-g(x): x \in R^{n}\right\}$     (19)

is called the conjugate function of g.

Let domg={x $\epsilon$Rⁿ:g(x)<∞} denote the effective domain of g. Suppose x₀ $\epsilon$domg, the subdifferential of g at x₀, denoted by ∂g(x₀), can be expressed as: ∂g(x₀):={y $\epsilon$Rⁿ: g(x)≥g(x₀)+〈x-x₀,y〉,∀x $\epsilon$Rⁿ}.

As a closed convex set in Rⁿ, the subdifferential ∂g(x₀) generalizes the derivative in the sense that g is diffrerntiable at x₀, if and only if ∂g(x₀) is reduced to a singleton which is exactly {g'(x₀)}.

DC programming is a programming problem dealing with DC functions. The DCA is a continuous approach based on local optimality and the duality of DC programming. The DCA aims to solve the DC program:

$\beta_{p}=\inf \left\{f(x)=g(x)-h(x): x \in R^{n}\right\}$  (20)

where, g and h are lower semi-continuous proper convex functions on Rⁿ. The dual program of formula (20) can be described as:

$\beta_{d}=\inf \left\{h^{*}(y)-g^{*}(y): y \in R^{n}\right\}$  (21)

Two sequences {x_k} and {y_k} are optimal solutions of primal and dual programs, respectively. The DCA consists of the two sequences such that the sequences {g(x^k)-h(x^k)} and {h^*(y^k)-g^*(y^k)} are decreasing. The two sequences are generated as follows: x^k+1 (resp.y^k) is a solution to the convex program (22) (resp. (23)) defined by

$\inf \left\{g(x)-h\left(x^{k}\right)-\left\langle x-x^{k}, y^{k}\right\rangle: x \in R^{n}\right\}$  (22)

$\inf \left\{h^{*}(y)-g^{*}\left(y^{k-1}\right)-\left\langle y-y^{k-1}, x^{k}\right\rangle: y \in R^{n}\right\}$  (23)

Problem (22) is obtained from primal DC program (20) by replacing h with its affine minimization defined by h_k(x):=h(x^k)+〈x-x^k, y^k〉 at a neighborhood of x^k. The solution set of problem (22) is ∂g^*(y^k). Similarly, problem (23) is obtained from the dual DC program (21) through the affine minimization of g^* defined by (g^*)_k(y):=g^*(y^k)+〈y-y^k,x^k+1〉 at a neighborhood of $y^{k}$, and the solution set of problem (23) is ∂h(x^k+1). Then, the DCA yields the next scheme:

$y^{k} \in \partial h\left(x^{k}\right) ; x^{k+1} \in \partial g^{*}\left(y^{k}\right)$   (24)

The complete introduction of DC programming and DCA is provided in [38-41].

To verify its empirical performance, the proposed Log-TC algorithm was compared with fixed point iterative method for low-rank tensor completion (FP-LRTC) [44] and fast low rank tensor completion algorithm (FaLRTC) [10] through numerical experiments. The FP-LRTC, which combines a fixed-point iterative method for solving the low n-rank tensor pursuit problem, and a continuation technique, has been applied to solve the relaxation model of the low-rank tensor completion problem. The FaLRTC utilizes a smoothing scheme to transform the original nonsmooth problem into a smooth problem, and was found more efficient than simple low rank tensor completion (SiLRTC), and high accuracy low rank tensor completion (HaLRTC).

This section is separated into two subsections. In Subsection 4.1, our algorithm was tested on the randomly generated low n-rank tensors. In Subsection 4.2, our algorithm was tested on color image recovery problem. All experiments were carried out in MATLAB R2016a on a computer running on Windows 10, with a 3.20GHz CPU and 4GB of memory. For fair comparison, each number is the average over 10 separate runs.

4.1 Synthetic data

First, our Log-TC algorithm was tested on several synthetic datasets for the tensor completion tasks. Following the method proposed by Rauhut et al. [45], a n-order tensor $\mathcal{J} \in R^{n_{1} \times n_{2} \times \cdots \times n_{N}}$ was created randomly with the n-rank (r₁, r₂, ..., r_N) of the form $\mathcal{T}:=\mathcal{M} \times_{1} U_{1} \times_{2} \cdots \times_{N} U_{N}$, where $\mathcal{M} \in R^{r_{1} \times r_{2} \times \cdots \times r_{N}}$ is the core tensor and $U_{i} \in R^{n_{i} \times r_{i}}(i=1,2, \cdots, N)$. Each entry of the core tensor $\mathcal{M}$ was sampled independently from a standard Gaussian distribution $\mathcal{N}(0,1)$. Then, m entries were sampled randomly from tensor $\mathcal{T}$ in a uniform manner. The sample ratio, denoted by SR, was defined as $\mathrm{SR}=m / I_{1} \cdots I_{N}$. The relative error (RelErr) of the recovered tensor defined as:

$\operatorname{Rel} \mathrm{Err}=\frac{\left\|x_{u l t}-\mathcal{T}\right\|_{F}}{\|\mathcal{T}\|_{F}}$  (62)

is used to estimate the closeness of $x_{u l t}$ to $\mathcal{T}$, where $x_{u l t}$ is the optimal solution to model (26) produced by the algorithms.

The ratio “fr” [44] between the degree of freedom (DOF) in a rank r matrix to the number of samples in the matrix completion was extended to tensor completion problem and denoted by FR:

$\mathrm{FR}=\frac{1}{N} \sum_{i=1}^{N} f r_{i}=\frac{1}{N} \sum_{i=1}^{N} \frac{r_{i}\left(I_{i}+\Pi_{i}-r_{i}\right)}{I_{1} \cdots I_{N} \times S R}$  (63)

where, $r_{i}=\operatorname{rank}\left(X_{(i)}\right)$; $\Pi_{i}=\prod_{k=1, k \neq i}^{N} I_{k}$. In the following results, “FR” is used to represent the hardness of the tensor completion. If FR is greater than 0.6, then the problems are “hard problems”; otherwise, the problems are “easy problems”. The parameters wee configured as: μ=1, θ_μ=1-sr, $\bar{\mu}=$$1 \times 10^{-8}$, τ=10, and λ=1.

The simulated performance of FP-LRTC, FaLRTC and Log-TC for random tensor completion problem, as measured by relative error and running time, are reported in Figure 2 and Table 1.

As shown in Figure 2, the relative error for tensor of size 60×60×60 with the n-rank was fixed at (6,6,6), and the sample rate changed from 10% to 80%. As it is known, when the n-rank is fixed, the complexity of the problem increases with the decrease of the SR. The plot illustrates that Log-TC is clearly superior to FP-LRTC and FaLRTC. The superiority is particularly evidence in the situation of SR=10%: Log-TC recovered the tensor with the relative error about 10^-4, while FaLRTC and FP-LRTC with 10⁰.

2.jpg

Figure 2. The relative error on tensors of size 60×60×60 with sample rate (SR) between 10% to 80% and fixed n-rank (6,6,6)

Table 1. Comparison of running time on tensors of size 60×60×60 with sample rate (SR) 70% and -rank (6,6,6)

Table 1 shows the performance of each algorithm with SR 70% in terms of the running time. Log-TC obviously outperformed the contrastive algorithms, under the same relative error. Even when the precision requirement is ultrahigh (10^-9), Log-TC can complete the trask in one minute.

Figure 3 shows the relative error for tensor of size 60×60×60, when the sample rate was fixed at SR=40% and n-rank r changed from 2 to 40. Like the previous situation, when sample rate is fixed, the complexity of the problem increases with the increase of the rank. As shown in Figure 3, Log-TC outperformed FP-LRTC and FaLRTC in all cases. Log-TC performed especially well at r=20: Log-TC recovered the tensors with the relative error about 10^-8, while FP-LRTC and FaLRTC with only about 10^-1. It can also be seen that, when r=24, our algorithm solved the problem with the relative error 10^-3, but the other two methods failed to effectively recover the tensors.

3.jpg

Figure 3. The relative error on tensors of size 60×60×60 n-rank (r,r,r) between 2 and 40 and fixed sample rate (SR=40%).

Next, the numerical results of the three algorithms on easy problems and hard problems were reported. Tables 2 and 3 list the numerical results of several cases in which the test tensors have different sizes, n-ranks and sample rates.

According to the results on easy problems (Table 2), the Log-TC outshined the other two algorithms in relative error and running time. FaLRTC is inferior to FP-LRTC and Log-TC in terms of relative error. FP-LRTC and Log-TC performed similarly on relative error, but Log-TC consumed the shorter time. Overall, Log-TC achieved the best performance on easy problems.

According to the results on hard problems (Table 3), Log-TC also realized better performance than the contrastive algorithms. With the increase of tensor order, the robustness of FaLRTC was getting worse. FP-LRTC and Log-TC were always capable of recovering the tensors efficiently, but the running time of FP-LRTC is about 2.5 times more than that of the Log-TC.

Table 2. Comparison results for easy problems

Table 3. Comparison results for hard problems

4.2 Image simulation

This subsection compared the three algorithms on the recovery of color images. Color images can be expressed as third order tensors. If the images are of low rank, the image recovery can be solved as a low n-rank tensor completion problem. Assuming that the images are well structured, FP-LRTC, FaLRTC and Log-TC were applied separately to solve the recovery problem.

The recovery effects of the Log-TC are illustrated in Figure 4 and Table 4. In Figure 4, (a1), (b1), (c1), and (d1) are original color images. Then, 80%, 60%, 40% and 20% entries were randomly removed, creating corrupted images (a2), (b2), (c2), and (d2), respectively. The images recovered by the Log-TC are (a3), (b3), (c3), and (d4), respectively. Only the effects of the Log-TC were provided, because the three algorithms achieved similar recovery accuracy. The recovery accuracy and running time of each algorithm are listed in Table 4.

4.jpg

Figure 4. Original images, corrupted images, and recovered image by Log-TC in the recovery of natural images

Table 4. Numerical results on the color images

SR	Log-TC		FP-LRTC		FaLRTC
	Relative Error	Time (s)	Relative Error	Time (s)	Relative Error	Time (s)
20%	9.14e-02	1300	9.62e-02	1681	8.31e-02	2137
40%	6.36e-02	598	8.45e-02	740	7.74e-02	1402
60%	4.89e-02	362	5.71e-02	298	2.57e-02	978
80%	3.40e-02	209	4.21e-02	302	3.91e-02	721

 

As shown in Figure 4, the Log-TC recovered the natural images with missing information. Table 4 shows that our algorithm achieved a recovery accuracy of about 10^-2, and a clear edge over the contrastive methods in running time.

5.jpg

Figure 5. Original images, corrupted images, recovered images by the Log-TC in the recovery of images with poor texture

Table 5. Numerical results on the color image with poor texture

SR	Log-TC		FP-LRTC		FaLRTC
	Relative Error	Time(s)	Relative Error	Time(s)	Relative Error	Time(s)
20%	1.17e-01	1835	5.72e-01	1920	4.53e-01	2341
40%	6.71e-02	857	8.45e-02	1103	7.47e-02	1780
60%	5.02e-02	493	2.71e-02	631	3.87e-02	1285
80%	3.37e-02	291	4.42e-02	329	2.57e-02	943

 

The original images in Figure 4 have a nice mixture of details, shading areas, flat regions, and textures. Next, the three algorithms were separately adopted to recover color images with poor texture. The recovery effect of our algorithm is displayed in Figure 5, and the numerical results of the three methods are compared in Table 5. The subgraphs in the middle column of Figure 5 were prepared in the same method as those in the middle column of Figure 4. It can be seen that, despite the poor texture of images, our algorithm still recovered the details of the images, and outperformed the other algorithms in running time.