A semi-proximal-based strictly contractive Peaceman-Rachford splitting method\footnote{This research is supported by a Project Funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions.
Yan Gu\thanks{School of Mathematical Sciences
Key Laboratory for NSLSCS of Jiangsu Province
Nanjing Normal University
Nanjing 210023
P.R. China.
The Peaceman-Rachford splitting method is very efficient for minimizing sum of two functions each depends on its variable, and the constraint is a linear equality. However, its convergence was not guaranteed without extra requirements. Very recently, He {et al. (SIAM J. Optim. 24: 1011 - 1040, 2014) proved the convergence of a strictly contractive Peaceman-Rachford splitting method by employing a suitable underdetermined relaxation factor. In this paper, we further extend the so-called strictly contractive Peaceman-Rachford splitting method by using two different relaxation factors, and to make the method more flexible, we introduce semi-proximal terms to the subproblems. We characterize the relation of these two factors, and show that one factor is always underdetermined while the other one is allowed to be larger than 1. Such a flexible conditions makes it possible to cover the Glowinski's ADMM whith larger stepsize. We show that the proposed modified strictly contractive Peaceman-Rachford splitting method is convergent and also prove $O(1/t)$ convergence rate in ergodic and nonergodic sense, respectively. The numerical tests on an extensive collection of problems demonstrate the efficiency of the proposed method.
Modified-S-PRSM-v1_0607.tex
\documentclass[10pt]{article}
\usepackage{epsfig}
\usepackage{epstopdf}
\usepackage[titletoc,title]{appendix}
\usepackage{graphicx}
\usepackage{latexsym, cite, bm, amsthm, amsmath}
\usepackage{contour,soul,color}
\usepackage[colorlinks=true]{hyperref}
\usepackage{booktabs,threeparttable} % Table Rule and Table notes
\usepackage[rotateright]{rotating} % Table rotating
\usepackage{tabularx,multirow} % Table
\usepackage{cases}
\usepackage{float}
\usepackage{algorithm}
\usepackage{algorithmic}
\usepackage{mathtools}
\graphicspath{{figure/}}
\usepackage{caption}
\usepackage{diagbox}
\usepackage{lscape}
\usepackage{amsfonts}\usepackage{epstopdf}
\hypersetup{urlcolor=blue, citecolor=blue,linkcolor=blue}

\ExecuteOptions{dvips} \marginparwidth 0pt \oddsidemargin 0.5 truecm
\evensidemargin 1.5 truecm \marginparsep 0pt \topmargin -55pt
\textheight 24.0 truecm \textwidth 16.0 truecm


%\renewcommand{\thefootnote}{\fnsymbol{footnote}}
\renewcommand\baselinestretch{1.25}
\numberwithin{equation}{section}
\newtheorem{definition}{Definition}[section]
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}{Lemma}[section]
\newtheorem{eg}{Example}[section]
\newtheorem{proposition}{Proposition}[section]
\newtheorem{remark}{Remark} [section]
\newtheorem{corollary}{Corollary}[section]
\newtheorem{ass}{Assumption}[section]

\def\u{{\mathbf{u}}}
\def\IR{\mathcal{R}}
\def\X{\mathcal{X}}
\def\Y{{\cal Y}}
\def\Z{{\cal Z}}
\def\la{\lambda}
\def \[{\begin{equation}}
\def \]{\end{equation}}
\def \dist{\hbox{dist}}
\def\R{{\mathbb R}}
\def\B{{\cal B}}
\def\G{{\cal G}}
\def\L{{\cal L}}
\def\W{{\cal W}}
\def\S{{\cal S}}
\def\D{{\cal D}}
\def\N{{\cal N}}
\def\V{{\cal V}}
\def\U{{\cal U}}
\def\ba{\bar}
\def\st{\mathrm{s. t.}}
\def\u{{\overline{\mathbf{u}}}}
\def\v{{\overline{\mathbf{v}}}}
\def\w{{\overline{\mathbf{w}}}}
\def\nn{\nonumber}
\def\Tsf{\mathsf{T}}
\def\Fsf{\mathsf{F}}
\def\tr{\mathrm{tr}}
\def\llang{\left\langle}
\def\rrang{\right\rangle}
\def\Sigt{\Sigma_{1}}
\def\Sigtt{\Sigma_{2}}
\def\SigT{\Sigma}
\def\Pone{S}
\def\Ptwo{T}
\def\P{P}
\def\tauhat{\widehat{\tau}}
\def\Ghat{\widehat{G}}
\newcommand{\rev}[1]{{\color{black}{#1}}}

\begin{document}


\title{\bf A semi-proximal-based strictly contractive Peaceman-Rachford splitting method\footnote{This research is supported by a Project Funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions.}}

\author{Yan Gu\thanks{School of Mathematical Sciences, Key Laboratory for NSLSCS of Jiangsu Province,
Nanjing Normal University, Nanjing 210023, P.R. China.},
\; \; Bo Jiang$^{\dag}$\thanks{Email: \texttt{jiangbo@njnu.edu.cn}. }
\/ \hskip 0.3cm  and\hskip 0.3cm
Deren Han$^{\dag}$\thanks{Email: \texttt{handeren@njnu.edu.cn}. }}
\maketitle

\begin{abstract}
The Peaceman-Rachford splitting method is very efficient for  minimizing sum of two functions each depends on its variable, and the constraint is a linear equality. However, its convergence was not guaranteed without extra requirements. Very recently, He {\it et al.}  (SIAM J. Optim. 24: 1011 - 1040, 2014) proved the convergence of a strictly contractive Peaceman-Rachford splitting method  by employing a suitable underdetermined relaxation factor. In this paper, we further extend the so-called strictly contractive Peaceman-Rachford splitting method by using two different relaxation factors, and to make the method more flexible, we introduce semi-proximal terms to the subproblems. We characterize the relation of these two factors, and show that one factor is always underdetermined while the other one is allowed to be larger than 1. Such a flexible conditions makes it possible to cover the Glowinski's ADMM whith larger stepsize. We show that the proposed modified strictly contractive Peaceman-Rachford splitting method is convergent and also prove $O(1/t)$ convergence rate in ergodic and nonergodic sense, respectively.   The numerical tests on an extensive collection of problems demonstrate the efficiency of the proposed method.
\end{abstract}
\bigskip\noindent
{\bf Key words}: Semi-proximal, strictly contractive, Peaceman-Rachford splitting method, convex minimization,   convergence rate.


\section{Introduction}
We consider the convex minimization problem with linear constraints and a separable objective function
\[\label{cp}
 \min \,  \theta_1(x)+\theta_2(y), \quad \st \quad   Ax+By=b, \; x \in \X, \; y \in \Y,
\]
where $\theta_1:\;\R^{n_1}\to \R$ and $\theta_2:\;\R^{n_2}\to \R$ are continuous closed convex (could be nonsmooth) functions; $A\in\R^{m\times {n_1}}$ and $B\in\R^{m\times {n_2}}$ are given matrices; $b\in\R^m$ is a given vector; $\X$ and $\Y$ are nonempty closed convex subsets of $\R^{n_1}$ and $\R^{n_2}$, respectively. Throughout, the solution set of (\ref{cp}) is assumed to be nonempty; and $\X$ and $\Y$ are assumed to be simple in the sense that it is easy to compute the projections under the Euclidean norm onto them (e.g., positive orthant, spheroidal or box areas).


Let  $\L_{\beta}(x,y,\la)$ be the augmented Lagrangian function for \eqref{cp} that defined by
\[\label{augL}
\L_{\beta}(x,y,\la):=\theta_1(x)+\theta_2(y)-\langle\la,Ax+By-b\rangle+\frac{\beta}{2}\|Ax+By-b\|^2,\]
in which $\la\in\R^m$ is the multiplier associated to the linear constraint and  $\beta>0$ is a penalty parameter.
Based on the classic Douglas-Rachford operator splitting method \cite{douglas1956numerical}, the alternating direction method of multipliers was proposed by Gabay and Mercier \cite{gabay1976dual}, Glowinski and Marrocco \cite{glowinski1975approximation} in the mid-1970s, which generates the iterative sequence via the following recursion:
\begin{subnumcases}{}\label{alx}
x^{k+1}=\arg\min_{x\in\X}\L_{\beta}(x,y^k,\la^k),\\
\label{aly}y^{k+1}=\arg\min_{y\in\Y}\L_{\beta}(x^{k+1},y,\la^k),\\
\label{all}\la^{k+1}=\la^k-\beta(Ax^{k+1}+By^{k+1}-b).
\end{subnumcases}
Based on another classic operator splitting method, i.e., Peaceman-Rachford operator splitting method \cite{peaceman1955numerical}, one can derive the following method for \eqref{cp}:
\begin{subnumcases}{}\label{alxp}
x^{k+1}=\arg\min_{x\in\X}\L_{\beta}(x,y^k,\la^k),\\
\label{allb}\la^{k+\frac{1}{2}}=\la^k-\beta(Ax^{k+1}+By^{k+1}-b),\\
\label{alyp}y^{k+1}=\arg\min_{y\in\Y}\L_{\beta}(x^{k+1},y,\la^{k+\frac{1}{2}}),\\
\label{allp}\la^{k+1}=\la^{k+\frac{1}{2}}-\beta(Ax^{k+1}+By^{k+1}-b).
\end{subnumcases}
While the global convergence of the alternating direction method of multipliers \eqref{alx}-\eqref{all} can be established under very mild conditions \cite{boyd2011distributed}, the convergence of the Peaceman-Rachford-based method \eqref{alxp}-\eqref{allp} can not be guaranteed without further conditions \cite{corman2014generalized}. Most recently, He et al. \cite{he2014strictly} propose a modification of \eqref{alxp}-\eqref{allp} by introducing a parameter $\alpha$ to the update scheme of the dual variable $\lambda$ in \eqref{allb} and
\eqref{allp}, yielding the following procedure:
\begin{subnumcases}{}\label{he1}
x^{k+1}=\arg\min_{x\in\X}\L_{\beta}(x,y^k,\la^k),\\
\label{he2}\la^{k+\frac{1}{2}}=\la^k-\alpha\beta(Ax^{k+1}+By^{k+1}-b),\\
\label{he3}y^{k+1}=\arg\min_{y\in\Y}\L_{\beta}(x^{k+1},y,\la^{k+\frac{1}{2}}),\\
\label{he4}\la^{k+1}=\la^{k+\frac{1}{2}}-\alpha\beta(Ax^{k+1}+By^{k+1}-b).
\end{subnumcases}
Note that when $\alpha=1$, \eqref{he1}-\eqref{he4} is exactly the same as \eqref{alxp}-\eqref{allp}. They explained the nonconvergence behavior of  \eqref{alxp}-\eqref{allp} from the contract perspective, i.e., the distance from the iterative point to the solution set is merely nonexpansive, but not contractive. The parameter $\alpha$ in \eqref{he1}-\eqref{he4} plays the essential role in forcing the strict contractiveness of the generated sequence. Under the condition that $\alpha\in (0,1)$, they proved the same sublinear convergence rate as that for ADMM \cite{he2012convergence}. Particularly, they showed that \eqref{he1}-\eqref{he4}  achieves an approximate solution of \eqref{cp} with the accuracy of $O(1/t)$ after $t$ iterations\footnote{As the work \cite{nesterov1983method, nesterov2007gradient} and many others, a worst-case $O(1/t)$ convergence rate means the accuracy to a solution under certain criteria is of the order $O(1/t)$ after $t$ iterations of an iterative scheme; or equivalently, it requires at most $O(1/\epsilon)$ iterations to achieve an approximate solution with an accuracy of $\epsilon$.}, both in the ergodic sense and the nonergodic sense.

Note that the parameter $\alpha$ plays different roles in \eqref{allb} and \eqref{allp}: the former only affects the update of the variable $y$ in \eqref{aly} while the latter is for the update of the dual variable $\lambda$. Hence, it is natural to choose different parameters in these two equalities. In this paper, we  give such a scheme by introducing a new parameter $\gamma$ in \eqref{allp}, i.e., the dual variable is updated by the following manner:
\[\label{alph}
\la^{k+1}=\la^{k+\frac{1}{2}}-\gamma\beta(Ax^{k+1}+By^{k+1}-b).\]
For convenience, we \rev{first} introduce the whole update scheme of the {\it modified strictly contractive semi-proximal Peaceman-Rachford splitting method} (sP-PRSM) as
\begin{subnumcases}{\label{equ:sP-PRSM}}
\label{equ:sP-PRSM1} x^{k+1}=\arg\min_{x\in\X}\ \L_{\beta}(x,y^k,\la^k) + \frac12 \|x - x^{k}\|_{\Pone}^2,\\
\label{equ:sP-PRSM2} \la^{k+\frac{1}{2}}=\la^k-\alpha\beta(Ax^{k+1}+By^{k}-b),\\
\label{equ:sP-PRSM3} y^{k+1}=\arg\min_{y\in\Y}\ \L_{\beta}(x^{k+1},y,\la^{k+\frac{1}{2}}) + \frac12 \|y - y^{k}\|_{\Ptwo}^2,\\
\label{equ:sP-PRSM4} \la^{k+1}=\la^{k+\frac{1}{2}}-\gamma\beta(Ax^{k+1}+By^{k+1}-b)
\end{subnumcases}
where $S$ and $T$ are two positive semi-definite matrices. In applications, by choosing different
matrices $S$ and $T$ customizing the problems' structures, we can obtain different efficient methods.

Our main contributions are  \rev{twofold.}

\begin{enumerate}
\item Motivated by the nice analysis techniques in \cite{he2014strictly} and \cite{xu2011class},  we proved that the sequence generated by sP-PRSM is  strictly contractive and thus convergent, under the requirement that
\[\label{equ:alpha:gamma}
\alpha\in [0,\,1),\quad \gamma\in \left(0,\,\frac{1 -\alpha  + \sqrt{ (1 + \alpha)^2 + 4(1 - \alpha^2)}}{2}\right).\]  Moreover, we proved that sP-PRSM is  $O(1/t)$ sublinearly convergent both in the ergodic and nonergodic sense. Note that the nonergodic convergence rate requires that $\gamma \in (0,1]$. We remark that the convergence of sP-PRSM \eqref{equ:sP-PRSM} can unify that of several existing splitting methods.
\begin{itemize}
\item Choosing $\alpha = 0,  \gamma \in (0, (1 + \sqrt{5})/2)$ and $T = 0$, $S = 0$, sP-PRSM \eqref{equ:sP-PRSM} reduces to the classical ADMM  while the convergence coincides with that of ADMM \cite{glowinski1984numerical};
\item Setting $\alpha = 0, \gamma \in (0, (1+ \sqrt{5})/2)$ and $T \succeq 0, S \succeq 0$, sP-PRSM \eqref{equ:sP-PRSM} covers the semi-proximal ADMM  considered in \cite{deng2012global, he2012convergence, fazel2013hankel, xu2011class} and the corresponding convergence results;
\item Setting $\alpha = \gamma \in (0,1)$, and $ T = 0$, $S = 0$, sP-PRSM \eqref{equ:sP-PRSM} reduces to the strictly contractive PRSM proposed in \cite{he2014strictly} and the convergence of the two methods is identical.
\end{itemize}

\item We added a proximal term to each of the two main subproblems in updating the $x$ and the $y$ variables.
In fact, for ADMM, Eckstein \cite{Eck1994} and He et al. \cite{HLHY2002} have already considered to add proximal terms to the subproblems for different purpose.  Recently, Fazel et al. \cite{fazel2013hankel}
proposed to allow $S$ and $T$ to be positive semi-definite, in contrast to the positive definite requirements in the classical algorithms, which makes the algorithm more flexible.
\end{enumerate}

The rest of this paper is organized as follows. In Sect.  \ref{section:preliminaries}, we give the optimality condition of \eqref{cp} by using the variational inequality and also list some assertions which will be used in later analysis. In Sect. \ref{section:convergence}, we first give the contraction analysis of  sP-PRSM \eqref{equ:sP-PRSM}, and then \rev{establish the global convergence.  We discuss the sublinear and linear convergence rate in Sects. \ref{section:sublinear} and \ref{section:linear}, respectively. } In Sect. \ref{section:numerical}, we test a variety of problems to show the efficiency of the proposed sP-PRSM \eqref{equ:sP-PRSM}. Finally, we make some conclusions in Sect. \ref{section:conclusions}.

On May 25, after attending the first author's thesis defense, Prof. He told us that they had also
%During writing the final version of this paper, we found that \cite{he2015step} also
considered to use different stepsizes $\alpha$ and $\gamma$ in \eqref{equ:sP-PRSM} with $T = 0$ and $S = 0$, and we got their manuscript on May 26 \cite{he2015step}. We found that they established the convergence when
\[ \label{equ:hmy}
  \gamma \in \left(0, \frac{1 + \sqrt{5}}{2}\right), \ \alpha \in (-1, 1),\   \alpha + \gamma >0, \  |\alpha| < 1 + \gamma - \gamma^2.
\]
Notice that \eqref{equ:hmy} covers the case when $\alpha <0$. If we restrict $\alpha \geq 0$ in \eqref{equ:hmy}, then the corresponding relation will be
\[ \label{equ:hmy2}
  \alpha \in [0,1), \quad \gamma \in \left(0, \frac{1+ \sqrt{5 - 4\alpha}}{2}\right).
\]
By some simple calculations, it is easy to see that the domain defined by \eqref{equ:hmy2} is a bit smaller than that defined by \eqref{equ:alpha:gamma}.  By some private communications with Prof. He, we learned that they also discovered the formula \eqref{equ:alpha:gamma}. However, to provide the intuitive understanding and the unified convergence analysis, they presented \eqref{equ:hmy}  in \cite{he2015step} which connects $\alpha$ and $1 - \gamma + \gamma^2$ directly.
\section{Preliminaries}\label{section:preliminaries}
%Let
%\[\label{m}
%M=\left(\begin{array}{cc} I_{n_2}&0\\ -\alpha\beta B& (\alpha+\gamma)I_m\end{array}\right)\]
%and
%\[\label{q0}
%Q_0=\left(\begin{array}{cc} 0&0\\ \beta B^TB&-\alpha B^T\\ - B& \frac{1}{\beta}I_m\end{array}\right).\]\
In this section, we give the optimality condition of \eqref{cp} and  some notations or relations which will be frequently used in our analysis. Let $\Omega = \X \times \Y \times \R^m$. Throughout this paper, we make the following assumption.
%\begin{ass}\label{ass:B}
%$B$ is of full column rank. ?????
%\end{ass}
\begin{ass}\label{ass:opt} Let $\Omega^*\subset \Omega$ be the set whose elements are the  optimal solutions  of \eqref{cp} and the associating dual solutions of \eqref{cp}. Throughout the paper, we assume that
$\Omega^*$ is non-empty.
\end{ass}

\subsection{Optimality condition of \eqref{cp}}

Owing to the convexity of $\theta_1(\cdot) $ and $\theta(\cdot)$,  there exist two positive semidefinite matrices  $\Sigt$ and $\Sigtt$ such that for any $x, x' \in \R^{n_1}$
\[ \label{equ:x:convex}
  \langle x - x',  \xi_x - \xi_{x}' \rangle  \geq \|x - x'\|_{\Sigt}^2,
\]
where $\xi_x \in \partial \theta_1(x)$, $\xi_x' \in \partial \theta_1(x')$, and for any $y, y' \in \R^{n_2}$
\[\label{equ:y:convex}
\langle y - y', \xi_y - \xi_y' \rangle \geq \|y - y' \|_{\Sigtt}^2,
\]
where $\xi_y \in \partial \theta_2(y)$, $\xi_y' \in \partial \theta_1(y')$.

Denote  $$u = \left(\begin{array}{c}x\\y\end{array}\right),\; v = \left(\begin{array}{c}y\\\lambda\end{array}\right),\; w = \left(\begin{array}{c}x\\ y\\\lambda\end{array}\right),\; \hbox{and} \; F(w) = \left(\begin{array}{c}
\xi_x -A^{\Tsf} \lambda \\ \xi_y -B^{\Tsf} \lambda \\ Ax + By - b \end{array} \right).$$  Due to the convexity of $\theta_1(\cdot)$ and $\theta_2(\cdot)$,  it is easy to show that the operator $F(\cdot)$ is monotonic. Specifically, for any $w, w' \in \Omega$, we have
\begin{equation}\label{equ:F:mono}
\langle w - w', F(w) - F(w') \rangle = \left \langle \left( \begin{array}{c} x - x' \\ y - y' \end{array} \right),  \left( \begin{array}{c} \xi_x - \xi_x' \\ \xi_y - \xi_y' \end{array} \right) \right\rangle \geq \|u - u'\|_{\SigT}^2,
\end{equation}
where $\SigT = \left( \begin{array}{cc} \Sigt & 0 \\ 0 & \Sigtt \end{array}\right)$
and  the inequality is due to \eqref{equ:x:convex} and \eqref{equ:y:convex}.
  As shown in \cite{he2014strictly}, we say that $w^* \in \Omega^*$ if  there holds that
\[\label{equ:opt:cond}
 \llang w - w^*, F(w^*) \rrang \geq 0, \quad \forall w \in \Omega.
 %\theta(u) - \theta(u^*) + (w - w^*)^{\Tsf} F(w^*)
\]
From Theorem 2.3.5 in \cite{facchinei2003finite} or Theorem 2.1 in  \cite{he2012convergence}, we can see that $\Omega^*$ is closed and convex, and it can be reformulated as
\[
  \Omega^* = \bigcap_{w \in \Omega} \{\tilde{w} \in \Omega \colon \llang w - \tilde{w},  F(w) \rrang \geq 0\}.\label{equ:opt:cond2}
  \]
  Let $\S$ denote the feasible set of \eqref{cp}, namely, $$\S = \{(x,y) : Ax + By = b, x \in \X, y \in \Y\}$$ and denote $\D = \S \times \R^m$. Let us restrict $w \in \Omega$ in \eqref{equ:opt:cond} and \eqref{equ:opt:cond2} to be $w \in \D \subset \Omega.$  By some suitable modification of the proof of Theorem 2.1 in \cite{he2012convergence}, it is easy to show that the optimal solution set $\Omega^*$ can be characterized as
  \[
\Omega^* = \bigcap_{w \in \D} \{\tilde{w} \in \Omega \colon  \llang w - \tilde{w},  F(w)  \rrang \geq 0\}.\label{equ:opt:cond2:D}
    \]
    Similarly, we say that $w^* \in \Omega^*$ if there holds that
    \[
 \llang w - w^*, F(w^*) \rrang \geq 0, \quad \forall w \in \D. \label{equ:opt:cond:D}
\]
   Note that the optimality condition \eqref{equ:opt:cond:D} will be frequently used in the convergence analysis. Owing to \eqref{equ:opt:cond:D},  we call $\tilde{w}$ to be an $\epsilon$ solution of \eqref{cp} if
    \[ \label{equ:epsilon_sol}
      \sup_{w \in \D}\, \llang  \tilde{w} - w, F(w)  \rrang\leq \epsilon.
      \]
\subsection{Some notations}
Given two  real matrices $C$ and $D$ of the same dimension, we define $\langle C, D \rangle = \tr(C^{\Tsf}D)$, where $\tr(\cdot)$ is the trace operator. We use $\|\cdot \|$ to denote the 2-norm of a vector. We denote  $\|z\|_G^2 = z^{\Tsf} Gz$  for $z \in \R^n$ and $G \in \R^{n \times n}$. For a real symmetric matrix ${\Pone}$, we mark ${\Pone}\succeq 0$ (${\Pone}\succ 0$) if ${\Pone}$ is positive semidefinite (positive definite). To make the analysis more elegant, we use $r^k = Ax^k + By^k - b$ for short.  Similarly, for any $w \in \D$, we denote $r(w) = Ax + By - b$. Obviously, there holds that $r(w) = 0$ for any $w \in \D$.  For ease of the analysis, we define the following matrices as

\begin{align} \label{equ:H}
  H =\frac{1}{\alpha + \gamma}\left( \begin{array}{cc} (\alpha + \gamma - \alpha \gamma) \beta B^{\Tsf} B & -\alpha B^{\Tsf} \\ -\alpha B & \frac{1}{\beta} I_m\end{array}\right) \quad
%  H_0 =\left( \begin{array}{cc} 0 & 0 \\ 0& H\end{array}\right),\nn
\end{align}
and
\[\label{equ:M}
M=\left(\begin{array}{cc} I_{n_2}&0\\ \alpha\beta B& (\alpha+\gamma)\beta I_m\end{array}\right).\]
We can easily verify that $H$ is positive semidefinite if $0 \le \alpha \leq 1$ and $\gamma > 0$, and
\begin{align}
M^{\Tsf} H M =
%{} & \frac{1}{\alpha + \gamma}\left(\begin{array}{cc} I_{n_2}& \alpha\beta B^{\Tsf} \\ 0& (\alpha+\gamma)\beta I_m\end{array}\right)
%\cdot\left( \begin{array}{cc} (\alpha + \gamma - \alpha \gamma) \beta B^{\Tsf} B & -\alpha B^{\Tsf} \\ -\alpha B & \frac{1}{\beta} I_m\end{array}\right)
%\cdot \left(\begin{array}{cc} I_{n_2}&0\\ \alpha\beta B& (\alpha+\gamma)\beta I_m\end{array}\right) \nn \\
%={} & \left(\begin{array}{cc} (1 - \alpha)\beta B^{\Tsf}B& 0 \\ -\alpha \beta B & I_m\end{array}\right)
%\cdot \left(\begin{array}{cc} I_{n_2}&0\\ \alpha\beta B& (\alpha+\gamma)\beta I_m\end{array}\right)
\left(\begin{array}{cc} (1 - \alpha)\beta B^{\Tsf}B& 0 \\ 0 & (\alpha+\gamma)\beta I_m\end{array}\right).  \label{equ:MtHM}
\end{align}
Denote $P = \left(\begin{array}{cc} {\Pone} & 0 \\ 0 & {\Ptwo}\end{array}\right)$  and define \begin{align} \label{equ:G}
  G \coloneqq \left(\begin{array}{cc} P & 0 \\ 0 & 0  \end{array}\right) + \left(\begin{array}{cc} 0 & 0 \\ 0 & H \end{array} \right) =
   \left( \begin{array}{ccc} {\Pone}  &  0 & 0 \\
  0 &  {\Ptwo} + \frac{\alpha + \gamma - \alpha \gamma}{\alpha + \gamma} \beta B^{\Tsf} B & -\frac{\alpha}{\alpha + \gamma} B^{\Tsf} \\
  0  & -\frac{\alpha}{\alpha + \gamma} B & \frac{1}{(\alpha + \gamma) \beta} I_m\end{array}\right)
%  H_0 =\left( \begin{array}{cc} 0 & 0 \\ 0& H\end{array}\right),\nn
\end{align}
and
\begin{align}\label{equ:Ghat}
\Ghat \coloneqq \left(\begin{array}{cc} \Sigma & 0 \\ 0 & 0  \end{array}\right) + G  =  \left( \begin{array}{ccc} {\Pone} + \Sigt   &  0 & 0 \\
  0 &  {\Ptwo} + \Sigtt + \frac{\alpha + \gamma - \alpha \gamma}{\alpha + \gamma} \beta B^{\Tsf} B & -\frac{\alpha}{\alpha + \gamma} B^{\Tsf} \\
  0  & -\frac{\alpha}{\alpha + \gamma} B & \frac{1}{(\alpha + \gamma) \beta} I_m\end{array}\right),
\end{align}
where the dimension of zero matrix can be easily identified from the context.
  For any $w, w' \in \Omega$, there hold that
\[ \label{equ:GH_norm}
\| w - w'\|_G^2 = \|u - u'\|_P^2 + \|v - v'\|_H^2
\]
and
\[\label{equ:GhatH_norm}
\| w - w'\|_{\Ghat}^2 = \|u - u'\|_{\Sigma}^2 + \|w - w'\|_G^2.
\]

Finally, we present some relations of the iterates. They will play a crucial role in the convergence analysis.  With the update scheme \eqref{equ:sP-PRSM2} and \eqref{equ:sP-PRSM4}, it is easy to have
\[ \label{equ:lak:lak1}
\la^k = \la^{k+1} + (\alpha + \gamma)\beta r^{k+1} + \alpha \beta B(y^k - y^{k+1}).
\]
Recalling the definition of $v$, with \eqref{equ:lak:lak1} and \eqref{equ:M}, we conclude that
\begin{align}\label{equ:vk:vk1}
  v^k - v^{k+1} ={} &\left(\begin{array}{c} y^k - y^{k+1} \\ \lambda^k - \lambda^{k+1} \end{array}\right)
%&  \left(\begin{array}{c} y^k - y^{k+1} \\ \alpha \beta B(y^{k} - y^{k+1}) + (\alpha + \gamma)\beta (Ax^{k+1} + B y^{k+1} - b) \end{array}\right)
= M \left(\begin{array}{c} y^k - y^{k+1} \\ r^{k+1} \end{array}\right).
\end{align}
With \eqref{equ:MtHM} and \eqref{equ:vk:vk1}, we thus have
\begin{align}
\|v^k - v^{k+1}\|_H^2 %={} & (v^k - v^{k+1})^{\Tsf} (M^{\Tsf} H M)(v^k - v^{k+1}) \nn \\
  ={} &  (1- \alpha) \beta  \|B(y^k - y^{k+1})\|^2 + (\alpha + \gamma) \beta \|r^{k+1}\|^2. \label{equ:lemma:vkvk1:b1}
\end{align}

\section{Convergence of sP-PRSM }\label{section:convergence}
In this section, we first show that the sequence $\{w_k\}$ generated by sP-PRSM \eqref{equ:sP-PRSM} is strictly contractive and thus can establish the convergence of the method. With the help of the contraction property, we further discuss the convergence rate in the ergodic and nonergodic sense.
\subsection{Contraction analysis}\label{subsection:contraction}
We first give the contraction result as follows, and then complete its proof at the end of this section. Note that the following result will play a key role in proving the convergence of \eqref{equ:sP-PRSM}.
%Now, we are ready to show the strictly contractive property of the sequence $\{w^k\}$.
\begin{theorem}\label{thm:wk:wk1:contraction}
Let the sequence $\{w^k\}$ be generated by sP-PRSM \eqref{equ:sP-PRSM}. If we choose $\alpha$ and $\gamma$ according to \eqref{equ:alpha:gamma}, then  for any $w^* \in W^*$, there holds that
\begin{align}\label{equ:wk:wk1:w*}
  &\left( \|w^k - w^*\|_{\Ghat}^2 + \rho_{\alpha,\gamma} \|r^k\|^2 + \eta_{\alpha,\gamma} \|y^{k} - y^{k-1}\|_{\Ptwo}^2\right)\nn \\
  -&\left( \|w^{k+1} - w^*\|_{\Ghat}^2 + \rho_{\alpha,\gamma} \|r^{k+1}\|^2 + \eta_{\alpha,\gamma} \|y^{k+1} - y^{k}\|_{\Ptwo}^2\right)   \geq{}  \widehat{\tau}_{\alpha, \gamma}\|w^{k} - w^{k+1}\|_{\Ghat}^2,
  \end{align}
where $\tauhat_{\alpha, \gamma} = \min(\frac12, \tau_{\alpha,\gamma})$, and $\rho_{\alpha, \gamma}, \eta_{\alpha, \gamma} \geq 0$ and $0< \tau_{\alpha, \gamma} < 1$ are some constants  which only depend on $\alpha$ and $\gamma$.
\end{theorem}

To prove Theorem \ref{thm:wk:wk1:contraction}, several lemmas  should be established first.  The following Lemma \ref{lemma:opt} is mainly based on the optimality conditions of \eqref{equ:sP-PRSM1}  and \eqref{equ:sP-PRSM2}.


\begin{lemma} \label{lemma:opt}
Let the sequence $\{w^k\}$ be generated by sP-PRSM \eqref{equ:sP-PRSM}. If $\alpha \geq 0$ and $\gamma >0$, then for any $w \in \D$, there holds that
\begin{align}\label{equ:lemma:opt:0}
(w - w^{k+1})^{\Tsf}  G (w^{k+1}-w^{k})
\ge{}  &(1 - \alpha - \gamma) \beta\left\|r^{k+1} \right\|^2  + (1 - \alpha)\beta \llang r^{k+1}, B(y^k - y^{k+1})\rrang \nn \\
 & + \llang w^{k+1} - w, F(w^{k+1}) \rrang.
\end{align}
%\begin{align}
%\left( \begin{array}{c} x - x^{k+1} \\ y - y^{k+1} \\ \lambda - \lambda^{k+1} \end{array} \right)^{\Tsf}
%\left(
%\begin{array}{c}
%(1 - \alpha - \gamma) \beta  A^{\Tsf} (r^{k+1}) + (1 - \alpha) \beta  A^{\Tsf} B (y^k -y^{k+1} ) \\
%(1 - \gamma) \beta  B^{\Tsf} (r^{k+1})\\
%- \frac{\alpha}{\alpha + \gamma} B(y^{k+1} - y^{k}) + \frac{1}{(\alpha + \gamma)\beta} (\la^{k+1} - \la^k)
%\end{array}
%\right) \geq 0
%\end{align
\end{lemma}
\begin{proof}
Consider any $\xi_x^{k+1} \in \partial \theta_1(x^{k+1})$ and $\xi_y^{k+1} \in \partial \theta_2(y^{k+1})$. We then obtain respectively from the optimality conditions for \eqref{equ:sP-PRSM1} and \eqref{equ:sP-PRSM3} that
\[
\llang x - x^{k+1},  {\Pone}(x^{k+1} - x^k) + \xi_x^{k+1}-A^{\Tsf}\la^k + \beta  A^{\Tsf} r^{k+1} + \beta A^{\Tsf}B (y^k - y^{k+1}) \rrang \geq 0, \quad \forall x \in \X\nn
\]
and
\[
\llang y - y^{k+1}, {\Ptwo}(y^{k+1} - y^k) +  \xi_y^{k+1}-B^{\Tsf}\la^{k + \frac12} + \beta  B^{\Tsf} r^{k+1} \rrang \geq 0, \quad \forall y \in \Y. \nn
\]
Substituting \eqref{equ:lak:lak1}  and  $\la^{k+\frac12} = \la^{k+1} + \gamma\beta r^{k+1}$ which follows from \eqref{equ:sP-PRSM4} into the above two inequalities, respectively, we obtain the new formulations of the optimality conditions for \eqref{equ:sP-PRSM1} and \eqref{equ:sP-PRSM3} as
\begin{align}
\left \langle x - x^{k+1}, {\Pone}(x^{k+1} - x^k)  + \xi_x^{k+1} -A^{\Tsf}\la^{k+1} + (1 - \alpha - \gamma) \beta  A^{\Tsf} r^{k+1} + (1 - \alpha) \beta  A^{\Tsf} B (y^k -y^{k+1} )\right\rangle & \geq 0, \nn \\
 \forall x \in \X. &\label{equ:opt:x}
\end{align}
and
\[
\llang y - y^{k+1}, {\Ptwo}(y^{k+1} - y^k) +   \xi_y^{k+1} -B^{\Tsf}\la^{k + 1} + (1 - \gamma) \beta  B^{\Tsf} r^{k+1}\rrang  \geq 0, \quad \forall y \in \Y.  \label{equ:opt:y}
\]
Rewriting \eqref{equ:lak:lak1} to be
\[
r^{k+1} - \frac{\alpha}{\alpha + \gamma} B(y^{k+1} - y^{k}) + \frac{1}{(\alpha + \gamma)\beta} (\la^{k+1} - \la^k) = 0. \label{equ:opt:la}
\]
Combing \eqref{equ:opt:x}, \eqref{equ:opt:y} and \eqref{equ:opt:la} in a suitable way, and recalling the definition of $F(w)$, we see that  for any $w \in \Omega$ there holds that
\begin{align}\label{equ:lemma:opt:d1}
&
 \llang w - w^{k+1},
\left(
\begin{array}{c}
{\Pone}(x^{k+1} - x^k) \\
{\Ptwo}(y^{k+1} - y^k) \\
0
\end{array}
\right)
+
\left(
\begin{array}{c}
0 \\
\alpha \beta B^{\Tsf}r^{k+1} + (1 - \alpha) \beta  B^{\Tsf} B (y^{k+1} - y^{k})\\
- \frac{\alpha}{\alpha + \gamma} B(y^{k+1} - y^{k}) + \frac{1}{(\alpha + \gamma)\beta} (\la^{k+1} - \la^k)
\end{array}
\right)
\rrang
\nn \\
& \geq
\left\langle w^{k+1} - w,
\left(
\begin{array}{c}
  A^{\Tsf} \\
  B^{\Tsf} \\
0
\end{array}
\right)
\left[(1 - \alpha - \gamma) \beta   r^{k+1} + (1 - \alpha) \beta   B (y^k -y^{k+1} ) \right]
\rrang
+ \llang
w^{k+1} - w,
F(w^{k+1})
\rrang.
%\left(
%\begin{array}{c}
%(1 - \alpha - \gamma) \beta  A^{\Tsf} r^{k+1} + (1 - \alpha) \beta  A^{\Tsf} B (y^k -y^{k+1} ) \\
%(1 - \alpha - \gamma) \beta  B^{\Tsf} r^{k+1} + (1 - \alpha) \beta  B^{\Tsf} B (y^k -y^{k+1} ) \\
%0
%\end{array}
%\right). \nn
\end{align}
With the relation \eqref{equ:opt:la} and the definition \eqref{equ:H} of $H$, by some easy calculations, we obtain that
\[\label{equ:lemma:opt:d2}
\left(\begin{array}{c}
\alpha \beta B^{\Tsf}r^{k+1} + (1 - \alpha) \beta  B^{\Tsf} B (y^{k+1} - y^{k})\\
- \frac{\alpha}{\alpha + \gamma} B(y^{k+1} - y^{k}) + \frac{1}{(\alpha + \gamma)\beta} (\la^{k+1} - \la^k)
\end{array}\right)
= H (v^{k+1} - v^k).
\]
Considering the definition of $r^{k+1}$ and $r(w)$, we have $(A, \,B, \, 0) (w^{k+1} - w)  = r^{k+1} - r(w)$. This, together with \eqref{equ:lemma:opt:d1} and \eqref{equ:lemma:opt:d2} indicates  that
\begin{align}\label{equ:lemma:opt:01}
(w - w^{k+1})^{\Tsf}  G (w^{k+1}-w^{k})
\ge{} &\llang r^{k+1} -r(w), (1 - \alpha - \gamma) \beta   r^{k+1} + (1 - \alpha) \beta   B (y^k -y^{k+1} ) \rrang \nn \\ &+ \llang w^{k+1} - w, F(w^{k+1}) \rrang.
\end{align}
Notice that  $r(w) = 0$ for any $w \in \D$,  we can immediately obtain \eqref{equ:lemma:opt:0}  from \eqref{equ:lemma:opt:01} since The proof is completed.
\end{proof}

%Further,
%\begin{align}
%& (w - w^{k+1})^{\Tsf}
%\left(
%\begin{array}{c}
%(1 - \alpha - \gamma) \beta  A^{\Tsf} (r^{k+1}) + (1 - \alpha) \beta  A^{\Tsf} B (y^k -y^{k+1} ) \\
%(1 - \gamma) \beta  B^{\Tsf} (r^{k+1})\\
%- \frac{\alpha}{\alpha + \gamma} B(y^{k+1} - y^{k}) + \frac{1}{(\alpha + \gamma)\beta} (\la^{k+1} - \la^k)
%\end{array}
%\right) \\
%& \geq \theta(u^{k+1}) - \theta(u^*) + (w^{k+1} - w)^{\Tsf}F(w^{k+1}) \geq 0.
%\end{align}


\begin{lemma} \label{lemma:vkvk1}
Let the sequence $\{w^k\}$ be generated by sP-PRSM \eqref{equ:sP-PRSM}. If $\alpha \geq 0$ and $\gamma >0$, then for any $w \in \D$, there holds that
\begin{align}\label{equ:lemma:vkvk1}
\|w^k - w\|_G^2 - \|w^{k+1} - w\|_G^2 \geq{} &\| u^k - u^{k+1} \|_P^2 +  (1- \alpha) \beta  \|B(y^k - y^{k+1})\|^2+ (2 - \alpha - \gamma) \beta\left\|r^{k+1} \right\|^2   \nn \\
&  + 2(1 - \alpha)\beta \llang r^{k+1}, B(y^k - y^{k+1}) \rrang  + 2\llang w^{k+1} - w, F(w^{k+1}) \rrang.
%\geq{} & \frac{(1- \alpha)^2}{1 + \alpha} \beta  \|B(y^k - y^{k+1})\|^2+ (2 - \alpha - \gamma) \beta\left\|r^{k+1} \right\|^2  \nn \\
%& + 2\frac{1 - \alpha}{1 + \alpha} (1 - \gamma)\beta \llang Ax^{k} + By^{k} -b, B(y^k - y^{k+1}) \rrang \label{equ:lemma:vkvk1:a0}
\end{align}
\end{lemma}
\begin{proof}
Using the identity
\[
\|a\|_G^2 - \|b\|_G^2 = \|a - b \|_G^2 + 2b^{\Tsf}G (a - b), \nn
\]
with $a = w - w^k$ and $b = w - w^{k+1}$, we have
\begin{align}
\|w^k - w\|_G^2 - \|w^{k+1} - w\|_G^2 = \|w^k - w^{k+1} \|_G^2 + 2 (w - w^{k+1})^{\Tsf} G (w^{k+1} - w^k).  \label{equ:lemma:vkvk1:a2}
\end{align}
Plugging \eqref{equ:lemma:opt:0} into \eqref{equ:lemma:vkvk1:a2} with \eqref{equ:GH_norm} and \eqref{equ:lemma:vkvk1:b1}, we have \eqref{equ:lemma:vkvk1}. The proof is completed.
\end{proof}

Denote
$K =\left( \begin{array}{cc} (1- \alpha) \beta B^{\Tsf} B & (1-\alpha) \beta B^{\Tsf} \\ (1-\alpha)\beta B &  (2 - \alpha - \gamma )\beta I_m\end{array}\right)$.  With  \eqref{equ:M} and \eqref{equ:vk:vk1}
, the inequality \eqref{equ:lemma:vkvk1} can be expressed as
 \begin{align}\label{equ:lemma:vkvk1K}
&\|w^k - w\|_G^2 - \|w^{k+1} - w\|_G^2  \nn \\
\geq{} & \| u^k - u^{k+1} \|_P^2  + (v^k - v^{k+1})^{\Tsf} M^{-\Tsf} K M^{-1} (v^k - v^{k+1}) + 2\llang w^{k+1} - w, F(w^{k+1}) \rrang.
\end{align}
Our task now is to estimate the first term in the righthand side of the inequality in \eqref{equ:lemma:vkvk1K}.

Consider the case where $0 \leq \alpha < 1$ and $0 < \gamma < 1$. It is easy to see that $K \succ 0$. Moreover, by some tedious calculations, we conclude that
 $$K \geq c \cdot M^{\Tsf}H M,$$ where
  $c = \frac{1 - \sqrt{1 - (\alpha + \gamma) ( 1 - \gamma)}}{\alpha + \gamma} \in (0,1).$ Thus we obtain from  \eqref{equ:lemma:vkvk1K} and \eqref{equ:GH_norm} that
\begin{align}
\|w^k - w\|_G^2 - \|w^{k+1} - w\|_G^2 \geq
c \|w^k - w^{k+1}\|_G^2 +2 \llang w^{k+1} - w, F(w^{k+1})\rrang,\label{equ:vkvk1:gamma_less_1}
\end{align}
which is sufficient to establish that the sequence $\{\|w^k - w^*\|_G^2\}$ is strictly contractive by choosing $w = w^*$ in \eqref{equ:vkvk1:gamma_less_1} and further to show the convergence of sP-PRSM \eqref{equ:sP-PRSM}.

Consider the case when $0< \alpha < 1$ but $\gamma \geq 1$. The matrix $K$ is not positive semidefinite any more. Instead of establishing the strictly contractive property of the sequence  $\{\|w^k - w^*\|_G^2\}$, we consider to prove the property of  another sequence related to the sequence  $\{\|w^k - w^*\|_G^2\}$.  To do so, we need to give a more careful estimation of the intersection term $\left\langle r^{k+1}, B(y^k - y^{k+1})\right\rangle$.

\begin{lemma} \label{lemma:AxBy-b}
Let the sequence $\{w^k\}$ be generated by sP-PRSM \eqref{equ:sP-PRSM}. If $\alpha \geq 0$ and $\gamma >0$, then for any $w \in \Omega$, there holds that
\begin{align}\label{equ:lemma:AxBy-b}
  \left\langle r^{k+1}, B(y^k - y^{k+1})\right\rangle  \geq{} & \frac{1 - \gamma}{1 + \alpha}  \left\langle r^{k}, B(y^k - y^{k+1}) \right\rangle  - \frac{\alpha}{1 + \alpha} \|B(y^k - y^{k+1})\|^2 \nn\\
& + \frac{1}{2(1 + \alpha)}\cdot \frac{1}{\beta}   \left(\|y^{k+1} - y^k\|_{\Ptwo}^2 - \|y^k - y^{k-1}\|_{\Ptwo}^2\right).
\end{align}
\end{lemma}
\begin{proof}
 Note that the optimality  condition \eqref{equ:opt:y} also holds with $k \coloneqq k-1$, then we have
\[\label{equ:lemma:AxBy-b:a0}
 \llang y - y^{k}, {\Ptwo}(y^k - y^{k-1}) + \xi_y^k -B^{\Tsf}\la^{k} + (1 - \gamma) \beta  B^{\Tsf} r^k \rrang \geq 0, \quad \forall y \in \Y,
\]
%  Consider the optimality condition \eqref{equ:opt:y}, letting $y$ in  be $y^k$, we obtain that
where $\xi_y^{k} \in \partial \theta_2(y^{k})$. Choosing $y$ to be $y^k$ and $y^{k+1}$ in \eqref{equ:opt:y}  and \eqref{equ:lemma:AxBy-b:a0} and then rearranging the obtained inequalities, respectively, we have that
\[\label{equ:lemma:AxBy-b:a1}
\llang B(y^k - y^{k+1}),  - \la^{k + 1} + (1 - \gamma) \beta   r^{k+1} \rrang \geq \| y^{k+1} - y^k \|_{\Ptwo}^2 + \llang y^{k+1} - y^k, \xi_y^{k+1} \rrang
\]
and
\[\label{equ:lemma:AxBy-b:a2}
\llang B( y^{k} - y^{k+1}),  \la^{k} - (1 - \gamma) \beta r^k \rrang \geq - \llang y^{k+1} - y^k, {\Ptwo}(y^k - y^{k-1})\rrang - \llang y^{k+1} - y^k, \xi_y^k \rrang.
\]
Summing \eqref{equ:lemma:AxBy-b:a1} and \eqref{equ:lemma:AxBy-b:a2} over the both sides and then using the relation \eqref{equ:lak:lak1}, we obtain that
\begin{align}\label{equ:lemma:AxBy-b:d1}
& (1 + \alpha) \beta \llang B(y^{k} - y^{k+1}), r^{k+1} \rrang  + \alpha \beta \|B(y^k - y^{k+1}) \|^2 - (1 - \gamma)\beta \llang B(y^k - y^{k+1}), r^k \rrang  \nn \\
\geq{} &  \| y^{k+1} - y^k \|_{\Ptwo}^2  + \llang y^{k+1} - y^k, \xi_y^{k+1} - \xi_y^k \rrang - \llang y^{k+1} - y^k, {\Ptwo}(y^k - y^{k-1}) \rrang  \nn \\
\geq{} & \frac12 ( \|y^{k+1} - y^k \|_{\Ptwo}^2 - \|y^k - y^{k-1} \|_{\Ptwo}^2 ),
\end{align}
where the second inequality is due to  \eqref{equ:y:convex} and the Cauchy-Schwarz inequality.
It is trivial to obtain  \eqref{equ:lemma:AxBy-b} from  \eqref{equ:lemma:AxBy-b:d1}. The proof is completed.
\end{proof}

With the help of Lemma \ref{lemma:AxBy-b}, we are ready to have the following theorem.
\begin{theorem}
Let the sequence $\{w^k\}$ be generated by sP-PRSM \eqref{equ:sP-PRSM}. If we choose $\alpha$ and $\gamma$ according to \eqref{equ:alpha:gamma} , then for any $w \in \D$, there holds that
\begin{align}\label{equ:wk:wk1}
  &\left( \|w^k - w\|_{G}^2 + \rho_{\alpha,\gamma} \|r^k\|^2 + \eta_{\alpha,\gamma} \|y^{k} - y^{k-1}\|_{\Ptwo}^2\right) - \left( \|w^{k+1} - w\|_{G}^2 + \rho_{\alpha,\gamma} \|r^{k+1}\|^2 + \eta_{\alpha,\gamma} \|y^{k+1} - y^{k}\|_{\Ptwo}^2\right)  \nn \\
  \geq{} & \tau_{\alpha, \gamma}\|w^{k} - w^{k+1}\|_G^2  + 2(w^{k+1} - w)^{\Tsf}F(w^{k+1}),
\end{align}
where $\rho_{\alpha,\gamma}, \eta_{\alpha,\gamma} \geq 0$ and $0< \tau_{\alpha, \gamma} <1$ are some constants which only depend on $\alpha$ and $\gamma$.
\end{theorem}
\begin{proof}
We consider three cases.

I). $0 < \gamma < 1$.  With  $\rho_{\alpha, \gamma} = \eta_{\alpha, \gamma} =  0$ and $\tau_{\alpha, \gamma} = \frac{1 - \sqrt{1 - (\alpha + \gamma) ( 1 - \gamma)}}{\alpha + \gamma}$, the inequality \eqref{equ:wk:wk1} coincides with \eqref{equ:vkvk1:gamma_less_1} which has been proved. 

Before we proceed, by combining \eqref{equ:lemma:AxBy-b} and \eqref{equ:lemma:vkvk1}, we derive the useful inequality as
\begin{align}
&\left(\|w^k - w\|_G^2  + \frac{1 - \alpha}{1 + \alpha} \|y^{k+1} - y^k \|_{\Ptwo}^2\right)- \left(\|w^{k+1} - w\|_G^2 + \frac{1 - \alpha}{1 + \alpha} \|y^k - y^{k-1}\|_{\Ptwo}^2 \right)  \nn\\
\geq{} & \|u^k - u^{k+1} \|_P^2 + \frac{(1- \alpha)^2}{1 + \alpha} \beta  \|B(y^k - y^{k+1})\|^2+ (2 - \alpha - \gamma) \beta\left\|r^{k+1} \right\|^2  \nn \\
&   + 2(1 - \gamma)\frac{1 - \alpha}{1 + \alpha} \beta \llang r^{k}, B(y^k - y^{k+1}) \rrang +  2\llang w^{k+1} - w, F(w^{k+1}) \rrang.
\label{equ:vkvk1:0}
\end{align}

II). $\gamma = 1$. It is not difficult to see from \eqref{equ:vkvk1:0}, \eqref{equ:GH_norm} and\eqref{equ:lemma:vkvk1:b1} that \eqref{equ:wk:wk1} holds with $\rho_{\alpha, \gamma} =0$ and $ \eta_{\alpha, \gamma} =  \tau_{\alpha, \gamma} = \frac{1 - \alpha}{1 + \alpha}$.

III). $1< \gamma < \frac{1 - \alpha + \sqrt{(1+ \alpha)^2 + 4(1 - \alpha^2)}}{2}$.
Let us choose any $\delta \in \left(\frac{\gamma - 1}{1 - \alpha}, \frac{1+ \alpha}{\gamma - 1} - \frac{1 + \alpha}{1-\alpha}\right)$. Note that this interval is well-defined due to the range of $\gamma$.  With the Cauchy-Schwarz inequality, we know that
\[
-2 \llang r^{k}, B(y^k - y^{k+1}) \rrang \geq - \delta \cdot\|r^{k}\|^2 - \frac{1}{\delta}\cdot\|B(y^k - y^{k+1})\|^2. \nn
\]
Plugging the above inequality into \eqref{equ:vkvk1:0}, we obtain that
\begin{align}
  &\left(\|w^k - w\|_G^2  + \frac{1 - \alpha}{1 + \alpha} \|y^{k+1} - y^k \|_{\Ptwo}^2\right)- \left(\|w^{k+1} - w\|_G^2 + \frac{1 - \alpha}{1 + \alpha} \|y^k - y^{k-1}\|_{\Ptwo}^2 \right)  \nn\\
  \geq{} & \|u^k - u^{k+1} \|_P^2 +   \frac{(1- \alpha)^2}{1 + \alpha} \left(1 - \frac{\gamma - 1}{1 - \alpha}\cdot \frac{1}{\delta} \right) \beta  \|B(y^k - y^{k+1})\|2
  + (\gamma - 1) \frac{1 - \alpha}{1 + \alpha}\left(\frac{1+ \alpha}{\gamma - 1} - \frac{1 + \alpha}{1-\alpha} -\delta \right) \left\|r^{k+1} \right\|^2  \nn \\
