Continuous Iteratively Reweighted Least Squares Algorithm for Solving Linear Models by Convex Relaxation

Xian Luo, Wanzhou Ye

Department of Mathematics, College of Science, Shanghai University, Shanghai, China.

DOI: 10.4236/apm.2019.96024                

Abstract

In this paper, we present continuous iteratively reweighted least squares algorithm (CIRLS) for solving the linear models problem by convex relaxation, and prove the convergence of this algorithm. Under some conditions, we give an error bound for the algorithm. In addition, the numerical result shows the efficiency of the algorithm.

Keywords

Linear Models, Continuous Iteratively Reweighted Least Squares, Convex Relaxation, Principal Component Analysis

1. Introduction

Low-dimensional linear models have broad applications in data analysis problems such as computer vision, pattern recognition, machine learning, and so on [1] [2] [3] . In most of these applications, the data set is noisy and contains numbers of outliers so that they are distributed in higher dimensional. Meanwhile, principal component analysis is the standard method [4] for finding a low-dimensional linear model. Mathematically, the problem can be present as

$\min {\displaystyle \underset{x\in \chi }{\sum }}{\Vert x-\Pi x\Vert }^2\mathrm{<mi>\; </mi>\; <mi>\; </mi>}\text{subject}\text{to}\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}\Pi \text{is an orthoprojector}\text{and}tr\Pi =d$ (1)

where $\chi$ is the data set consisting of N points in ${ℝ}^D$, $x\in {ℝ}^D$, a target dimension $d\in \left\{\mathrm{1,2,}\cdots \mathrm{,}D-1\right\}$, ${\Vert \cdot \Vert }^2$ denotes the $l_2$ norm on vectors, and tr refers to the trace.

Unfortunately, if the data set contains a large number of noise in the inliers and a substantial number of outliers, these nonidealities points can interfere with the linear models. To guard the subspace estimation procedure against outliers, statistics have proposed to replace the $l_2$ norm [5] [6] [7] with $l_1$ norm that is less sensitive to outliers. This idea leads the following optimization problem

$\min {\displaystyle \underset{x\in \chi }{\sum }}\Vert x-\Pi x\Vert \mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}\text{subject}\text{to}\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}\Pi \text{is an orthoprojector}\text{and}tr\Pi =d$ (2)

The optimization problem (2) is not convex, and we have no right to expect that the problem is tractable [8] [9] . Wright [10] proved that most matrices can be efficiently and exactly recovered from most error sign-and-support patterns by solving a simple convex program, for which they give a fast and provably convergent algorithm. Later, Candes [11] presented that under some suitable assumptions, it is possible to recover the subspace by solving a very convenient convex program called Principal Component Pursuit. Base on convex optimization, Lerman [12] proposed to use a relaxation of the set of orthogonal projectors to reach the convex formulation, and give the linear model as follows

$\min {\displaystyle \underset{x\in \chi }{\sum }}\Vert x-Px\Vert \mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}\text{subjectto}\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>0}\preccurlyeq P\preccurlyeq I\mathrm{<mi>\; </mi>\; <mi>\; </mi>}\text{and}\mathrm{<mi>\; </mi>\; <mi>\; </mi>}trP=d$ (3)

where the matrix P is the relaxation of orthoprojector $\Pi$, whose eigenvalues lie in the interval $\left[\mathrm{0,1}\right]$ form a convex set. The curly inequality $\preccurlyeq$ denotes the semidefinite order: for symmetric matrices A and B, we write $A\preccurlyeq B$ if and only if $B-A$ is positive semidefinite, the problem (3) is called REAPER.

To obtain a d-dimensional linear model from a minimizer of REAPER, we need to consider the auxiliary problem

$\min {\displaystyle \underset{x\in \chi }{\sum }}{\Vert P_{\star }-\Pi \Vert }_{S_1}\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}\text{subjectto}\mathrm{<mi>\; </mi>\; <mi>\; </mi>}\Pi \text{is an orthoprojector}\text{and}tr\Pi =d$ (4)

where $P_{\star }$ is an optimal point of REAPER, ${\Vert \cdot \Vert }_{S_1}$ denotes Schatten 1-norm, in other words, the orthoprojector ${\Pi }_{\star }$ is closest to $P_{\star }$ in the Schatten 1-norm, and the range of ${\Pi }_{\star }$ is the linear model we want. Fortunately, Lerman has given the error bound between the d-dimensional subspace with the d-dimensional orthoprojector ${\Pi }_{\star }$ in Theorem 2.1 [12] .

In this paper, we improve that the algorithm calls continuous iteratively reweighted least squares algorithm for solving REAPER (3), and under a weaker assumption on the data set, we can prove the algorithm is convergent. In the experiment part, we compare the algorithm with the IRLS algorithm [12] .

The rest of the paper is organized as follows. In Section 2, we develop the CIRLS algorithm for solving problem (3). We present a detail convergence analysis for the CIRLS algorithm in Section 3. An efficient numerical is reported in Section 4. Finally, we conclude this paper in Section 5.

2. The CIRLS Algorithm

We give the CIRLS algorithm for solving optimization problem (3), and the algorithm is summarized as Algorithm 1.

In the next section, we prove the sequence ${\left\{P^{\left(k\right)}\right\}}_{k\ge 1}$ is convergent to $P_{\delta }$. Furthermore, we also present the $P_{\delta }$ satisfied the bound with the optimal point $P_{\star }$ of REAPER.

3. Convergence of CIRLS Algorithm

In this section, we prove that the sequence ${\left\{P^{\left(k\right)}\right\}}_{k\ge 1}$ generated by Algorithm 1 is convergent to $P_{\delta }$ and we provide the $P_{\delta }$ satisfied the bound with the optimal point $P_{\star }$ of REAPER. Firstly, we start from the Lemma prepare for the proof of the following theorem.

