Fast Singular Value Shrinkage with Chebyshev Polynomial Approximation Based on Signal Sparsity

TL;DR

This paper introduces a Chebyshev polynomial approximation method to efficiently perform singular value thresholding, significantly reducing computational costs in high-dimensional matrix rank minimization tasks like image processing.

Contribution

It presents a novel CPA-based approach that approximates singular value manipulation without explicit SVD, leveraging signal sparsity for faster computations.

Findings

Reduces computation time in matrix rank minimization tasks.

Maintains high approximation accuracy in image processing applications.

Demonstrates effectiveness through experimental results.

Abstract

We propose an approximation method for thresholding of singular values using Chebyshev polynomial approximation (CPA). Many signal processing problems require iterative application of singular value decomposition (SVD) for minimizing the rank of a given data matrix with other cost functions and/or constraints, which is called matrix rank minimization. In matrix rank minimization, singular values of a matrix are shrunk by hard-thresholding, soft-thresholding, or weighted soft-thresholding. However, the computational cost of SVD is generally too expensive to handle high dimensional signals such as images; hence, in this case, matrix rank minimization requires enormous computation time. In this paper, we leverage CPA to (approximately) manipulate singular values without computing singular values and vectors. The thresholding of singular values is expressed by a multiplication of certain matrices, which is derived from a characteristic of CPA. The multiplication is also efficiently computed using the sparsity of signals. As a result, the computational cost is significantly reduced. Experimental results suggest the effectiveness of our method through several image processing applications based on matrix rank minimization with nuclear norm relaxation in terms of computation time and approximation precision.

Tables

Table 1. TABLE I: Total and Average Computation Time and RMSEs (Image Inpainting Shown in Fig.  6 )

Bricks
Approximation order	5	10	15	20	EVD-based method
Total computation time (s)	401.85	373.91	347.53	348.38	464.81
Average computation time of singular value soft-shrinkage (s)	1.35	1.64	1.81	2.00	2.48
RMSE between our method and SVD-based method ($\times$10-3)	9.92	9.83	9.77	9.73	-

Table 2. TABLE II: Total and Average Computation Time and RMSEs (Background Modeling Shown in Fig.  7 )

Approximation order	5	10	15	20	EVD-based method
Total computation time (s)	1200.21	1127.60	1150.27	2192.24	6457.53
Average computation time of singular value soft-shrinkage (s)	20.95	21.21	21.79	22.20	30.77
RMSE between our method and EVD-based method ($\times$10-3)	56.27	23.89	10.86	3.71	-

Table 3. TABLE III: Total and Average Computation Time and RMSEs of Proposed and Existing Methods.

Proposed method (see Section V-E for the explanation)
Used transformation methods of CPA-based method	DWT	Block DCT	DCT
Total computation time (s)	347.53	344.03	380.23
Average computation time of singular value soft-shrinkage (s)	1.81	1.06	1.33
RMSE between our method and SVD-based method ($\times$10-3)	9.77	3.82	3.81

Table 4. TABLE IV: Comparisons of DCT with Adaptive and Fixed k 𝑘 k . Recall that k 1 = ( 0.5 × 10 − 4 ) ​ n 2 subscript 𝑘 1 0.5 superscript 10 4 superscript 𝑛 2 k_{1}=(0.5\!\times\!10^{-4})n^{2} , k 2 = ( 1.5 × 10 − 4 ) ​ n 2 subscript 𝑘 2 1.5 superscript 10 4 superscript 𝑛 2 k_{2}=(1.5\!\times\!10^{-4})n^{2} , k 3 = ( 0.5 × 10 − 3 ) ​ n 2 subscript 𝑘 3 0.5 superscript 10 3 superscript 𝑛 2 k_{3}=(0.5\!\times\!10^{-3})n^{2} , and k 4 = n 2 subscript 𝑘 4 superscript 𝑛 2 k_{4}=n^{2} .

Methods	DCT with adaptive $k$	DCT with fixed $k$
Used $k$ for $𝒦(\left|\underset{¯}{𝚽}\right|,k)$	Defined in (35)	$k_1$	$k_2$	$k_3$	$k_4$
Total computation time (s)	358.16	Not converged	416.52	380.23	765.84
Average computation time of singular value soft-shrinkage (s)	1.26	1.25	1.26	1.33	5.92
RMSE between our method and SVD-based method ($\times$10-3)	3.81	-	3.82	3.81	3.81

Equations

$$\hat{h}(x):=\frac{1}{2}{c}_0+\sum _{k=1}^{\alpha -1}{c}_k{\psi }_k(x),$$

$$\psi_{k}(x):=\cos\bigl{(}k\arccos(x)\bigr{)}.$$

$$\begin{array}{rl}& {\psi }_k(x)=2x{\psi }_{k-1}(x)-{\psi }_{k-2}(x),\\ & {\psi }_0(x)=1,{\psi }_1(x)=x.\end{array}$$

$$c_k:=\frac{2}{\alpha }\sum _{l=1}^{\alpha }\cos (k\theta (l))h(\cos \theta (l)),$$

$$\mathcal{H}(\mathbf{A}):=\mathbf{P}\,\mathrm{diag}\bigl{(}h(\lambda^{\mathrm{A}}_{1}),\ldots,h(\lambda^{\mathrm{A}}_{n})\bigr{)}\mathbf{P}^{-1},$$

$$\hat{\mathcal{H}}(\mathbf{A}):=\frac{1}{2}{\hat{c}}_0\mathbf{Id}+\sum _{k=1}^{\alpha -1}{\hat{c}}_k{\Psi }_k(\hat{\mathbf{A}}),$$

$$\hat{\mathbf{A}}:=\frac{2}{{\lambda }_{\max }^{\mathrm{A}}}\mathbf{A}-\mathbf{Id},$$

$${\Psi }_k(\hat{\mathbf{A}})$$

$$=\mathbf{P}{\Psi }_k\left(\frac{2}{{\lambda }_{\max }^{\mathrm{A}}}{\mathbf{\Lambda }}_{\mathrm{A}}-\mathbf{Id}\right){\mathbf{P}}^{-1}$$

$$=\mathbf{P}\mathrm{diag}(\cos k{\theta }_1,\dots ,\cos k{\theta }_n){\mathbf{P}}^{-1},$$

$$\begin{array}{rl}& {\Psi }_k(\hat{\mathbf{A}})=2\hat{\mathbf{A}}{\Psi }_{k-1}(\hat{\mathbf{A}})-{\Psi }_{k-2}(\hat{\mathbf{A}}),\\ & {\Psi }_0(\hat{\mathbf{A}})=\mathbf{Id},{\Psi }_1(\hat{\mathbf{A}})=\hat{\mathbf{A}}.\end{array}$$

$$\widehat{c}_{k}=\frac{2}{\alpha}\sum^{\alpha}_{l=1}\cos(k\theta(l))\,h\!\hskip 1.13809pt\Bigl{(}\frac{\lambda^{\mathrm{A}}_{\mathrm{max}}}{2}(\cos\theta(l)+1)\Bigr{)}.$$

$$\widehat{\mathcal{H}}(\mathbf{A})=\mathbf{P}\,\mathrm{diag}\bigl{(}\widehat{h}(\lambda^{\mathrm{A}}_{1}),\ldots,\widehat{h}(\lambda^{\mathrm{A}}_{n})\bigr{)}\mathbf{P}^{-1}.$$

$$\mathbf{\Sigma }=\left[\begin{array}{ccc}{\sigma }_1& & \mathbf{O}\\ & \ddots & \\ & & {\sigma }_n\\ \mathbf{O}& & \end{array}\right].$$

$$\mathcal{G}(\mathbf{B}):=\mathbf{U}\left[\begin{array}{ccc}g({\sigma }_1)& & \mathbf{O}\\ & \ddots & \\ & & g({\sigma }_n)\\ \mathbf{O}& & \end{array}\right]{\mathbf{V}}^{\top },$$

$$\mathcal{G}(\mathbf{B})=\mathbf{B}\mathcal{H}({\mathbf{B}}^{\top }\mathbf{B}),$$

$$\mathcal{G}(\mathbf{B})$$