& - \delta(\gamma - 1) \frac{1 - \alpha}{1 + \alpha}\left( \|r^{k}\|^2 -  \|r^{k+1}\|^2\right) +  2\llang w^{k+1} - w, F(w^{k+1}) \rrang,\nn
%\label{equ:lemma:vkvk1:a1}
\end{align}
which with \eqref{equ:GH_norm} and \eqref{equ:lemma:vkvk1:b1} implies that \eqref{equ:wk:wk1} holds with $\rho_{\alpha, \gamma} = \delta(\gamma - 1) \frac{1 - \alpha}{1 + \alpha}$, $\eta_{\alpha, \gamma} = \frac{1 - \alpha}{1 + \alpha}$ and
\[
  \tau_{\alpha, \gamma} =  \frac{1- \alpha}{1 + \alpha} \cdot \min\left\{\left(1 - \frac{\gamma - 1}{1 - \alpha}\cdot \frac{1}{\delta} \right),\frac{\gamma - 1}{\alpha + \gamma}\left(\frac{1+ \alpha}{\gamma - 1} - \frac{1 + \alpha}{1-\alpha} -\delta \right) \right\}. \label{equ:tau:alpha:gamma}
\]
The proof is completed.
\end{proof}

Finally, we end this subsection by finishing the proof of Theorem \ref{thm:wk:wk1:contraction}.
%\begin{proof}