Lemma 1 ( [12] , Theorem 4.1) Assume that the set $\chi$ of observations does not lie in the union of two strict subspaces of ${ℝ}^D$. Then the iterates of IRLS Algorithm with $\epsilon =0$ converge to a point $P_{\delta }$ that satisfies the constraints of the REAPER problem. Moreover, the objective value at $P_{\delta }$ satisfies the bound

${\displaystyle \underset{x\in \chi }{\sum }}\Vert x-P_{\delta }x\Vert -{\displaystyle \underset{x\in \chi }{\sum }}\Vert x-P_{\star }x\Vert \le \frac{1}{2}\delta \left|\chi \right|\mathrm{,}$ (5)

where $P_{\star }$ is an optimal point of REAPER.

In lemma 1, under the assumption that the set $\chi$ of observations does not lie in the union of two strict subspaces of ${ℝ}^D$, they consider the convergence of the IRLS algorithm. However, verify whether a data set satisfies this assumption requires amounts of computation in theory, so we give a assumption that is easier to verify in following theorem.

Theorem 1 Assume that the set $\chi$ of observations satisfies $span\left\{\chi \right\}={ℝ}^D$, if $P_{\delta }$ is the limit point of the sequence ${\left\{P^{\left(k\right)}\right\}}_{k\ge 1}$ generated by $CIRL{S}_{\delta }$ algorithm, then $P_{\delta }$ is an optimal point of the following optimization problem

$\begin{array}{l}\min \left({\displaystyle \underset{\left\{x\mathrm{<mo>\; :\; </mo>}\Vert x-Px\Vert \ge \delta \right\}}{\sum }}\Vert x-Px\Vert +\frac{1}{2}{\displaystyle \underset{\left\{x\mathrm{<mo>\; :\; </mo>}\Vert x-Px\Vert <\delta \right\}}{\sum }}\left\{\frac{{\Vert x-Px\Vert }^2}{\delta }+\delta \right\}\right)\\ \text{subjectto}\mathrm{<mi>\; </mi>\; <mi>\; </mi>0}\preccurlyeq P\preccurlyeq I\mathrm{<mi>\; </mi>\; <mi>\; </mi>}\text{and}\mathrm{<mi>\; </mi>\; <mi>\; </mi>}trP\mathrm{<mo>\; =\; </mo>}d\mathrm{,}\end{array}$ (6)

and satisfies

${\displaystyle \underset{x\in \chi }{\sum }}\Vert x-P_{\delta }x\Vert -{\displaystyle \underset{x\in \chi }{\sum }}\Vert x-P_{\star }x\Vert \le \frac{1}{2}\delta \left|\chi \right|\mathrm{,}$ (7)

where $P_{\star }$ is an optimal point of REAPER(3).

Proof. Firstly, we consider the optimization model with Algoritnm 1,

$\begin{array}{l}\min \left({\displaystyle \underset{\left\{x\mathrm{<mo>\; :\; </mo>}\Vert x-Px\Vert \ge \delta \right\}}{\sum }}\Vert x-Px\Vert +\frac{1}{2}{\displaystyle \underset{\left\{x\mathrm{<mo>\; :\; </mo>}\Vert x-Px\Vert <\delta \right\}}{\sum }}\left\{\frac{{\Vert x-Px\Vert }^2}{\delta }+\delta \right\}\right)\\ \text{subjectto}\mathrm{<mi>\; </mi>\; <mi>\; </mi>0}\preccurlyeq P\preccurlyeq I\mathrm{<mi>\; </mi>\; <mi>\; </mi>}\text{and}\mathrm{<mi>\; </mi>\; <mi>\; </mi>}trP=d\mathrm{,}\end{array}$ (8)

then, we define iterative point $P^{\left(k+1\right)}$ of the optimization problem

$P^{\left(k+1\right)}\mathrm{<mo>\; :\; </mo>}=\arg \min {\displaystyle \underset{x\in \chi }{\sum }}{\omega }_{{\delta }_k\mathrm{,}x}{\Vert x-Px\Vert }^2\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}\text{subjectto}\mathrm{<mi>\; </mi>\; <mi>\; </mi>0}\preccurlyeq P\preccurlyeq I\mathrm{<mi>\; </mi>\; <mi>\; </mi>}\text{and}\mathrm{<mi>\; </mi>\; <mi>\; </mi>}trP=d\mathrm{.}$ (9)

where ${\omega }_{{\delta }_k,x}:=\frac{1}{\max \left\{{\delta }_k,\Vert x-P^{\left(k\right)}x\Vert \right\}}$.

For convenience, let $Q=I-P$, and the optimization model(9) convert into

$Q^{\left(k+1\right)}\mathrm{<mo>\; :\; </mo>}=\underset{\begin{array}{c}0\preccurlyeq Q\preccurlyeq I\\ trQ=D-d\end{array}}{\arg \min }{\displaystyle \underset{x\in \chi }{\sum }}{\omega }_{{\delta }_k\mathrm{,}x}{\Vert Qx\Vert }^2\mathrm{.}$ (10)

where ${\omega }_{{\delta }_k\mathrm{,}x}\mathrm{<mo>\; :\; </mo>}=\frac{1}{\max \left\{{\delta }_k\mathrm{,}\Vert Q^{\left(k\right)}x\Vert \right\}}$.

Similarly, we convert the optimization model (8) into