$$=\mathbf{U\Sigma }{\mathbf{V}}^{\top }\mathbf{V}\mathrm{diag}\left(\frac{g({\sigma }_1)}{{\sigma }_1},\dots ,\frac{g({\sigma }_n)}{{\sigma }_n}\right){\mathbf{V}}^{\top }$$

$$=\mathbf{BV}\mathrm{diag}\left(\frac{g({\sigma }_1)}{{\sigma }_1},\dots ,\frac{g({\sigma }_n)}{{\sigma }_n}\right){\mathbf{V}}^{\top }.$$

$$\mathcal{G}(\mathbf{B})$$

$$\displaystyle=\mathbf{B}\mathbf{V}\,\mathrm{diag}\bigl{(}h(\sigma^{2}_{1}),\ldots,h(\sigma^{2}_{n})\bigr{)}\mathbf{V}^{\top}$$

$$=\mathbf{B}\mathcal{H}({\mathbf{B}}^{\top }\mathbf{B}).$$

$$\hat{\mathbf{x}}=\mathbf{B}\mathcal{H}({\mathbf{B}}^{\top }\mathbf{B})\mathbf{x},$$

$$\mathcal{G}(\mathbf{B})=\mathbf{B}{\mathbf{T}}^{\top }\mathcal{H}(\mathbf{T}{\mathbf{B}}^{\top }\mathbf{B}{\mathbf{T}}^{\top })\mathbf{T}=\mathbf{B}{\mathbf{T}}^{\top }\mathcal{H}(\mathbf{\Phi })\mathbf{T},$$

$$\widehat{\mathcal{H}}(\mathbf{\Phi})=\mathbf{T}\mathbf{V}\,\mathrm{diag}\bigl{(}\widehat{h}(\sigma^{2}_{1}),\ldots,\widehat{h}(\sigma^{2}_{n})\bigr{)}\mathbf{V}^{\top}\mathbf{T}^{\top}.$$