\noindent{\bf Proof of Theorem \ref{thm:wk:wk1:contraction}}.
Consider any $w^* \in \Omega^*$, recalling the monotonicity \eqref{equ:F:mono}
of  $F(w)$,  we thus obtain from the optimality \eqref{equ:opt:cond}  of $w^*$ that
$$
 \llang w^{k+1} - w^*, F(w^{k+1}) \rrang \geq \llang w^{k+1} - w^*, F(w^*) \rrang + \|u^{k+1} - u^* \|_{\SigT}^2\geq \|u^{k+1} - u^* \|_{\Sigma}^2.
$$
With the above statement, choosing $w = w^*$ in \eqref{equ:wk:wk1}, we  can easily  obtain \eqref{equ:wk:wk1:w*}.
\begin{align}\label{equ:wk:wk1:*}
  &\left( \|w^k - w^*\|_{G}^2 + \rho_{\alpha,\gamma} \|r^k\|^2 + \eta_{\alpha,\gamma} \|y^{k} - y^{k-1}\|_{\Ptwo}^2\right) - \left( \|w^{k+1} - w^*\|_{G}^2 + \rho_{\alpha,\gamma} \|r^{k+1}\|^2 + \eta_{\alpha,\gamma} \|y^{k+1} - y^{k}\|_{\Ptwo}^2\right)  \nn \\
  \geq{} & \tau_{\alpha, \gamma}\|w^{k} - w^{k+1}\|_G^2+  2\|u^{k+1} - u^* \|_{\Sigma}^2,
