Adaptive Singular Value Thresholding

TL;DR

This paper introduces an adaptive thresholding method for low-rank matrix recovery that dynamically adjusts the threshold during iterations, leading to improved performance over traditional fixed-threshold approaches.

Contribution

The novel adaptive singular value thresholding (ASVT) method adjusts thresholds iteratively, enhancing low-rank recovery under affine constraints compared to existing fixed-threshold methods.

Findings

Better recovery performance demonstrated in simulations

Adaptive threshold decreases during iterations

Outperforms fixed-threshold methods

Abstract

In this paper, we propose an Adaptive Singular Value Thresholding (ASVT) for low rank recovery under affine constraints. Unlike previous iterative methods that the threshold level is independent of the iteration number, in our proposed method, the threshold in adaptively decreases during iterations. The simulation results reveal that we get better performance with this thresholding strategy.

Figures

Tables

Table 1. TABLE I: Comparison of the proposed algorithm and the SVT algorithm.

			ASVT		SVT [12]
Size	Rank	$\frac{m}{n_1{n}_2}$	$\#$iters	RE	$\#$iters	RE
$500\times 500$	10	0.15	75	9.60e-4	200	5.02e-3
	50	0.4	64	9.62e-4	200	4.38e-2
	100	0.5	138	9.64e-4	200	3.28e-1
$1000\times 1000$	10	0.15	38	9.22e-4	62	1.49e-3
	50	0.4	29	9.21e-4	70	1.60e-3
	100	0.5	37	9.88e-4	103	1.89e-3
$1500\times 1500$	10	0.15	23	8.17e-4	49	1.28e-3
	50	0.4	22	8.69e-4	53	1.32e-3
	100	0.5	26	8.79e-4	68	1.51e-3
$2000\times 2000$	10	0.1	49	9.24e-4	54	1.31e-3
	50	0.3	27	8.05e-4	56	1.35e-3
	100	0.4	28	9.98e-4	68	1.53e-3
$2500\times 2500$	10	0.1	38	9.60e-4	47	1.65e-3
	50	0.25	33	8.36e-4	56	1.37e-3
	100	0.3	37	8.96e-4	77	1.61e-3
$3000\times 3000$	10	0.05	29	9.29e-4	65	1.53e-3
	50	0.25	29	9.33e-4	50	1.33e-3
	100	0.3	32	8.37e-4	66	1.51e-3

Equations

$$\begin{array}{rl}\underset{\mathbf{X}}{\min }& rank(\mathbf{X}),\\ s.t.:& \mathcal{A}(\mathbf{X})=\mathbf{b},\end{array}$$

$$\mathbf{A}vec(\mathbf{X})=\mathbf{v}$$

$$\begin{array}{rl}\underset{\mathbf{X}}{\min }& \Vert \sigma (\mathbf{X}){\Vert }_0,\\ s.t.:& \mathcal{A}(\mathbf{X})=\mathbf{b},\end{array}$$

$$\begin{array}{rl}\underset{\mathbf{X}}{\min }& rank(\mathbf{X}),\\ s.t.:& [\mathbf{X}{]}_{i,j}=[{\mathbf{X}}_{\ast }{]}_{i,j},\forall (i,j)\in \Omega ,\end{array}$$

$$\begin{array}{rl}\underset{\mathbf{X}}{\min }& \Vert \sigma (\mathbf{X}){\Vert }_1,\\ s.t.:& \mathcal{A}(\mathbf{X})=\mathbf{b}.\end{array}$$