$\begin{array}{l}\min F_{\delta }\left(Q\right)=\min \left({\displaystyle \underset{\left\{x\mathrm{<mo>\; :\; </mo>}\Vert Qx\Vert \ge \delta \right\}}{\sum }}\Vert Qx\Vert +\frac{1}{2}{\displaystyle \underset{\left\{x\mathrm{<mo>\; :\; </mo>}\Vert Qx\Vert <\delta \right\}}{\sum }}\left(\frac{{\Vert Qx\Vert }^2}{\delta }+\delta \right)\right)\\ \text{subjectto}\mathrm{<mi>\; </mi>\; <mi>\; </mi>0}\preccurlyeq Q\preccurlyeq I\mathrm{<mi>\; </mi>}\text{and}\mathrm{<mi>\; </mi>\; <mi>\; </mi>}trQ=D-d\mathrm{.}\end{array}$ (11)

Next we prove the convergence of $P^{\left(k\right)}$, consider the Huber-like function

$H_{{\delta }_k}\left(x,y\right)=\{\begin{array}{l}\frac{1}{2}\left(\frac{x^2}{{\delta }_k}+{\delta }_k\right)0\le y<{\delta }_k\\ \frac{1}{2}\left(\frac{x^2}{y}+y\right)y\ge {\delta }_k\end{array}$ (12)

since $\frac{1}{2}\left(\frac{x^2}{y}+y\right)\ge \left|x\right|$, for any $y>0$, then

$H_{{\delta }_k}\left(x,y\right)\ge H_{{\delta }_k}\left(x,\left|x\right|\right).$ (13)

holds for any $y>0$, $x\in R$. We introduce the convex function

$\begin{array}{c}{F}_{{\delta }_k}\left(Q\right):={\displaystyle \underset{x\in \chi }{\sum }}{H}_{{\delta }_k}\left(\Vert Qx\Vert ,\Vert Qx\Vert \right)\\ ={\displaystyle \underset{\left\{x\in \chi :\Vert Qx\Vert \ge {\delta }_k\right\}}{\sum }}\Vert Qx\Vert +\frac{1}{2}{\displaystyle \underset{\left\{x\in \chi :\Vert Qx\Vert <{\delta }_k\right\}}{\sum }}\left(\frac{{\Vert Qx\Vert }^2}{{\delta }_k}+{\delta }_k\right),\end{array}$ (14)

note that F is continuously differentiable at each matrix Q, and the gradient is

$\begin{array}{c}\nabla F_{{\delta }_k}\left(Q\right)={\displaystyle \underset{\left\{x\in \chi :\Vert Qx\Vert \ge {\delta }_k\right\}}{\sum }}\frac{Qx{x}^{\text{T}}}{\Vert Qx\Vert }+{\displaystyle \underset{\left\{x\in \chi :\Vert Qx\Vert <{\delta }_k\right\}}{\sum }}\frac{Qx{x}^{\text{T}}}{{\delta }_k}\\ ={\displaystyle \underset{x\in \chi }{\sum }}\frac{Qx{x}^{\text{T}}}{\max \left\{{\delta }_k,\Vert Qx\Vert \right\}}.\end{array}$ (15)

in addition, we introduce the function

$\begin{array}{l}{G}_{{\delta }_k}\left(Q,Q^{\left(k\right)}\right)\\ :={\displaystyle \underset{x\in \chi }{\sum }}{H}_{{\delta }_k}\left(\Vert Qx\Vert ,\Vert Q^{\left(k\right)}x\Vert \right)\\ ={\displaystyle \underset{\left\{x\in \chi :\Vert Q^{\left(k\right)}x\Vert \ge {\delta }_k\right\}}{\sum }}\frac{1}{2}\left(\frac{{\Vert Qx\Vert }^2}{\Vert Q^{\left(k\right)}x\Vert }+\Vert Q^{\left(k\right)}x\Vert \right)+{\displaystyle \underset{\left\{x\in \chi :\Vert Q^{\left(k\right)}x\Vert <{\delta }_k\right\}}{\sum }}\frac{1}{2}\left(\frac{{\Vert Qx\Vert }^2}{{\delta }_k}+{\delta }_k\right),\end{array}$ (16)

then

$\nabla G_{{\delta }_k}\left(Q\mathrm{,}{Q}^{\left(k\right)}\right)={\displaystyle \underset{x\in \chi }{\sum }}\frac{Qx{x}^{\text{T}}}{\max \left\{{\delta }_k\mathrm{,}\Vert Q^{\left(k\right)}x\Vert \right\}}\mathrm{.}$ (17)

By the definition of $G_{{\delta }_k}\left(Q\mathrm{,}{Q}^{\left(k\right)}\right)$ we know that

$G_{{\delta }_k}\left(Q^{\left(k\right)},Q^{\left(k\right)}\right)={\displaystyle \underset{x\in \chi }{\sum }}{H}_{{\delta }_k}\left(\Vert Q^{\left(k\right)}x\Vert ,\Vert Q^{\left(k\right)}x\Vert \right)=F_{{\delta }_k}\left(Q^{\left(k\right)}\right),$ (18)

and

$\nabla G_{{\delta }_k}\left(Q^{\left(k\right)},Q^{\left(k\right)}\right)={\displaystyle \underset{x\in \chi }{\sum }}\frac{Q^{\left(k\right)}x{x}^{\text{T}}}{\max \left\{{\delta }_k,\Vert Q^{\left(k\right)}x\Vert \right\}}=\nabla F_{{\delta }_k}\left(Q^{\left(k\right)}\right),$ (19)

it is obvious that $G_{{\delta }_k}\left(\mathrm{<mi>\; </mi>}\cdot \mathrm{<mi>\; </mi>,}{Q}^{\left(k\right)}\right)$ is a smooth quadratic function, we may relate $G_{{\delta }_k}\left(\mathrm{<mi>\; </mi>}\cdot \mathrm{<mi>\; </mi>,}{Q}^{\left(k\right)}\right)$ through the expansion in $Q^{\left(k\right)}$ as follows

$\begin{array}{c}{G}_{{\delta }_k}\left(Q,Q^{\left(k\right)}\right)=F_{{\delta }_k}\left(Q^{\left(k\right)}\right)+\langle Q-Q^{\left(k\right)},\nabla F_{{\delta }_k}\left(Q^{\left(k\right)}\right)\rangle \\ +\frac{1}{2}\langle Q-Q^{\left(k\right)},C_k\left(Q^{\left(k\right)}\right)\left(Q-Q^{\left(k\right)}\right)\rangle ,\end{array}$ (20)

where $C_k\left(Q^{\left(k\right)}\right)={\displaystyle \underset{x\in \chi }{\sum }}\frac{x{x}^{\text{T}}}{\max \left\{{\delta }_k\mathrm{,}\Vert Q^{\left(k\right)}x\Vert \right\}}$ is the Hessian matrix of $G_{{\delta }_k}\left(Q,Q^{\left(k\right)}\right)$.

By the definition of $G_{{\delta }_k}\left(Q,Q^{\left(k\right)}\right)$ we know that $F_{{\delta }_k}\left(Q\right)=G_{{\delta }_k}\left(Q,Q\right)$, combines with (13) we have

$\begin{array}{c}{F}_{{\delta }_k}\left(Q\right)={\displaystyle \underset{x\in \chi }{\sum }}{H}_{{\delta }_k}\left(\Vert Qx\Vert ,\Vert Qx\Vert \right)\\ \le {\displaystyle \underset{x\in \chi }{\sum }}{H}_{{\delta }_k}\left(\Vert Qx\Vert ,\Vert Q^{\left(k\right)}x\Vert \right)\\ =G_{{\delta }_k}\left(Q,Q^{\left(k\right)}\right),\end{array}$ (21)

and since the optimization model(10) is equivalent to the model

$Q^{\left(k+1\right)}\mathrm{<mo>\; :\; </mo>}=\underset{\begin{array}{c}0\preccurlyeq Q\preccurlyeq I\\ trP=D-d\end{array}}{\arg \min }{G}_{{\delta }_k}\left(Q\mathrm{,}{Q}^{\left(k\right)}\right)\mathrm{,}$ (22)

then we have the monotonicity property

$F_{{\delta }_k}\left(Q^{\left(k+1\right)}\right)\le G_{{\delta }_k}\left(Q^{\left(k+1\right)}\mathrm{,}{Q}^{\left(k\right)}\right)\le G_{{\delta }_k}\left(Q^{\left(k\right)}\mathrm{,}{Q}^{\left(k\right)}\right)=F_{{\delta }_k}\left(Q^{\left(k\right)}\right)\mathrm{.}$ (23)

We note that when $x\ge 0$, we have

$H_{{\delta }_k}\left(x,x\right)=\{\begin{array}{l}xx\ge {\delta }_k\\ \frac{1}{2}\left(\frac{x^2}{{\delta }_k}+{\delta }_k\right)0\le x<{\delta }_k\end{array}$ (24)

$H_{{\delta }_{k+1}}\left(x,x\right)=\{\begin{array}{l}xx\ge {\delta }_{k+1}\\ \frac{1}{2}\left(\frac{x^2}{{\delta }_{k+1}}+{\delta }_{k+1}\right)0\le x<{\delta }_{k+1}\end{array}$ (25)

since ${\delta }_{k+1}<{\delta }_k$, on the one hand, when $x\ge {\delta }_{k+1}$, $H_{{\delta }_{k+1}}\left(x,x\right)=x$ holds, and $H_{{\delta }_k}\left(x\mathrm{,}x\right)\ge x$ holds for any $x\in R$, therefore the inequality $H_{{\delta }_{k+1}}\left(x\mathrm{,}x\right)\le H_{{\delta }_k}\left(x\mathrm{,}x\right)$ holds for $x\ge {\delta }_{k+1}$. On the other hand, when $x<{\delta }_{k+1}$,

$H_{{\delta }_{k+1}}\left(x,x\right)=\frac{1}{2}\left(\frac{x^2}{{\delta }_{k+1}}+{\delta }_{k+1}\right)$, and $H_{{\delta }_k}\left(x,x\right)=\frac{1}{2}\left(\frac{x^2}{{\delta }_k}+{\delta }_k\right)$ holds, then we have

$\begin{array}{c}{H}_{{\delta }_k}\left(x,x\right)-H_{{\delta }_{k+1}}\left(x,x\right)=\frac{1}{2}\left(\frac{x^2}{{\delta }_k}-\frac{x^2}{{\delta }_{k+1}}+{\delta }_k-{\delta }_{k+1}\right)\\ =\frac{1}{2}\left({\delta }_k-{\delta }_{k+1}\right)\left(1-\frac{x^2}{{\delta }_k{\delta }_{k+1}}\right).\end{array}$ (26)

since ${\delta }_k>{\delta }_{k+1}$, and $x<{\delta }_{k+1}$, $H_{{\delta }_k}\left(x,x\right)-H_{{\delta }_{k+1}}\left(x,x\right)>0$ holds for all $x\ge 0$, therefore

$\begin{array}{c}{F}_{{\delta }_{k+1}}\left(Q^{\left(k+1\right)}\right)={\displaystyle \underset{x\in \chi }{\sum }}{H}_{{\delta }_{k+1}}\left(\Vert Q^{\left(k+1\right)}x\Vert ,\Vert Q^{\left(k+1\right)}x\Vert \right)\\ \le {\displaystyle \underset{x\in \chi }{\sum }}{H}_{{\delta }_k}\left(\Vert Q^{\left(k+1\right)}x\Vert ,\Vert Q^{\left(k+1\right)}x\Vert \right)\\ =F_{{\delta }_k}\left(Q^{\left(k+1\right)}\right),\end{array}$ (27)

and according to (23), we have the following result

$F_{{\delta }_{k+1}}\left(Q^{\left(k+1\right)}\right)\le F_{{\delta }_k}\left(Q^{\left(k\right)}\right)\mathrm{.}$ (28)

Since

$Q^{\left(k+1\right)}=\underset{\begin{array}{c}0\preccurlyeq Q\preccurlyeq I\\ trQ=D-d\end{array}}{\arg \min }{G}_{{\delta }_k}\left(Q\mathrm{,}{Q}^{\left(k\right)}\right)\mathrm{,}$ (29)

combined with the convex optimization variational inequalities with some constraints, we have

$0\le \langle Q-Q^{\left(k+1\right)}\mathrm{,}\nabla G_{{\delta }_k}\left(Q^{\left(k+1\right)}\mathrm{,}{Q}^{\left(k\right)}\right)\rangle \mathrm{<mi>\; </mi>\; <mi>\; </mi>}\forall 0\preccurlyeq Q\preccurlyeq I\mathrm{,\; <mi>\; </mi>}trQ=D-d\mathrm{,}$ (30)

base on (20), we have the equation

$\nabla G_{{\delta }_k}\left(Q,Q^{\left(k\right)}\right)=\nabla F_{{\delta }_k}\left(Q^{\left(k\right)}\right)+C_k\left(Q^{\left(k\right)}\right)\left(Q-Q^{\left(k\right)}\right)\mathrm{,}$ (31)

then we have

$0\le \langle Q^{\left(k\right)}-Q^{\left(k+1\right)}\mathrm{,}\nabla F_{{\delta }_k}\left(Q^{\left(k\right)}\right)+C_k\left(Q^{\left(k\right)}\right)\left(Q^{\left(k+1\right)}-Q^{\left(k\right)}\right)\rangle \mathrm{,}$ (32)

therefore

$\begin{array}{l}{G}_{{\delta }_k}\left(Q^{\left(k+1\right)}\mathrm{,}{Q}^{\left(k\right)}\right)\\ =F_{{\delta }_k}\left(Q^{\left(k\right)}\right)+\langle Q^{\left(k+1\right)}-Q^{\left(k\right)}\mathrm{,}\nabla F_{{\delta }_k}\left(Q^{\left(k\right)}\right)\rangle \\ +\frac{1}{2}\langle Q^{\left(k+1\right)}-Q^{\left(k\right)}\mathrm{,}{C}_k\left(Q^{\left(k\right)}\right)\left(Q^{\left(k+1\right)}-Q^{\left(k\right)}\right)\rangle \\ \le F_{{\delta }_k}\left(Q^{\left(k\right)}\right)-\langle Q^{\left(k+1\right)}-Q^{\left(k\right)}\mathrm{,}{C}_k\left(Q^{\left(k\right)}\right)\left(Q^{\left(k+1\right)}-Q^{\left(k\right)}\right)\rangle \end{array}$

$\begin{array}{l}+\frac{1}{2}\langle Q^{\left(k+1\right)}-Q^{\left(k\right)}\mathrm{,}{C}_{{\delta }_k}\left(Q^{\left(k\right)}\right)\left(Q^{\left(k+1\right)}-Q^{\left(k\right)}\right)\rangle \\ \le F_{{\delta }_k}\left(Q^{\left(k\right)}\right)-\frac{1}{2}\langle Q^{\left(k+1\right)}-Q^{\left(k\right)}\mathrm{,}{C}_k\left(Q^{\left(k\right)}\right)\left(Q^{\left(k+1\right)}-Q^{\left(k\right)}\right)\rangle \mathrm{,}\end{array}$ (33)

and according to(23), then we get the inequality

$F_{{\delta }_k}\left(Q^{\left(k+1\right)}\right)\le F_{{\delta }_k}\left(Q^{\left(k\right)}\right)-\frac{1}{2}\langle Q^{\left(k+1\right)}-Q^{\left(k\right)}\mathrm{,}{C}_k\left(Q^{\left(k\right)}\right)\left(Q^{\left(k+1\right)}-Q^{\left(k\right)}\right)\rangle \mathrm{,}$ (34)

add to $F_{{\delta }_{k+1}}\left(Q^{\left(k+1\right)}\right)\le F_{{\delta }_k}\left(Q^{\left(k+1\right)}\right)$, then we have

$\frac{1}{2}\langle Q^{\left(k+1\right)}-Q^{\left(k\right)}\mathrm{,}{C}_{{\delta }_k}\left(Q^{\left(k\right)}\right)\left(Q^{\left(k+1\right)}-Q^{\left(k\right)}\right)\rangle \le F_{{\delta }_k}\left(Q^{\left(k\right)}\right)-F_{{\delta }_{k+1}}\left(Q^{\left(k+1\right)}\right)\mathrm{.}$ (35)

Let $m:=\underset{k\ge 1}{\max }\left\{{\delta }_k,\underset{x\in \chi }{\max }\left\{\Vert x\Vert \right\}\right\}$, since $Q^{\left(k+1\right)}$, $Q^{\left(k\right)}$ and $C_k\left(Q^{\left(k\right)}\right)$ are all symmetric matrix, then the inequality

$\begin{array}{l}\langle Q^{\left(k+1\right)}-Q^{\left(k\right)}\mathrm{,}{C}_k\left(Q^{\left(k\right)}\right)\left(Q^{\left(k+1\right)}-Q^{\left(k\right)}\right)\rangle \\ \ge tr\left({\left(Q^{\left(k+1\right)}-Q^{\left(k\right)}\right)}^2{C}_k\left(Q^{\left(k\right)}\right)\right)\mathrm{,}\end{array}$ (36)

holds, in addition, $0\preccurlyeq Q^{\left(k+1\right)}\preccurlyeq I$, $0\preccurlyeq Q^{\left(k\right)}\preccurlyeq I$, thus ${\left(Q^{\left(k+1\right)}-Q^{\left(k\right)}\right)}^2$ is positive semidefinite matrix, and $\Vert Q^{\left(k\right)}x\Vert \le \Vert Q^{\left(k\right)}\Vert \cdot \Vert x\Vert \le \Vert x\Vert$, then we have $\max \left\{{\delta }_k\mathrm{,}\Vert Q^{\left(k\right)}x\Vert \right\}\le \max \left\{{\delta }_k\mathrm{,}\Vert x\Vert \right\}\le m$, therefore

$C_k\left(Q^{\left(k\right)}\right)={\displaystyle \underset{x\in \chi }{\sum }}\frac{x{x}^{\text{T}}}{\max \left\{{\delta }_k\mathrm{,}\Vert Q^{\left(k\right)}x\Vert \right\}}\succcurlyeq \frac{1}{m}{\displaystyle \underset{x\in \chi }{\sum }}x{x}^{\text{T}}\succcurlyeq \frac{1}{m}{\lambda }_{\min }\left({\displaystyle \underset{x\in \chi }{\sum }}x{x}^{\text{T}}\right)I\mathrm{,}$ (37)

so we have

${\left(Q^{\left(k+1\right)}-Q^{\left(k\right)}\right)}^2{C}_k\left(Q^{\left(k\right)}\right)\succcurlyeq \frac{1}{m}{\lambda }_{\min }\left({\displaystyle \underset{x\in \chi }{\sum }}x{x}^{\text{T}}\right){\left(Q^{\left(k+1\right)}-Q^{\left(k\right)}\right)}^2\mathrm{,}$ (38)

thus we can get the inequality as follows

$tr{\left(Q^{\left(k+1\right)}-Q^{\left(k\right)}\right)}^2{C}_k\left(Q^{\left(k\right)}\right)\ge \frac{1}{m}{\lambda }_{\min }\left({\displaystyle \underset{x\in \chi }{\sum }}x{x}^{\text{T}}\right){\Vert Q^{\left(k+1\right)}-Q^{\left(k\right)}\Vert }_F^2\mathrm{.}$ (39)

Let $\mu =\frac{1}{2m}{\lambda }_{\min }\left({\displaystyle \underset{x\in \chi }{\sum }}x{x}^{\text{T}}\right)$, then $\mu >0$, and we have

${\Vert Q^{\left(k+1\right)}-Q^{\left(k\right)}\Vert }_F^2\le \frac{1}{\mu }\left(F_{{\delta }_k}\left(Q^{\left(k\right)}\right)-F_{{\delta }_{k+1}}\left(Q^{\left(k+1\right)}\right)\right)\mathrm{,}$ (40)

since $F_{{\delta }_k}\left(Q\right)\ge 0$, so we have

${\displaystyle \underset{k=1}{\overset{\infty }{\sum }}}{\Vert Q^{\left(k+1\right)}-Q^{\left(k\right)}\Vert }_F^2\le \frac{1}{\mu }{F}_{{\delta }_1}\left(Q^{\left(1\right)}\right)<+\infty ,$ (41)

therefore, we get the following limits

$\underset{k\to \infty }{\lim }{\Vert Q^{\left(k+1\right)}-Q^{\left(k\right)}\Vert }_F=0.$ (42)

Since the set $\left\{Q\mathrm{<mo>\; :\; </mo>0}\preccurlyeq Q\preccurlyeq I\mathrm{,}trQ=D-d\right\}$ is bounded and closed convex set, and $Q^{\left(k+1\right)}\in \left\{Q\mathrm{<mo>\; :\; </mo>0}\preccurlyeq Q\preccurlyeq I\mathrm{,}trQ=D-d\mathrm{,}k=\mathrm{1,2,}\cdots \right\}$, then the sequence ${\left\{Q^{\left(k\right)}\right\}}_{k\ge 1}$ exist convergent subsequences, we suppose ${\left\{Q^{\left(k_i\right)}\right\}}_{i\ge 1}$ is one of the subsequences and

$\tilde{Q}=\underset{i\to +\infty }{\lim }\left\{Q^{\left(k_i\right)}\right\}\mathrm{.}$ (43)

Since for any $i\ge 1$, we have

$\begin{array}{c}{\Vert Q^{\left(k_i+1\right)}-\tilde{Q}\Vert }_F={\Vert Q^{\left(k_i+1\right)}-Q^{\left(k_i\right)}+Q^{\left(k_i\right)}-\tilde{Q}\Vert }_F\\ \le {\Vert Q^{\left(k_i+1\right)}-Q^{\left(k_i\right)}\Vert }_F+{\Vert Q^{\left(k_i\right)}-\tilde{Q}\Vert }_F\mathrm{,}\end{array}$ (44)

and

$\underset{k\to \infty }{\lim }{\Vert Q^{\left(k+1\right)}-Q^{\left(k\right)}\Vert }_F=0,$ (45)

so we have

$\underset{i\to +\infty }{\lim }{\Vert Q^{\left(k_i+1\right)}-\tilde{Q}\Vert }_F=0,$ (46)

in other words

$\underset{i\to +\infty }{\lim }{Q}^{\left(k_i+1\right)}=\tilde{Q}.$ (47)

According to the definition of $Q^{\left(k_i+1\right)}$ and (10), for any $Q\mathrm{<mo>\; :\; </mo>0}\preccurlyeq Q\preccurlyeq I\mathrm{,}trQ=D-d$, we have

$0\le \langle Q-Q^{\left(k_i+1\right)}\mathrm{,}\nabla F_{{\delta }_{k_i}}\left(Q^{\left(k_i\right)}\right)+C_{k_i}\left(Q^{\left(k_i\right)}\right)\left(Q^{\left(k_i+1\right)}-Q^{\left(k_i\right)}\right)\rangle \mathrm{.}$ (48)

Since

$\begin{array}{c}\underset{i\to +\infty }{\lim }\nabla F_{{\delta }_{k_i}}\left(Q^{\left(k_i\right)}\right)=\underset{i\to +\infty }{\lim }{\displaystyle \underset{x\in \chi }{\sum }}\frac{Q^{\left(k_i\right)}x{x}^{\text{T}}}{\max \left\{{\delta }_{k_i},\Vert Q^{\left(k_i\right)}x\Vert \right\}}\\ ={\displaystyle \underset{x\in \chi }{\sum }}\underset{i\to +\infty }{\lim }\frac{Q^{\left(k_i\right)}x{x}^{\text{T}}}{\max \left\{{\delta }_{k_i},\Vert Q^{\left(k_i\right)}x\Vert \right\}}\\ ={\displaystyle \underset{x\in \chi }{\sum }}\frac{\tilde{Q}x{x}^{\text{T}}}{\max \left\{\delta ,\Vert \tilde{Q}x\Vert \right\}}\\ =\nabla F_{\delta }\left(\tilde{Q}\right)\end{array}$ (49)

and

$\begin{array}{l}{\Vert C_{k_i}\left(Q^{\left(k_i\right)}\right)\left(Q^{\left(k_i+1\right)}-Q^{\left(k_i\right)}\right)\Vert }_F\\ ={\Vert {\displaystyle \underset{x\in \chi }{\sum }}\frac{x{x}^{\text{T}}}{\max \left\{{\delta }_{k_i},\Vert Q^{\left(k_i\right)}x\Vert \right\}}\left(Q^{\left(k_i+1\right)}-Q^{\left(k_i\right)}\right)\Vert }_F\\ \le {\displaystyle \underset{x\in \chi }{\sum }}\frac{{\Vert x{x}^{\text{T}}\Vert }_F}{\max \left\{{\delta }_{k_i},\Vert Q^{\left(k_i\right)}x\Vert \right\}}\cdot {\Vert Q^{\left(k_i+1\right)}-Q^{\left(k_i\right)}\Vert }_F\\ \le \left({\displaystyle \underset{x\in \chi }{\sum }}\frac{{\Vert x{x}^{\text{T}}\Vert }_F}{{\delta }_{k_i}}\right)\cdot {\Vert Q^{\left(k_i+1\right)}-Q^{\left(k_i\right)}\Vert }_F\\ \le \frac{1}{\delta }\left({\displaystyle \underset{x\in \chi }{\sum }}{\Vert x{x}^{\text{T}}\Vert }_F\right)\cdot {\Vert Q^{\left(k_i+1\right)}-Q^{\left(k_i\right)}\Vert }_F\end{array}$ (50)

as well as ${\Vert Q^{\left(k_i+1\right)}-Q^{\left(k_i\right)}\Vert }_F\to 0\left(i\to +\infty \right)$, so we have

$\underset{i\to +\infty }{\lim }{\Vert C_{k_i}\left(Q^{\left(k_i\right)}\right)\left(Q^{\left(k_i+1\right)}-Q^{\left(k_i\right)}\right)\Vert }_F=0,$ (51)

thus

$\underset{i\to +\infty }{\lim }{C}_{k_i}\left(Q^{\left(k_i\right)}\right)\left(Q^{\left(k_i+1\right)}-Q^{\left(k_i\right)}\right)=0.$ (52)

Taking the limit at both ends of the inequality (48) and using the continuity of the inner product, for any $Q\mathrm{<mo>\; :\; </mo>0}\preccurlyeq Q\preccurlyeq I\mathrm{,}trQ=D-d$, we have

$0\le \langle Q-\tilde{Q}\mathrm{,}\nabla F_{\delta }\left(\tilde{Q}\right)\rangle \mathrm{,}$ (53)

the variational inequalities demonstrate that

$\tilde{Q}=\arg \min F_{\delta }\left(Q\right)\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}\text{subjectto}\mathrm{<mi>\; </mi>\; <mi>\; </mi>0}\le Q\le I\mathrm{,\; <mi>\; </mi>\; <mi>\; </mi>}trQ=D-d\mathrm{,}$ (54)

it means that the limit points $Q_{\delta }$ of any convergent subsequence generated by ${\left\{Q^{\left(k\right)}\right\}}_{k\ge 1}$ is an optimal point of the following optimization problem

$Q_{\delta }=\arg \min F_{\delta }\left(Q\right)\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}\text{subjectto}\mathrm{<mi>\; </mi>\; <mi>\; </mi>0}\le Q\le I,\mathrm{<mi>\; </mi>\; <mi>\; </mi>}trQ=D-d.$ (55)

Now we set $P_{\delta }=I-Q_{\delta }$ and define $F_0\left(P\right)\mathrm{<mo>\; :\; </mo>}={\displaystyle \underset{x\in \chi }{\sum }}\Vert \left(I-P\right)x\Vert$ with $0\le P\le I$ and $trP=d$. And we define $F_{\delta }\left(P\right):={\displaystyle \underset{x\in \chi }{\sum }}{H}_{\delta }\left(\Vert \left(I-P\right)x\Vert ,\Vert \left(I-P\right)x\Vert \right)$, and $P_{\star }:=\arg \min F_0\left(P\right)$ with respect to the feasible set $0\le P\le I$, and $trP=d$, then we have

$P_{\delta }=\arg \min F_{\delta }\left(P\right),\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}\text{subjectto}\mathrm{<mi>\; </mi>\; <mi>\; </mi>0}\le P\le I\mathrm{<mi>\; </mi>\; <mi>\; </mi>},trQ=d.$ (56)