\end{align}
Adding the term $\|u^k - u^*\|_{\Sigma}^2 - \|u^{k+1} - u^*\|_{\Sigma}^2$ on both sides of \eqref{equ:wk:wk1:*}, we have
\begin{align}\label{equ:wk:wk1:a:*}
  &\left(\|u^k - u^*\|_{\Sigma}^2 +  \|w^k - w^*\|_{G}^2 + \rho_{\alpha,\gamma} \|r^k\|^2 + \eta_{\alpha,\gamma} \|y^{k} - y^{k-1}\|_{\Ptwo}^2\right) \nn \\
  -{}&  \left( \|u^{k+1} - u^*\|_{\Sigma}^2 + \|w^{k+1} - w^*\|_{G}^2 + \rho_{\alpha,\gamma} \|r^{k+1}\|^2 + \eta_{\alpha,\gamma} \|y^{k+1} - y^{k}\|_{\Ptwo}^2\right)  \nn \\
  \geq{} & \tau_{\alpha, \gamma}\|w^{k} - w^{k+1}\|_G^2 +  \|u^{k+1} - u^* \|_{\Sigma}^2 + \|u^k - u^*\|_{\Sigma}^2 \nn \\
  \geq{} & \tau_{\alpha, \gamma}\|w^{k} - w^{k+1}\|_G^2 + \frac12 \|u^k - u^{k+1}\|_{\Sigma}^2,