$$\left\{\begin{array}[]{l}\mathbf{X}_{k}=\mathcal{D}(\mathbf{Y}_{k-1},\tau)\\ {}\\ \mathbf{Y}_{k}=\mathbf{Y}_{k-1}+\delta_{k}\mathbf{\mathcal{A}}(\mathbf{X_{*}}-\mathbf{X}_{k})\end{array}\right.$$

$$\mathcal{D}(\mathbf{A},\tau )=\mathbf{U}diag(\{\max ({\sigma }_i-\tau ,0){\}}_{i=1}^{\min (n_1,n_2)}){\mathbf{V}}^T.$$

$$f_{\delta }(x)=\exp \left(-\frac{x^2}{2{\delta }^2}\right).$$

$$rank(\mathbf{X})\approx min(n_{1},n_{2})-\sum_{i=1}^{min(n_{1},n_{2})}f_{\delta}\big{(}\sigma_{i}(\mathbf{X})\big{)}$$

$${\tau }_k=b\exp (-ak),$$

$$\left\{\begin{array}[]{l}\mathbf{X}_{k}=\mathcal{T}(\mathbf{Y}_{k-1},\tau_{k})\\ {}\\ \mathbf{Y}_{k}=\mathbf{Y}_{k-1}+\delta_{k}\mathbf{\mathcal{A}}^{*}(\mathbf{b}-\mathbf{\mathcal{A}}\mathbf{X}_{k})\end{array}\right.$$

$$\mathcal{T}(\mathbf{X},\tau_{k})\buildrel\Delta\over{=}\mathbf{U}\big{[}\mbox{if~{}}\sigma_{i,j}<\tau_{k}:0,\mbox{~{}else:~{}}\sigma_{i,j}\big{]}_{i=1,j=1}^{n_{1},n_{2}}\mathbf{V}^{T}.$$

$$RE(\mathbf{X},\hat{\mathbf{X}})=\frac{\Vert \hat{\mathbf{X}}-\mathbf{X}{\Vert }_F}{\Vert \mathbf{X}{\Vert }_F}$$

Full text

I Introduction

Finding a minimum rank matrix subject to affine constraints is the heart of many engineering applications such as minimal realization theory [1], minimum order controller design [2, 3, 4], collaborative filtering [5], and machine learning [6, 7, 8]. This problem is usually known as Affine Rank Minimization (ARM). Mathematically speaking, we want to solve the following minimization problem:

[TABLE]

where $\mathbf{X}\in\mathbb{R}^{n_{1}\times n_{2}}$ is the unknown matrix, $\mathbf{\mathcal{A}}:\mathbb{R}^{n_{1}\times n_{2}}\to\mathbb{R}^{m}$ denotes the linear mapping, and $\mathbf{b}$ is the measurement vector. The affine constraint $\mathbf{\mathcal{A}}(\mathbf{X})=\mathbf{b}$ can be simplified as

[TABLE]

where $\mathbf{A}\in\mathbb{R}^{m\times n_{1}n_{2}}$ is the matrix representation of $\mathbf{\mathcal{A}}$ and $vec(.)$ is the standard vectorization operator.

Since rank is the number of non-zero Singular Values (SVs), the minimization problem (1) is the same as

[TABLE]

where $\mathbf{(\sigma)}(\mathbf{X})$ is the vector of SVs of $\mathbf{X}$.

An important special case of ARM problem introduced in (1) is the matrix completion problem:

[TABLE]

where $\mathbf{X}$ is the unknown matrix, $[\mathbf{X}]_{i,j}$ is the element of matrix $\mathbf{X}$ at the intersection of $i$-th row and $j$-th column, and $\mathbf{X_{*}}$ is the low rank matrix that we would like to recover but we only know a small subset of all entries of $\mathbf{X_{*}}$, that is, $\Omega\subseteq\{1,2,\cdots,n_{1}\}\times\{1,2,\cdots,n_{2}\}$.

There are about two main approaches to solve the ARM problem or better saying approximate the non-convex minimization problem (2).

In the first approach, the rank function is approximated by the nuclear norm. The nuclear norm of a matrix is defined as $\|\mathbf{X}\|_{*}=\sum_{i=1}^{r}\sigma_{i}(\mathbf{X})$ where $\sigma_{i}$ is the $i$-th largest SV of $\mathbf{X}$. Therefore, we can write the nuclear norm as $\|\mathbf{X}\|_{*}=\|\sigma(\mathbf{X})\|_{1}$ and finally, the nuclear norm minimization will be equivalent to an $\ell_{1}$-norm minimization as follows [9, 10]:

[TABLE]

Above nuclear norm minimization can be efficiently solved using some convex optimization techniques such as a semi-definite program [11]. For the purpose of matrix completion, a Singular Value Thresholding (SVT) algorithm is proposed to solve the nuclear norm minimization [12]. It is an iterative algorithm that for a fixed value of $\tau>0$ and positive step-sizes $\{\delta_{k}\}$, starts with $\mathbf{Y}_{0}=\mathbf{0}$ and repeats

[TABLE]

until a stopping criterion is satisfied. $\mathcal{D}(.,.)$ is called the singular value shrinkage operator. For a matrix $\mathbf{A}\in\mathbb{R}^{n_{1}\times n_{2}}$ with Singular Value Decomposition (SVD) of $\mathbf{A}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{T}$, we have

[TABLE]

In the second approach, the $\ell_{0}$-norm used in the rank function is approximated by families of smoothing $\ell_{0}$-norm functions [13]. For example, a class of Gaussian smoothing functions are defined as:

[TABLE]

As $\delta\to 0$, $f_{\delta}(.)$ gives better approximation of the Kronecker delta function. In this approach, the rank of $\mathbf{X}\in\mathbb{R}^{n_{1}\times n_{2}}$ is approximated as follows:

[TABLE]

and the final optimization problem is solved using the Gradient Projection (GP) algorithm [14].

In this paper, we consider the ARM problem (1) and propose an adaptive singular value thresholding method. Here, the threshold level, unlike the threshold level in the SVT algorithm is not fixed.

II The Proposed Algorithm

We aim to recover a low rank matrix under some linear constraints. As stated in Section I, this problem is known as ARM. It should be noted that we don’t know the rank of the primary matrix. Since the low rank matrix has sparsity in SV domain, we propose to threshold the SVs. But the main difference between our algorithm and previous ones is that the threshold level is adaptively changed during iterations. The idea of adaptive thresholding is used in [15] for recovery of the original signal from its random samples. As suggested in [15], we choose threshold level of $k$-th iteration as follows:

[TABLE]

where $a$ and $b$ are two constants. The proposed algorithm is as follows:

[TABLE]

where $\mathbf{\mathcal{A}}^{*}(.)$ is the adjoint of $\mathbf{\mathcal{A}}(.)$. For a matrix $\mathbf{X}\in\mathbb{R}^{n_{1}\times n_{2}}$ with SVD $\mathbf{X}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{T}$, $\mathcal{T}(\mathbf{X},\tau_{k})$ is defined as follows:

[TABLE]

We call our algorithm Adaptive Singular Value Thresholding (ASVT) which is shown in Algorithm 1.

III Simulation Results

The numerical experiment results are presented in this section. All simulations are done by MATLAB R2015a on Intel(R) Core(TM) i7-5960X @ 3GHz with 32GB-RAM.

We consider matrix completion as a special case of linear constraints. The random matrices are generated according to the scheme proposed in [12]. An $r$-rank matrix $\mathbf{X}$ of size $n_{1}\times n_{2}$ is generated as $\mathbf{X}=\mathbf{A}\mathbf{B}$ where $\mathbf{A}\in\mathbb{R}^{n_{1}\times r}$ and $\mathbf{B}\in\mathbb{R}^{r\times n_{2}}$ are sampled independently form a standard Gaussian distribution ($\mathcal{N}(0,1)$). The measurements or the subset of observed entries $\Omega$ is sampled uniformly at random of size $m$. As [12], we choose constant step-sizes as $\delta=1$.

To evaluate the performance of the algorithms, we report the relative error of the reconstruction defined as follows:

[TABLE]

where $\mathbf{X}$ and $\mathbf{\hat{X}}$ are original and reconstructed matrices, respectively.

III-A Comparison with the SVT Algorithm

We generate matrices of size $n_{1}=n_{2}=100,300,500,1000,\mbox{and},2000$ with different ranks and different fraction of observations $\left({\frac{m}{n_{1}n_{2}}}\right)$ and provide a comparison between our algorithm with the SVT algorithm proposed in [12]. We have downloaded the SVT MATLAB codes from http://svt.stanford.edu/code.html.

According to the results of Table I, our proposed method has better performance in terms of relative error even with a fewer number of iterations. For example, when the size of the matrix is $3000\times 3000$ and its rank is 100 and only 0.3 of its entries are observed, our algorithm recovers it after 29 iterations while the SVT algorithm does this after 65 iterations. Moreover, the relative error of our algorithm and the SVT algorithm is $9.29\times 10^{-4}$ and $1.53\times 10^{-3}$, respectively.

III-B Effect of Parameters

In this Subsection, we investigate the effect of damping factor in thresholding operator ($\alpha$) and step size ($\delta$). For this simulation, we generate $300\times 300$ random matrices of rank $10$ where only $1/3$ of their entries are observed. We test the proposed method with different values of $\alpha$ and $\delta$. The simulation results are shown in Figs. 1 and 2. According to the results of Figs. 1 and 2, the run time of the algorithm decrease as $\alpha$ and $\delta$ increases.

III-C Phase Transition Plot

To better show the performance of the proposed method, we plot the phase transition between $\frac{m}{n_{1}n_{2}}$ and $\frac{d_{r}}{m}$ where $d_{r}=r(n_{1}+n_{2}-r)$ is the degree of freedom.

In this simulation, we generate random matrix $\mathbf{M}$ of size $80\times 80$. Each experiment is repeated 100 times and a matrix is declared to be successfully recovered if the relative error between the recovered matrix and the original one is less than $10^{-3}$. The phase transition is charted in Fig. 3. The gray color of each cell indicates the probability of success recovery during simulations. For example, the white color indicates that the proposed algorithm can successfully recover the original matrix while the black indicates that it fails.

IV Conclusion

In this paper, we proposed an adaptive singular value thresholding for recovery of low rank matrices under affine constraints. The word adaptive means that the threshold level is adaptively changed during the iterations. Recently, the SVT algorithm is proposed for the matrix completion which is an important case of low rank matrix recovery under affine constraints. We also compared our algorithm with the SVT algorithm and the simulation results show that we can achieve improvement with adaptive thresholding. We show that the relative error of the proposed method is less than the relative error of the SVT algorithm. Moreover, our algorithm does this task in a fewer number of iterations.