According to the lemma 3.1, we have

$0\le F_0\left(P_{\delta }\right)-F_0\left(P_{\star }\right)\le \frac{1}{2}\delta \left|\chi \right|\mathrm{,}$ (57)

where $\left|\chi \right|$ denotes the number of elements of $\chi$.

4. Numerical Experiments

In this section, we present a numerical experiment to show the efficiency of the CIRLS algorithm for solving problem (3). We compare the performance of our algorithm with IRLS on the data generated from the following model. In the test, we randomly choose Nin inliers sampled from the d-dimensional Multivariate Normal distribution $N\left(\mathrm{0,}{\Pi }_L\right)$ on subspace L and add Nout outliers sampled from a uniform distribution on ${\left[\mathrm{0,1}\right]}^D$, we also add a Gaussian noise $N\left(\mathrm{0,}{\epsilon }^{2}I\right)$. The experiment is performed in R and all the experimental results were averaged over 10 independent trials.

The parameters were set the same as [12] that $\epsilon ={10}^{-15}$, $\epsilon =0.01$ and

$\delta ={10}^{-10}$, we choose $\eta$ from $\left[\mathrm{0,1}\right]$ and we set $\eta =\frac{1}{2}$ in general. All the

other parameters of the two algorithms were set to be the same, Nout means the number of outliers and Nin is the number of inliers, D is the ambient dimension and d is the subspace dimension we want. From the model above, we get the data sets with values in the table below, then we calculate the iterative number