\end{align}
which with \eqref{equ:Ghat} and \eqref{equ:GhatH_norm} leads to \eqref{equ:wk:wk1:w*}.
The proof is completed.

  \subsection{Global convergence}

 \begin{theorem}
    Let the sequence $\{w^k\}$ be generated by sP-PRSM \eqref{equ:sP-PRSM}.  If we choose $\alpha$ and $\gamma$ according to \eqref{equ:alpha:gamma} , then $\{w^{k}\}$ converges
 to an optimal solution of \eqref{cp}.
  \end{theorem}
\begin{proof}
We divide the proof into three steps.

(i) To show that the sequences $\{w^{k}\}$ is  bounded.  It is trivial to know from \eqref{equ:wk:wk1:w*} that $\lim_{k \rightarrow \infty} \|w^k - w^{k+1}\|_{\Ghat}^2 = 0$, which with \eqref{equ:GH_norm}, \eqref{equ:GhatH_norm} and \eqref{equ:lemma:vkvk1:b1} indicates that
\[ \label{equ:r:limit}
\lim_{k \rightarrow \infty} \|x^k - x^{k+1}\|_{\Pone + \Sigt} = \lim_{k \rightarrow \infty} \|y^k - y^{k+1}\|_{\Ptwo + \Sigtt} = \lim_{k \rightarrow \infty} \|r^k\| = \lim_{k \rightarrow \infty} \|B(y^k - y^{k+1})\| = 0.
\]
 By some simple calculations, there holds that there exists a constant $c_0= \frac{2(\alpha + \gamma) - \alpha \gamma - \alpha \sqrt{\gamma^2 + 4(\alpha + \gamma)}}{2(\alpha + \gamma)} \in (0,1]$ such that
\[  \label{equ:Ghat:a}
\Ghat \succeq c_0\left(\begin{array}{ccc} S + \Sigt  & 0 & 0 \\ 0 & T + \Sigtt + \beta B^{\Tsf}B & 0 \\ 0 & 0 & \frac{1}{\beta} I_m  \end{array}\right).
\]
Besides, it is straightforward to see from \eqref{equ:wk:wk1:w*} that $\|w^k - w^*\|_{\Ghat}^2$ is bounded. This together with \eqref{equ:Ghat:a} shows that
\[
  \|x^k - x^*\|_{\Pone + \Sigt},  \quad \|y^k - y^*\|_{T + \Sigtt + \beta B^{\Tsf}B}, \quad \|\lambda^k - \lambda^* \| \nn
\]
are all bounded.  Obviously, we claim that $\{\lambda^k\}$ is bounded. Moreover,   we see that $\{y^k\}$ is bounded as $\Ptwo + \Sigtt + \beta B^{\Tsf}B \succ 0$. Notice that $\|Ax^k - Ax^*\| = \| r^k  + B(y^k - y^*) \| \leq \|r^k \| + \|B(y^k - y^*)\|$, with \eqref{equ:r:limit}, this implies  that $\|x^k - x^*\|_{\beta A^{\Tsf}A}$ is bounded.  Recalling that $ \Pone + \Sigt + \beta A^{\Tsf}A \succ 0$ and $\|x^k - x^*\|_{\Pone + \Sigt}$ is bounded, it is safe to say that $\{x^k\}$ is also bounded.

(ii) To show that any cluster point  of the sequence $\{w^k\}$ is an optimal solution of \eqref{cp}. Let $\{w^{k_i}\}$ be a subsequence of the sequence $\{w^k\}$ and
$\lim\limits_{k_i \rightarrow \infty} w^{k_i}  = w^{\infty}$. Since the graphs of $\partial \theta_1(\cdot)$ and $\partial \theta_2(\cdot)$ are both closed, taking the limit with respect $k_i \rightarrow \infty$ on both sides of  \eqref{equ:lemma:opt:0} and using  \eqref{equ:r:limit}, we have that
\[
  (w - w^{\infty})^{\Tsf} F(w^{\infty}) \geq 0, \quad \forall w \in \D, \nn
\]
which means that $w^{\infty}$ is an optimal solution of \eqref{cp}.

(iii) To show that the sequence $\{w^k\}$ has only one cluster point.  We first replace $x^*$ with $x^{\infty}$ in the analysis of Steps (i) and (ii).  It follows from $ w^{k_i}  \rightarrow  w^{\infty}$ that $
\|w^{k_i} - w^{\infty}\|_{\Ghat}^2 + \rho_{\alpha,\gamma} \|r^{k_i}\|^2 + \eta_{\alpha,\gamma} \|y^{k_i} - y^{k_i-1}\|_{\Ptwo}^2 \rightarrow 0$.  Owing to the monotonicity of the sequence $\|w^k - w^{\infty}\|_{\Ghat}^2 + \rho_{\alpha,\gamma} \|r^k\|^2 + \eta_{\alpha,\gamma} \|y^{k} - y^{k-1}\|_{\Ptwo}^2$, we can see that
\[
\lim_{k \rightarrow \infty} \|w^k - w^{\infty}\|_{\Ghat}^2 + \rho_{\alpha,\gamma} \|r^k\|^2 + \eta_{\alpha,\gamma} \|y^{k} - y^{k-1}\|_{\Ptwo}^2 = 0. \nn
\]
This together with \eqref{equ:Ghat:a} shows that
\[ \label{equ:lim:proof:yishengtanxi}
  \lim_{k \rightarrow \infty} \|x^k - x^{\infty}\|_{\Pone + \Sigt}  = \lim_{k \rightarrow \infty}  \|y^k - y^{\infty}\|_{T + \Sigtt + \beta B^{\Tsf}B} = \lim_{k \rightarrow \infty}  \quad \|\lambda^k - \lambda^{\infty} \|  = 0.
\]
Using again the inequality $\|Ax^k - Ax^{\infty}\| = \| r^k  + B(y^k - y^{\infty}) \| \leq \|r^k \| + \|B(y^k - y^{\infty})\|$ and $\lim_{k \rightarrow \infty} \|r_k\| = 0$,  we have
\[\label{equ:lim:proof:yishengtanxi:2}
\lim_{k \rightarrow \infty} \|A(x^k - x^{\infty})\| = 0.
\]
Combing \eqref{equ:lim:proof:yishengtanxi} and \eqref{equ:lim:proof:yishengtanxi:2}, and using that
$\Pone + \Sigt + \beta A^{\Tsf}A \succ 0$ and $\Ptwo + \Sigtt + \beta B^{\Tsf}B \succ 0$, we immediately have
\[
 \lim_{k \rightarrow \infty} w^k = w^{\infty}. \nn
\]
The proof is completed.
\end{proof}
\section{Sublinear convergence of sP-PRSM}\label{section:sublinear}
The rate of convergence of an algorithm can help us have a deeper understanding of the algorithm. Thus, in this section, we establish the sublinear rate of convergence of sP-PRSM, in ergodic sense and nonergodic sense, respectively.
\subsection{Convergence rate in the ergodic sense}
We now give the sublinear rate of convergence of sP-PRSM in the ergodic sense, which is very easy due to the key inequality \eqref{equ:wk:wk1}.
\medskip

\begin{theorem}\label{thm:ergodic}
Let the sequence $\{w^k\}$ be generated by sP-PRSM \eqref{equ:sP-PRSM}. For any integer $t>0$, define
\[
  \bar{w}_t = \frac{1}{t} \sum_{k = 1}^{t} w^{k+1}. \nn
  \]
  Then for any $w \in \D$,  we have that
\[\label{thm:ergodic:00}
\sup_{w \in \Omega} \, \llang \bar{w}_t - w, F(w) \rrang \le \frac{1}{2t} \left( \|w^1 - w^0\|_{G}^2 + \rho_{\alpha,\gamma} \|r^1\|^2 + \eta_{\alpha,\gamma} \|y^{1} - y^{0}\|_{\Ptwo}^2\right).
\]
\end{theorem}
\begin{proof}
It follows from \eqref{equ:wk:wk1} and $(w^{k+1} - w)^{\Tsf}F(w^{k+1}) \geq (w^{k+1} - w)^{\Tsf}F(w)$ following from \eqref{equ:F:mono} that
\begin{align}
  \llang w^{k+1} - w, F(w) \rrang
\leq {}& \frac12 \left( \|w^k - w\|_{G}^2 + \rho_{\alpha,\gamma} \|r^k\|^2 + \eta_{\alpha,\gamma} \|y^{k} - y^{k-1}\|_{\Ptwo}^2\right)  \nn \\
  {}& - \frac12 \left( \|w^{k+1} - w\|_{G}^2 + \rho_{\alpha,\gamma} \|r^{k+1}\|^2 + \eta_{\alpha,\gamma} \|y^{k+1} - y^{k}\|_{\Ptwo}^2\right) \nn
\end{align}
For all $k = 1, \ldots, t$, summing the above inequality from both sides and noting the notation of $\bar{w}_t$, we derive that
\[\label{thm:ergodic:a0}
  t \llang \bar{w}_t - w, F(w) \rrang \leq \frac12 \left( \|w^1 - w^0\|_{G}^2 + \rho_{\alpha,\gamma} \|r^1\|^2 + \eta_{\alpha,\gamma} \|y^{1} - y^{0}\|_{\Ptwo}^2\right) .
  \]
  Since \eqref{thm:ergodic:a0} holds for any $w \in \D$, we can easily have  \eqref{thm:ergodic:00}. The proof is completed.
\end{proof}
\medskip

\noindent{\bf Remark}: Consider the case when the  set $\D$ is compact, \eqref{thm:ergodic:00} implies the $O(1/t)$ ergodic convergence of sP-PRSM \eqref{equ:sP-PRSM}.

\subsection{Convergence rate in the nonergodic sense}
\rev{
  \begin{lemma}
    Let the sequence $\{w^k\}$ be generated by sP-PRSM \eqref{equ:sP-PRSM}.  If $0 \leq \alpha \leq 1$ and $\gamma >0$, then $w^{k+1} \in \Omega^*$, namely, $w^{k+1}$ is one optimal solution of \eqref{cp},  if $$\|w^k - w^{k+1}\|_{G} = 0.$$
  \end{lemma}
  \begin{proof}
  Following from \eqref{equ:GH_norm}  and \eqref{equ:lemma:vkvk1:b1},  the statement $\|w^k - w^{k+1}\|_G = 0$ indicates that $r^{k+1} = 0$. Noting that $G$ is always positive semidefinite, it is trivial to  have $G(w^k - w^{k+1}) = 0.$  Plugging the above two assertions into \eqref{equ:lemma:opt:0}, we see that  $ \llang w - w^{k+1}, F(w^{k+1}) \rrang \geq 0$  holds for any $w \in \D$. Following the optimality condition \eqref{equ:opt:cond:D},  we immediately see that  $w^{k+1} \in \Omega^*$. The proof is completed.
  \end{proof}
  Based on this lemma,  we can use  $\|w^k - w^{k+1}\|_G$ to measure the accuracy of $w^{k+1}$. Similarly to the definition \eqref{equ:epsilon_sol} of $\epsilon$-solution of \eqref{cp}, we also call  $w^{k+1}$ as an $\epsilon$-solution of \eqref{cp} when $\|w^k - w^{k+1}\|_G \leq \epsilon.$  Together with \eqref{equ:GH_norm} and \eqref{equ:lemma:vkvk1:b1}, this also provide a practical stopping condition for sP-PRSM \eqref{equ:sP-PRSM} as
  \[\label{equ:stop0}
   \max\{ \|x^{k+1} - x^k\|_{\Pone}, \|y^{k+1} - y^k \|_{\Ptwo},  \|B(y^{k+1} - y^k)\|, \|r^{k+1}\|\} \leq \mathrm{tol},
  \]
  where $\mathrm{tol}$ is some tolerance. Note that in \eqref{equ:stop0},  we may pay some more price to compute  $\|x^{k+1} - x^k\|_{\Pone}, \|y^{k+1} - y^k \|_{\Ptwo}$. Thus, when the spectral norm of $S$ or $T$ is bounded, we can also simply stop the method when
  \[\label{equ:stop}
   \max\{\|B(y^{k+1} - y^k)\|, \|r^{k+1}\|\} \leq \mathrm{tol}.
  \]
}

\begin{lemma}\label{lemma:decrease}
Let the sequence $\{w^k\}$ be generated by sP-PRSM \eqref{equ:sP-PRSM}.  If we choose $\alpha$ and $\gamma$ according to \eqref{equ:alpha:gamma}, then there holds that
\[ \label{equ:lemma:decrease:00}
\|w^{k+1} - w^k \|_G^2 - \|w^{k+2} - w^{k+1} \|_G^2 \geq \frac{1 - \gamma}{\alpha + \gamma} \|(w^{k+1}-w^{k}) - (w^{k+2} - w^{k+1} )\|_G^2.
    \]
\end{lemma}
\begin{proof}
%It follows from \eqref{equ:lemma:opt:0} that
%\begin{align}
 % &(v - v^{k+1})^{\Tsf}  H (v^{k+1}-v^{k}) - [\theta(u^{k+1}) - \theta(u) + (w^{k+1} - w)^{\Tsf}F(w^{k+1})] \nn \\
