A multi-stage convex relaxation approach to noisy structured low-rank matrix recovery

TL;DR

This paper introduces a multi-stage convex relaxation method for noisy structured low-rank matrix recovery, providing theoretical guarantees and demonstrating its effectiveness through numerical experiments.

Contribution

It proposes a novel multi-stage convex relaxation approach based on an exact penalty formulation, with proven convergence and error reduction guarantees.

Findings

The method achieves geometric convergence of the error sequence.

The approach outperforms standard nuclear norm relaxation in experiments.

Theoretical analysis confirms error reduction and rank approximation bounds.

Abstract

This paper concerns with a noisy structured low-rank matrix recovery problem which can be modeled as a structured rank minimization problem. We reformulate this problem as a mathematical program with a generalized complementarity constraint (MPGCC), and show that its penalty version, yielded by moving the generalized complementarity constraint to the objective, has the same global optimal solution set as the MPGCC does whenever the penalty parameter is over a threshold. Then, by solving the exact penalty problem in an alternating way, we obtain a multi-stage convex relaxation approach. We provide theoretical guarantees for our approach under a mild restricted eigenvalue condition, by quantifying the reduction of the error and approximate rank bounds of the first stage convex relaxation (which is exactly the nuclear norm relaxation) in the subsequent stages and establishing the geometric convergence of the error sequence in a statistical sense. Numerical experiments are conducted for some structured low-rank matrix recovery examples to confirm our theoretical findings.

Tables

Table 1. Table 1: Reduction rate of the error bounds of the first stage in the 2nd- 5 5 5 th stage

	${\rho }_1$	$c$	0	0.1	0.3	0.5	0.7	0.9
${\varphi }_1$	$[\frac{0.29\alpha }{{\sigma }_r\left(\overline{X}\right)},\frac{\alpha }{{\sigma }_r\left(\overline{X}\right)})$	$\Xi \left({\tilde{\gamma }}_1\right)/\Xi \left({\gamma }_0\right)$	0.819	0.766	0.658	0.547	0.433	0.316
		$\Xi \left({\tilde{\gamma }}_2\right)/\Xi \left({\gamma }_0\right)$	0.818	0.763	0.652	0.537	0.420	0.302
		$\Xi \left({\tilde{\gamma }}_3\right)/\Xi \left({\gamma }_0\right)$	0.818	0.763	0.651	0.536	0.420	0.302
		$\Xi \left({\tilde{\gamma }}_4\right)/\Xi \left({\gamma }_0\right)$	0.818	0.763	0.651	0.536	0.420	0.302
${\varphi }_2$	$[\frac{0.24\alpha }{{\sigma }_r\left(\overline{X}\right)},\frac{4.42\alpha }{{\sigma }_r\left(\overline{X}\right)}]$	$\Xi \left({\tilde{\gamma }}_1\right)/\Xi \left({\gamma }_0\right)$	0.975	0.969	0.955	0.934	0.905	0.856
		$\Xi \left({\tilde{\gamma }}_2\right)/\Xi \left({\gamma }_0\right)$	0.967	0.958	0.931	0.888	0.816	0.689
		$\Xi \left({\tilde{\gamma }}_3\right)/\Xi \left({\gamma }_0\right)$	0.965	0.954	0.920	0.760	0.752	0.572
		$\Xi \left({\tilde{\gamma }}_4\right)/\Xi \left({\gamma }_0\right)$	0.965	0.953	0.915	0.744	0.714	0.516
		$\Xi \left(0\right)/\Xi \left({\gamma }_0\right)$	0.817	0.759	0.644	0.528	0.413	0.297

Table 2. Table 2: Performance for low-rank correlation matrix recovery problems with n = 1000 𝑛 1000 n=1000

			Algorithm 3.1
$r$​	​$\begin{array}{c}\mathrm{off}-\\ \mathrm{diag}\end{array}$​	​eigr ​	The first stage		The first two stages		Final result
			relerr(rank)	​iter​	relerr(rank)	​iter​	$k$	relerr(rank)	iter
	0	1.19	5.94e-1(1000)	89	1.59e-1(6)	1100	5	1.56e-1(5)	4458
	0	2.86	4.43e-1(1000)	101	1.51e-1(5)	629	4	1.52e-1(5)	1805
	0	4.36	3.57e-1(1000)	98	1.49e-1(6)	769	5	1.54e-1(5)	2006
5	100	1.17	5.81e-1(1000)	88	1.53e-1(6)	894	5	1.50e-1(5)	4098
	100	2.79	4.48e-1(1000)	100	1.47e-1(5)	697	4	1.48e-1(5)	1734
	100	4.23	3.60e-1(1000)	98	1.48e-1(6)	832	5	1.49e-1(5)	2004
	0	1.36	4.16e-1(1000)	67	1.48e-1(10)	608	4	1.43e-1(10)	1987
	0	3.52	3.49e-1(1000)	70	1.42e-1(10)	464	4	1.41e-1(10)	1152
	0	6.39	2.95e-1(1000)	59	1.34e-1(10)	373	4	1.37e-1(10)	852
10	100	1.42	3.97e-1(1000)	67	1.46e-1(10)	562	4	1.42e-1(10)	1934
	100	3.31	3.43e-1(1000)	70	1.40e-1(10)	470	4	1.40e-1(10)	1218
	100	6.35	2.87e-1(1000)	66	1.38e-1(10)	344	5	1.42e-1(10)	818

Table 3. Table 3: Performance for low-rank covariance matrix recovery problems with n = 1000 𝑛 1000 n=1000

			Algorithm 3.1
$r$	​$\begin{array}{c}(\mathrm{diag},\\ \mathrm{offdiag})\end{array}$​	​eigr​	The first stage		The first two stages		Final result
			​relerr(rank)​	iter	​relerr(rank)​	iter	$k$	​relerr(rank)​	iter
	$(200,0)$	1.18	4.80e-1(36)	787	2.42e-1(7)	1309	4	2.25e-1(5)	1996
	$(200,0)$	4.59	3.19e-1(32)	676	1.92e-1(8)	1214	4	1.82e-1(5)	1800
5	$(0,200)$	1.20	4.86e-1(36)	275	2.42e-1(7)	451	4	2.21e-1(5)	996
	$(0,200)$	4.26	3.24e-1(33)	349	1.98e-1(7)	523	4	1.90e-1(5)	1016
	$(100,100)$	1.21	4.74e-1(36)	1038	2.42e-1(6)	1734	4	2.24e-1(5)	2746
	$(100,100)$	4.07	3.33e-1(33)	879	1.92e-1(7)	1354	4	1.80e-1(5)	2155
	$(200,0)$	5.33	2.20e-1(53)	253	1.56e-1(13)	409	3	1.54e-1(13)	554
	$(200,0)$	7.72	1.90e-1(48)	258	1.46e-1(13)	393	3	1.45e-1(13)	557
13	$(0,200)$	5.17	2.20e-1(53)	177	1.58e-1(13)	321	3	1.55e-1(13)	466
	$(0,200)$	9.01	1.78e-1(45)	172	1.44e-1(13)	323	3	1.43e-1(13)	484
	$(100,100)$	4.58	2.30e-1(54)	256	1.59e-1(13)	402	3	1.56e-1(13)	548
	$(100,100)$	8.11	1.85e-1(48)	221	1.47e-1(13)	369	3	1.46e-1(13)	530

Equations

$$\min_{X\in\mathbb{R}^{n_{1}\times n_{2}}}\!\big{\{}{\rm rank}(X)\!:\ \|\mathcal{A}X-b\|\leq\delta,\,X\in\Omega\big{\}},$$

$$\mathrm{int}(\mathrm{dom}\varphi )\supseteq [0,1],1>t^{\ast }:={\mathrm{argmin}}_{0\le t\le 1}\varphi (t),\varphi (t^{\ast })=0\mathrm{and}{\varphi }_{-}^{\prime }(1)<+\infty .$$