Table 1. The iterative number of the CIRLS and IRLS for different dimension.

$n\left(\text{iterative}\right)$ with the two algorithms through the data sets and the given parameters. The results are shown in Table 1.

We focus on the convergence speed of the two algorithms. Table 1 reports the numerical results of the two algorithms for different space dimension. From the result, we can see that in different dimension, the CIRLS algorithm performs better than IRLS algorithm in convergent efficiency.

5. Conclusion

In this paper, we propose an efficient continuous iteratively reweighted least squares algorithm for solving REAPER problem, and we prove the convergence of the algorithm. In addition, we present a bound between the convergent limit and the optimal point of REAPER problem. Moreover, in the experiment part, we compare the algorithm with the IRLS algorithm to show that our algorithm is convergent and performs better than IRLS algorithm in the rate of convergence.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Torre, F.D. and Black, M.J. (2001) Robust Principal Component Analysis for Computer Vision. Proceedings of Eighth IEEE International Conference on Computer Vision, Vancouver, 7-14 July 2001, 362-369. https://doi.org/10.1109/ICCV.2001.937541
[2]	Wanger, A., Wright, J. and Ganesh, A. (2009) Towards a Practical Face Recongnition System: Robust Registration Tillumination by Sparse Represention. 2009 IEEE Conference on Computer Vision and Pattern Recognition, Miami, FL, 20-25 June 2009, 20-25. https://doi.org/10.1109/CVPR.2009.5206654
[3]	Ho, J., Yang, M., Lim, J., Lee, K. and Kriegman, D. (2003) Clustering Appearances of Objects under Varying Illumination Conditions. 2003 IEEE Computer Vision and Pattern Recognition, Madison, WI, 18-20 June 2003, 11-18.
[4]	Jolliffe, I.T. (2002) Principal Component Analysis. Journal of Marketing Research, 25, 513.
[5]	Watson, G.A. (2002) On the Gauss-Newton Method for l1 Orthogonal Distance Regression. IMA Journal of Numerical Analysis, 22, 345-357. https://doi.org/10.1093/imanum/22.3.345
[6]	Ding, C., Zhou, D., He, X. and Zha, H. (2006) R1-PCA: Rotational Invariant L1-Norm Principal Component Analysis for Robust Subspace Factorization. The 23rd International Conference on Machine Learning, 23, 281-288. https://doi.org/10.1145/1143844.1143880
[7]	Kwak, N. (2008) Principal Component Analysis Based on L1-Norm Maximization. IEEE Transactions on Pattern Analysis and Machine Intelligence, Pittsburgh, PA, 25-29 June 2006, 1672-1680. https://doi.org/10.1109/TPAMI.2008.114
[8]	Maronna, R.A., Martin, D.R. and Yohai, V.J. (2006) Robust Statistics. Wiley Series in Probability and Statistics. Wiley, Chichester. https://doi.org/10.1002/0470010940
[9]	Zhang, T., Szlam, A. and Lerman, G. (2009) Median K-Flats for Hybrid Linear Modeling with Many Outliers. IEEE International Conference on Computer Vision Workshops, 12, 234-241.
[10]	Wright, J., Peng, Y.G. and Ma, Y. (2009) Robust Principal Component Analysis: Exact Recovery of Corrupted Low-Rank Matrices by Convex Optimization. 2009 23rd Annual Conference on Neural Information Processing Systems, Vancouver, Canada, 7-10 December 2009, 2080-2088.
[11]	Candes, E.J., Li, X., Ma, Y. and Wright, J. (2011) Robust Principal Component Analysis? Journal of the ACM, 58, Article No. 11. https://doi.org/10.1145/1970392.1970395
[12]	Lerman, G., McCoy, M.B., Tropp, J.T. and Zhang, T. (2015) Robust Computation of Linear Models by Convex Relaxation. Foundations of Computational Mathematics, 15, 363-410. https://doi.org/10.1007/s10208-014-9221-0