%\ge{}& (1 - \alpha - \gamma) \beta \langle r^{k+1}, r^{k+1} - r\rangle   + (1 - \alpha)\beta \langle B(y^k - y^{k+1}), r^{k+1} - r\rangle . \label{equ:lemma:opt:nonergodic}
%\end{align}
Note that \eqref{equ:lemma:opt:01} also holds with $k \coloneqq k +1$, then we have
\begin{align}
(w - w^{k+2})^{\Tsf}  G (w^{k+2}-w^{k+1})
\ge{} & \llang r^{k+2} -r(w), (1 - \alpha - \gamma) \beta   r^{k+2} + (1 - \alpha) \beta   B (y^{k+1} -y^{k+2}) \rrang  \nn \\ &+ \llang w^{k+2} - w, F(w^{k+2}) \rrang.  \label{equ:lemma:opt:nonergodic:2}
\end{align}
%\begin{align}
 % &(v - v^{k+2})^{\Tsf}  H (v^{k+2}-v^{k}) - [\theta(u^{k+2}) - \theta(u) + (w^{k+2} - w)^{\Tsf}F(w^{k+2})] \nn \\
%\ge{}& (1 - \alpha - \gamma) \beta \langle r^{k+2}, r^{k+2} - r\rangle   + (1 - \alpha)\beta \langle B(y^k - y^{k+2}), r^{k+2} - r\rangle . \label{equ:lemma:opt:nonergodic:2}
%\end{align}
Choosing $w$ to be $w^{k+2}$ and $w^{k+1}$, respectively, in \eqref{equ:lemma:opt:01} and \eqref{equ:lemma:opt:nonergodic:2} leads to
\begin{align}
 (w^{k+2} - w^{k+1})^{\Tsf}  G (w^{k+1}-w^{k})
\ge{}& \llang r^{k+1} -r^{k+2}, (1 - \alpha - \gamma) \beta   r^{k+1} + (1 - \alpha) \beta   B (y^{k} -y^{k+1}) \rrang \nn \\
 &+ \llang w^{k+1} - w^{k+2}, F(w^{k+1}) \rrang. \label{equ:lemma:opt:nonergodic:b0}
\end{align}
and
\begin{align}
 (w^{k+1} - w^{k+2})^{\Tsf}  G (w^{k+2}-w^{k+1})
\ge{} &\llang r^{k+2} -r^{k+1}, (1 - \alpha - \gamma) \beta   r^{k+2} + (1 - \alpha) \beta   B (y^{k+1} -y^{k+2}) \rrang\nn \\
&+ \llang w^{k+2} - w^{k+1}, F(w^{k+2}) \rrang. \label{equ:lemma:opt:nonergodic:b1}
\end{align}
Adding \eqref{equ:lemma:opt:nonergodic:b0} and \eqref{equ:lemma:opt:nonergodic:b1} and noting
$\llang w^{k+2} - w^{k+1}, F(w^{k+2}) - F(w^{k+1}) \rrang \geq 0$, following from \eqref{equ:F:mono},   we obtain that
\begin{align}
  &(w^{k+2} - w^{k+1})^{\Tsf}  G \left[(w^{k+1}-w^{k}) - (w^{k+2} - w^{k+1}) \right] \nn \\
  \ge{}& (1 - \alpha - \gamma) \beta \| r^{k+1} - r^{k+2}\|^2  + (1 - \alpha)\beta \left \langle B\left[(y^k - y^{k+1}) - (y^{k+2} - y^{k+1})\right], r^{k+1} - r^{k+2}\right\rangle . \label{equ:lemma:opt:nonergodic:b3}
\end{align}
Following the deriving process of \eqref{equ:lemma:vkvk1:b1}, we have that
\begin{align}
&  \|(w^{k+1}-w^{k}) - (w^{k+2} - w^{k+1} )\|_G^2 \nn \\
={} & (1-\alpha)\beta \left\|B[ (y^{k+1} - y^k) - (y^{k+2} - y^{k+1})]\right\|^2 + (\alpha + \gamma) \beta \|r^{k+1} - r^{k+2}\|^2.\label{equ:lemma:opt:nonergodic:b4}
\end{align}
Thus we conclude that
\begin{align}
  & \|w^{k+1} - w^k \|_H^2 - \|w^{k+2} - w^{k+1} \|_H^2 \nn \\
  ={}& 2(w^{k+2} - w^{k+1})^{\Tsf}  G \left[(w^{k+1}-w^{k}) - (w^{k+2} - w^{k+1}) \right] + \|(w^{k+1}-w^{k}) - (w^{k+2} - w^{k+1} )\|_H^2 \nn \\
  \ge{}& (2 - \alpha - \gamma) \beta \| r^{k+1} - r^{k+2}\|^2  + 2(1 - \alpha)\beta \left \langle B\left[(y^k - y^{k+1}) - (y^{k+2} - y^{k+1})\right], r^{k+1} - r^{k+2}\right\rangle \nn \\
  & + (1-\alpha)\beta \left\|B[ (y^{k+1} - y^k) - (y^{k+2} - y^{k+1})]\right\|^2 \nn \\
  \geq{} &(1 -\gamma) \beta \| r^{k+1} - r^{k+2}\|^2 \geq \frac{1 - \gamma}{\alpha + \gamma} \|(w^{k+1}-w^{k}) - (w^{k+2} - w^{k+1} )\|_G^2,\nn
 % \label{equ:lemma:opt:nonergodic:b5}
\end{align}
where the first inequality is due to \eqref{equ:lemma:opt:nonergodic:b3} and \eqref{equ:lemma:opt:nonergodic:b4}, the second inequality is trivial and the last inequality owes to \eqref{equ:lemma:opt:nonergodic:b4}.
The proof is completed.
\end{proof}
\begin{theorem}\label{theorem:nonergodic}
Let the sequence $\{v^k\}$ be generated by sP-PRSM \eqref{equ:sP-PRSM}.  If $\gamma \in (0,1]$, then there holds that
  \[
    \|w^{t+1} - w^{t} \|_G^2 \leq \frac{1}{\tau_{\alpha,\gamma}}\cdot\frac{1}{t} \cdot \left(\|w^1 - w^0\|_G^2 + \rho_{\alpha,\gamma} \|r^1\|^2 + \eta_{\alpha,\gamma} \|y^1 - y^0\|_{\Ptwo}^2 \right). \label{theorem:nonergodic:00}
    \]
\end{theorem}
\begin{proof}
  Summing \eqref{equ:wk:wk1:*}over $k = 1, \ldots, t$ leads to
  \[\label{theorem:nonergodic:b0}
    \tau_{\alpha,\gamma} \cdot   \sum_{k=1}^t \|w^{k+1} - w^{k} \|_G^2 \leq \|w^1 - w^*\|_G^2 + \rho_{\alpha,\gamma} \|r^1\|^2 + \eta_{\alpha, \gamma} \|y^1 - y^0\|_{\Ptwo}^2.
    \]
Observe that $\gamma \in (0, 1]$, it follows from \eqref{equ:lemma:decrease:00} that $\|w^{k+1} - w^k\|_G^2 \geq \|w^{k+2} - w^{k+1}\|_G^2$, which together with \eqref{theorem:nonergodic:b0} can easily suggest \eqref{theorem:nonergodic:00}. The proof is completed.
\end{proof}


\section{Linear convergence of sP-PRSM \eqref{equ:sP-PRSM}}\label{section:linear}
For simplicity, we only consider the case when $T = 0$ and $S = 0$.  In this case,  we always have that
\[ \label{equ:reduce:G}
G = \left( \begin{array}{cc} 0 & 0 \\ 0 & H\end{array} \right) \quad \mbox{and} \quad \|w - w'\|_G^2 = \|v - v'\|_H^2.
\]
\begin{ass} \label{ass:Y}
The set $\Y = \R^{n_2}$.
\end{ass}


\begin{ass}\label{ass:gradient_lips}
The gradient $\nabla \theta_2(\cdot)$ is Lipschtiz continuous on $\R^{n_2}$  with $L>0$. Thus for any $y_1, y_2 \in \R^{n_2} $ there holds that
\[ \label{equ:gradient_lips}
\|  \nabla \theta_2(y_1) - \nabla \theta_2(y_2) \| \leq L  \| y_1 - y_2\|
\]
and
\[ \label{equ:gradient_lips:2}
\llang y_1 - y_2, \nabla \theta_2(y_1) - \nabla \theta_2(y_2) \rrang \geq
\frac{1}{L} \| \nabla \theta_2(y_1) - \nabla \theta_2(y_2)\|^2.
\]
\end{ass}


\begin{ass}\label{ass:strong_convex}
The function $\theta_2(\cdot)$ is strongly convex on $\R^{n_2}$ with constant $\mu > 0$. Thus for any $y_1, y_2 \in \R^{n_2} $ there holds that
\[\label{equ:strong_convex:2}
\llang y_1 - y_2, \nabla \theta_2(y_1) - \nabla \theta_2(y_2) \rrang \geq \mu \|y_1 - y_2\|^2.
\]
\end{ass}

%Following from Theorem 2.1.11 in \cite{nesterov2004introductory}, we have that

Since $\Y = \R^{n_2}$, it follows from \eqref{equ:opt:y} that
\[\label{equ:theta_yk1}
\nabla \theta_2(y^{k+1}) = B^{\Tsf} \lambda^{k+1} - (1 - \gamma) \beta B^{\Tsf} r^{k+1}.
\]
Similarly, we obtain from \eqref{equ:opt:cond:D} that
\[ \label{equ:theta_y_opt}
\nabla \theta_2(y^*) = B^{\Tsf} \lambda^*.
\]

\begin{theorem}\label{thm:linear}
Let the sequence $\{w^k\}$ be generated by sP-PRSM \eqref{equ:sP-PRSM}.  Under Assumptions  \ref{ass:Y}, \ref{ass:strong_convex}, \ref{ass:gradient_lips}, there must exists some constant $c >0$ such that
\[
\|v^k - v^*\|_{H}^2 + \rho_{\alpha,\gamma} \|r^k\|^2 \geq (1 + c) \cdot \left( \|v^{k+1} - v^*\|_{H}^2 + \rho_{\alpha,\gamma}\|r^{k+1}\|^2 \right). \nn
\]
\end{theorem}
\begin{proof}
Letting $v = v^*$ in \eqref{equ:wk:wk1}, with \eqref{equ:reduce:G},  we can obtain
\begin{align}\label{equ:theorem:linear_convergence:00}
 & \left( \|v^k - v^*\|_{H}^2 + \rho_{\alpha,\gamma} \|r^k\|^2  \right)-\left( \|v^{k+1} - v^*\|_{H}^2 + \rho_{\alpha,\gamma}\|r^{k+1}\|^2 \right)  \nn \\
\geq{}  & \tau_{\alpha, \gamma}\|v^{k} - v^{k+1}\|_H^2 +2 \llang w^{k+1} - w^*, F(w^{k+1})\rrang.
\end{align}
Next we will estimate the two terms of the righthand side of \eqref{equ:theorem:linear_convergence:00}, respectively. Firstly, with \eqref{equ:lemma:vkvk1:b1}, we see that
\[ \label{equ:thm:linear:vkvk1}
\|v^{k} - v^{k+1}\|_H^2 \geq (\alpha + \gamma)\beta \|r^{k+1}\|^2.
\]
It follows from \eqref{equ:F:mono},  \eqref{equ:opt:cond} and the convexity of $\theta_1(\cdot)$ that
\begin{align}\label{equ:thm:linear:a1}
      \llang w^{k+1} - w^*,  F(w^{k+1}) \rrang  \geq{} &
      \llang
      y^{k+1} - y^*,  \nabla \theta_2(y^{k+1}) - \nabla \theta_2(y^*)
      \rrang
   \end{align}
For any $0\leq t \leq 1$, consider the convex combination of  \eqref{equ:strong_convex:2} and \eqref{equ:gradient_lips:2},  we know from \eqref{equ:thm:linear:a1}, \eqref{equ:theta_yk1} and  \eqref{equ:theta_y_opt} that there holds that
 \begin{align}\label{equ:thm:linear:a2}
      \llang w^{k+1} - w^*,  F(w^{k+1}) \rrang  \geq{} &
        \mu (1 - t) \|y^{k+1} - y^* \|^2 +
\frac{\la_{\min} (BB^{\Tsf})}{L} t \cdot  \| (\lambda^{k+1} - \lambda^*) - (1 - \gamma) \beta r^{k+1}\|^2.
      \end{align}
      By applying the inequality
      \[
      \|a + b \|^2 \geq \left(1 - \frac{1}{\rho}\right) \|a\|^2 + ( 1 - \rho) \|b\|^2, \quad \rho>1 \nn
      \]
      with $a = \lambda^{k+1} - \lambda^*$, $b = -(1-\gamma) r^{k+1}$ and $\tau>1$,
      we have
      \[ \label{equ:thm:linear:m1}
      \|  (\lambda^{k+1} - \lambda^*) - (1 - \gamma)  r^{k+1}\|^2 \geq \left(1 - \frac{1}{\rho}\right) \|\lambda^{k+1} - \lambda^*\|^2 + (1 - \rho) |\gamma - 1| \beta \cdot \| r^{k+1}\|^2.
      \]
  Plugging \eqref{equ:thm:linear:m1} into     \eqref{equ:thm:linear:a2} , we further have
  \begin{align}\label{equ:thm:linear:m2}
      \llang w^{k+1} - w^*,  F(w^{k+1}) \rrang  \geq{} &
  c_1 \left(  \beta \|B\|^2 \cdot \| y^{k+1} - y^*\|^2 +  \frac{1}{\beta} \|\lambda^{k+1} - \lambda^* \|^2 \right) + c_2 \beta \|r^{k+1}\|^2,
  \end{align}
        where $c_1 = \min\left(\frac{\mu}{\beta \|B\|^2} (1 - t) , \left(1 - \rho^{-1}\right)\frac{\la_{\min} (BB^{\Tsf})}{L} \beta t\right)$, $c_2 = (1-\rho)|\gamma - 1|\frac{\la_{\min} (BB^{\Tsf})}{L}  t$.
        With the definition \eqref{equ:H} of $H$, we can easily see that there exits a positive constant $c_{\alpha,\gamma} = \frac{2(\alpha + \gamma) - \alpha \gamma - \alpha \sqrt{\gamma^2 + 4(\alpha + \gamma)}}{2(\alpha + \gamma)(1 - \alpha)} $ such that