$$\psi(t):=\!\left\{\!\begin{array}[]{cl}\phi(t)&\textrm{if}\ t\in[0,1],\\ +\infty&\textrm{otherwise}.\end{array}\right.$$

$$\displaystyle\partial\psi^{*}(t)=\big{[}(\psi^{*})_{-}^{\prime}(t),(\psi^{*})_{+}^{\prime}(t)\big{]}\subset[0,1]\quad\forall t\in\mathbb{R},\qquad\qquad$$

$$({\psi }^{\ast }{)}_{+}^{\prime }(t_1)\le ({\psi }^{\ast }{)}_{-}^{\prime }(t)\le ({\psi }^{\ast }{)}_{+}^{\prime }(t)\le ({\psi }^{\ast }{)}_{-}^{\prime }(t_2)\forall t_1<t<t_2.$$

$${\vartheta }_{+}(k):=\underset{0<\mathrm{rank}(X)\le k}{\sup }\frac{⟨X,{\mathcal{A}}^{\ast }\mathcal{A}(X)⟩}{\Vert X{\Vert }_F^2}\mathrm{and}{\vartheta }_{-}(k):=\underset{0<\mathrm{rank}(X)\le k}{\inf }\frac{⟨X,{\mathcal{A}}^{\ast }\mathcal{A}(X)⟩}{\Vert X{\Vert }_F^2},$$

$$\phi(1)\,{\rm rank}(X)=\min_{W\in\mathbb{R}^{n_{1}\times n_{2}}}\!\Big{\{}{\textstyle\sum_{i=1}^{n_{1}}}\phi(\sigma_{i}(W))\!:\ \|X\|_{*}\!-\!\langle W,X\rangle=0,\,\|W\|\leq 1\Big{\}}.$$

$$\Vert X{\Vert }_{\ast }=⟨W,X⟩\le ⟨\sigma (W),\sigma (X)⟩\le \Vert \sigma (X){\Vert }_1=\Vert X{\Vert }_{\ast },$$

$$\underset{X,W\in {\mathbb{R}}^{n_1\times n_2}}{\min }{\sum }_{i=1}^{n_1}\varphi ({\sigma }_i(W))$$

$$s.t.\Vert \mathcal{A}(X)-b\Vert \le \delta ,X\in \Omega ,$$

$$\Vert X{\Vert }_{\ast }-⟨W,X⟩=0,\Vert W\Vert \le 1$$

$$\underset{X,W\in {\mathbb{R}}^{n_1\times n_2}}{\min }{\sum }_{i=1}^{n_1}\varphi ({\sigma }_i(W))+\rho (\Vert X{\Vert }_{\ast }-⟨W,X⟩)$$

$$s.t.\Vert \mathcal{A}(X)-b\Vert \le \delta ,X\in \Omega ,\Vert W\Vert \le 1.$$

$$\{1\}={\textstyle\mathop{\arg\min}_{t\in[0,1]}}\big{\{}\phi(t)+\rho\sigma_{i}(X)(1-t)\big{\}}.$$

$$r\varphi (1)\ge {\sum }_{i=1}^{n_1}\varphi ({\sigma }_i(W^{\ast }))+\rho (\Vert X^{\ast }{\Vert }_{\ast }-⟨W^{\ast },X^{\ast }⟩).$$

$${\sum }_{i=1}^{n_1}\varphi ({\sigma }_i(W^{\ast }))+\rho (\Vert X^{\ast }{\Vert }_{\ast }-⟨W^{\ast },X^{\ast }⟩)$$

$$\displaystyle\geq{\textstyle\sum_{i=1}^{r}}\big{[}\phi(\sigma_{i}(W^{*}))+\rho\sigma_{i}(X^{*})(1-\sigma_{i}(W^{*}))\big{]}$$

$$\displaystyle\geq{\textstyle\sum_{i=1}^{r}}\min_{t\in[0,1]}\big{[}\phi(t)\!+\!\rho\sigma_{i}(X^{*})(1\!-\!t)\big{]}=r\phi(1),$$

$${\sum }_{i=1}^{n_1}\varphi ({\sigma }_i(W^{\ast }))+\rho (\Vert X^{\ast }{\Vert }_{\ast }-⟨W^{\ast },X^{\ast }⟩)$$

$$\displaystyle=\!{\textstyle\sum_{i=1}^{r}}\min_{t\in[0,1]}\big{[}\phi(t)\!+\!\rho\sigma_{i}(X^{*})(1\!-\!t)\big{]}=r\phi(1).$$

$$X^{k}\in\mathop{\arg\min}_{X\in\mathbb{R}^{n_{1}\times n_{2}}}\big{\{}\|X\|_{*}\!-\!\langle W^{k-1},X\rangle\!:\ \|\mathcal{A}(X)\!-b\|\leq\delta,\ X\in\Omega\big{\}}.$$

$$W^{k}\in\mathop{\arg\min}_{W\in\mathbb{R}^{n_{1}\times n_{2}}}\big{\{}{\textstyle\sum_{i=1}^{n_{1}}}\phi(\sigma_{i}(W))\!-\rho_{k}\langle W,X^{k}\rangle\!:\ \|W\|\leq 1\big{\}}.$$

$$\underset{z\in {\mathbb{R}}^{n_1}}{\min }\left\{{\sum }_{i=1}^{n_1}\psi (z_i)-\rho ⟨z,\sigma (X^k)⟩\right\},$$

$$W^{k}:=U^{k}[{\rm Diag}\big{(}w_{1}^{k},w_{2}^{k},\ldots,w_{n_{1}}^{k})\ \ 0](V^{k})^{\mathbb{T}}\ \ {\rm with}\ \ w_{i}^{k}\in\partial\psi^{*}(\rho_{k}\sigma_{i}(X^{k})).$$

$$\psi_{1}^{*}(s)=\left\{\!\begin{array}[]{cl}s\!-\!1&\textrm{if}\ s>1;\\ 0&\textrm{if}\ s\leq 1\end{array}\right.\ {\rm and}\ \ \partial\psi_{1}^{*}(s)=\left\{\begin{array}[]{cl}\{1\}&\textrm{if}\ s>1;\\ \,[0,1]&\textrm{if}\ s=1;\\ \{0\}&\textrm{if}\ s<1.\\ \end{array}\right.$$

$$\partial\psi_{2}^{*}(s)=\left\{\!\begin{array}[]{cl}\{1\}&\textrm{if}\ s\geq\epsilon^{\frac{1}{q-1}}-1;\\ \{1+\epsilon-(s\!+\!1)^{q-1}\}&\textrm{if}\ (1\!+\!\epsilon)^{\frac{1}{q-1}}\!-\!1<s<\!\epsilon^{\frac{1}{q-1}}\!-\!1;\\ \{0\}&\textrm{if}\ s\leq(1\!+\!\epsilon)^{\frac{1}{q-1}}-1.\end{array}\right.$$

$$\frac{{\vartheta }_{+}(r)}{{\vartheta }_{-}(2r+2r)}\le \frac{1+{\delta }_{4r}}{1-{\delta }_{4r}}<1+\frac{2\sqrt{2}-2}{2-\sqrt{2}}<1+2\times 0.843^2.$$

$$\max \left(\frac{{\vartheta }_{+}(r/2)}{{\vartheta }_{-}(2r+r)},\frac{{\vartheta }_{+}((r-1)/2)}{{\vartheta }_{-}(2r+r-1)}\right)\le \frac{1+{\delta }_{r/2}}{1-{\delta }_{3r}}\le \frac{1+{\delta }_{3r}}{1-{\delta }_{3r}}<1+\frac{4\sqrt{5}-8}{5-2\sqrt{5}}.$$

$${\gamma }_{k-1}:=\frac{\Vert {\mathcal{P}}_{\mathcal{T}}(W^{k-1}-{\overline{U}}_1{\overline{V}}_1^{\mathbb{T}}){\Vert }_F}{\sqrt{2r}(1-\Vert {\mathcal{P}}_{{\mathcal{T}}^{\perp }}(W^{k-1})\Vert )}.$$

$$\big{\|}X^{k}\!-\!\overline{X}\big{\|}_{F}\leq\Xi(\gamma_{k-1})\ \ {\rm and}\ \ \big{\|}\mathcal{P}_{\mathcal{T}^{\perp}}(X^{k})\big{\|}_{*}\leq\Gamma(\gamma_{k-1}),$$

$$\Xi (t):=\frac{2\delta \sqrt{{\vartheta }_{+}(2r+s)}}{{\vartheta }_{-}(2r+s)}\cdot \frac{1}{1-ct}\sqrt{1+\frac{r{t}^2}{2s}}\mathrm{and}\Gamma (t):=\frac{2\delta \sqrt{{\vartheta }_{+}(2r+s)}}{{\vartheta }_{-}(2r+s)}\cdot \frac{\sqrt{2r}t}{1-ct}.$$

$$\big{\|}X^{1}\!-\!\overline{X}\big{\|}_{F}\leq\Xi(\gamma_{0})=\Xi(1/\sqrt{2})\ \ {\rm and}\ \ \big{\|}\mathcal{P}_{\mathcal{T}^{\perp}}(X^{1})\big{\|}_{*}\leq\Gamma(\gamma_{0})=\Gamma(1/\sqrt{2}).$$

$$\displaystyle\Xi(\gamma_{0})=\frac{\sqrt{\vartheta_{+}(2.5r)}\sqrt{6}\delta}{\big{(}1\!-\!{c}/{\sqrt{2}}\big{)}\vartheta_{-}(2.5r)}\leq\frac{\sqrt{1+\delta_{3r}}\sqrt{6}\delta}{\big{(}1\!-\!{c}/{\sqrt{2}}\big{)}(1-\delta_{3r})}<\frac{3\delta\sqrt{1+\delta_{3r}}}{1\!-\!\delta_{3r}(1\!+\!{\sqrt{5}}/{2})}.$$

$${\tilde{\gamma }}_k:=\frac{\sqrt{r}(1-{\tilde{b}}_k)+(\sqrt{2}{\tilde{a}}_k+1){\tilde{\beta }}_k}{\sqrt{2r}(1-{\tilde{a}}_k)(1-{\tilde{\beta }}_k^2)}\mathrm{with}{\tilde{a}}_k=({\psi }^{\ast }{)}_{+}^{\prime }\left[{\rho }_k\Xi ({\tilde{\gamma }}_{k-1})\right],$$

Full text

Keywords Structured rank minimization; MPGCC; Exact penalty; Convex relaxation

Mathematics Subject Classification(2010). 90C27, 90C33, 49M20

1 Introduction

The task of noisy structured low-rank matrix recovery is to find a low-rank matrix with a certain structure consistent with some noisy linear measurements. Let $\overline{X}$ be the target matrix to be recovered and $b=\mathcal{A}\overline{X}+\xi$ be the noisy measurement vector, where $\mathcal{A}\!:\mathbb{R}^{n_{1}\times n_{2}}\to\mathbb{R}^{m}$ is the sampling operator and $\xi\in\mathbb{R}^{m}$ is the noisy vector with $\|\xi\|\leq\delta$ for some $\delta>0$. The noisy structured low-rank matrix recovery problem can be modeled as

[TABLE]

where $\Omega\subseteq\mathbb{R}^{n_{1}\times n_{2}}$ is a compact convex set to represent the structure of $\overline{X}$. Without loss of generality, we assume that $\mathcal{A}X\!:=(\langle A_{1},X\rangle,\ldots,\langle A_{m},X\rangle)^{\mathbb{T}}$ for $X\in\mathbb{R}^{n_{1}\times n_{2}}$, where $A_{1},\ldots,A_{m}$ are the given matrices in $\mathbb{R}^{n_{1}\times n_{2}}$. Such a structured rank minimization problem has wide applications in system identification and control [10, 12], signal and image processing [16, 6], machine learning [33], multi-dimensional scaling in statistics [28], finance [27], quantum tomography [14], and so on. For instance, one is often led to seek a low-rank Hankel matrix in system identification and control, a low-rank correlation matrix in finance and a low-rank density matrix in quantum tomography.

Due to the combinatorial property of the rank function, problem (1) is generally NP-hard. One popular way to deal with NP-hard problems is to use the convex relaxation technique, which typically yields a desirable local optimal or at least a feasible solution via solving a single or a sequence of numerically tractable convex optimization problems. Fazel [10] initiated the research for the nuclear norm relaxation method, motivated by the fact that the nuclear norm is the convex envelope of the rank function in the unit ball on the spectral norm. In the past ten years, this relaxation method has received much attention from many fields such as information, computer science, statistics, optimization, and so on (see, e.g., [4, 14, 30, 19, 20, 25, 35]), and it has been shown that a single nuclear norm minimization problem can recover the target matrix $\overline{X}$ under a certain restricted isometry property (RIP) of $\mathcal{A}$ when $\delta\!=0$ [30] or can yield a solution satisfying a certain error bound when $\delta>0$ [5]. For its recoverability and error bounds under other conditions, the interested readers may refer to the literature [9, 25, 31] and references therein.

Most of the existing low-rank matrix optimization models are focused on the case that $\Omega=\mathbb{R}^{n_{1}\times n_{2}}$. When the structure on the target matrix is known, it is reasonable to consider the rank minimization problem (1) with $\Omega$ indicating the available information. However, the (hard) constraint $X\in\Omega$ often contradicts the role of the nuclear norm in promoting a low-rank solution. For example, when $\Omega$ consists of the set of correlation matrices, the nuclear norm relaxation method for (1) may fail in generating a low-rank solution since the nuclear norm becomes a constant in the set $\Omega$. In addition, although some error bounds have been established for the nuclear norm relaxation method in the noisy setting [5, 25, 26], they are minimax-optimal up to a logarithmic factor of the dimension [26], instead of a constant factor like the $l_{1}$-norm relaxation method for sparse regression [29]. These two considerations motivate us to seek more efficient convex relaxations.

1.1 Our main contribution

The main contribution of this work is the introduction of a multi-stage convex relaxation approach via an equivalent Lipschitz optimization reformulation. This approach can efficiently reduce the error bounds obtained from the nuclear norm convex relaxation. More specifically, we reformulate problem (1) as an equivalent MPGCC by using a variational characterization of the rank function and verify that its penalized version, yielded by moving the generalized complementarity constraint to the objective, has the same global optimal solution set as the MPGCC does once the penalty parameter is over a threshold. This exact penalty problem not only has a convex feasible set but also possesses a Lipschitz objective function with a bilinear structure, which offers a favorable Lipschitz reformulation for problem (1). To the best of our knowledge, this is the first equivalent Lipschitz characterization for low-rank matrix optimization problems. Then with this reformulation, we propose a multi-stage convex relaxation approach by solving the exact penalty problem in an alternating way. In particular, under a restricted eigenvalue condition weaker than the RIP condition used in [5, 23], we quantify the reduction of the error and approximate bounds derived from the first stage nuclear norm convex relaxation in the subsequent stages, and establish the geometric convergence of the error sequence in a statistical sense. Among others, the latter entails an upper estimation for the stage number of the convex relaxations to make the estimation error to reach the statistical error level. The analysis shows that the error and approximate rank bounds of the nuclear norm relaxation are reduced most in the second stage and the reduction rate is at least $40\%$ for those problems with a relatively worse restricted eigenvalue property, and the reduction becomes less as the number of stages increases and can be ignored after the fifth stage.

1.2 Related works

The idea of using the multi-stage convex relaxation for low-rank optimization problems is not new. In order to improve the solution quality of the nuclear norm relaxation method, some researchers pay their attention to nonconvex surrogates of low-rank optimization problems. Since seeking a global optimal solution of a nonconvex surrogate problem is almost as difficult as solving a low-rank optimization problem itself, they relax nonconvex surrogates into a sequence of simple matrix optimization problems, and develop the reweighted minimization methods (see [11, 24, 21]). In contrast to our multi-stage convex relaxation approach, such sequential convex relaxation methods are designed by solving a sequence of convex relaxation problems of nonconvex surrogates instead of the equivalent reformulation. We also notice that the theoretical analysis in [23] for the reweighted trace norm minimization method [11] depends on the special property of the log-determinant function, which is not applicable to general low-rank optimization problems, and the theoretical guarantees in [21] were established only for the noiseless recovery problem.

Additionally, some researchers have reformulated low-rank optimization problems as smooth nonconvex problems with the help of low-rank decomposition of matrices in the attempt to achieve a desirable solution by solving the smooth nonconvex problems in an alternating way (actually by solving a sequence of simple convex matrix optimization problems); see, e.g., [32, 18]. This class of convex relaxation methods has a theoretical guarantee, but is not applicable to those problems with hard constraints such as (1).

Finally, it is worthwhile to point out that our multi-stage convex relaxation approach is highly relevant to the one proposed by Zhang [39] for sparse regularization problems and the rank-corrected procedure for the matrix completion problem with fixed coefficients [22]. The former is designed via solving a sequence of convex relaxation problems for the nonconvex surrogates of the zero-norm regularization problem. Since the singular values vectors are involved in low-rank matrix recovery, the analysis technique in [39] is not applicable to our multi-stage convex approach to problem (1). In particular, for low-rank matrix recovery, it is not clear whether the error sequence yielded by the multi-stage convex relaxation approach shrinks geometrically or not in a statistical sense, and if it does, under what conditions. We will answer these questions affirmatively in Section 4. The rank-corrected procedure [22] is actually a two-stage convex relaxation approach in which the first-stage is to find a reasonably good initial estimator and the second-stage is to solve the rank-corrected problem. This procedure has already been applied to nonlinear dimensionality reduction problems [8] and tensor completion problems [3]. However, when the rank of the true matrix is unknown, the rank-corrected problem in [22] needs to be constructed heuristically with the knowledge of the initial estimator, while each subproblem in our multi-stage convex relaxation approach stems from the global exact penalty of the equivalent MPGCC. In addition, the analytical technique used in [22] is more reliant on concentration inequalities in probability analysis, whereas our analysis is deterministic and relies on the restricted eigenvalue property of linear operators.

1.3 Notation

We stipulate $n_{1}\leq n_{2}$ and let $\mathbb{R}^{n_{1}\times n_{2}}$ be the vector space of all $n_{1}\!\times\!n_{2}$ real matrices endowed with the trace inner product $\langle\cdot,\cdot\rangle$ and its induced norm $\|\cdot\|_{F}$. For $X\in\mathbb{R}^{n_{1}\times n_{2}}$, we denote $\sigma(X)\in\mathbb{R}^{n_{1}}$ by the singular value vector of $X$ with entries arranged in a non-increasing order, and $\|X\|_{*}$ and $\|X\|$ by the nuclear norm and the spectral norm of $X$, respectively. Let $\mathbb{O}^{n\times\kappa}$ be the set in $\mathbb{R}^{n\times\kappa}$ consisting of all matrices whose columns are of unit length and are mutually orthogonal to each other, and denote $\mathbb{O}^{n}=\mathbb{O}^{n\times n}$. Let $e$ be the vector of all ones whose dimension is known from the context.

Let $\Phi$ be the family of closed proper convex functions $\phi\!:\mathbb{R}\to(-\infty,+\infty]$ satisfying

[TABLE]

For each $\phi\in\Phi$, let $\psi:\mathbb{R}\to\mathbb{R}\cup\{+\infty\}$ be the associated proper closed convex function:

[TABLE]

Then from convex analysis [34] we know that the conjugate $\psi^{*}$ of $\psi$ has the properties:

[TABLE]

We also need the eigenvalues of $\mathcal{A}^{*}\mathcal{A}$ restricted to a set of low-rank matrices, where $\mathcal{A}^{*}$ denotes the adjoint of $\mathcal{A}$. To this end, for a given positive integer $k$, we define

[TABLE]

which can be viewed as the largest and the smallest rank $k$-restricted eigenvalues of $\mathcal{A}^{*}\mathcal{A}$, respectively.

2 Exact penalty for an equivalent reformulation

First of all, we shall provide an equivalent reformulation of the rank minimization problem (1) with the help of the following variational characterization of the rank function.

Lemma 2.1

Let $\phi\in\Phi$. Then, for any given $X\in\mathbb{R}^{n_{1}\times n_{2}}$, it holds that

[TABLE]

Proof: Fix an arbitrary matrix $X\in\mathbb{R}^{n_{1}\times n_{2}}$. Write $\kappa={\rm rank}(X)$ and assume that $X$ has the SVD as $U[{\rm Diag}(\sigma(X))\ \ 0]V^{\mathbb{T}}$, where $U\!=[U_{1}\ \ U_{2}]\in\mathbb{O}^{n_{1}}$ and $V=[V_{1}\ \ V_{2}]\in\mathbb{O}^{n_{2}}$ with $U_{1}\in\mathbb{O}^{n_{1}\times\kappa}$ and $V_{1}\in\mathbb{O}^{n_{2}\times\kappa}$. Let $W\in\mathbb{R}^{n_{1}\times n_{2}}$ be an arbitrary feasible point of the minimization problem (6). It then follows from [17, Equation (3.3.25)] that

[TABLE]

which implies that $\sum_{i=1}^{n_{1}}(1-\!\sigma_{i}(W))\sigma_{i}(X)=0$. Along with $\sigma_{i}(W)\!\in[0,1]$ for $i=1,\ldots,n_{1}$, we obtain that $\sigma_{i}(W)=1$ if $\sigma_{i}(X)\neq 0$, and $\phi(\sigma_{i}(W))\geq 0$ if $\sigma_{i}(X)=0$. Consequently, $\sum_{i=1}^{n_{1}}\phi(\sigma_{i}(W))\geq\phi(1)\,{\rm rank}(X),$ i.e., $\phi(1)\,{\rm rank}(X)$ is a lower bound for the optimal value of (6). Clearly, $W^{*}\!=U_{1}V_{1}^{\mathbb{T}}+t^{*}U_{2}[{\rm Diag}(e)\ \ 0]V_{2}^{\mathbb{T}}$ with $t^{*}$ defined in (2) is feasible to (6) with the objective value being $\phi(1)\,{\rm rank}(X)$. This shows that $W^{*}$ is an optimal solution of the minimization problem (6) with the optimal value equal to $\phi(1)\,{\rm rank}(X)$. $\Box$

Recall that $\phi(1)>0$ for $\phi\in\Phi$. By Lemma 2.1, we readily have the following result.

Proposition 2.1

Let $\phi\in\Phi$. Then, the rank minimization problem (1) is equivalent to

[TABLE]

in the following sense: if $X^{*}=U^{*}[{\rm Diag}(\sigma(X^{*}))\ \ 0](V^{*})^{\mathbb{T}}$ is a global optimal solution of (1), where $U^{*}=[U_{1}^{*}\ \ U_{2}^{*}]\in\mathbb{O}^{n_{1}}$ and $V^{*}=[V_{1}^{*}\ \ V_{2}^{*}]\in\mathbb{O}^{n_{2}}$ with $U_{1}^{*}\in\mathbb{O}^{n_{1}\times r}$ and $V_{1}^{*}\in\mathbb{O}^{n_{2}\times r}$ for $r={\rm rank}(X^{*})$, then $(X^{*},U_{1}^{*}(V_{1}^{*})^{\mathbb{T}}\!+\!t^{*}U_{2}^{*}[{\rm Diag}(e)\ \ 0](V_{2}^{*})^{\mathbb{T}})$ is globally optimal to (2.1); and conversely, if $(X^{*},W^{*})$ is a global optimal solution to (2.1), then $X^{*}$ is globally optimal to (1).

The constraints $\|X\|_{*}\!-\!\langle W,X\rangle=0$ and $\|W\|\leq 1$ involve a complementarity relation that, for the positive semidefinite (PSD) rank minimization problem, is exactly the PSD cone complementarity relation. In view of this, we call problem (2.1) an MPGCC.

Due to the presence of the nonconvex constraint $\|X\|_{*}\!-\langle W,X\rangle\!=0,$ the MPGCC (2.1) is as difficult as the original problem (1). Nevertheless, it provides us a new view to tackle the difficult rank minimization problem (1). Since numerically it is usually more convenient to handle nonconvex objective functions than to handle nonconvex constraints, we are motivated to investigate the following penalized problem of the MPGCC (2.1):

[TABLE]

Next we shall verify that (2) is an exact penalty version for (2.1) in the sense that there exists a constant $\overline{\rho}>0$ such that the global optimal solution set of (2) associated to any $\rho>\overline{\rho}$ coincides with that of (2.1). To the best of our knowledge, there are only a few works devoted to mathematical programs with matrix cone complementarity constraints [7, 37], which mainly focus on the optimality conditions, but not the exact penalty conditions.

Theorem 2.1

Let the optimal value of (1) be $r>0$. Then, there exists a constant $\alpha>0$ such that $\sigma_{r}(X)\geq\alpha$ for all $X\in\mathcal{F}$, where $\mathcal{F}$ is the feasible set of (1), and for $\phi\in\Phi$ the global optimal solution set of (2) associated to any $\rho>\frac{\phi_{-}^{\prime}(1)}{\alpha}$ is the same as that of (2.1).

Proof: Suppose that there exists a sequence $\{X^{k}\}\subset\mathcal{F}$ such that $\sigma_{r}(X^{k})\rightarrow 0$. Notice that $\{X^{k}\}$ is bounded since $\mathcal{F}$ is bounded. Let $\widehat{X}$ be an accumulation point of $\{X^{k}\}$. By the closedness of $\mathcal{F}$ and the continuity of $\sigma_{r}(\cdot)$, $\widehat{X}\in\mathcal{F}$ and $\sigma_{r}(\widehat{X})=0$, so that ${\rm rank}(\widehat{X})\leq r-1$. This is a contradiction which establishes the existence of $\alpha$.

We denote by $\mathcal{S}$ and $\mathcal{S}^{*}$ the feasible set and the global optimal solution set of (2.1), respectively. For any given $\rho>0$, let $\mathcal{S}_{\rho}$ and $\mathcal{S}_{\rho}^{*}$ be the feasible set and the global optimal solution set of the corresponding penalty problem (2), respectively. By the first part of our conclusion, there exists a constant $\alpha>0$ such that $\sigma_{r}(X)\geq\alpha$ for all $X\in\mathcal{F}$. Let $\rho$ be an arbitrary constant with $\rho>{\phi_{-}^{\prime}(1)}/{\alpha}$. Then, for any $X\in\mathcal{F}$ and each $i\in\{1,\ldots,r\}$,

[TABLE]

First we verify that each $(X^{*},W^{*})\in\mathcal{S}_{\rho}^{*}$ satisfies $\|X^{*}\|_{*}-\langle W^{*},X^{*}\rangle=0$ and ${\rm rank}(X^{*})=r$. Indeed, since $\mathcal{S}^{*}\subset\mathcal{S}\subset\mathcal{S}_{\rho}$ and $r\phi(1)$ is the optimal value of problem (2.1), it holds that

[TABLE]

In addition, from [17, Equation (3.3.25)], it follows that

[TABLE]

where the second inequality is by the nonnegativity of $\phi(\sigma_{i}(W^{*}))$ and $\sigma_{i}(X^{*})(1-\sigma_{i}(W^{*}))$ for all $i$, and the last one is due to (9). Together with (10), we obtain that

[TABLE]

This, along with (9), implies that $\sigma_{i}(W^{*})=1$ for $i=1,\ldots,r$. Substituting $\sigma_{i}(W^{*})=1$ for $i=1,\ldots,r$ into the last equation and using the nonnegativity of $\phi$ in $[0,1]$, we deduce that $\sum_{i=r+1}^{n_{1}}\phi(\sigma_{i}(W^{*}))=0$ and $\|X^{*}\|_{*}=\langle W^{*},X^{*}\rangle=\langle\sigma(X^{*}),\sigma(W^{*})\rangle$. This means that $\sigma_{i}(W^{*})=t^{*}$ for $i=r+\!1,\ldots,n_{1}$ and ${\rm rank}(X^{*})=r$, where $t^{*}<1$ is defined in (2). Then, $\mathcal{S}_{\rho}^{*}\subset\mathcal{S}$ and $\sum_{i=1}^{n_{1}}\phi(\sigma_{i}(W^{*}))=r\phi(1)$ for $(X^{*},W^{*})\in\mathcal{S}_{\rho}^{*}$. Since the global optimal value of (2.1) is $r\phi(1)$, we have $\mathcal{S}_{\rho}^{*}\subseteq\mathcal{S}^{*}$. For the reverse inclusion, let $(X^{*},W^{*})$ be an arbitrary point from $\mathcal{S}^{*}$. Then $(X^{*},W^{*})\in\mathcal{S}_{\rho}$ and $\sum_{i=1}^{n_{1}}\phi(\sigma_{i}(W^{*}))=r\phi(1)$. While the optimal value of (2) is $r\phi(1)$ by the last equation. Thus, we get $S^{*}\subseteq\mathcal{S}^{*}_{\rho}$. $\Box$

Theorem 2.1 extends the exact penalty result of [1, Theorem 3.3] for the zero-norm minimization problem to the matrix setting, and further develops the exact penalty result of the rank-constrained minimization problems in [2, Theorem 3.1]. Observe that the objective function of problem (2) is globally Lipschitz continuous over its feasible set. Combining Theorem 2.1 with Proposition 2.1, we conclude that the rank minimization problem (1) is equivalent to the Lipschitz optimization problem (2).

3 A multi-stage convex relaxation approach

In the last section, we prove that the rank minimization problem (1) is equivalent to a single penalty problem (2). This penalty problem depends on the parameter $\alpha$, the lower bound for the $r$th largest singular value of all $X\in\mathcal{F}$, which can be difficult to estimate. This means that a sequence of penalty problems of the form (2) with non-decreasing $\rho$ should be solved so as to target achieving a global optimal solution of (1). The penalty problem (2) associated to a given $\rho>0$ is not globally solvable due to the nonconvexity of the objective function. However, it becomes a nuclear semi-norm minimization with respect to $X$ if the variable $W$ is fixed and has a closed form solution of $W$ (as will be shown later) if the variable $X$ is fixed. This motivates us to propose a multi-stage convex relaxation approach to (1) by solving a single penalty problem (2) in an alternating way.

Algorithm 3.1

(A multi-stage convex relaxation approach)

(S.0)

Choose a function $\phi\in\Phi$. Let $W^{0}:=0$ and set $k:=1$.

(S.1)

Solve the following nuclear semi-norm minimization problem

X^{k}\in\mathop{\arg\min}_{X\in\mathbb{R}^{n_{1}\times n_{2}}}\big{\{}\|X\|_{*}\!-\!\langle W^{k-1},X\rangle\!:\ \|\mathcal{A}(X)\!-b\|\leq\delta,\ X\in\Omega\big{\}}.

(11)

If $k=1$, select a suitable $\rho_{1}>0$ and go to Step (S.3); else go to Step (S.2).

(S.2)

Select a suitable ratio factor $\mu_{k}\geq 1$ and set $\rho_{k}:=\mu_{k}\rho_{k-1}$.

(S.3)

Solve the following minimization problem

W^{k}\in\mathop{\arg\min}_{W\in\mathbb{R}^{n_{1}\times n_{2}}}\big{\{}{\textstyle\sum_{i=1}^{n_{1}}}\phi(\sigma_{i}(W))\!-\rho_{k}\langle W,X^{k}\rangle\!:\ \|W\|\leq 1\big{\}}.

(12)

(S.4)

Let $k\leftarrow k+1$, and then go to Step (S.1).

The subproblem (11) corresponds to the penalty problem (2) associated to $\rho_{k-1}$ with the variable $W$ fixed to $W^{k-1}$. Since the set $\Omega$ is assumed to be compact, its solution $X^{k}$ is well defined. Assume that $X^{k}$ has the SVD as $U^{k}[{\rm Diag}(\sigma(X^{k}))\ \ 0](V^{k})^{\mathbb{T}}$. By [17, Eq.(3.3.25)], it is easy to check that if $z^{*}$ is optimal to the convex minimization

[TABLE]

then $Z^{*}=U^{k}[{\rm Diag}(z^{*})\ \ 0](V^{k})^{\mathbb{T}}$ is a global optimal solution to (12); and conversely if $W^{*}$ is globally optimal to (12), then $\sigma(W^{*})$ is an optimal solution to (13). Write

[TABLE]

Then, together with [34, Theorem 23.5], it follows that such $W^{k}$ is optimal to the subproblem (12). This means that the main computational work of Algorithm 3.1 consists of solving a sequence of subproblems (11). Unless otherwise stated, in the sequel we choose $w_{i}^{k}=w_{j}^{k}$ when $\sigma_{i}(X^{k})=\sigma_{j}(X^{k})$, which ensures that $1\geq w_{1}^{k}\geq\cdots\geq w_{n_{1}}^{k}\geq 0$.

Since $\|W^{k-1}\|\leq 1$, the function $\|\cdot\|_{*}-\langle W^{k-1},\cdot\rangle$ defines a semi-norm over $\mathbb{R}^{n_{1}\times n_{2}}$. So, the subproblem (11) is a nuclear semi-norm minimization problem. When $k=1$, it reduces to the nuclear norm minimization problem, i.e., the first stage of Algorithm 3.1 is exactly the nuclear norm convex relaxation. It should be emphasized that Algorithm 3.1 is different from the reweighted trace norm minimization method [11, 23] and the iterative reweighted algorithm [21]. The former is proposed from the primal and dual viewpoint by solving an equivalent Lipschitz reformulation in an alternating way, whereas the latter is proposed from the primal viewpoint by relaxing a smooth nonconvex surrogate of (1).

To close this section, we illustrate the choice of $w_{i}^{k}$ in formula (14) with two $\phi\in\Phi$.

Example 3.1

Let $\phi_{1}(t)=t$ for $t\in\mathbb{R}$. Clearly, $\phi_{1}\in\Phi$ with $t^{*}=0$. Moreover, for the function $\psi_{1}$ defined by (3) with $\phi_{1}$, an elementary calculation yields that

[TABLE]

Thus, one may choose w_{i}^{k}\!=\!\left\{\begin{array}[]{ll}1&\textrm{if}\ \sigma_{i}(X^{k})\geq\frac{1}{\rho_{k}};\\ 0&\textrm{otherwise}\end{array}\right. for the matrix $W^{k}$ in formula (14).

Example 3.2

Let $\phi_{2}(t)=-t-\frac{q-1}{q}(1-t+\epsilon)^{\frac{q}{q-1}}+\epsilon+\frac{q-1}{q}$ for $t\in(-\infty,1+\epsilon)$ with $0<q<1$, where $\epsilon\in(0,1)$ is a constant. One can check that $\phi_{2}\in\Phi$ with $t^{*}=\epsilon$. For the function $\psi_{2}$ defined by the equation (3) with $\phi_{2}$, an elementary calculation yields that

[TABLE]

Hence, one may take w_{i}^{k}=\min\big{[}1+\epsilon-(\rho_{k}\sigma_{i}(X^{k})\!+\!1)^{q-1},1\big{]} for the matrix $W^{k}$ in (14),

Remark 3.1

A constant $\epsilon\in(0,1)$ is introduced in the function $\phi_{2}$ so as to ensure that $(\phi_{2})_{-}^{\prime}(1)<+\infty$, and then problem (2) is a global exact penalization of (1). Thus, once $(\widehat{X},\widehat{W})$ yielded by Algorithm 3.1 satisfies $\|X\|_{*}\!-\!\langle X,W\rangle=0$, $\widehat{X}$ is at least a local minimum of the problem (1) since each feasible solution of (1) is locally optimal.

4 Theoretical guarantees of Algorithm 3.1

In this section, we shall provide the theoretical guarantees of Algorithm 3.1 under a mild condition for the restricted eigenvalues of $\mathcal{A}^{*}\mathcal{A}$, which is stated as follows.

Assumption 1

There exist a constant $c\in[0,\sqrt{2})$ and an integer $s\in[1,\frac{n-2r}{2}]$ such that $\frac{\vartheta_{+}(s)}{\vartheta_{-}(2r+2s)}\leq 1+\frac{2c^{2}s}{r}$, where $\vartheta_{+}(\cdot)$ and $\vartheta_{-}(\cdot)$ are the functions defined by (5).

Assumption 1 requires the restricted eigenvalue ratio of $\mathcal{A}^{*}\mathcal{A}$ to grow sublinearly in $s$. This condition, extending the sparse eigenvalue condition used for the analysis of sparse regularization (see [38, 39]), is weaker than the RIP condition $\delta_{4r}<\sqrt{2}-1$ used in [5] for $n\geq 4r$, where $\delta_{kr}$ is the $kr$-restricted isometry constant of $\mathcal{A}$ defined as in [5]. Indeed, from the definitions of $\vartheta_{+}(\cdot)$ and $\vartheta_{-}(\cdot)$, it is immediate to have that

[TABLE]

This shows that $c=0.843$ is such that $\frac{\vartheta_{+}(s)}{\vartheta_{-}(2r+2s)}\leq 1+\frac{2c^{2}s}{r}$ for $s=r$. In addition, this condition is also weaker than the RIP condition $\delta_{3r}<2\sqrt{5}\!-\!4$ used in [23] for $n\geq 3r$, where $r$ is an even number or $r$ is an odd number greater than $11$. To see this, assume that $\delta_{3r}<2\sqrt{5}\!-\!4$, and $r$ is an even number or is an odd number greater than $11$. Then,

[TABLE]

So, $c=1.34$ and $1.403$ are respectively such that $\frac{\vartheta_{+}(s)}{\vartheta_{-}(2r+2s)}\leq 1+\frac{2c^{2}s}{r}$ for $s=\frac{r}{2}$ and $\frac{r-1}{2}$.

In the sequel, we let $\overline{X}$ have the SVD as $\overline{U}[{\rm Diag}(\sigma(\overline{X}))\ \ 0]\overline{V}^{\mathbb{T}}$, where $\overline{U}=[\overline{U}_{1}\ \ \overline{U}_{2}]\in\mathbb{O}^{n_{1}}$ and $\overline{V}=[\overline{V}_{1}\ \ \overline{V}_{2}]\in\mathbb{O}^{n_{2}}$ with $\overline{U}_{1}\in\mathbb{O}^{n_{1}\times r}$ and $\overline{V}_{1}\in\mathbb{O}^{n_{2}\times r}$ for $r={\rm rank}(\overline{X})$, and write $\mathcal{T}=\mathcal{T}(\overline{X})$ where $\mathcal{T}(\overline{X})$ is the tangent space at $\overline{X}$ associated to the rank constraint ${\rm rank}(X)\leq r$ (see equation (28) for its definition). For convenience, for $k=1,2,\ldots$, let

[TABLE]

The proofs of all the results in the subsequent subsections are given in Appendix C.

4.1 Error and approximate rank bounds

Under Assumption 1, when $\gamma_{k-1}\in[0,1/c)$ for some $k\geq 1$, we can establish the following error bound and approximate rank bound for the solution $X^{k}$ of the $k$th subproblem.

Proposition 4.1

Suppose that Assumption 1 holds and $0\leq\gamma_{k-1}<{1}/{c}$ for some $k\!\geq 1$. Then

[TABLE]

where $\Xi\!:[0,1/c)\to\mathbb{R}_{+}$ and $\Gamma\!:[0,1/c)\to\mathbb{R}_{+}$ are the increasing functions defined by

[TABLE]

Remark 4.1

(a)* Since $\|\mathcal{P}_{\mathcal{T}^{\perp}}(X^{k})\|_{*}\!=0$ implies that ${\rm rank}(X^{k})\!\leq 2r$, it is reasonable to view $\|\mathcal{P}_{\mathcal{T}^{\perp}}(X^{k})\|_{*}$ as a measure for the approximate rank of $X^{k}$. So, the second inequality in (18) provides an approximate rank bound for $X^{k}$. The error and approximate rank bounds in (18) consist of two parts: one part is the statistical error $\Xi(0)=\frac{2\delta\sqrt{\vartheta_{+}(2r\!+\!s)}}{\vartheta_{-}(2r\!+\!s)}$ from the noise and the operator $\mathcal{A}$, and the other part is the estimation error from $\gamma_{k-1}$.*

(b)* Since $W^{0}=0$, we have $\gamma_{0}=\frac{1}{\sqrt{2r}}\|\overline{U}_{1}\overline{V}_{1}^{\mathbb{T}}\|_{F}=\!\frac{1}{\sqrt{2}}<\frac{1}{c}$. Hence, under Assumption 1, the error and approximate rank bounds of the nuclear norm convex relaxation are*

[TABLE]

Moreover, if Assumption 1 is satisfied with $s=r/2$ and $c<\sqrt{2}-\frac{2(1-\delta_{3r}(1+{\sqrt{5}}/{2}))}{\sqrt{{3}}(1-\delta_{3r})}$ for $\delta_{3r}<2\sqrt{5}\!-\!4$, then the error bound $\Xi(\gamma_{0})$ is tighter than the bound $\frac{3\delta\sqrt{1+\delta_{3r}}}{1-\delta_{3r}(1+{\sqrt{5}}/{2})}$ given by **[23, Theorem III.1]** with $C_{1,1}=1$ for the nuclear norm relaxation because

[TABLE]

Remark 4.1 (b) says that under Assumption 1 the solution $X^{1}$ of the first stage convex relaxation has the error and approximate rank bounds as in (19). However, it is not clear whether $X^{k}\ (k\geq 2)$ has such error and approximate rank bounds or not. The following theorem states that if in addition $\sigma_{r}(\overline{X})>2\Xi(\gamma_{0})$ and $\rho_{1}$ and $\mu_{k}$ are appropriately chosen, all $X^{k}\ (k\geq 2)$ have the bounds as in (18), and more importantly, their error and approximate rank bounds are, respectively, smaller than those of $X^{1}$. To achieve this result, we need the sequence $\{\widetilde{\gamma}_{k}\}_{k\geq 1}$, which is defined recursively with $\widetilde{\gamma}_{0}=\gamma_{0}$ as

[TABLE]

Theorem 4.1

Suppose that Assumption 1 holds and $\sigma_{r}(\overline{X})>2\Xi(\gamma_{0})$. If the parameters $\rho_{1}$ and $\mu_{k}$ are respectively chosen such that $\widetilde{a}_{1}<\!\frac{(\widetilde{b}_{1}-\widetilde{\beta}_{1}^{2})\sqrt{r}-\widetilde{\beta}_{1}}{(1-\widetilde{\beta}_{1}^{2})\sqrt{r}+\sqrt{2}\widetilde{\beta}_{1}}$ and \mu_{k}\in\big{[}1,\frac{\Xi(\widetilde{\gamma}_{k-2})}{\Xi(\widetilde{\gamma}_{k-1})}\big{]}, then all $X^{k}\ (k\geq 1)$ satisfy the inequalities in (18) and for $k\geq 2$ it also holds that

[TABLE]

Remark 4.2

(a)* Theorem 4.1 shows that under Assumption 1 and $\sigma_{r}(\overline{X})\!>\!2\Xi(\gamma_{0})$, if $\rho_{1}$ and $\mu_{k}$ are chosen appropriately, then the error and approximate rank bounds of $X^{k}\ (k\geq 2)$ improve those of $X^{1}$ at least by $1-\frac{\Xi(\widetilde{\gamma}_{k-1})}{\Xi(\gamma_{0})}$ and $1\!-\!\frac{\Gamma(\widetilde{\gamma}_{k-1})}{\Gamma(\gamma_{0})}$, respectively.*

(b)* The choice of $\rho_{1}$ depends on $\Xi(\gamma_{0})$. For instance, take the function $\phi_{1}$ in Example 3.1. If $\sigma_{r}(\overline{X})=\alpha\Xi(\gamma_{0})$ for $\alpha\geq 2.5$, then by virtue of the definitions of $\widetilde{a}_{1},\widetilde{b}_{1}$ and $\widetilde{\beta}_{1}$ and equation (15) it is easy to check that $\widetilde{a}_{1}=0,\widetilde{b}_{1}=1$ and $\widetilde{\beta}_{1}\in[0,0.6)$, and consequently $(\widetilde{b}_{1}-\widetilde{\beta}_{1}^{2})\sqrt{r}-\widetilde{\beta}_{1}>0$. This means that \big{(}\frac{1}{(\alpha-1)\Xi(\gamma_{0})},\frac{1}{\Xi(\gamma_{0})}\big{)} is the range of choice for $\rho_{1}$. For numerical computations, one may estimate $r$ and $\sigma_{r}(\overline{X})$ with the help of $\sigma(X^{1})$.*

To close this subsection, we illustrate the ratios $\frac{\Xi(\widetilde{\gamma}_{k-1})}{\Xi(\gamma_{0})}$ and $\frac{\Gamma(\widetilde{\gamma}_{k-1})}{\Gamma(\gamma_{0})}$ by the functions $\phi_{1}$ and $\phi_{2}$ with $q=1/2$ and $\epsilon=10^{-3}$. For this purpose, we suppose that Assumption 1 holds with $r=10,s\!=\!r/2$ and $\sigma_{r}(\overline{X})=\alpha\Xi(\gamma_{0})$ for $\alpha\!\geq\!4.5$. Then, for those $c$ in the first row of Table 1, one may compute the ratios $\frac{\Xi(\widetilde{\gamma}_{k-1})}{\Xi(\gamma_{0})}$ and $\frac{\Gamma(\widetilde{\gamma}_{k-1})}{\Gamma(\gamma_{0})}$ as those in the last six columns of Table 1 with $\rho_{1}$ chosen as the middle point of the interval and $\mu_{k}\equiv 1$. We see that the error bound of the first stage is reduced most in the second stage, and as the number of stages increases, the reduction becomes less. For Algorithm 3.1 with $\phi_{1}$, the reduction is close to the limit ${\Xi(0)}/{\Xi(\widetilde{\gamma}_{0})}$ when $k=5$, but for Algorithm 3.1 with $\phi_{2}$, there is a little room for the reduction especially for those $\mathcal{A}^{*}\mathcal{A}$ with $c\geq 0.5$.

4.2 Geometric convergence

Generally speaking, because of the presence of the noise, it is impossible for the error sequence $\{\|X^{k}\!-\!\overline{X}\|_{F}\}_{k\geq 1}$ to decrease and then converge geometrically. However, one may establish its geometric convergence in a statistical sense as in the following theorem.

Theorem 4.2

Suppose that Assumption 1 holds and $\sigma_{r}(\overline{X})\!>\max(2,\sqrt{2}\!+\!\alpha)\Xi(\gamma_{0})$ for $\alpha=\frac{1+\sqrt{2}\widetilde{a}_{1}}{(1-\widetilde{a}_{1})(1\!-\!\widetilde{\beta}_{1}^{2})\sqrt{r+4s}}$. If $\rho_{1}$ and $\mu_{k}$ are chosen as in Theorem 4.1, then for $k\geq 1$,

[TABLE]

Remark 4.3

(a)* The requirement $\sigma_{r}(\overline{X})\!>\max(2,\sqrt{2}\!+\!\alpha)\Xi(\gamma_{0})$ in Theorem 4.2 is bit stronger than $\sigma_{r}(\overline{X})\!>2\Xi(\gamma_{0})$. Take $\phi_{1}$ for example. When $\sigma_{r}(\overline{X})\geq 2.4\Xi(\gamma_{0})$, this requirement is automatically satisfied. Also, now we have that $\varrho:=\frac{\alpha\,\Xi(\gamma_{0})}{\sigma_{r}(\overline{X})\!-\!\sqrt{2}\Xi(\gamma_{0})}\leq 0.76$.*

(b)* The first term of the sum on the right hand side of (21) represents the statistical error arising from the noise and the sampling operator $\mathcal{A}$, and the second term is the estimation error related to the multi-stage convex relaxation. Clearly, the statistical error is of a certain order of $\Xi(0)$. Thus, to guarantee that the second term is less than the statistical error, at most $\overline{k}$ stage convex relaxations are required, where*

[TABLE]

Take $\varrho=0.7$ for example. When $s=r$, one can calculate that $\overline{k}\leq 2$ if $c=0.3$, and $\overline{k}\leq 4$ if $c=0.7$. This means that, for those $\mathcal{A}^{*}\mathcal{A}$ with a worse restricted eigenvalue condition, more than two stage convex relaxations are needed to yield a satisfactory solution.

For the analysis in the previous two subsections, the condition $\sigma_{r}(\overline{X})\!\geq\alpha\Xi(\gamma_{0})$ for a certain $\alpha>2$ is required for the decreasing of the error and approximate rank bounds of the first stage convex relaxation and the contraction of the error sequence. Such a condition is necessary for the low-rank recovery since, when the smallest nonzero singular is mistaken as a zero, the additional singular vectors will yield a large error. In fact, in the geometric convergence analysis of sparse vector optimization (see [39]), the error bound of the first stage was implicitly assumed not to be too large. In addition, we observe that the structure information of $\overline{X}$ does not lend any help to the low-rank matrix recovery in terms of convergence rates. However, when the true matrix has a certain structure, it is necessary to incorporate such structure information into model (1). Otherwise, the solution $X^{k}$ yielded by the multi-stage convex relaxation may not satisfy the structure constraint, and then it is impossible to control the error of $X^{k}$ to the true matrix $\overline{X}$.

Finally, we point out that when the components $\xi_{1},\xi_{2},\ldots,\xi_{m}$ of the noisy vector $\xi$ are independent (but not necessarily identically distributed) sub-Gaussians, i.e., there exists a constant $\sigma\geq 0$ such that $\mathbb{E}[e^{t\xi_{i}}]\leq e^{\sigma^{2}t^{2}/2}$ holds for all $i$ and any $t\in\mathbb{R}$, by Lemma 7 in Appendix C, the conclusions of Theorems 4.1 and 4.2 hold with $\delta=\sqrt{m}\sigma$ with probability at least $1-\exp(1-\frac{c_{1}m}{4})$ for an absolute constant $c_{1}>0$. For the random $\mathcal{A}$, the following result is immediate by [5, Theorem 2.3] and the first inequality in (16).

Theorem 4.3

Fix $\overline{\delta}\in(0,1/2)$ and let $\mathcal{A}$ be a random measurement ensemble obeying the following conditions: for any given $X\in\mathbb{R}^{n_{1}\times n_{2}}$ and any fixed $0<t<1$,

[TABLE]

for fixed constants $C,c_{2}>0$ (which may depend on $t$). If $m\geq 3C(n_{1}\!+\!n_{2}\!+\!1)r$ with $C>\frac{\log(36\sqrt{2}/\overline{\delta})}{c_{2}}$, then Assumption 1 holds for $s=r/2$ and $c=\sqrt{\frac{2\overline{\delta}}{1-\overline{\delta}}}$ with probability exceeding $1\!-\!2\exp(-dm)$ where $d\!=c_{2}-\frac{\log(36\sqrt{2}/\overline{\delta})}{C}$. Consequently, when $0\leq\!\gamma_{k-1}<\!1/c$, the bounds in (18) holds with probability at least $1\!-\!2\exp(-dm)$ for such random measurements.

As remarked after [5, Theorem 2.3], the condition in (22) holds when $\mathcal{A}$ is a Gaussian random measurement ensemble (i.e., $A_{1},\ldots,A_{m}$ are independent from each other and each $A_{i}$ contains i.i.d. entries $\mathcal{N}(0,1/m)$); or when each entry of each $A_{i}$ has i.i.d. entries that are equally likely to take $\frac{1}{\!\sqrt{m}}$ or $-\frac{1}{\!\sqrt{m}}$; or when $\mathcal{A}$ is a random projection (see [30]).

5 Numerical experiments

In this section, we shall test the theoretical results in Section 4 by applying Algorithm 3.1 to some low-rank matrix recovery problems, including matrix sensing and matrix completion problems. During the testing, we chose $\phi_{2}$ with $q=1/2$ and $\epsilon=10^{-3}$ for the function $\phi$ in Algorithm 3.1. Although Table 1 shows that Algorithm 3.1 with $\phi_{1}$ reduces the error faster than Algorithm 3.1 with $\phi_{2}$ does, our preliminary testing indicates that the latter has a little better performance in reducing the relative error. This accounts for choosing $\phi_{2}$ instead of $\phi_{1}$ for our numerical testing. In addition, we always chose $\rho_{1}=10/\|X^{1}\|$ and $\mu_{k}={5}/{4}\ (k\geq 1)$. All the results in the subsequent subsections were run on the Windows system with an Intel(R) Core(TM) i3-2120 CPU 3.30GHz.

5.1 Low-rank matrix sensing problems

We tested the performance of Algorithm 3.1 with some matrix sensing problems, for which some entries are known exactly. Specifically, we assumed that ${\bf 5}$ entries of the true $\overline{X}\in\mathbb{R}^{n_{1}\times n_{2}}$ of rank $r$ are known exactly, and generated $\overline{X}$ in the following command

     XR = randn(n1,r);  XL = randn(n2,r); Xbar = XR*XL’.

We successively generated the matrices $A_{1},\ldots,A_{m}\in\mathbb{R}^{n_{1}\times n_{2}}$ with i.i.d. standard normal entries to formulate the sampling operator $\mathcal{A}$. Such $\mathcal{A}$ satisfies the RIP property with a high probability by [30], which means that the restricted eigenvalues of $\mathcal{A}^{*}\mathcal{A}$ can satisfy Assumption 1 with a high probability from the discussions after Assumption 1. Then, we successively generated the standard Gaussian noises $\xi_{1},\ldots,\xi_{m}$ to formulate $b$ by

[TABLE]

Take $\delta=0.1\|b\|$ and \Omega=\big{\{}X\in\mathbb{R}^{n_{1}\times n_{2}}\ |\ \mathcal{B}X=d,\,\|X\|\leq R\big{\}} for $R=2\|\overline{X}\|$, where $\mathcal{B}X:=(X_{ij})$ for $(i,j)\in\Upsilon_{\rm fix}$ with $\Upsilon_{\rm fix}$ being the index set of known entries, and $d\in\mathbb{R}^{|\Upsilon_{\rm fix}|}$ is the vector consisting of $\overline{X}_{ij}$ for $(i,j)\in\Upsilon_{\rm fix}$. Let $\mathbb{I}_{S}(\cdot)$ denote the indicator function over a set $S$. The subproblem (11) in Algorithm 3.1 now has the form

[TABLE]

where $\mathcal{R}\!:=\{z\in\mathbb{R}^{m}\ |\ \|z\|\leq\delta\}$ and $\Lambda\!:=\{X\in\mathbb{R}^{n_{1}\times n_{2}}\ |\ \|X\|\leq R\}$. The dual of (5.1) is

[TABLE]

During the testing, we solved the subproblems of the form (5.1) until the primal and dual relative infeasibility is less than $10^{-6}$ and the difference between the primal objective value and the dual one is less than $10^{-5}$, with the powerful Schur-complement based semi-proximal ADMM (alternating direction method of multipliers) [15] for problem (5.1).

5.1.1 Performance of Algorithm 3.1 in different stages

We generated randomly a matrix sensing problem with some entries known as above with $n_{1}=n_{2}=100$, $r=6$ and $m=2328$ to test the performance of Algorithm 3.1 in different stages. Figure 1 plots the relative error of Algorithm 3.1 in the first fifteen stages. We see that Algorithm 3.1 reduces the relative error of the nuclear norm relaxation method most in the second stage, and after the third stage the reduction becomes insignificant. This performance coincides with the analysis results shown as in Table 1.

5.1.2 Performance of Algorithm 3.1 with different samples

We generated randomly a matrix sensing problem with some entries known as above with $n_{1}=n_{2}=100$, $r=5$ and $m=\nu r(2n\!-\!r)$ for $\nu\in\{1.0,1.1,\ldots,3.0\}$ to test the performance of Algorithm 3.1 under different samples. Figure 2 depicts the relative error curves and the rank curves of the first stage convex relaxation and the first five stages convex relaxation, respectively. We see that the relative errors of the first stage convex relaxation and the first five stages convex relaxation decrease as the number of samples increases, but the relative error of the latter is always smaller than that of the former. Moreover, the first five stages convex relaxation reduces those of the first stage convex relaxation at least $25\%$ for $\nu\in[1.0,3.0]$, and the reduction becomes less as the number of samples increases. In particular, the rank of $X^{1}$ is higher than that of $\overline{X}$ even for $30\%$ sampling ratio, but the rank of $X^{5}$ equals that of $\overline{X}$ even for $12\%$ sampling ratio.

5.2 Low-rank PSD matrix completion problems

We applied Algorithm 3.1 to two classes of low-rank PSD matrix completion problems. Although the sampling operators for such problems do not satisfy the RIP property, it is possible for the restricted eigenvalues of $\mathcal{A}^{*}\mathcal{A}$ to satisfy Assumption 1. For these problems, \Omega=\big{\{}X\in\mathbb{S}_{+}^{n}\ |\ \mathcal{E}_{1}(X)=g_{1},\,\mathcal{E}_{2}(X)\leq g_{2}\big{\}} where $\mathcal{E}_{1}\!:\mathbb{S}^{n}\to\mathbb{R}^{l_{1}}$ and $\mathcal{E}_{2}\!:\mathbb{S}^{n}\to\mathbb{R}^{l_{2}}$ are the linear operators, and $g_{1}\in\mathbb{R}^{l_{1}}$ and $g_{2}\in\mathbb{R}^{l_{2}}$ are the given vectors. For this case, the subproblem (11) in Algorithm 3.1 now takes the form of

[TABLE]

After an elementary calculation, the dual problem of (5.2) has the following form

[TABLE]

During the testing, we solved the subproblems of the form (5.2) until the primal and dual relative infeasibility is less than $10^{-6}$ and the difference between the primal objective value and the dual one is less than $10^{-5}$ with the Schur-complement based semi-proximal ADMM [15] for problem (5.2), and stopped Algorithm 3.1 at the $k$th iteration once

[TABLE]

where ${\rm rank}(X^{k})$ is the number of nonzero singular values of $X^{k}$ less than $10^{-10}\cdot\sigma_{\rm max}(X^{k})$.

5.2.1 Low-rank correlation matrix completion problems

A correlation matrix is a real symmetric PSD matrix with all diagonals being $1$. We generated the true correlation matrix $\overline{X}\in\mathbb{S}_{+}^{n}$ of rank $r$ in the following command:

L = randn(n,r);  W = weight*L(:,1:1);  L(:,1:1) = W; G = L*L’;
M = diag(1./sqrt(diag(G)))*G*diag(1./sqrt(diag(G))); Xbar = (M+M’)/2.

In this way, one can control the ratio of the largest eigenvalue and the smallest nonzero eigenvalue of $\overline{X}$ by weight. We assume that some off-diagonal entries of $\overline{X}$ are known. Thus, $\mathcal{E}_{1}(X)=\left(\begin{matrix}{\rm diag}(X)\\ \mathcal{B}(X)\end{matrix}\right)$ for $X\in\mathbb{S}^{n}$, $g_{1}=\left(\begin{matrix}e\\ d\end{matrix}\right)$, $\mathcal{E}_{2}\equiv 0$ and $g_{2}=0$, where the operator $\mathcal{B}\!:\mathbb{S}^{n}\to\mathbb{R}^{|\Upsilon_{\rm fix}|}$ and the vector $d\in\mathbb{R}^{|\Upsilon_{\rm fix}|}$ are defined as in Subsection 5.1. The noise vector $\xi$ and the observation vector $b$ are generated in the same way as in (23).

Table 2 reports the numerical results of Algorithm 3.1 for some examples generated randomly. The information of $\overline{X}$ is reported in the first three columns, where the second column lists the number of known off-diagonal entries for $\overline{X}$, and the third column gives the ratio of the largest eigenvalue of $\overline{X}$ to the smallest nonzero eigenvalue of $\overline{X}$. For each test example, we sampled partial unknown off-diagonal entries uniformly at random to formulate the operator $\mathcal{A}$, where the sample ratio is ${\bf 1.92\%}$ for ${\rm rank}(\overline{X})=5$ and ${\bf 4.32\%}$ for ${\rm rank}(\overline{X})=10$. The fourth and the fifth columns report the results of the first stage convex relaxation and the first two stages convex relaxation, respectively, and the sixth column reports the final result, where relerr(rank) means the relative error and the rank of solutions, and iter is the total number of iterations required by the Schur-complement based semi-proximal ADMM for the corresponding convex relaxation.

We see that the solution given by the trace-norm relaxation method has a high relative error and a full rank, while the two-stage convex relaxation reduces the relative error of the trace-norm relaxation at least $50\%$ for all test problems. Although the two-stage convex relaxation may yield the desirable relative error for all the test problems, the ranks of some problems (for example, the third and the fourth) are higher than that of $\overline{X}$. With the number of stages increasing, Algorithm 3.1 yields the same rank as that of $\overline{X}$. This indicates that for the problems with suitable sample ratios, the two-stage convex relaxation is enough; while for the problem with very low sample ratios, more than two stages convex relaxation is needed. In addition, since all the constraints to define the set $\Omega$ are of the hard type, some of the relative errors in the sixth column are little higher than those in the fifth column.

5.2.2 Low-rank covariance matrix completion problems

We generated the true covariance matrix $\overline{X}\in\mathbb{S}_{+}^{n}$ of rank $r$ in the following command:

   L = randn(n,r)/sqrt(sqrt(n));  W = weight*L(:,1:1);
   L(:,1:1) = W;  G = L*L’; Xbar = (G+G’)/2.

In this case, $\mathcal{E}_{1}=\mathcal{B}$ and $g_{1}=d$ where $\mathcal{B}\!:\mathbb{S}^{n}\to\mathbb{R}^{|\Upsilon_{\rm fix}|}$ and $d\in\mathbb{R}^{|\Upsilon_{\rm fix}|}$ are defined as in Subsection 5.1, $\mathcal{E}_{2}(X):=(X_{ii})$ for $(i,i)\in\Upsilon_{\rm diag}$ with $\Upsilon_{\rm diag}$ being the index set of unknown diagonal entries of $\overline{X}$, and $g_{2}\in\mathbb{R}^{|\Upsilon_{\rm diag}|}$ is the vector consisting of the upper bounds for unknown diagonal entries of $\overline{X}$. We set $g_{2}=(1+0.01\textbf{rand}(1,1))\|\overline{X}\|_{\infty}\textbf{ones}(|\Upsilon_{\rm diag}|,1)$.

Table 3 reports the numerical results of Algorithm 3.1 for some problems generated randomly. The information of the true covariance matrix $\overline{X}$ is reported in the first two columns, where the second column lists the number of known diagonal and off-diagonal entries of $\overline{X}$, and the third column reports the ratio of the largest eigenvalue of $\overline{X}$ to the smallest nonzero eigenvalue of $\overline{X}$. For each test example, we sampled the upper triangular entries uniformly at random to formulate the sampling operator $\mathcal{A}$, where the sample ratio is ${\bf 1.91\%}$ for ${\rm rank}(\overline{X})=5$ and ${\bf 5.72\%}$ for ${\rm rank}(\overline{X})=13$. The fourth and the fifth columns report the results of the first stage convex relaxation and the first two stages convex relaxation, respectively, and the last one lists the final results.

We see that the solution yielded by the trace-norm relaxation method has a high relative error and rank, and the solution given by the first two stages convex relaxation has the desirable relative error but its rank is still higher than that of the true matrix for those problems with low sample ratios. This shows that for these difficult problems, more than two stages convex relaxation is required. For the problems with with $1.91\%$ sample ratio, the two-stage convex relaxation reduces the error bounds of the trace-norm relaxation method at least $38\%$, while for those problems with $5.32\%$ sample ratio, the reduction rate is over $19\%$. In addition, combining with the results in Table 2, we see that the performance of our multi-stage convex relaxation has no direct like with the ratio of the largest eigenvalue and the smallest nonzero eigenvalue of the true $\overline{X}$.

6 Conclusions

We have proposed a multi-stage convex relaxation approach to the structured rank minimization problem (1) by solving the exact penalty problem of its equivalent MPGCC in an alternating way. It turned out that this approach not only has favorable theoretical guarantees but also reduces the error of the nuclear norm relaxation method effectively. There are several topics worthwhile to pursue, such as to develop the fast and effective algorithms for seeking the solution of subproblems, to establish the theoretical guarantee for the case where the subproblems are solved inexactly, and apply this approach to other classes of low-rank optimization problems, say, low-rank plus sparse problems.

Acknowledgements  This work is supported by the National Natural Science Foundation of China under project No.11571120 and the Natural Science Foundation of Guangdong Province under project No. 2015A030313214. The authors would like to thank Dr Xudong Li on help us improve our computing codes.