$$\underline{{\Phi }_{ij}}=\{\begin{array}{ll}{\Phi }_{ij}& \text{if}|{\Phi }_{ij}|\ge \epsilon ,\\ 0& \text{otherwise},\end{array}$$

$$\hat{\mathcal{H}}(\mathbf{\Phi })\approx \hat{\mathcal{H}}\left(\underline{\mathbf{\Phi }}\right).$$

$$\mathcal{G}(\mathbf{B})$$

$$=\mathbf{B}\mathcal{H}({\mathbf{B}}^{\top }\mathbf{B})$$

$$=\mathbf{B}{\mathbf{T}}^{\top }\mathcal{H}(\mathbf{\Phi })\mathbf{T}$$

$$\approx \mathbf{B}{\mathbf{T}}^{\top }\hat{\mathcal{H}}(\mathbf{\Phi })\mathbf{T}$$

$$\approx \mathbf{B}{\mathbf{T}}^{\top }\hat{\mathcal{H}}\left(\underline{\mathbf{\Phi }}\right)\mathbf{T}=\hat{\mathcal{G}}(\mathbf{B}),$$

$$\begin{array}{rl}& {\Psi }_k(\hat{\mathbf{B}})=2\hat{\mathbf{B}}{\Psi }_{k-1}(\hat{\mathbf{B}})-{\Psi }_{k-2}(\hat{\mathbf{B}}),\\ & {\Psi }_0(\hat{\mathbf{B}})=\mathbf{Id},{\Psi }_1(\hat{\mathbf{B}})=\hat{\mathbf{B}},\end{array}$$

$$\hat{\mathbf{B}}:=\frac{2}{{\Lambda }_{\max }}\underline{\mathbf{\Phi }}-\mathbf{Id},$$

$$h_{\mathrm{hard}}(x):=\{\begin{array}{ll}1& \text{if}x>{\tau }_{\mathrm{hard}},\\ 0& \text{otherwise},\end{array}$$

Full text

I Introduction

The low-rank structure inherent in various signals has been widely exploited in many signal processing applications, such as matrix and tensor completion [1, 2, 3, 4], image decomposition [5, 6], photometric stereo [7, 8], image alignment [9, 10], colorization [11], inpainting [12, 13], background modeling [14, 15, 16, 17, 18], color artifact removal [19], cognitive radio [20], and voice separation [21]. In such applications, the low-rank structure is incorporated into a minimization problem involving the rank function or its continuous relaxation. The problem is solved using some iterative algorithms, with which the thresholding of singular values is usually required at each iteration. We refer to this methodology as matrix rank minimization.

There are two representative approaches of matrix rank minimization. One is the exact method. It is an ideal formulation, but the resulting problem is very difficult to solve due to the non-convexity and combinatorial nature of the rank function. The other is the nuclear norm relaxation [22]. Since the nuclear norm, the sum of the singular values of a matrix, is the tightest convex relaxation of the rank function, we can efficiently solve the resulting problem via convex optimization techniques. Weighted nuclear norm relaxation [23, 24] has recently been proposed as a non-convex but continuous approximation of the rank function.

Essentially, both of the above methods require the thresholding of singular values, which we call singular value shrinkage, at each iteration of certain optimization methods (Fig. 1). That is, most methods for matrix rank minimization must carry out singular value decomposition (SVD) many times. This is a serious problem in terms of computational cost when we handle large matrices, even with high-spec computers.

Several methods have been proposed to tackle this issue [25, 26, 27]. The basic concept of [25, 26] is to approximately compute partial singular values and/or vectors. These methods can drastically reduce the computation time of singular value shrinkage but would not be suitable for the matrix rank minimization. Since the number of singular values above a threshold is not identified without the full decomposition, many singular values above a threshold are reduced to zero in each iteration. As a result, large approximation errors are produced, which results in an unstable convergence in the matrix rank minimization. With the other method [27], singular value shrinkage is carried out by computing neither singular values nor vectors, but the reduction in the computation time is still limited. This is because the method requires a complete orthogonal decomposition [28] and the calculation of the inverse of a large matrix. It also leads to large approximation errors, i.e., an unstable convergence in the matrix rank minimization. We consider a method similar to that by Cai and Osher [27]: We only need a “processed” matrix with thresholded singular values.

In this paper, we propose a fast singular value shrinkage method for reducing the computational cost in the matrix rank minimization of high dimensional matrices. Note that the proposed method computes neither singular values nor vectors during the process of singular value shrinkage, similar to the method by Cai and Osher [27]. Furthermore, our method maintains computational precision to lead matrix rank minimization algorithms to stable convergence. The two key tools of our method are described as follows.

• •

Chebyshev polynomial approximation (CPA) [29, 30, 31]: This tool is often used for designing filters in signal processing [32, 33] and is a key tool for reducing computational cost. The applications of CPA have been studied by Saad et al.[34, 35, 36, 37]. With the applications by Saad et al., CPA is used to calculate a vector after being transformed by a matrix with singular value shrinkage. That is, it requires the iterative multiplications of a matrix and vector to derive the Chebyshev polynomials. The concept of the applications has recently been used for improving the performance of image filtering methods such as bilateral filter, non-local means, and BM3D [38, 39, 40]. In contrast, we propose a method to obtain a matrix whose singular values are processed by using CPA. Since CPA results in truncation errors, such as ripples in the lower-order approximations, we also investigate the designs of thresholding functions and appropriate approximation order for reducing approximation errors.

• •

Sparsity of signals: By using CPA, our method can represent singular value shrinkage as a multiplication of matrices. Since the multiplication can be computed efficiently when the matrices are sparse, our method exploits the inherent sparsity of signals in their frequency domain for further acceleration.

Since matrix rank minimization plays a central role in various signal processing tasks, our method offers many promising applications. For this study, we validated the proposed method by using two image processing applications: image inpainting [12] and background modeling[14, 15, 16, 17, 18]. In these applications, target problems are formulated as convex optimization problems involving the nuclear norm so that they can be efficiently solved using the alternating direction method of multiplier (ADMM) [41] with our method.

Although the ADMM is widely known as a robust method for computation errors in each iteration, optimization methods (including the ADMM) with the other fast singular value shrinkage methods [25, 26, 27] do not converge well due to their large approximation errors. In contrast, our CPA-based singular value shrinkage method leads optimization methods to stable convergence. We validated this advantage experimentally by comparing our method with the other fast singular value shrinkage methods in several image processing tasks and a synthetic data.

The preliminary version of this study, without using signal sparsity, analysis of our method, and new applications, has previously been published [42].

This paper is organized as follows. Section II defines notations and preliminaries. We discuss our CPA-based singular value shrinkage method, which is the main contribution in this paper, in Section III. We discuss an approximation order of CPA for reducing the size of approximation errors in Section IV and verification of our method through applications in Section V. Finally, we conclude the paper in Section VI.

II Notations and Preliminaries

II-A Notations

Bold-face capital and small letters indicate a matrix and a vector, respectively. Superscript $\cdot^{\top}$ is the transpose of a matrix and a vector, and superscript $\cdot^{-1}$ is the inverse of a non-singular matrix. The matrices $\mathbf{Id}$ and $\mathbf{O}$ are the identity matrix and null matrix, respectively. The vector $\mathbbm{1}_{n}:=[\underbrace{1,\ldots,1}_{n}]^{\top}$. The $\ell_{p}$ norm for $p\geq 1$ is defined as $\|\mathbf{x}\|_{p}\!:=\!(\sum^{n}_{i=1}|x_{i}|^{p})^{\frac{1}{p}}~{}(\forall\mathbf{x}\!\in\!\mathbb{R}^{n})$. We also use CPA as follows.

II-B Chebyshev Polynomial Approximation

Let $h(x)$ and $\widehat{h}(x)$ be a real-valued function defined on the interval $x\!\in\![-1,1]$ and its approximated function by using CPA, respectively. Chebyshev polynomial approximation [29, 30, 31] gives an approximate solution of $h(x)$ by using the truncated Chebyshev series:

[TABLE]

where $c_{k}$ and $\alpha$ denote a Chebyshev coefficient (described later) and an approximation order, respectively. Additionally, $\psi_{k}(\cdot)$ denotes the $k$-th order Chebyshev polynomials of the first kind, defined as

[TABLE]

It can also be computed using the stable recurrence relation:

[TABLE]

The initial condition is defined as $\psi_{0}(x)$ and $\psi_{1}(x)$. Since the polynomials consist of cosine functions, the value of $\psi_{k}(x)$ is bounded between $-1$ and $1$ for $x\!\in\![-1,1]$. By using $\psi_{k}(x)$ and the orthogonality of the cosine function, $c_{k}$ is calculated as

[TABLE]

where $\displaystyle\theta(l)\!:=\!\frac{\pi(l-\frac{1}{2})}{\alpha}$.

III Singular Value Shrinkage using Chebyshev Polynomial Approximation by Exploiting Sparsity

We discuss singular value shrinkage using CPA. First, the CPA of a matrix form, which can approximately shrink the eigenvalues of a matrix (eigenvalue shrinkage), is indicated then extended to the singular one.

III-A Chebyshev Polynomial Approximation for Matrix

Let $\mathbf{A}\!\in\!\mathbb{R}^{n\times n}$ be a full rank matrix and $\mathbf{A}\!=\!\mathbf{P}\mathbf{\Lambda}_{\mathrm{A}}\mathbf{P}^{-1}$ be its eigendecomposition (EVD), where $\mathbf{P}\!\in\!\mathbb{R}^{n\times n}$ is the matrix composed of eigenvectors and $\mathbf{\Lambda}_{\mathrm{A}}\!=\!\mathrm{diag}(\lambda^{\mathrm{A}}_{1},\ldots,\lambda^{\mathrm{A}}_{i},\ldots,\lambda^{\mathrm{A}}_{n})$ is the diagonal matrix with the corresponding eigenvalues. We assume that the eigenvalues are bounded between [math] and $\lambda^{\mathrm{A}}_{\mathrm{max}}$, where $\lambda^{\mathrm{A}}_{\mathrm{max}}\!>\!1$. Hence, the eigenvalues of $\mathbf{A}$ are shrunk as

[TABLE]

where $\mathcal{H}(\cdot)$ is the eigenvalue shrinkage function, and $h(x)$ is the filter kernel defined in $x\!\in\![0,\lambda^{\mathrm{A}}_{\mathrm{max}}]$. In this subsection, we consider the approximated solution of (5) using the CPA.

The CPA of the matrix form [30, 34, 36, 37] gives an approximated solution of the eigenvalue shrinkage function $\mathcal{H}(\cdot)$ by using truncated Chebyshev series as

[TABLE]

where $\widehat{c}_{k}$ and $\mathit{\Psi}_{k}(\widehat{\mathbf{A}})$ are Chebyshev coefficients and Chebyshev polynomials, respectively, which are defined later. Additionally, $\widehat{\mathbf{A}}$ is the eigenvalue-shifted matrix given by

[TABLE]

whose eigenvalues are obviously within $[-1,1]$. Thanks to (7), the $k$-th order Chebyshev polynomial of $\widehat{\mathbf{A}}$ is computed as

[TABLE]

where $\theta_{i}:=\arccos\left(\frac{2}{\lambda^{\mathrm{A}}_{\mathrm{max}}}\lambda^{\mathrm{A}}_{i}-1\right)$. Similarly to (3), the Chebyshev polynomials are obtained using the recurrence relation:

[TABLE]

Recall that $\mathit{\Psi}_{k}(\widehat{\mathbf{A}})$ is defined only in the interval $[-1,1]$. Therefore, the range of the filter kernel is modified by deriving $\widehat{c}_{k}$ as

[TABLE]

The term $h\left(\lambda^{\mathrm{A}}_{\mathrm{max}}/2\,(\cos\theta(l)+1)\right)$ returns the shifted range back to the original range $[0,\lambda^{\mathrm{A}}_{\mathrm{max}}]$. From (10), $\widehat{\mathcal{H}}(\mathbf{A})$ can also be represented using $\widehat{h}(\lambda^{\mathrm{A}}_{i})$ as

[TABLE]

The function $\widehat{\mathcal{H}}(\cdot)$, which is referred to as the CPA-based eigenvalue shrinkage function, results in approximate eigenvalue shrinkage. The CPA-based eigenvalue shrinkage actually computes neither eigenvalues nor vectors thanks to the recurrence relation (9).

III-B CPA-based Singular Value Shrinkage

Let $\mathbf{B}\!\in\!\mathbb{R}^{m\times n}~{}(m\!>\!n)$ be a rectangular matrix and $\mathbf{B}\!=\!\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}$ be its singular value decomposition, where $\mathbf{U}\!\in\!\mathbb{R}^{m\times m}$ and $\mathbf{V}\!\in\!\mathbb{R}^{n\times n}$ are orthogonal matrices. The $\mathbf{\Sigma}\!\in\!\mathbb{R}^{m\times n}$ is the singular value matrix represented as

[TABLE]

Without loss of generality, we can assume $\sigma_{1}\geq\ldots\geq\sigma_{n}$. The singular values of $\mathbf{B}$ are shrunk with the singular value shrinkage function $\mathcal{G}(\cdot)$ as

[TABLE]

where $g(\cdot)$ is an arbitrary function.

The eigenvalue shrinkage in (5) can be extended to $\mathcal{G}(\mathbf{B})$ in (13) as [35]

[TABLE]

where $h(x)\!:=\!g(\sqrt{x})/\sqrt{x}$ in $\mathcal{H}(\cdot)$ in (5). Equation (14) is derived as follows. First, (13) can be expanded as

[TABLE]

When the eigenvalue matrix of $\mathbf{B}^{\!\top}\!\mathbf{B}$ is defined as $\mathbf{\Lambda}_{\mathrm{B}^{\!\top}\!\mathrm{B}}\!=\!\mathrm{diag}(\lambda^{\mathrm{B}^{\!\top}\!\mathrm{B}}_{1},\ldots,\lambda^{\mathrm{B}^{\!\top}\!\mathrm{B}}_{n})$, it is obviously represented using the singular values of $\mathbf{B}$ as $\lambda^{\mathrm{B}^{\!\top}\!\mathrm{B}}_{i}\!=\!\sigma^{2}_{i}$. Consequently, (15) is equally calculated using $\mathcal{H}(\cdot)$ in (5) as

[TABLE]

Note that [35] aims to calculate a vector represented as

[TABLE]

where $\mathbf{x}\in\mathbb{R}^{n}$ and $\widehat{\mathbf{x}}\in\mathbb{R}^{m}$ are the input and output vectors, respectively. The CPA is applied to $\mathcal{H}(\mathbf{B}^{\!\top}\!\mathbf{B})\mathbf{x}$ to quickly derive $\widehat{\mathbf{x}}$ in [35]. In contrast, our method is focused on deriving the matrix in (16) itself. When the matrix, whose singular values are shrunk using CPA, is represented as $\mathbf{B}\mathcal{H}(\mathbf{B}^{\!\top}\!\mathbf{B})$, the explicit SVD of $\mathbf{B}$ can be avoided. However, deriving the matrix, not the vectors, usually requires enormous computation time because of multiplication of dense matrices. To accelerate the calculation, the sparseness of a matrix is exploited with our method, as indicated below.

Assume that $\mathbf{B}$ is a matrix composed of an inherently sparse signal. Let $\mathbf{T}\!\in\!\mathbb{R}^{n\times n}$ be an arbitrary orthogonal matrix that efficiently sparsifies $\mathbf{B}$, e.g., $\mathbf{T}$ is considered as the discrete Fourier transform [43], discrete cosine transform (DCT) [44], and discrete wavelet transform (DWT) [45]. Note that the matrix $\mathbf{T}$ is the forward transform, i.e., when let $\mathbf{y}\in\mathbb{R}^{n}$ be a column vector, the forward transform is represented as $\mathbf{Ty}$. From the above, (14) is further rewritten with $\mathbf{T}$ as

[TABLE]

where $\mathbf{\Phi}:=\mathbf{T}\mathbf{B}^{\!\top}\!\mathbf{B}\mathbf{T}^{\top}$ for simplicity. With the CPA-based eigenvalue shrinkage function $\widehat{\mathcal{H}}(\cdot)$ in (6) and (11), $\mathcal{H}(\mathbf{\Phi})$ in (18) is efficiently approximated as

[TABLE]

The form of (19) enables us to use the sparsity of a signal in its frequency domain.

For further enhancing the sparsity of $\mathbf{\Phi}$, its components are thresholded as

[TABLE]

where $\Phi_{ij}$ and $\underline{\Phi_{ij}}$ are the $i$-th row and $j$-th column of $\mathbf{\Phi}$ and its truncated coefficient, and $\varepsilon\!\in\!\mathbb{R}$ is an arbitrary small value. We show that this truncation has little effect on the performance of our method and provides recommended settings of $\varepsilon$ in Section V. As a result, the CPA-based eigenvalue shrinkage of $\widehat{\mathcal{H}}(\mathbf{\Phi})$ is approximately given by

[TABLE]

In summary, the singular value shrinkage of $\mathbf{B}$ is approximately represented as

[TABLE]

where the function $\widehat{\mathcal{G}}(\cdot)$ is the CPA-based singular value shrinkage function. It can be calculated with the recurrence relation as

[TABLE]

where

[TABLE]

in which $\Lambda_{\mathrm{max}}\!=\!\lambda^{\underline{\mathbf{\Phi}}}_{\mathrm{max}}$ is the maximum eigenvalue of $\underline{\mathbf{\Phi}}$. The pseudocode of the CPA-based singular value shrinkage is indicated in Algorithm 1.

III-C Computational Complexity of CPA-based Singular Value Shrinkage

We now discuss the computational complexity of our method. Assume that matrices $\widehat{\mathbf{B}}\in\mathbb{R}^{n\times n}$ and $\mathit{\Psi}_{k}(\widehat{\mathbf{B}})\in\mathbb{R}^{n\times n}$ have $M$ and $M_{k}$ nonzero elements, respectively. The maximum number of multiplications of nonzero elements required to calculate $\widehat{\mathbf{B}}\mathit{\Psi}_{k}(\widehat{\mathbf{B}})$ is represented as $MM_{k}$ in the case of a sparse matrix. The computational complexity of line 9 in Algorithm 1 can be represented as $\mathcal{O}(\sum^{\alpha-2}_{k=1}MM_{k})$ due to the multiplication $\widehat{\mathbf{B}}\mathit{\Psi}_{k-1}(\widehat{\mathbf{B}})$. At line 10 in Algorithm 1, the computation takes $\mathcal{O}(\sum^{\alpha-1}_{k=2}M_{k})$ from the multiplication $c_{k}\mathit{\Psi}_{k}(\widehat{\mathbf{B}})$. That is, the total computational complexity is represented as $\mathcal{O}((M+1)\max_{k}\{M_{k}\})$, where $\max_{k}\{M_{k}\}$ represents the maximum value among $M_{k}$. From the above, when $\max_{k}\{M_{k}\}$ becomes small, the computational cost is also reduced.

For low computational complexity, the matrix should be constructed so as not to increase the number of its nonzero elements as much as possible in the multiplication of matrices.

IV Shrinkage Functions and Approximation Order

In this section, we discuss suitable approximation orders for shrinkage functions approximated by CPA, which has small truncation errors. Additionally, we argue that CPA is a reasonable choice for our method among a variety of polynomial approximation methods.

As an introduction, we consider the shrinkage function shown in Fig. 2(a). Let $h_{\mathrm{hard}}(x)$ be the hard shrinkage response defined as

[TABLE]

where $\tau_{\mathrm{hard}}$ is an arbitrary real value and $x\!\in\![0,1]$. Chebyshev polynomial approximation gives an approximate response of $h_{\mathrm{hard}}(x)$ in (1). As in Fig. 2(b), the approximated response $\widehat{h}_{\mathrm{hard}}(x)$ has ripples, which is widely known in digital filter design [46, 47, 32, 48, 49]. Therefore, studying the design of appropriate shrinkage responses and approximation orders is an important topic, even for our method.

IV-A Approximation Order

Possible shrinkage responses handled with our method can be expressed as the following generic form:

[TABLE]

where $w(x)$ is a weight function, and $\tau$ and $\rho$ are arbitrary thresholding values. The choices of $w(x)$ and $\tau$ determine the characteristics of (26) as follows:

• •

$\displaystyle h\left(x;0,\tau_{\mathrm{hard}}\right)$  : Hard-shrinkage.

• •

$\displaystyle h\left(x;\frac{w(x)}{\rho},\frac{w(x)}{\rho}\right)$  : Weighted soft-shrinkage.

• •

$\displaystyle h\left(x;\frac{1}{\rho},\frac{1}{\rho}\right)$  : Soft-shrinkage111Soft-shrinkage is widely known as $g(x)\!:=\!\max(x\!-\!1/\rho,0)$ in which $\max(x_{1},x_{2})$ is an operator choosing the greater one out of $x_{1}$ and $x_{2}$. However, we call $h(x;1/\rho,1/\rho)$ soft-shrinkage because it is finally transformed into $g(x)\!:=\!\max(x\!-\!1/\rho,0)$ in (22)..

Note that we defined $g(\sqrt{x})=\sqrt{x}h(x)$ in (14), where $h(x)$ is an arbitrary shrinkage function for the eigenvalues of $\mathbf{B}^{\!\top}\!\mathbf{B}$.

As a result, $g(x)$ becomes the hard-shrinkage, weighted soft-shrinkage, or soft-shrinkage functions when $h(x)$ is set as above.

These choices among the shrinkages are shown in Fig. 3. It is clear that hard-shrinkage has a sharp transition band (see Fig. 3(a)); therefore, CPA, which is computed as a linear combination of cosine functions, may not approximate it well. In contrast, one can expect that a response that has a smooth transition band is suitable for CPA. To verify this numerically, the approximated responses were compared among hard-shrinkage, weighted soft-shrinkage, and soft-shrinkage. In this experiment, $\rho\!=\!0.002$ was used. For weighted soft-shrinkage, $w(x)\!=\!0.5$ and $1.5$ were used. Additionally, the thresholding value for hard-shrinkage was set to $\tau_{\mathrm{hard}}\!=\!500$.

Figure 3 also shows the approximated shrinkage responses of (26) for various shrinkage conditions. Hard-shrinkage yields larger errors than soft ones. Empirically, hard-shrinkage requires more than the 50th-oder approximation. In contrast, soft-shrinkages only require 10–20th-order approximations. To be more specific, $\alpha\!=\!20$ is recommended for a small weight shown in Fig. 3(b), whereas $\alpha\!=\!10$ for the soft-shrinkage response shown in Fig. 3(d).

IV-B Suitability of CPA

There are many polynomial approximations. Even among them, minimax polynomial approximation [50, 46, 47, 51, 52] and least squares approximation [53] are well known as the best approximation in the sense of the minimization of the infinity norm and the least squares error w.r.t the difference between an exact and approximated responses, respectively. To derive polynomial coefficients, their optimization requires a minimization of $\ell_{p}$ norm represented as

[TABLE]

where $h^{\prime}(x)$ is an approximated shrinkage response with the above two polynomial approximations. Clearly, (27) requires the exact response $h(x)$ for $x\!\in\!\mathbb{R}$. When $h(x)$ is precisely represented using many sampling points, $h^{\prime}(x)$ exhibits good performance. However, computational complexity becomes high when many sampling points are used, especially in the case of least squares approximation. Let $\mathbf{c}_{\alpha}\!:=\![c_{0},c_{1},\ldots,c_{\alpha-1}]^{\top}$ be the column vector of coefficients for the polynomial approximation. That is, the approximated shrinkage response can be calculated as $h^{\prime}(x)\!=\!c_{0}+c_{1}x+c_{2}x^{2}+\ldots+c_{\alpha-1}x^{\alpha-1}$. Additionally, let $\mathbf{x}\!:=\![x_{1},\ldots,x_{n}]^{\top}$ and $\mathbf{h}\!:=\![h(x_{1}),\ldots,h(x_{n})]^{\top}$ be the column vectors composed of real values, respectively. The Vandermonde matrix $\mathbf{\Upsilon}\!\in\!\mathbb{R}^{n\times\alpha}$ is defined as

[TABLE]

From the above definitions, coefficients of the least squares approximation are calculated as $\mathbf{c}_{\alpha}\!=\!\mathbf{\Upsilon}^{+}\mathbf{h}$, where $\cdot^{+}$ is the pseudo inverse of a matrix. The calculation requires high computational cost when $n$ and/or $\alpha$ are large. In contrast, CPA only requires the inner product of $[h(\cos\theta_{1}),\ldots,h(\cos\theta_{\alpha})]^{\top}$ and $[\cos k\theta_{1},\ldots,\cos k\theta_{\alpha}]^{\top}$ to derive coefficients of polynomials from (4), where $\theta_{i}\!\in\![0,\pi]$. Additionally, CPA performs better approximation than other optimization methods. To verify the exellent approximation, CPA was compared with minimax approximation and least squares approximation, as shown in Fig. 4. As can be seen, CPA and least squares approximations have a similar oscillation pattern. Chebyshev polynomial approximation sufficiently attenuates ripples, compared with the other methods in the stopband, as shown in the differences of Fig. 4.

V Applications

We compared our CPA-based singular value shrinkage method with the exact and approximate singular value shrinkage methods. Specifically, we applied our method to two applications using nuclear norm relaxation, i.e., inpainting of texture images and background subtraction of videos. Additionally, we compared our CPA-based method with the existing methods, i.e., the exact partial singular value decomposition (PSVD) based method and fast singular value shrinkage methods [25, 26, 27], in Section V-F. The computation time and approximation precision were indicated for the comparisons.

V-A Experimental Conditions

The applications were implemented with MATLAB R2015b and run on a 3.2-GHz Intel Xeon E5-2667 processor with 512-GB RAM. We compared our method with the SVD-based naive method (denoted as SVD-based method) in (13) and EVD-based methods in (15) with respect to approximation precision and computation time. Both SVD and EVD-based methods are exact singular value soft-shrinkage methods. The EVD-based method222The EVD-based method could lead to loss of computational precision compared with the SVD-based one. Though the errors may affect the performance of applications, we did not encounter such a problem in the experiments described in this paper. is usually faster than the SVD-based method and is widely used in many applications. Therefore, the computation time of only the EVD-based method is indicated for the results of the exact methods. With the SVD-based method, the SVD of an arbitrary matrix $\mathbf{X}\!\in\!\mathbb{R}^{m\times n}$ is first performed, then the obtained singular values are shrunk as $\max(\sigma_{i}(\mathbf{X})\!-\!1/\rho,0)$, where $\sigma_{i}(\mathbf{X})$ indicates the $i$-th largest singular value of $\mathbf{X}$. The EVD-based method uses the relation between singular value shrinkage and eigenvalue shrinkage: the EVD of $\mathbf{X}^{\!\top}\!\mathbf{X}$ is first computed, then the obtained eigenvalues are shrunk as \max\bigl{(}\sqrt{\lambda^{\mathrm{X}^{\!\top}\!\mathrm{X}}_{i}}-1/\rho,0\bigr{)}/\sqrt{\lambda^{\mathrm{X}^{\!\top}\!\mathrm{X}}_{i}}, where $\lambda_{i}^{\mathbf{X}^{\!\top}\!\mathbf{X}}$ denotes the $i$-th largest eigenvalue of $\mathbf{X}^{\!\top}\!\mathbf{X}$, to derive singular value shrinkage. Also, the shrinkage function $h(\cdot)$ in (14) is defined as the soft-shrinkage case given by $h\left(x;\frac{1}{\rho},\frac{1}{\rho}\right)$ from (26) for our method. The DWT [45] was used in (22) to sparsify the signals. We used Haar wavelet transform as the DWT. In the DWT, one level transform was performed and all high frequency components were set to [math]. The selection of a transform method naturally affects the computation time of our method. Therefore, we indicate the effect of the selection in Section V-E. To indicate the approximation precision, root mean squared error (RMSE) was used, which was computed using the results of our method and those of the SVD/EVD-based methods. Furthermore, the computation times of all the methods are shown, and the average computation times of the CPA-based/exact singular value shrinkage in each iteration are also indicated. In all applications, we used the 5th, 10th, 15th, and 20th-order approximations. We also used the following optimization tools to solve the above applications.

V-B Optimization Tools

V-B1 Proximity Operator

Let $\boldsymbol{\Gamma}_{0}(\mathbb{R}^{N})$ be the set of all proper lower semicontinuous convex functions333A function $f:\mathbb{R}^{N}\!\rightarrow\mathbb{R}\cup\{\infty\}$ is called proper lower semicontinuous convex if $\mbox{dom}(f)\!:=\!\{\mathbf{x}\!\in\!\mathbb{R}^{N}|~{}f(\mathbf{x})\!<\!\infty\}\!\neq\!\emptyset$, $\mbox{lev}_{\leq a}(f)\!:=\!\{\mathbf{x}\!\in\!\mathbb{R}^{N}|~{}f(\mathbf{x})\!\leq\!a\}$ is closed in $\forall a\!\in\!\mathbb{R}$, and $f(\eta\mathbf{x}+(1-\eta)\mathbf{y})\!\leq\!\eta f(\mathbf{x})+(1-\eta)f(\mathbf{y})$ in $\forall\mathbf{x},\mathbf{y}\!\in\!\mathbb{R}^{N}$ and $\forall\eta\!\in\!(0,1)$, respectively. over $\mathbb{R}^{N}$. The proximity operator [54] of a function $f\!\in\!\boldsymbol{\Gamma}_{0}(\mathbb{R}^{N})$ of index $\gamma\!>\!0$ is defined as

[TABLE]

The proximity operator plays a central role in the optimization of applications, as discussed in this section. When function $f$ is defined as the nuclear norm, i.e., $\mbox{prox}_{\gamma\|\cdot\|_{*}}$, the proximity operator can be calculated by singular value shrinkage with the thresholding parameter $\gamma$ [2]. Therefore, our CPA-based method is applied to the operator in the case of the nuclear norm.

V-B2 Alternating Direction Method of Multipliers

The ADMM [41] is an algorithm for solving a convex optimization problem represented as

[TABLE]

where $f\!\in\!\boldsymbol{\Gamma}_{0}(\mathbb{R}^{n_{1}})$, $g\!\in\!\boldsymbol{\Gamma}_{0}(\mathbb{R}^{n_{2}})$ and $\mathbf{K}\!\in\!\mathbb{R}^{n_{2}\times n_{1}}$. For arbitrary $\mathbf{z}_{0}$, $\mathbf{p}_{0}\!\in\!\mathbb{R}^{n_{2}}$, and $\rho\!>\!0$, the ADMM algorithm is given by

[TABLE]

We recall a convergence analysis of the ADMM by Eskstein and Bertsekas [41].

Fact 1** (Convergence of the ADMM [41])**

Consider Prob. (30). Assume that $\mathbf{K}^{\!\top}\!\mathbf{K}$ is invertible and that a saddle point of its unaugmented Lagrangian $L_{0}(\mathbf{x},\mathbf{z},\mathbf{u}^{\prime})\!:=\!f(\mathbf{x})+g(\mathbf{z})-\langle\mathbf{u}^{\prime},\mathbf{K}\mathbf{x}\!-\!\mathbf{z}\rangle$ exists, where $\mathbf{u}^{\prime}\!:=\!\rho\mathbf{u}$. Then the sequence $(\mathbf{x}_{t})_{(t\geq 1)}$ generated using (31) converges to a solution of Prob. (30).

We used the ADMM algorithm to practically solve the following applications. In all applications, the stopping criterion444For example, in (39), which is indicated in Appendix A-A, the criterion is evaluated using $\|\mathbf{l}_{t+1}-\mathbf{l}_{t}\|_{2}/\|\mathbf{l}_{t+1}\|_{2}$. in the ADMM algorithm was set to $1.0\!\times\!10^{-4}$.

V-C Texture Image Inpainting [12, 13]

The objective with this application is to recover a missing region (as shown in the later Fig. 6(b)).

Let $\mathbf{L}$ and $\mathbf{I}\!\in\!\mathbb{R}^{m\times n}$ be a texture image and a given image with missing regions, respectively. Then, let $\Omega$ and $\overline{\Omega}$ be observed and missing regions and $P_{\Omega}(\cdot)$ and $P_{\overline{\Omega}}(\cdot)$ be linear operators extracting pixels in their regions. From the notations, the missing region is represented as $\mathbf{M}\!=\!P_{\overline{\Omega}}(\mathbf{L})$. The pixels surrounding $\mathbf{M}$ with the size of five pixels, as shown in Fig. 5, are defined as $\mathbf{M}_{\partial}$. Let $\mathbf{T}_{1}\!\in\!\mathbb{R}^{m\times m}$ and $\mathbf{T}_{2}\!\in\!\mathbb{R}^{n\times n}$ be the DCT matrices in the horizontal and vertical matrix directions, i.e., these matrices transform an image to its frequency domain. Since a regular texture image is basically sparse in its frequency domain, it can be represented as $\mathbf{L}\!=\!\mathbf{T}_{1}\mathbf{S}\mathbf{T}_{2}^{\top}$, where $\mathbf{S}$ is the coefficients on the frequency domain of $\mathbf{L}$. Additionally, the set of a normalized dynamic range constraint is defined as $\mathcal{D}\!:=\!\{\mathbf{x}\!:=\![x_{i}]_{i=1}^{mn}|~{}x_{i}\!\in\![0,1]\}$. When $\mathbf{L}$ and $\mathbf{S}$ are assumed to be low rank and sparse, the reconstruction problem can approximately be solved using the nuclear norm555The nuclear norm of $\mathbf{X}\!\in\!\mathbb{R}^{m\times n}$ is defined as $\|\mathbf{X}\|_{*}:=\sum^{K}_{i=1}\sigma_{i}(\mathbf{X})$, where $i\!\in\!\{1,2,\ldots,K\}(K:=\min(m,n))$. and the $\ell_{1}$ norm as

[TABLE]

where a positive real value $\eta$ is a regularization parameter, $\mathrm{ave}(\cdot)$ calculates the arithmetic average, and $\mathrm{vec}(\cdot)$ is the operator vectorizing a matrix. The average pixel value on the recovered region is assumed to be identical to that of its surrounding pixel values. As can be seen, (32) is composed of nuclear norm relaxation so that our method can be used for its efficient calculation. Hereafter, we discuss the validation of our method by applying the ADMM to (32) to obtain the optimal solution. In Appendix A-A, (32) is converted to the form to which the ADMM is applicable.

Eight-bit color images Bricks and Office widows666The images are available at http://www.mayang.com/textures/., as shown in Figs. 6(a) and (f), were used for the application, where the size of each color component of the images was 1920 $\times$ 2560. The pixel values of each color component of the images were in the range from 0 to 1. In the application, each color component was inpainted separately. The observed image $\mathbf{I}\!\in\!\mathbb{R}^{2560\times 1920}$ with the missing regions were defined, as shown in Fig. 6(b) and (g)777In the experiment, the images are transposed to “portrait”.. The number of missing pixels in Fig. 6(b) was $\mathbf{M}\!\in\!\mathbb{R}^{300\times 300}$ and that in Fig. 6(g) was $\mathbf{M}\!\in\!\mathbb{R}^{180\times 1250}$. In the ADMM applicable form (see Appendix A-A), the column vectors $\mathbf{l}_{0}$ and $\mathbf{u}_{0}$ were initialized by all-one vectors. Additionally, for $\mathrm{prox}_{1/\rho\|\cdot\|_{*}}$ and $\mathrm{prox}_{\eta/\rho\|\cdot\|_{1}}$, the thresholding parameters $(1/\rho,\eta/\rho)$ were set to $(6,0.1)$ in Bricks and $(5,0.1)$ in Office windows, where the parameters were determined for the fast and stable convergence of the optimization. For fast computation, parallel processing888The MATLAB function parfor, which is contained in the parallel computing toolbox, was used for the parallel computing only in the image inpainting method. was performed in the application of the color components.

Figure 6 shows the results of image inpainting with the 20th-order approximation. The resulting image recovered using our method was practically equivalent to that with the SVD-based method by comparing Figs. 6(c), (d), (h), and (i). In Figs. 6(e) and (j), it is clear that the exact and approximated solutions had little differences, which visually indicates the high approximation precision of our method.

Table I lists the computation time and RMSE comparisons. Our method was faster than the EVD-based method while maintaining reconstruction performance. Regarding the total computation times of Bricks and Office windows, our method with the 5th-order and 10th-order approximations was slower than with the 15th-order approximations. This is because our method with the low-order approximations did not converge well due to the low approximate precision.

V-D Background Modeling of Video[14, 15, 16, 17, 18]

The objective with this application is to divide a video sequence into background and object sequences (as shown in Fig. 7(a) and (b)).

Let $\mathbf{I}^{(i)}\!\in\!\mathbb{R}^{m\times n}$ be the $i$-th frame of a video sequence. The sequence is rearranged into a matrix $\mathbf{I}\!\in\!\mathbb{R}^{mn\times K}$ as

[TABLE]

Then, let $\mathbf{L}$ and $\mathbf{S}\!\in\!\mathbb{R}^{mn\times K}$ be the background sequence and sequence of moving objects of a video. In $\mathbf{L}$, pixel values corresponding to $\mathbf{S}$ are zero and vice versa. The background and moving objects can be assumed to be low rank and sparse; hence, the background modeling is solved as the following convex optimization problem:

[TABLE]

Problem (34) is also solved using the ADMM.

In the background modeling, an eight-bit grayscale video Laboratory999This video was recorded with our video camera. It was downsampled and transformed into grayscale for the experiment. was used. The $\mathbf{I}^{(i)}\!\in\!\mathbb{R}^{360\times 240}$ was the $i$-th frame of the video in $i\!\in\!\{1,2,\ldots,5000\}$; hence, the matrix of the sequence was $\mathbf{I}\!\in\!\mathbb{R}^{86400\times 5000}$. The pixel values of the video were in the range from 0 to 1. In the ADMM applicable form (see Appendix A-B), the column vectors $\mathbf{l}_{0}$, $\mathbf{s}_{0}$, and $\mathbf{u}_{0}$ were initialized by all-one vectors.

Additionally, the thresholding parameters $(1/\rho,\eta/\rho)$ were set to $(480,0.12)$, where the parameters were determined for the fast and stable convergence of the optimization. For comparison, low rank and sparse components in the 160th frame are shown in Fig. 7 with the 20th-order approximation.

Our method effectively decomposed the video sequences to low rank and sparse sequences, as shown in Figs. 7(a) and (b). They are almost equivalent to those with the EVD-based method; hence, the difference between the low rank images are not displayed even though the difference is amplified (bottom left of Fig. 7).

Table II summarizes the computation times and RMSEs between the background modeling of our method and that of the EVD-based method101010The SVD-based method ran out of memory in our machine so that the results of the EVD-based method were used for the RMSEs.. Our method was sufficiently faster than the EVD-based method in all approximation orders.

This is because the reduction in computational complexity by thresholding the transformed coefficients described in (21) is effective for our method to have low computational complexity while retaining high approximation precision. However, our method with the 5th-order approximations took more time than that with 10th-order approximation because it did not converge well due to its low approximate precision.

V-E Effects of Selections: Transform Matrix and Thresholding Value

The selections of a transform matrix $\mathbf{T}$ in (22) and a thresholding value $\varepsilon$ in (20) affect computation time and size of approximation error. In this subsection, we indicate these effects experimentally by using the image inpainting method for Bricks. In all experiments, the 15th-order approximation was used for our method.

V-E1 Effects of Selected Transform Matrix

We compared the DWT with the DCT and the block diagonal forms of the DCT (block DCT) [32] whose block size was $8\times 8$ for indicating the differences among chosen transform matrices $\mathbf{T}$ in (22). The threshold values in (20) were fixed to $\varepsilon\!=\!250$ for the DCTs. The other experimental conditions were the same as those discussed in Section V-C.

The results of the proposed method in Table III show the performance comparisons. Our method with the DWT is as fast as that with the block DCT in the total computation time, though the singular value shrinkage with our method with the DWT takes more time than the others. This is because the maximum iteration of our method with the DWT is only 83, whereas those of our method with the block DCT and DCT are 101 and 97, respectively. Therefore, the DWT leads our method to be stable convergence. However, our method with the DWT indicates a higher RMSE than the others since all the high frequency components were removed in it. The fact is not fatal problem because the images derived by using our method with the DWT were very similar to that of the exact one as shown in the previous section. Additionally, the block DCT can substantially sparsify an image compared to the DWT and DCT. Therefore, our method with the block DCT is faster than the others. In spite of fast computation, the results with the block DCT shows an RMSE as low as that of the DCT.

V-E2 Effects of Thresholding Value

When $\varepsilon$ is an excessively large value, our method becomes fast but matrix rank minimization cannot converge well due to the errors w.r.t. the reduction in the number of components. For achieving the fast computation and stable convergence, we present a recommended guideline on the thresholding values.

Let $E_{t}$ be an error at the $t$-th iteration of the ADMM, i.e., $E_{t}\!=\!\|\mathbf{l}_{t}-\mathbf{l}_{t-1}\|/\|\mathbf{l}_{t}\|$ from Appendix A-A and $\mathcal{K}(\mathbf{X},k)$ be a function that returns the $k$-th largest element in a matrix $\mathbf{X}$, where $\mathcal{K}(\mathbf{X},k)$ is used as the threshold. Basically, when $E_{t}$ is a large value, a small $k$ does not have any problem to decrease the error. In contrast, when $E_{t}$ is a small value, i.e., the optimization almost converges, $k$ should be large for a stable convergence. For this purpose, we recommend the thresholding percentage for $\mathcal{K}(|\underline{\mathbf{\Phi}}|,k)$ as

[TABLE]

where $E_{\ell}\!:=\!0.3\!\times\!10^{-1}$ and $E_{\mathrm{m}}\!:=\!6\!\times\!10^{-4}$. This was determined experimentally. Recall that the size of $\underline{\mathbf{\Phi}}$ is $n\times n$. To verify (35), it was compared with fixed threshold. Four values $k_{1}=(0.5\!\times\!10^{-4})n^{2}$, $k_{2}=(1.5\!\times\!10^{-4})n^{2}$, $k_{3}=(0.5\!\times\!10^{-3})n^{2}$, and $k_{4}=n^{2}$ were used for this comparison, where $k_{4}$ means that all components of a matrix are retained. Additionally, the DCT was exploited for the sparsifying matrix.

The results are listed in Table IV. The adaptive method was faster than the fixed method. Our method with $k_{1}$ was the fastest but it did not converge well for the large size of approximation errors, where $k_{1}$ is considered as the limitation of the thresholding percentage in our method. In addition, the singular value shrinkage of our CPA-based method by thresholding the matrix components was about five times faster than that by maintaining those.

V-F Comparison with Existing Methods

As previously mentioned, there are several fast singular value shrinkage methods [25, 26, 27]. To illustrate the advantage of our method, we compared it with the singular value shrinkage by using the exact PSVD, fast randomized singular value shrinkage (FRSVS) [25], singular value shrinkage by using the Nyström method (NSVS) [26], and the fast singular value shrinkage without the exact SVD (FSVS) [27]. The experiments were conducted using the image inpainting method for Bricks and a synthetic data. The synthetic data is constructed as a block diagonal matrix whose number of the main diagonal blocks is equal to its rank. Let $\mathbf{D}\in\mathbb{R}^{1000\times 1000}$ be the matrix form of the synthetic data with rank $n$, and this is defined as $\mathbf{D}:=\mathbf{J}-\mathrm{blkdiag}(\mathbf{D}_{\mathrm{s}},n)$, where $\mathbf{J}:=\mathbbm{1}_{1000}\mathbbm{1}^{\top}_{1000}$, $\mathbf{D}_{\mathrm{s}}:=0.5\mathbbm{1}_{1000/n}\mathbbm{1}^{\top}_{1000/n}$ and $\mathrm{blkdiag}(\mathbf{D}_{\mathrm{s}},n):=\mathrm{diag}(\underbrace{\mathbf{D}_{\mathrm{s}},\mathbf{D}_{\mathrm{s}},\ldots,\mathbf{D}_{\mathrm{s}}}_{n})$. In the experiment, $n=$10, 100, 200, and 500 were used, and $x\%$ of elements in $\mathbf{D}$ were randomly replaced with zero, where $x\in\{1,10,20\}$. To restore the corrupted data, the optimization problem in (32) was solved without using the average term, i.e., \mathrm{ave}\bigl{(}\mathrm{vec}(\mathbf{M})\bigr{)}=\mathrm{ave}\bigl{(}\mathrm{vec}(\mathbf{M}_{\partial})\bigr{)}, because the matrices $\mathbf{M}$ and $\mathbf{M}_{\partial}$ can hardly be defined for the random corruption. The function svdsechon111111Available at https://www.mathworks.com/matlabcentral/fileexchange/47132-fast-svd-and-pca was used for carrying out the exact PSVD, and the 500 largest singular values were calculated for Bricks. All preferences of the FRSVS [25] was determined in the original code121212Available at http://thohkaistackr.wixsite.com/page/projectfrsvt provided by the authors. In the FRSVS, the 550 largest singular values were approximately derived for Bricks, whose number of singular values was experimentally determined for carrying out precise and fast singular value shrinkage. In addition, the 200 largest singular values were derived for the synthetic data in the exact PSVD and the FRSVS. The partial singular values derived in the exact PSVD and the FRSVS were soft-thresholded: The $i$th partial singular value $\sigma^{\mathrm{p}}_{i}$ is shrunk to $\max(\sigma^{\mathrm{p}}_{i}-1/\rho,0)$. Since the Nyström method requires a square matrix to derive partial eigenvalues, it was applied to $\mathbf{X}^{\!\top}\!\mathbf{X}$, as in the EVD-based method. The 500 and 200 largest eigenvalues were calculated for the Bricks and the synthetic data, respectively. Those eigenvalues were then shrunk in the same way as with the EVD-based method, whose number of calculated eigenvalues was experimentally determined from the same reason as the FRSVS. All preferences used in the FSVS were directly used as suggested in [27].

The results for Bricks are indicated in Table III and Fig. 8. Note that the experiments of the existing methods were stopped at the 80th iterations since these methods did not converge. The concept of the FSVS is similar to our method, but it requires longer computation time as shown in Table III and does not converge. For the FRSVS and NSVS, although their average computation times are slightly less than ours, they result in much larger errors. This would be because many singular values or eigenvalues above the threshold $1/\rho$ were reduced to zero, so that the exact PSVD, FRSVS, and NSVS produced the large errors in each iteration leading to unstable convergence. In contrast, our method is stable and does not affect the convergence of the optimization method because the CPA-based method can shrink the entire singular values.

Figure 9 shows the comparison of errors in the case of the synthetic data in several conditions according to the reduction rate of data elements and the matrix rank. The optimization methods using the existing methods do not converge well when the matrix rank is 500. This is because many singular values above $1/\rho$ are discarded, i.e., enormous computation errors are produced in each iteration. From the results, the matrix rank of target data should be estimated beforehand, and then the numbers of partial singular values and vectors should be estimated to be larger than the matrix rank, in order to make the optimization method converged. In contrast, our method can lead the optimization method to stable convergence. It certainly generates some approximation errors, but it can process all singular values, which means that most singular values above $1/\rho$ are remained.

VI Conclusion

We proposed a fast thresholding method of singular values without computing singular values and vectors. The key tool of the proposed method is CPA. From CPA characteristics, singular value shrinkage could be computed by a multiplication of matrices. The proposed method was further accelerated using the sparsity of a signal, where the frequency transform was used for obtaining sparse coefficients. Moreover, we studied the approximation order for reducing the size of approximation errors. The experimental results revealed that our method was much faster than the exact methods with high approximation precision in the case of a large data size. In addition, our method can lead the optimization method to be stable convergence in comparison of the existing fast singular value shrinkage methods because of its approximation precision.

Appendix A ADMM Applicable Forms

A-A Texture Image Inpainting

Let $\mathbf{i}\!:=\!\mathrm{vec}(\mathbf{I})$, $\mathbf{l}\!:=\!\mathrm{vec}(\mathbf{L})$, $\mathbf{m}\!:=\!\mathrm{vec}(\mathbf{M})$, and $\mathbf{m}_{\partial}\!:=\!\mathrm{vec}(\mathbf{M}_{\partial})$. The 2-D DCT matrix is represented as $\mathbf{\Psi}$, and the matrix form of $P_{\Omega}$ and $P_{\overline{\Omega}}$ are defined as $\mathbf{\Omega}$ and $\overline{\mathbf{\Omega}}$. In addition, the indicator functions of the sets $\mathcal{I}\!:=\!\{\mathbf{x}\!\in\!\mathbb{R}^{mn}|~{}\mathbf{x}\!=\!\mathbf{\Omega}\mathbf{i}\}$, $\mathcal{M}\!:=\!\{\mathbf{x}\!\in\!\mathbb{R}^{N}\ |\ \mathrm{ave}(\mathbf{x})\!=\!\mathrm{ave}(\mathbf{m}_{\partial})\}$ and $\mathcal{D}$ are denoted as $\iota_{\mathcal{I}}$, $\iota_{\mathcal{M}}$, and $\iota_{\mathcal{D}}$, respectively, where $N$ is the size of $\mathrm{vec}(\mathbf{M})$. By using the above definitions, (32) is redefined as

[TABLE]

Let the vector $\mathbf{z}\!\in\!\mathbb{R}^{5mn}$ be

[TABLE]

Finally, (36) is represented as

[TABLE]

Equation (38) can be applied to the ADMM algorithm in (31) which is indicated in Appendix V-B. Let $\mathbf{u}_{0}\!:=\![(\mathbf{u}^{(1)}_{0})^{\top},(\mathbf{u}^{(2)}_{0})^{\top},(\mathbf{u}^{(3)}_{0})^{\top},(\mathbf{u}^{(4)}_{0})^{\top},(\mathbf{u}^{(5)}_{0})^{\top}]^{\top}$ be an arbitrary auxiliary vector, where $\mathbf{u}^{(i)}_{0}\!\in\!\mathbb{R}^{mn}$ in $i\!=\!1,\ldots,5$. Applying the ADMM to (38) yields the following algorithm:

[TABLE]

where the update of $\mathrm{prox}_{1/\rho\|\cdot\|_{*}}$ in ($\ast$) of (39) can be computed by singular value shrinkage, which is performed using our CPA-based method. Additionally, the update of $\mathrm{prox}_{\eta/\rho\|\cdot\|_{1}}$ is approximated by soft-thresholding, i.e., $\mathrm{sgn}(X_{ij})\max(|X_{ij}|\!-\!\eta/\rho,0)$, where $X_{ij}$ is the entry of an arbitrary matrix $\mathbf{X}$ and $\mathrm{sgn}(\cdot)$ is the sign function. In the following applications, the same calculation is used for the updates of $\mathrm{prox}_{1/\rho\|\cdot\|_{*}}$ and $\mathrm{prox}_{\eta/\rho\|\cdot\|_{1}}$. In (39), $\Pi_{\mathcal{I}}(\cdot)$, $\Pi_{\mathcal{D}}(\cdot)$, and $\Pi_{\mathcal{M}}(\cdot)$ are the metric projections onto $\mathcal{I}$, $\mathcal{D}$, and $\mathcal{M}$, respectively. Practically, $\Pi_{\mathcal{I}}(\cdot)$ is given by maintaining the assigned pixels, and $\Pi_{\mathcal{D}}(\cdot)$ is calculated by pushing the entries outside $[0,1]$ into [math] or $1$ (the nearest is chosen). Additionally, the auxiliary value is calculated as the difference between the average value of $\mathbf{m}_{\partial}$ on $\mathbf{i}$ and the average value of $\overline{\mathbf{\Omega}}\mathbf{l}_{t+1}+\mathbf{u}_{t}^{(5)}$ on the recovered region. The $\Pi_{\mathcal{M}}(\cdot)$ is derived by adding the auxiliary value to $\overline{\mathbf{\Omega}}\mathbf{l}_{t+1}+\mathbf{u}_{t}^{(5)}$ on the recovered region.

A-B Background Modeling

Let $\mathbf{i}\!:=\!\mathrm{vec}(\mathbf{I})$, $\mathbf{l}\!:=\!\mathrm{vec}(\mathbf{L})$, and $\mathbf{s}\!:=\!\mathrm{vec}(\mathbf{S})$. The indicator function of the set $\mathcal{I}\!:=\!\{\mathbf{x}\!\in\!\mathbb{R}^{mnK}|~{}\mathbf{x}\!=\!\mathbf{i}\}$ is defined as $\iota_{\mathcal{I}}$. By using the above definitions, (34) is rewritten as

[TABLE]

When an auxiliary vector $\mathbf{z}$ is represented as

[TABLE]

Problem (40) is further rewritten as

[TABLE]

Let $\mathbf{u}_{0}\!:=\![(\mathbf{u}^{(1)}_{0})^{\top},(\mathbf{u}^{(2)}_{0})^{\top},(\mathbf{u}^{(3)}_{0})^{\top}]^{\top}$ be an auxiliary vector for the ADMM, where $\mathbf{u}^{(i)}_{0}\!\in\!\mathbb{R}^{mnK}$ in $i\!=\!1,2,3$. Applying the ADMM to (42) yields the following algorithm:

[TABLE]

where the update of $\mathbf{z}^{(1)}$ in ($\ast$) of (43) is calculated using our CPA-based method. In (43), $\Pi_{\mathcal{I}}(\cdot)$ is the metric projection onto $\mathcal{I}$, which is given by maintaining the observed pixel values of the original sequences.