\[ \label{equ:thm:linear:m3}
 \beta \|B\|^2 \cdot \| y^{k+1} - y^*\|^2 +  \frac{1}{\beta} \|\lambda^{k+1} - \lambda^* \|^2 \geq c_{\alpha,\gamma} \cdot \|v^{k+1} - v^* \|_H^2.
\]
Combing \eqref{equ:thm:linear:vkvk1}, \eqref{equ:thm:linear:m2} and \eqref{equ:thm:linear:m3},  we see that
\begin{align}
\mbox{RHS of \eqref{equ:theorem:linear_convergence:00}} \geq{} &  2(c_1 \cdot c_{\alpha,\gamma})  \|v^{k+1} - v^* \|_H^2 + (2c_2 + \alpha + \gamma) \beta \|r^{k+1}\|^2 \nn \\
\geq{} & c \left( \|v^{k+1} - v^*\|_{H}^2 + \rho_{\alpha,\gamma}\|r^{k+1}\|^2 \right),
\end{align}
where $c = \min\left(2c_1c_{\alpha, \gamma},  \frac{2c_2 + \alpha + \gamma}{\rho_{\alpha,\gamma}}\beta\right).$ The proof is completed.
\end{proof}
\medskip

\noindent{ \textbf{Remark}.} Let us consider as a special case, namely,  $0 \leq \alpha <1, \gamma = 1$. Here, $\rho_{\alpha,\gamma} = 0$ and    \eqref{equ:thm:linear:m1} reduces to
 $\|(\lambda^{k+1} - \lambda^*) - (1 - \gamma)  r^{k+1}\|^2 =\|\lambda^{k+1} - \lambda^*\|^2$.
 Hence  $c = 2c_1c_3$. Letting $t = \frac{\mu L}{\mu L + \beta^2 \|B\|^2 \lambda_{\min}(BB^{\Tsf})}$, we have that $$c_1 =\frac{\mu \beta \lambda_{\min}(BB^{\Tsf}) }{\mu L + \beta^2 \|B\|^2 \lambda_{\min}(BB^{\Tsf})}  \leq \frac{1}{2 \kappa_{B}  \sqrt{\kappa_{\theta_2}}}$$
 where $\kappa_{B} = \sqrt{\lambda(BB^{\Tsf})/\lambda_{\min}(BB^{\Tsf})}$, $\kappa_{\theta_2} = L/{\mu}$ and the equality attains when $\beta = \sqrt{\frac{\mu L}{\|B\|^2 \lambda_{\min}(BB^{\Tsf})}}$. The  largest $c$ is
 \[
 c_{\max} =  \frac{c_{\alpha,\gamma}}{ \kappa_{B}  \sqrt{\kappa_{\theta_2}}}.
 \]
 which reduces to (3.16) in \cite{deng2012global} when $\alpha = 0$.



  \section{Numerical results}\label{section:numerical}
%\input{numResults.tex}
In this section, we demonstrate the  potential efficiency of our method sP-PRSM \eqref{equ:sP-PRSM} by solving the $l_1$-regularized least square problem:
\[\label{prob:lasso}
 \min_{x \in \R^n} \ \frac{1}{2}\|Ax-b\|_2^{2}+ \mu\|x\|_1, \quad \mu >0,
 \]
where $A \in \R^{m \times n}$ is the data matrix, $m$ is the number of the data points, $n$ is the number of features, $x \in \R^n$ is the vector of feature coefficients to be estimated and $\|x\|_1:=\sum\nolimits_{i=1}^{m}|x_i|$, $b \in \R^m$ is the observation vector and $\mu \in \R$ is the regularization parameter.
Given  $n$, $m = n/2$ and a  $p$-sparse vector $\bar{x}\in \R^m$ ($p$ is the number of nonzero elements in $x$ over $m$), the Matlab codes for generating  the data of \eqref{prob:lasso} are given as
\begin{verbatim}
       xbar = sprandn(n,1,p);  D = randn(m,n);
       D = D*spdiags(1./sqrt(sum(D.^2))',0,n,n); % normalize columns
       b = D*xbar + sqrt(0.001)*randn(m,1);
       mu = 0.1* norm(D'*b, 'inf').
\end{verbatim}


By introducing an auxiliary variable $y\in \R^{n}$, we reformulate the problem \eqref{prob:lasso} as
\[\label{prob:lasso2}
\min_{ x\in \R^{n},\  y\in \R^{n}} \  \ \frac{1}{2}\|Ax- b\|_2^{2}+ \mu\|y\|_1, \quad \st \quad  x-y=0.
\]
We consider to apply the sP-PRSM \eqref{equ:sP-PRSM} to solve \eqref{prob:lasso2}. In our implementation, we always choose $\Ptwo = 0$ and $\beta = 1$. Starting from $x^0 = y^0 = 0$ and $\lambda^0 = 0$,  with some suitably chosen  proximal matrix  $\Pone$,  the iterative scheme is given as
\[\label{lasso:qp_prsm}
\left\{
\begin{aligned}
x^{k+1}&= (A^{T}A+\beta I + S)^{-1}(A^{T}b+\beta y^{k}+\la^{k} + Sx^k) ,\\
\la^{k+\frac{1}{2}}&=\la^k-\alpha\beta(x^{k+1}-y^{k}),\\
y^{k+1}&=\S_{\mu/\beta}(x^{k+1}-\la^{k+\frac{1}{2}}/\beta),\\
\la^{k+1}&=\la^{k+\frac{1}{2}}-\gamma\beta(x^{k+1}-y^{k+1}),
\end{aligned}
\right.\]
where the shrinkage operator is defined as $\S_{\nu}(y)_i= \mathrm{sgn}(y_i) \cdot \max\{|y_i| - \nu, 0\}$ and $\mathrm{sgn}(y_i) = 1$ if $y_i \geq 0$ and $\mathrm{sgn}(y_i) = -1$ if $y_i < 0$.
\rev{According to \eqref{equ:stop}}, we stop the iterative scheme when
\[
  \max\{\|x^{k+1} - y^{k+1}\|, \|y^{k+1} - y^k\|\} \leq 10^{-6}
\]
or the iterative counter $k \geq 1000$.

%\subsection{$S = 0$}
Considering that the role of $\alpha$ and $\gamma$ in the sP-PRSM \eqref{equ:sP-PRSM} are equal, we can easily extend the convergence domain to
\[\label{equ:alpha:gamma2}
\alpha\in \left[0,\,\frac{1 + \sqrt{5}}{2} \right),\quad \gamma\in \left(0,\,\frac{1 -\alpha  + \sqrt{ (1 + \alpha)^2 + 4(1 - \alpha^2)}}{2}\right).\]
Note that \cite{he2015step} also establish the convergence of  the sP-PRSM \eqref{equ:sP-PRSM} with $\Pone = 0$ and $\Ptwo = 0$ when
\[ \label{equ:alpha:gamma:he2}
\alpha \in (-1, 0), \gamma \in (0, \alpha^2 - \alpha - 1),\  \mbox{or} \ \alpha \in (0,1), \gamma \in (\alpha^2 - \alpha - 1, 0).
\]
The union domain defined by \eqref{equ:alpha:gamma2} and \eqref{equ:alpha:gamma:he2} will be considered below.
We generate the mesh grid with respect to $\alpha$ and $\gamma$ by equally dividing  the intervals $[-1, 0]$ and  $[0,  1.618]$ into 10 parts, respectively.  We set $n = 2000$, $p = 0.2$. The corresponding scatter diagram is depicted in Figure \ref{figure:1}.  For consideration of space, only the detailed results for \eqref{equ:alpha:gamma2} are shown in Table \ref{results:1}.  From this table, we see that when $\alpha + \gamma = 1.618$ and $\min(\alpha, \gamma) \geq 0.485$, the sP-PRSM \eqref{equ:sP-PRSM} always perform  best; when $\alpha$ and $\gamma$ satisfy \eqref{equ:alpha:gamma:he2}, the results are not so good.

~\\
\begin{figure}[!htp]
\centering
\includegraphics[width=4.0in]{iter_lasso.eps}
\caption{$n = 2000$, $p=0.2$, the intervals $[-1,0]$ and $[0, 1.618]$ are both equally  divided into 10 parts. Note that the `iter'-axis and `alpha'-axis are both in reverse order.} \label{figure:1}
\end{figure}


\begin{table}[!htp]
\centering
\caption{The number of iterations taken by  sP-PRSM \eqref{equ:sP-PRSM}  with different $\alpha$ and $\gamma$ for solving \eqref{prob:lasso} with $n = 2000$, $p = 0.2$,  `-' means that the corresponding $\alpha$ and $\gamma$ do not satisfy \eqref{equ:alpha:gamma}.} \label{results:1}
\small
\begin{tabular}{cccccccccccc}\hline
\diagbox{$\gamma$}{$\alpha$}& 0.000& 0.162& 0.324& 0.485& 0.647& 0.809& 0.971& 1.133& 1.294& 1.456& 1.618\\ \hline
0.000& -- & 384 & 189 & 124 & 92 & 76 & 66 & 61 & 56 & 62 & 65 \\
0.162& 384 & 189 & 124 & 92 & 73 & 63 & 56 & 52 & 48 & 51 & -- \\
0.324& 189 & 124 & 92 & 73 & 61 & 55 & 49 & 45 & 43 & 60 & -- \\
0.485& 124 & 92 & 73 & 61 & 53 & 48 & 44 & 40 & 60 & -- & -- \\
0.647& 92 & 73 & 61 & 53 & 48 & 44 & 40 & 61 & 246 & -- & -- \\
0.809& 75 & 63 & 54 & 48 & 44 & 40 & 61 & 246 & -- & -- & -- \\
0.971& 64 & 55 & 48 & 43 & 40 & 60 & 246 & -- & -- & -- & -- \\
1.133& 60 & 52 & 45 & 40 & 60 & 243 & -- & -- & -- & -- & -- \\
1.294& 56 & 48 & 43 & 61 & 242 & -- & -- & -- & -- & -- & -- \\
1.456& 60 & 50 & 62 & -- & -- & -- & -- & -- & -- & -- & -- \\
1.618& 63 & -- & -- & -- & -- & -- & -- & -- & -- & -- & -- \\ \hline
\end{tabular}
\end{table}


Based on the above observation, we below focus on the choice of $\alpha = 0.618$ and $\gamma = 1$. For simplicity, we call the corresponding sP-PRSM \eqref{equ:sP-PRSM} as sP-PRSM*.  Recall that the classical ADMM always take the dual stepsize as 1 or 1.618 (named by ADMM-1, ADMM-2), covered by our sP-PRSM by taking $\alpha = 0, \gamma = 1$ or $\alpha = 0, \gamma = 1.618$, we next compare these three methods for solving \eqref{prob:lasso}.


We first consider the sP-PRSM \eqref{equ:sP-PRSM} without proximal terms, namely, $S = 0$. The comparison results are shown in Table \ref{table:s0}.  Comparing with ADMM-1, sP-PRSM* can always reduce the number of iterations in 40\%. While as compared to ADMM-2, sP-PRSM* can reduce the number of iterations at about 40\%  for the hard case when the sparsity level $p$ is 0.2 or 0.3, and lead to about 2-7\% reduction for the easy case when $p$ is 0.1.
\begin{table}[!htp]
\centering
\caption{Comparison among the number of iterations taken by ADMM-1, ADMM-2 and sP-PRSM* for solving \eqref{prob:lasso} with $\Pone = 0$} \label{table:s0}
\begin{tabular}{ccccccccc}\hline
$n$ &  $p$ &&  ADMM-1 & ADMM-2 &sP-PRSM* && ratio1(\%) & ratio2(\%) \\ \hline
2000 &  0.1 &&   65  &   43 &   40  && 61.5 & 93.0  \\
2000 &  0.2 &&   63  &   63 &   39  && 61.9 & 61.9  \\
2000 &  0.3 &&   85  &   93 &   53  && 62.4 & 57.0  \\
4000 &  0.1 &&   68  &   43 &   42  && 61.8 & 97.7  \\
4000 &  0.2 &&   67  &   62 &   42  && 62.7 & 67.7  \\
4000 &  0.3 &&   80  &   88 &   49  && 61.3 & 55.7  \\
8000 &  0.1 &&   70  &   44 &   43  && 61.4 & 97.7  \\
8000 &  0.2 &&   66  &   62 &   41  && 62.1 & 66.1  \\
8000 &  0.3 &&   69  &   79 &   42  && 60.9 & 53.2  \\
\hline
\end{tabular}
\end{table}

We now consider  the sP-PRSM \eqref{equ:sP-PRSM} with semi-positive proximal matrix $S$, which takes the form as
 \[\label{equ:s_semi}
 S = \beta\left( (\xi  - 1)I  - \frac{1}{2\beta}A^{\Tsf}A\right) \ \mbox{with}  \ \xi = 1.01 \lambda_{\max} ( I + \frac{1}{2\beta} A^{\Tsf}A).
\]
Now we name ADMM-1 and ADMM-2 as sPADMM-1 and sPADMM-2, respectively. The results are shown in Table \ref{table:s_pos}. Note that in this case, sPADMM-1 performs better than sPADMM-2, while sP-PRSM* is slightly better than sPADMM-1.
\begin{table}[!htp]
\centering
\caption{Comparison among the number of iterations taken by sPADMM-1, sPADMM-2 and sP-PRSM for solving \eqref{prob:lasso} with semi-positive $S$}\label{table:s_pos}
\begin{tabular}{ccccccccc}\hline
$n$ &  $p$ &&  sPADMM-1 & sPADMM-2 &sP-PRSM && ratio1(\%) & ratio2(\%) \\ \hline
2000 &  0.1 &&  236  &  245 &  231  && 97.9 & 94.3  \\
2000 &  0.2 &&  401  &  408 &  385  && 96.0 & 94.4  \\
2000 &  0.3 &&  574  &  579 &  549  && 95.6 & 94.8  \\
4000 &  0.1 &&  195  &  201 &  189  && 96.9 & 94.0  \\
4000 &  0.2 &&  395  &  402 &  380  && 96.2 & 94.5  \\
4000 &  0.3 &&  542  &  548 &  518  && 95.6 & 94.5  \\
8000 &  0.1 &&  210  &  220 &  206  && 98.1 & 93.6  \\
8000 &  0.2 &&  392  &  399 &  377  && 96.2 & 94.5  \\
8000 &  0.3 &&  491  &  498 &  471  && 95.9 & 94.6  \\ \hline
\end{tabular}
\end{table}

Finally, we  follow the way in \cite{li2014majorized} to generate  an indefinite proximal matrix $S$ as
\[\label{equ:s_inf}
S = \beta(\xi  - 1)I  - A^{\Tsf}A, \ \mbox{with}\
\xi = \lambda_{\max} ( I + \frac{1}{\beta} A^{\Tsf}A).
\]
Now we name ADMM-1 and ADMM-2 as iPADMM-1 and iPADMM-2, respectively.  We list the results in Table \ref{figure:s_inf}.  From this table, we can see that sP-PRSM* performs best among the three methods. On average, it can lead to 5\% reduction in the number of iterations as compared to iPADMM-1 and 10\% reduction as compared to iPADMM-2.  We also remark that for our test instances, no matter choosing the  semi-definite or indefinite $S$, ADMM-1 is always better than ADMM-2. Besides, for each method, the version with indefinite $S$ in \eqref{equ:s_inf} can always bring out 40-50\% improvement in the number of iterations as compared to the version with semi-definite $S$ in \eqref{equ:s_semi}. This coincides with the observation in \cite{li2014majorized}.

\begin{table}[!htp]
\centering
\caption{Comparison among the number of iterations taken by iPADMM-1, iPADMM-2 and sP-PRSM for solving \eqref{prob:lasso} with indefinite $S$}\label{figure:s_inf}
\begin{tabular}{ccccccccc}\hline
$n$ &  $p$ &&  iPADMM-1 & iPADMM-2 &sP-PRSM && ratio1(\%) & ratio2(\%) \\ \hline
2000 &  0.1 &&  124  &  138 &  122  && 98.4 & 88.4  \\
2000 &  0.2 &&  231  &  240 &  216  && 93.5 & 90.0  \\
2000 &  0.3 &&  342  &  348 &  315  && 92.1 & 90.5  \\
4000 &  0.1 &&  105  &  116 &  101  && 96.2 & 87.1  \\
4000 &  0.2 &&  228  &  237 &  213  && 93.4 & 89.9  \\
4000 &  0.3 &&  322  &  328 &  297  && 92.2 & 90.5  \\
8000 &  0.1 &&  107  &  122 &  106  && 99.1 & 86.9  \\
8000 &  0.2 &&  226  &  235 &  211  && 93.4 & 89.8  \\
8000 &  0.3 &&  288  &  296 &  267  && 92.7 & 90.2  \\
\hline
\end{tabular}
\end{table}

\medskip

From Table 4 we can observe that we can use indefinite matrices in the subproblems \eqref{equ:sP-PRSM1} and 
\eqref{equ:sP-PRSM3}, which can further improve the efficiency of the algorithm. In fact, under certain conditions on the problems' data and suitable requirements on the matrices, we can theoretically establish its convergence. Since the analysis is very similar to those in this paper, we do not describe them here.   
\section{Conclusions}\label{section:conclusions}

In this paper, we proposed a modification of the Peaceman-Rachford splitting method by introducing two different parameters in updating the dual variable, and by introducing semi-proximal terms to the subproblems in updating the primal variables. We established the relationship between the two parameters under which we proved the global convergence of the algorithm. We also analyzed the sublinear rate convergence under ergodic and non-ergodic senses, respectively. Under further conditions of one of the objective functions, we proved the linear convergence of the algorithm. Finally, we reported extensive numerical results, indicating the efficiency of the proposed algorithm.

Note that the parameters  $\alpha$ and $\gamma$ are essential to the efficiency of the algorithm, which should be variable along with the iteration. Allowing the parameter $\alpha$ and $\gamma$  varying with the process of the iterate may give us the freedom of choosing them in a self-adaptive manner.
Suitable updating rules are among our future research tasks.

On the other hand, considering that the main task in the algorithm is to solve the $x$-optimization problem (e.g., \eqref{he1}) or the $y$-optimization problem (e.g., \eqref{he3}), solving them in an  inexact  manner may improve the efficiency of the algorithm. An approximate version of the proposed sP-PRSM with practical accuracy criteria is also our future research topic.
\small
%\bibliographystyle{siam}
%\bibliography{SPRSM}


\begin{thebibliography}{10}

\bibitem{boyd2011distributed}
{\sc S.~Boyd, N.~Parikh, E.~Chu, B.~Peleato, and J.~Eckstein}, {\em Distributed
  optimization and statistical learning via the alternating direction method of
  multipliers}, Foundations and Trends{\textregistered} in Machine Learning, 3
  (2011), pp.~1--122.

\bibitem{corman2014generalized}
{\sc E.~Corman and X.~Yuan}, {\em A generalized proximal point algorithm and
  its convergence rate}, SIAM Journal on Optimization, 24 (2014),
  pp.~1614--1638.

\bibitem{deng2012global}
{\sc W.~Deng and W.~Yin}, {\em On the global and linear convergence of the
  generalized alternating direction method of multipliers}, tech. report, DTIC
  Document, 2012.

\bibitem{douglas1956numerical}
{\sc J.~Douglas and H.~Rachford}, {\em On the numerical solution of heat
  conduction problems in two and three space variables}, Transactions of the
  American mathematical Society,  (1956), pp.~421--439.

\bibitem{Eck1994} {\sc J. Eckstein and M. Fukushima},
{\em Some reformulation and applications of the alternating directions method
of multipliers}. In: Hager, W.W. et al., eds., Large Scale Optimization:
State of the Art, pp. 115-134. Kluwer Academic Publishers, 1994.

\bibitem{facchinei2003finite}
{\sc F.~Facchinei and J.-S. Pang}, {\em Finite-dimensional variational
  inequalities and complementarity problems}, vol.~1, Springer, 2003.

\bibitem{fazel2013hankel}
{\sc M.~Fazel, T.~K. Pong, D.~Sun, and P.~Tseng}, {\em Hankel matrix rank
  minimization with applications to system identification and realization},
  SIAM Journal on Matrix Analysis and Applications, 34 (2013), pp.~946--977.

\bibitem{gabay1976dual}
{\sc D.~Gabay and B.~Mercier}, {\em A dual algorithm for the solution of
  nonlinear variational problems via finite element approximation}, Computers
  \& Mathematics with Applications, 2 (1976), pp.~17--40.

\bibitem{glowinski1984numerical}
{\sc R.~Glowinski}, {\em Numerical Methods for Nonlinear Variational Problems},
  Springer, 1984.

\bibitem{glowinski1975approximation}
{\sc R.~Glowinski and A.~Marroco}, {\em Sur l'approximation, par
  {\'e}l{\'e}ments finis d'ordre un, et la r{\'e}solution, par
  p{\'e}nalisation-dualit{\'e} d'une classe de probl{\`e}mes de dirichlet non
  lin{\'e}aires}, ESAIM: Mathematical Modelling and Numerical
  Analysis-Mod{\'e}lisation Math{\'e}matique et Analyse Num{\'e}rique, 9
  (1975), pp.~41--76.

\bibitem{HLHY2002} {\sc B.S. He, L.Z. Liao, D.R. Han, H. Yang}, {\em A new inexact alternating directions method for
monotone variational inequalities}, Mathematical Programming, 92 (2002), pp.~103--118.

\bibitem{he2014strictly}
{\sc B.~He, H.~Liu, Z.~Wang, and X.~Yuan}, {\em A strictly contractive
  {P}eaceman--{R}achford splitting method for convex programming}, SIAM Journal
  on Optimization, 24 (2014), pp.~1011--1040.

\bibitem{he2015step}
{\sc B.~He, F.~Ma, and X.~Yuan}, {\em On the step size of symmetric alternating
  directions method of multipliers}, tech. report, Optimization online, 2015.

\bibitem{he2012convergence}
{\sc B.~He and X.~Yuan}, {\em On the ${O}(1/n)$ convergence rate of the
  {D}ouglas-{R}achford alternating direction method}, SIAM Journal on Numerical
  Analysis, 50 (2012), pp.~700--709.

\bibitem{li2014majorized}
{\sc M.~Li, D.~Sun, and K.-C. Toh}, {\em A majorized admm with indefinite
  proximal terms for linearly constrained convex composite optimization}, arXiv
  preprint arXiv:1412.1911,  (2014).

\bibitem{nesterov1983method}
{\sc Y.~Nesterov}, {\em A method of solving a convex programming problem with
  convergence rate ${O}(1/k^2)$}, in Soviet Mathematics Doklady, vol.~27, 1983,
  pp.~372--376.

\bibitem{nesterov2007gradient}
\leavevmode\vrule height 2pt depth -1.6pt width 23pt, {\em Gradient methods for
  minimizing composite objective function}, 2007.

\bibitem{peaceman1955numerical}
{\sc D.~W. Peaceman and H.~H. Rachford, Jr}, {\em The numerical solution of
  parabolic and elliptic differential equations}, Journal of the Society for
  Industrial \& Applied Mathematics, 3 (1955), pp.~28--41.

\bibitem{xu2011class}
{\sc M.~Xu and T.~Wu}, {\em A class of linearized proximal alternating
  direction methods}, Journal of Optimization Theory and Applications, 151
  (2011), pp.~321--337.

\end{thebibliography}

\end{document}



1506.02221.extracted_bib
{boyd2011distributed {\sc S.~Boyd, N.~Parikh, E.~Chu, B.~Peleato, and J.~Eckstein, {\em Distributed optimization and statistical learning via the alternating direction method of multipliers, Foundations and Trends{\textregistered in Machine Learning, 3 (2011), pp.~1--122.
1506.02221.extracted_bib
{corman2014generalized {\sc E.~Corman and X.~Yuan, {\em A generalized proximal point algorithm and its convergence rate, SIAM Journal on Optimization, 24 (2014), pp.~1614--1638.
1506.02221.extracted_bib
{deng2012global {\sc W.~Deng and W.~Yin, {\em On the global and linear convergence of the generalized alternating direction method of multipliers, tech. report, DTIC Document, 2012.
1506.02221.extracted_bib
{douglas1956numerical {\sc J.~Douglas and H.~Rachford, {\em On the numerical solution of heat conduction problems in two and three space variables, Transactions of the American mathematical Society, (1956), pp.~421--439.
1506.02221.extracted_bib
{Eck1994 {\sc J. Eckstein and M. Fukushima, {\em Some reformulation and applications of the alternating directions method of multipliers. In: Hager, W.W. et al., eds., Large Scale Optimization: State of the Art, pp. 115-134. Kluwer Academic Publishers, 1994.
1506.02221.extracted_bib
{facchinei2003finite {\sc F.~Facchinei and J.-S. Pang, {\em Finite-dimensional variational inequalities and complementarity problems, vol.~1, Springer, 2003.
1506.02221.extracted_bib
{fazel2013hankel {\sc M.~Fazel, T.~K. Pong, D.~Sun, and P.~Tseng, {\em Hankel matrix rank minimization with applications to system identification and realization, SIAM Journal on Matrix Analysis and Applications, 34 (2013), pp.~946--977.
1506.02221.extracted_bib
{gabay1976dual {\sc D.~Gabay and B.~Mercier, {\em A dual algorithm for the solution of nonlinear variational problems via finite element approximation, Computers \& Mathematics with Applications, 2 (1976), pp.~17--40.
1506.02221.extracted_bib
{glowinski1984numerical {\sc R.~Glowinski, {\em Numerical Methods for Nonlinear Variational Problems, Springer, 1984.
1506.02221.extracted_bib
{glowinski1975approximation {\sc R.~Glowinski and A.~Marroco, {\em Sur l'approximation, par {\'el{\'ements finis d'ordre un, et la r{\'esolution, par p{\'enalisation-dualit{\'e d'une classe de probl{\`emes de dirichlet non lin{\'eaires, ESAIM: Mathematical Modelling and Numerical Analysis-Mod{\'elisation Math{\'ematique et Analyse Num{\'erique, 9 (1975), pp.~41--76.
1506.02221.extracted_bib
{HLHY2002 {\sc B.S. He, L.Z. Liao, D.R. Han, H. Yang, {\em A new inexact alternating directions method for monotone variational inequalities, Mathematical Programming, 92 (2002), pp.~103--118.
1506.02221.extracted_bib
{he2014strictly {\sc B.~He, H.~Liu, Z.~Wang, and X.~Yuan, {\em A strictly contractive {Peaceman--{Rachford splitting method for convex programming, SIAM Journal on Optimization, 24 (2014), pp.~1011--1040.
1506.02221.extracted_bib
{he2015step {\sc B.~He, F.~Ma, and X.~Yuan, {\em On the step size of symmetric alternating directions method of multipliers, tech. report, Optimization online, 2015.
1506.02221.extracted_bib
{he2012convergence {\sc B.~He and X.~Yuan, {\em On the ${O(1/n)$ convergence rate of the {Douglas-{Rachford alternating direction method, SIAM Journal on Numerical Analysis, 50 (2012), pp.~700--709.
1506.02221.extracted_bib
{li2014majorized {\sc M.~Li, D.~Sun, and K.-C. Toh, {\em A majorized admm with indefinite proximal terms for linearly constrained convex composite optimization, arXiv preprint arXiv:1412.1911, (2014).
1506.02221.extracted_bib
{nesterov1983method {\sc Y.~Nesterov, {\em A method of solving a convex programming problem with convergence rate ${O(1/k^2)$, in Soviet Mathematics Doklady, vol.~27, 1983, pp.~372--376.
1506.02221.extracted_bib
{nesterov2007gradient \leavevmode\vrule height 2pt depth -1.6pt width 23pt, {\em Gradient methods for minimizing composite objective function, 2007.
1506.02221.extracted_bib
{peaceman1955numerical {\sc D.~W. Peaceman and H.~H. Rachford, Jr, {\em The numerical solution of parabolic and elliptic differential equations, Journal of the Society for Industrial \& Applied Mathematics, 3 (1955), pp.~28--41.
1506.02221.extracted_bib
{xu2011class {\sc M.~Xu and T.~Wu, {\em A class of linearized proximal alternating direction methods, Journal of Optimization Theory and Applications, 151 (2011), pp.~321--337.