\section{Introduction}
In this paper, we consider the problem
\begin{equation}\label{eq: prob_temp}
\min_{x\in\mathcal{X}} f(x) + g(x) + h(Ax),
\end{equation}
where $f, g\colon \mathcal{X} \to \mathbb{R}\cup \{+\infty\}$ and $h\colon \mathcal{Y} \to \mathbb{R}\cup\{+\infty\}$ are proper, lower semicontinuous, convex functions, $A\colon \mathcal{X} \to \mathcal{Y}$ is a linear operator.
$\mathcal{X}$ and $\mathcal{Y}$ are Euclidean spaces such that $\mathcal{X} = \prod_{i=1}^n \mathcal{X}_i$, and $\mathcal{Y}=\prod_{j=1}^m \mathcal{Y}_j$.
Moreover, $f$ is assumed to have coordinatewise Lipschitz continuous gradients and $g, h$ admit easily computable proximal operators.



\begin{table*}[t]
  \centering
  \begin{tabular}{| l | l | c | c | c | c |}
 \hline
      									& \shortstack{Step sizes\\ with dense data} 				& \shortstack{per iteration \\ cost} 			& \shortstack{block-wise \\ Lipschitz}  &\shortstack{probability \\ law} & \shortstack{Efficient \\ implementation} \\\hline
    \cite{chambolle2018stochastic} 			& \cellcolor{lightgray}{$n\tau_i \sigma \| A_i \|^2 < 1$} 		& $m$ 								& N/A 						& \cellcolor{lightgray}{arbitrary}	& \cellcolor{lightgray}{direct$^{\dagger}$} \\ \hline
        \cite{fercoq2019coordinate} 			& \centering{$n^2\tau_i\sigma\|A_i\|^2 < 1$} 						& \cellcolor{lightgray}{$\vert J(i) \vert^{*}$} 		& \cellcolor{lightgray}{Yes} & uniform & \cellcolor{lightgray}{\shortstack{direct or duplication}}  \\ \hline
                \cite{latafat2019new} 				& $n^2\tau_i\sigma\|A_i\|^2 < 1$ 						& \cellcolor{lightgray}{$\vert J(i) \vert^{*}$} 		& No 				& \cellcolor{lightgray}{arbitrary} & duplication \\ \hline
    PURE-CD 							& \cellcolor{lightgray}{$n\tau_i\sigma\|A_i\|^2 < 1$} 			& \cellcolor{lightgray}{$\vert J(i) \vert^{*}$} 		& \cellcolor{lightgray}{Yes} & \cellcolor{lightgray}{arbitrary} & \cellcolor{lightgray}{direct} \\
    \hline
  \end{tabular}
  \caption{\small{Comparison of primal-dual coordinate descent methods. Note that we only compare here the most related methods to ours and include a comprehensive review of other existing methods with comparison to PURE-CD in Section~\ref{sec: lit}. 
  In the last column, we refer to the way one needs to implement the algorithm, for it to be efficient in both sparse and dense settings.
  $^{*}J(i)$ is defined in~\eqref{eq: def_ji}. $^{\dagger}$SPDHG only has implementation for dense setting and not for sparse. The concept of duplication for PDCD is described in \cite{fercoq2019coordinate}.}}
  \label{tab:1}
\end{table*}

Problem~\eqref{eq: prob_temp} is a general template that covers many problems in different fields, such as regularized empirical risk minimization~\cite{shalev2013stochastic,zhang2017stochastic}, optimization with large number of constraints~\cite{patrascu2017nonasymptotic,fercoq2019almost}, and total variation (TV) regularized problems~\cite{chambolle2018stochastic,fercoq2019coordinate}.

The classic choice for solving problem~\eqref{eq: prob_temp} is to use primal-dual methods~\cite{chambolle2011first,vu2013splitting,condat2013primal}.
These methods utilize the proximal operators for $g, h^\ast$ and gradient of the differentiable component $f$.
Randomized versions that we refer to as primal dual coordinate descent (PDCD), are proposed in several works~\cite{zhang2017stochastic,dang2014randomized,gao2019randomized,fercoq2019coordinate,chambolle2018stochastic,latafat2019new}.

First advantage of coordinate-based methods is that they access to blocks of $A$ and update a subset of variables, resulting in cheap per iteration costs.
Moreover, they utilize larger step sizes depending on the properties of the problem in selected blocks.

Existing PDCD methods fail to retain both these advantages, as sparsity of $A$ varies.
In particular, methods that have cheap per-iteration costs with sparse $A$~\cite{fercoq2019coordinate,latafat2019new}, are restricted to use small step sizes with dense $A$.
On the other hand, methods that can use large step sizes with dense $A$~\cite{chambolle2018stochastic}, have high per-iteration costs with sparse $A$.

\textbf{Contributions.} In this paper, we identify random extrapolation as the key to design a method that combines the benefits of the methods in two camps and propose the primal-dual method with random extrapolation and coordinate descent (PURE-CD).
PURE-CD exhibits the advantages of~\cite{fercoq2019coordinate,latafat2019new} in the sparse setting and the advantages of~\cite{chambolle2018stochastic} in the dense setting simultaneously, achieving the best of both worlds.
As PURE-CD has the favorable properties in both ends of the spectrum, it has the best performance in the regime in between: moderately sparse data.

In addition to adapting to the sparsity of $A$, we prove that PURE-CD also adapts to unknown structures in the problem, and obtains linear rate of convergence, without any modifications in the step sizes.
Our linear convergence results apply to strongly convex-strongly concave problems, linear programs, and problems with piecewise linear quadratic functions, involving Lasso, support vector machines and linearly constrained problems with piecewise linear-quadratic objectives.
In the general convex case, we prove that the iterates of PURE-CD converges almost surely to a solution of problem~\eqref{eq: prob_temp}.
Moreover, we show that in this case, the ergodic sequence obtains the optimal $\mathcal{O}(1/k)$ sublinear rate of convergence.



\section{Preliminaries}
\subsection{Notation}
For a positive definite matrix $V$, we denote $\| x \|_V^2 = \langle x, Vx \rangle$.
We define the distance of a point $x$ from a set $\mathcal{X}$ as $\dist(x, \mathcal{X}) = \min_{u\in\mathcal{X}} \| u - x \|^2$.
Given an index $i\in\{1,\ldots, n\}$, the corresponding coordinate of the gradient vector is $\nabla _i f(x)$.
Graph of mapping $F$ is denoted by $\mathrm{gra~} F$.
Recall that $\mathcal{X} = \prod_{i=1}^n \mathcal{X}_i$, and $\mathcal{Y}=\prod_{j=1}^m \mathcal{Y}_j$ and $\mathcal{Z}=\mathcal{X}\times\mathcal{Y}$.
For $u\in\mathcal{X}_i$, $U_i(u)\in\mathcal{X}$ is such that each element of $U_i(u)$ is $0$, except the block $i$ which contains $u$.
We denote the indicator function of a set $\mathcal{X}$ as $\delta_\mathcal{X}$.

Proximal operator with a positive definite $V$ is defined as
\begin{equation}
\prox_{V, g}(x) = \arg\min_{u} g(u) + \frac{1}{2} \| u - x \|^2_{V^{-1}}.\notag
\end{equation}
We will need the following notation for the sparse setting,
\begin{equation}\label{eq: def_ji}
\begin{aligned}
J(i) = \{ j\in\{1, \dots, m\}\colon A_{j, i} \neq 0 \} \\
I(j) = \{ i\in\{1, \dots, n\}\colon A_{j, i} \neq 0 \}.
\end{aligned}
\end{equation}
Given a matrix $A$ and $i \in\{1,\ldots,n\}$, $J(i)$ denotes the row indices that correspond to nonzero values in the column indexed by $i$.
Similarly, with $j\in\{1,\ldots,m\}$, $I(j)$ gives the column indices corresponding to nonzero values in the row indexed by $j$.

Moreover, given positive probabilities $(p_i)_{1 \leq i \leq n}$, we define
\begin{equation}
\pi_j = \sum_{i\in I(j)} p_i.
\end{equation}
In the simple case of $p_i = 1/n$, it is easy to see that $n \pi_j$ corresponds to number of nonzeros in the row indexed by $j$.

At iteration $k$, the algorithm randomly picks an index $i_{k+1}\in\{1,\ldots, n\}$. 
To govern the selection rule, we define the probability matrix $P=\diag(p_1, \dots, p_n)$, where $p_i = \mathbb{P}(i_{k+1} = i)$, and $\underline p = \min_i p_i$. 
We define as $\mathcal{F}_k$ the filtration generated by the random indices $\{i_1, \dots, i_k\}$.

Denoting $z=(x,y)$, we define the functions
\begin{align}
&D_p(x_{k+1}; z) = f(x_{k+1}) + g(x_{k+1}) - f(x) - g(x) \notag\\
&\qquad\qquad\qquad\qquad\qquad\qquad+ \langle A^\top y, x_{k+1} - x \rangle,\notag\\
&D_d(\bar{y}_{k+1}; z) = h^\ast(\bar{y}_{k+1}) - h^\ast(y) - \langle Ax, \bar{y}_{k+1} - y \rangle.\notag
\end{align}

\subsection{Optimality}
Problem~\eqref{eq: prob_temp} has the following saddle point formulation
\begin{align}
\min_{x\in\mathcal{X}} \max_{y\in\mathcal{Y}} f(x) + g(x) + \langle Ax, y \rangle - h^\ast(y).\notag
\end{align}
Karush-Kuhn-Tucker (KKT) conditions state that the vector $z_\star = (x_\star, y_\star)$ is a primal-dual solution of the problem when
\begin{align}\label{eq: kkt}
0 \in \begin{bmatrix} \nabla f(x_\star) + \partial g(x_\star) + A^\top y_\star \\
Ax_\star - \partial h^\ast(y_\star) 
\end{bmatrix} =: F(z_\star).
\end{align}
We call $\mathcal{Z}_\star$ the set of such solutions.


\subsection{Metric subregularity}\label{sec: ms}
We utilize the metric subregularity assumption for proving linear convergence.
This assumption has been used in primal-dual optimization literature for both deterministic~\cite{liang2016convergence} and randomized algorithms~\cite{latafat2019new,alacaoglu2019convergence}.
\begin{definition}\label{eq: def_ms}
A set valued mapping $F\colon \mathcal{X} \rightrightarrows \mathcal{Y}$ is metrically subregular at $\bar{x}$ for $\bar{y}$, with $(\bar{x}, \bar{y}) \in \mathrm{gra~} F$, if there exists $\eta > 0$ with a neighborhood of regularity $\mathcal{N}(\bar{x})$ such that
\begin{equation}
\dist(x, F^{-1} \bar{y}) \leq \eta \dist(\bar{y}, Fx),~~~~\forall x \in \mathcal{N}(\bar{x}).\notag
\end{equation}
\end{definition}
We will be interested in the metric subregularity of KKT operator $F$, given in~\eqref{eq: kkt}, for $0$. Intuitively speaking, as $0 \in F(z_\star),\forall z^\star \in \mathcal{Z}_\star$, metric subregularity of $F$ for $0$ essentially gives us a way to characterize the behavior of the iterates around the solution set.

Even though Definition~\ref{eq: def_ms} looks daunting, fortunately, one does not need to check it for a given problem.
Metric subregularity is well-studied in the literature and it is known to be satisfied in the following cases:
\begin{example}\label{ex: ms_examples}~\\
$\triangleright$ If $f+g$ and $h^\ast$ are strongly convex, Definition~\ref{eq: def_ms} holds with $\mathcal{N}(\bar{x}) = \mathbb{R}^d$~\citep[Lemma IV.2]{latafat2019new}. \\
$\triangleright$ If $f, g, h$ are piecewise linear quadratic (PLQ) functions, Definition~\ref{eq: def_ms} holds with any bounded neighborhood $\mathcal{N}(\bar{x})$~\citep[Lemma IV.4]{latafat2019new}.
\end{example}
PLQ functions include $\ell_1$ norm, hinge loss, indicator of polyhedral sets. Therefore third bullet point apply to Lasso, support vector machines, elastic net, and linearly constrained problems with PLQ loss functions~\cite{latafat2019new}.

We now state our main assumptions which are standard in the literature~\cite{fercoq2019coordinate,chambolle2018stochastic,latafat2019new,bauschke2011convex}:
\begin{assumption}\label{asmp: asmp1}~\\
$\triangleright$ $f$, $g$ and $h$ are proper, lower semicontinuous, convex. \\[1mm]
$\triangleright$ $g$ is separable, i.e., $g(x) = \sum_{i=1}^n g_i(x_i)$, and $f$ has coordinatewise Lipschitz gradients such that $\forall x \in \mathcal{X}, \forall u \in \mathcal{X}_i$,
\begin{equation}
f(x + U_i(u)) \leq f(x) + \langle \nabla_i f(x), u \rangle + \frac{\beta_i}{2} \| u \|^2.
\end{equation}
$\triangleright$ Set of solutions to problem~\eqref{eq: prob_temp}, defined in~\eqref{eq: kkt} is nonempty. \\[1mm]
$\triangleright$ Slater's condition holds, which states that $0 \in \mathrm{ri}(\dom h - A\dom g)$ where $\mathrm{ri}$ denotes the relative interior.
\end{assumption}

\section{Algorithm}\label{sec: alg}
In this section, we will sketch the main ideas behind our algorithm.
Primal-dual method\footnote{This method is also known as V\~u-Condat algorithm.}, due to~\cite{chambolle2011first,condat2013primal,vu2013splitting} reads as
\begin{equation}\label{eq: vc_alg}
\begin{aligned}
&\bar{x}_{k+1} = \prox_{\tau, g} \left(\bar{x}_k - \tau\left( \nabla f(\bar x_k) + A^\top \bar{y}_k \right) \right)  \\
&\bar{y}_{k+1} = \prox_{\sigma, h^\ast} \left( \bar{y}_k + \sigma A (2\bar{x}_{k+1} - \bar{x}_k)\right).
\end{aligned}
\end{equation}
The main intuition behind PDCD methods proposed by~\cite{zhang2017stochastic,fercoq2019coordinate,chambolle2018stochastic} is to incorporate coordinate based updates.
Among these methods,~\cite{zhang2017stochastic} specializes in strongly convex-strongly concave problems, whereas the other other ones apply to more general classes of problems.

A closely related approach focused on the following interpretation of primal-dual method~\eqref{eq: vc_alg} which is named as TriPD in~\citep[Algorithm 1]{latafat2019new}
\begin{equation}\label{eq: tripd}
\begin{aligned}
&\bar{y}_{k+1} = \prox_{\sigma, h^\ast}\left( \hat{y}_k + \sigma A \bar{x}_k \right)  \\
&\bar{x}_{k+1} = \prox_{\tau, g}\left(\bar{x}_k - \tau\left( \nabla f(\bar{x}_k) + A^\top \bar{y}_{k+1} \right) \right)  \\
&\hat{y}_{k+1} = \bar{y}_{k+1} + \sigma A(\bar{x}_{k+1} - \bar{x}_k). 
\end{aligned}
\end{equation}
We notice that by moving the $\bar{y}_{k+1}$ update in TriPD to take place after $\hat{y}_{k+1}$ update, one obtains~\eqref{eq: vc_alg}.

As observed in~\cite{latafat2019new}, this particular interpretation of primal-dual method is useful for randomization.
TriPD-BC as proposed in~\cite{latafat2019new} iterates as
\begin{align}
&\bar{y}_{k+1} = \prox_{\sigma, h^\ast}\left( {y}_k + \sigma A {x}_k \right) \notag \\
&\bar{x}_{k+1} = \prox_{\tau, g}\left({x}_k - \tau\left( \nabla f(x_k) + A^\top \bar{y}_{k+1} \right) \right) \notag \\
&\hat{y}_{k+1} = \bar{y}_{k+1} + \sigma A(\bar{x}_{k+1} - x_k) \notag \\
&\text{Draw an index } i_{k+1}\in\{1,\dots,n\} \text{ randomly.} \notag \\
&x_{k+1}^{i_{k+1}} = \bar{x}_{k+1}^{i_{k+1}},~~~x_{k+1}^j = x_k^j, \forall j\neq i_{k+1}\notag\\
&y_{k+1}^j = \hat{y}_{k+1}^j, \forall j \in J(i_{k+1}),~~~y_{k+1}^j = y_k^j, \forall j\not\in J(i_{k+1}). \notag
\end{align}
One immediate limitation of TriPD-BC is that to update $y_{k+1}$, one needs to know $\bar{x}_{k+1}$, whereas only $\bar{x}_{k+1}^{i_{k+1}}$ is needed to update $x_{k+1}$.
As also discussed in~\cite{latafat2019new}, this scheme is suitable when $A$ has special structure such as sparsity.
When $A$ is dense, one needs to update all elements of $y_{k+1}$ and $\hat{y}_{k+1}$, in which case one needs to compute both $\bar{y}_{k+1}$ and $\bar{x}_{k+1}$ which has the same cost as a deterministic algorithm.

In the dense setting, for an efficient implementation, one can use duplication of dual variables as described in~\cite{fercoq2019coordinate}. 
However, in this case one is restricted to use small step sizes as discussed in~\cite{fercoq2019coordinate}. Compared to SPDHG in~\cite{chambolle2018stochastic}, the step sizes can be $n$ times worse, deteriorating the performance of the method considerably in the dense setting.

On the other hand, the drawback of~SPDHG is that it needs to update all dual variables at every iteration, whereas the methods in~\cite{fercoq2019coordinate,latafat2019new} update only a subset of dual variables depending on the sparsity of $A$. When the dual dimension is high, per iteration cost of~\cite{chambolle2018stochastic} becomes prohibitive.

Our idea, inspired by~\cite{chambolle2018stochastic}, to make TriPD-BC efficient for dense setting is to use $x_{k+1}$ rather than $\bar{x}_{k+1}$ in the update of $\hat{y}_{k+1}$.
Although simple to state, this modification makes $\hat{y}_{k+1}$ random, rendering the analysis of~\cite{latafat2019new} and other analyses based on monotone operator theory not applicable.

This leads to our algorithm, primal-dual method with random extrapolation and coordinate descent (PURE-CD). Our method uses large step sizes as in~\cite{chambolle2018stochastic} in the dense setting, while staying efficient in terms of per iteration costs in the sparse setting as in~\cite{fercoq2019coordinate,latafat2019new}; leading to the first general PDCD algorithm that obtains favorable properties in both sparse and dense settings.

\begin{algorithm}[h]
\caption{Primal-dual method with random extrapolation and coordinate descent (PURE-CD)}
\label{alg:stripd}
\begin{algorithmic}[1]
    \STATE { \textbf{Input:}} Diagonal matrices $\theta, \tau, \sigma > 0$, chosen according to~\eqref{eq: theta_choice},~\eqref{eq: ss_choice}.
    \FOR{$k = 0,1\ldots $} 
	\STATE $\bar{y}_{k+1} = \prox_{\sigma, h^\ast}\left(y_k + \sigma Ax_k\right)$
         \STATE $\bar{x}_{k+1} = \prox_{\tau, g}\left(x_k - \tau \left(\nabla f(x_k) + A^\top \bar{y}_{k+1}\right)\right)$
        \STATE Draw $i_{k+1} \in \{ 1, \dots, n \}$ with $\mathbb{P}(i_{k+1} = i) = p_i$
         \STATE $x_{k+1}^{i_{k+1}} = \bar{x}_{k+1}^{i_{k+1}}$ 
         \STATE $x_{k+1}^{j} = x_k^j, \forall j \neq i_{k+1}$
         \STATE $y_{k+1}^{j} = \bar{y}_{k+1}^j + \sigma_j \theta_j(A (x_{k+1} - x_k))_j, \forall j \in J(i_{k+1})$ $y_{k+1}^j = y_k^j, \forall j\not\in J(i_{k+1})$
	\ENDFOR
\end{algorithmic}	
\end{algorithm}


\section{Convergence Analysis}
In this section, we  analyze the behavior of Algorithm~\ref{alg:stripd} under various assumptions.
We first start with a lemma analyzing one iteration of the algorithm.
\vspace{0.05cm}
\begin{lemma}\label{lem: lem1}
Under Assumption~\ref{asmp: asmp1}, let $\theta = \diag(\theta_1, \ldots, \theta_m)$ and $\pi = \diag(\pi_1, \ldots, \pi_m)$ be chosen as
\begin{equation}
\theta_j = \frac{\pi_j}{\underline p}, \text{ where } \pi_j = \sum_{i \in I(j)} p_i, \text{ and } \underline p = \min_i p_i.\label{eq: theta_choice}
\end{equation}
We define the functions, given $z=(x, y)$,
\begin{align}
&V(z) = \frac{\underline p}{2} \| x \|^2_{\tau^{-1}P^{-1}} + \frac{\underline p}{2} \| y \|^2_{\sigma^{-1}\pi^{-1}},\notag\\
&\tilde{V}(z) = \frac{\underline p}{2} \| x \|^2_{C(\tau)} + \frac{\underline p}{2} \| y \|^2_{\sigma^{-1}},\notag
\end{align}
where 
$
C(\tau)_i = \frac{2p_i}{\underline p \tau_i} - \frac{1}{\tau_i} - p_i \sum_{j=1}^m \pi_j^{-1}\sigma_j \theta_j^2 A_{j, i}^2 - \frac{\beta_i p_i}{\underline p}$.

Then, for the iterates of Algorithm~\ref{alg:stripd}, $\forall z=(x, y) \in \mathcal{Z}$, it holds that:
\begin{multline}
\mathbb{E}_k \left[ D_p(x_{k+1}; z) \right] +\underline p D_d(\bar{y}_{k+1}; z) + \mathbb{E}_k \left[ V(z_{k+1} - z) \right] \\ \leq (1-\underline p)D_p (x_k; z) + V(z_k - z) - \tilde{V}(\bar{z}_{k+1} - z_k).\notag
\end{multline}
\end{lemma}
The main technical challenge in the proof of the lemma, compared to the corresponding results in~\cite{latafat2019new} and~\cite{chambolle2018stochastic} is handling stochasticity in both variables $x_{k+1}, y_{k+1}$ (and also $\hat{y}_{k+1}$ for~\cite{latafat2019new}).
Using coordinatewise Lipschitz constants of $f$ with arbitrary sampling also requires an intricate analysis.

The result of Lemma~\ref{lem: lem1} is promising for deriving convergence results for Algorithm~\ref{alg:stripd}. 
As $D_p(x_{k+1}; z_\star) \geq 0$, $D_d(\bar{y}_{k+1}; z_\star)\geq0$ and when step sizes are chosen such that $\tilde{V}$ is nonnegative, Lemma~\ref{lem: lem1} describes a stochastic monotonicity property.
In particular, it shows that $D_p(x_{k+1}; z_\star) + V(z_{k+1} - z_\star)$ which measures the distance to solution in a Bregman distance sense, is monotonically nonincreasing in expectation.

\subsection{Almost sure convergence}\label{sec: as}
Almost sure convergence is a fundamental property for randomized methods describing the limiting behavior of the iterates in different realization of the algorithm.
The following theorem states that the iterates of Algorithm~\ref{alg:stripd} converge almost surely to a point in the solution set.
\begin{theorem}\label{thm: thm1}
Let Assumption~\ref{asmp: asmp1} hold, $\theta$ and $\pi$ be as in~\Cref{lem: lem1}, and the step sizes $\tau, \sigma$ satisfy
\begin{equation}
\tau_i < \frac{2p_i - \underline{p}}{\beta_i p_i + \underline p^{-1}p_i\sum_{j=1}^m \pi_j \sigma_j A_{j, i}^2}.\label{eq: ss_choice}
\end{equation}
The iterates $z_k$ are produced by~\Cref{alg:stripd}. 
Then, almost surely, there exist $z_\star = (x_\star, y_\star) \in \mathcal{Z}_\star$ such that $z_k \to z_\star$.
\end{theorem}

We analyze the step size rule~\eqref{eq: ss_choice} in Theorem~\ref{thm: thm1} and compare with existing efficient methods in dense and sparse settings.
\begin{remark}\label{rm: rm1}
~\\
$\triangleright$ Let $A$ be dense, with all its elements being nonzero, $p_i = 1/n$ and $f(\cdot)=0$, then the step size rule reduces to
\begin{align}
\tau_i < \frac{1}{n \sigma \| A_i \|^2},\notag
\end{align}
which is the step size rule of SPDHG~\cite{chambolle2018stochastic,alacaoglu2019convergence}, which is shown to be favorable in the dense setting. In contrast, step size rules of~\cite{fercoq2019coordinate,latafat2019new} are $n$ times worse due to duplication, in this case.\\[2mm]
$\triangleright$ Let $A$ be such that it contains one nonzero element per row, and we use $p_i = \frac{1}{n}$, which results in $\pi_j = \frac{1}{n}$. Then,
\begin{align}
\tau_i < \frac{1}{\beta_i + \sum_{j=1}^m \sigma_j A_{j, i}^2},\notag
\end{align}
which is the step size rule of Vu-Condat-CD~\cite{fercoq2019coordinate}, upon using the definition of $J(i)$ from~\eqref{eq: def_ji}.
Similarly, Algorithm~\ref{alg:stripd} updates $1$ dual coordinate and $1$ primal coordinate, in this case. In contrast, SPDHG~\cite{chambolle2018stochastic} updates $m$ dual coordinates, resulting in $m$ times higher per iteration cost.
\end{remark}

We note that the step size of TriPD-BC~\cite{latafat2019new} depends on global Lipschitz constant of $f$ rather than coordinatewise Lipschitz constants.
In fact, using coordinatewise Lipschitz constants is very important one of the most important reasons of the success of coordinate descent methods, as it results in larger step sizes~\cite{nesterov2012efficiency,richtarik2014iteration,fercoq2015accelerated}.

The takeaway from Remark~\ref{rm: rm1} is that Algorithm~\ref{alg:stripd} recovers the characteristics of the best performing methods in fully dense and fully sparse settings.
Moreover, as it is the only method that has the desirable dependencies in both cases, it has the best properties in the moderate sparse cases.
We also validate this observation in numerical experiments. 

\subsection{Linear convergence}
Linear convergence of primal-dual methods in practice is a widely observed phenomenon~\cite{chambolle2011first,liang2016convergence}.
We show that Algorithm~\ref{alg:stripd} also shares this property and obtains linear convergence under metric subregularity, without any modification on the algorithm.

We define the Bregman-type projection onto the solution set 
\begin{equation}\label{eq: def_zkstar}
z_k^\star = \arg\min_{u\in\mathcal{Z}_\star} D_p(x_k; u) + V(z_k - u).
\end{equation}
We now show that $z_k^\star$ is well-defined under our assumptions.
First, the solution set is convex and closed. 
Second, $D_p(x_{k}; u) \geq 0$ for all $u\in\mathcal{Z}_\star$ and it is also lower semicontinuous. 
Third, we remark that $V(z_k-u)$ is a squared norm (see Lemma~\ref{lem: lem1}), thus coercive, therefore the sum is coercive and lower semicontinuous over $\mathcal{Z}_\star$.
Hence, $z_k^\star$ exists.

The definition of $z_k^\star$ in~\eqref{eq: def_zkstar} is more involved compared to the corresponding quantity in~\cite{latafat2019new}.
This is in fact due to us using coordinatewise Lipschitz constants in our step sizes, rather than the global Lipschitz constant in~\cite{latafat2019new}.
\begin{assumption}\label{asmp: asmp2}~\\
KKT operator $F$ is metrically subregular at all $z_\star \in \mathcal{Z}_\star$ for $0$, and $\bar{z}_k\in\mathcal{N}(z_\star), \forall z_\star, \forall k$.
\end{assumption}

\begin{theorem}\label{thm: thm2}
Let~\Cref{asmp: asmp1} and~\ref{asmp: asmp2} hold.
Let $\theta$ and the step sizes $\tau, \sigma$ be chosen according to~\eqref{eq: theta_choice} and~\eqref{eq: ss_choice}.
Moreover, $z_k^\star = (x_k^\star, y_k^\star)$ is as defined in~\eqref{eq: def_zkstar}.
Then, for $z_k$ generated by Algorithm~\ref{alg:stripd}, it follows that
\begin{multline}
\mathbb{E} \left[\frac{\underline p}{2} \| x_k - x_k^\star \|^2_{\tau^{-1} P^{-1}} + \frac{\underline p}{2} \| y_k - y_k^\star \|^2_{\sigma^{-1}\pi^{-1}} \right] \\ 
\leq (1-\rho)^k \Delta_0,\notag
\end{multline}
where $\rho = \min\left( \underline p, \frac{C_{2, \tilde V}}{C_{V, 2} ((2+2c)+(1+c)(\eta \| H-M\|+\beta))^2} \right)$, $\Delta_0 = D_p(x_0; z_0^\star) + V(z_0 - z_0^\star)$, $\bar \beta$ is the global Lipschitz constant of $f$, \\$C_{2, \tilde V} = \frac{\underline{p}}{2} \min \left\{ \min_i C(\tau)_i, \min_j \sigma_j^{-1} \right\}$, \\
$C_{V, 2}=\frac{1}{2} \max\left\{ \max_i \frac{1}{\tau_i}, \max_j \frac{1}{\sigma_j} \right\}$, \\
$C_{2,V} = \sqrt{\frac{2}{\underline p\min_i\{\tau_i^{-1}p_i^{-1}\}}}$, $c=C_{2,V}\sqrt{\| A\|/2}$, and
\begin{equation}
H = \begin{bmatrix} \tau^{-1} & A^\top \\ 0 & \sigma^{-1} \end{bmatrix}, ~~~~ M = \begin{bmatrix} 0 & A^\top \\ -A & 0 \end{bmatrix}.\notag
\end{equation}
\end{theorem}


The first remark about Theorem~\ref{thm: thm2} is that since metric subregularity constant $\eta$ is not required in the algorithm, the step sizes to achieve linear convergence are the same step sizes as~\eqref{eq: ss_choice}.
Therefore, PURE-CD adapts to structures on the problem, without any need to modify the algorithm, and attains linear rate of convergence.
This supports the well-known observation that primal-dual algorithms converge linearly on most problems, with standard step sizes in~\eqref{eq: ss_choice}.

In particular, a direct corollary of our theorem is that for problems listed in Example~\ref{ex: ms_examples}, PURE-CD obtains linear rate of convergence.
For the first case in Example~\ref{ex: ms_examples}, our result applies directly since the neighborhood of subregularity $\mathcal{N}(z_\star)$ is the whole space.
For the second case, we have to assume additionally that $\bar{z}_k$ is contained in a compact set, since the $\mathcal{N}(z_\star)$ is not the whole space, and is bounded.
A sufficient assumption for this is when the domains of $g$ and $h^\ast$ are compact.
We note that compactness is only required for this result in our paper. 
This is common to other results for PDCD methods with metric subregularity~\cite{latafat2019new,alacaoglu2019convergence}.
The issue, as explained in~\cite{alacaoglu2019convergence}, stems from a fundamental limitation of the existing analyses of PDCD methods.

Many results in the literature for linear convergence only applies to the first case in Example~\ref{ex: ms_examples}, when $g, h^\ast$ are strongly convex~\cite{zhang2017stochastic,chambolle2018stochastic}.
Moreover, these results require setting step sizes depending on strong convexity constants of $g, h^\ast$, therefore not applicable when strong convexity is absent.
Our result applies to more general problems and it uses step sizes independent of these constants.
Our algorithm can be directly applied to any problem satisfying Assumption~\ref{asmp: asmp1} and fast convergence will occur provably, if the selected problem is in Example~\ref{ex: ms_examples}.

Compared with the linear convergence rate in~\cite{latafat2019new} for TriPD-BC, our result have a similar contraction factor, however, due to larger step sizes (see~\Cref{rm: rm1}), the rate comes with a better constant.


\subsection{Ergodic rates}\label{sec: ergodic}
In this section, we study Algorithm~\ref{alg:stripd} in the general case, under Assumption~\ref{asmp: asmp1}, and show the optimal $\mathcal{O}(1/k)$ convergence rate on the ergodic sequence.
The quantity of interest is the primal-dual gap function~\cite{chambolle2011first}
\begin{align}
G(\bar{x}, \bar{y}) &= \sup_{z=(x, y)\in\mathcal{Z}} f(\bar{x})+g(\bar{x}) + \langle A\bar{x}, y \rangle - h^\ast(y) \notag \\
&-f(x) - g(x) - \langle Ax, \bar{y} \rangle + h^\ast(\bar{y}).\label{eq: def_pdgap}
\end{align}
A related quantity is the restricted gap function (see~\cite{chambolle2011first}) for any set $\mathcal{C}\subset\mathcal{Z}$
\begin{align}
G_\mathcal{C}(\bar{x}, \bar{y}) &= \sup_{z\in\mathcal{C}} f(\bar{x})+g(\bar{x}) + \langle A\bar{x}, y \rangle - h^\ast(y) \notag \\
&-f(x) - g(x) - \langle Ax, \bar{y} \rangle + h^\ast(\bar{y}).
\end{align}

Due to randomization in PDCD, we are interested in the expected primal-dual gap, denoted as $\mathbb{E}\left[ G_{\mathcal{C}}(\bar{x}, \bar{y}) \right]$.
As noted by~\citet{dang2014randomized}, it is technically challenging to prove rates for this quantity as it is the expectation of supremum.
Recently,~\cite{alacaoglu2019convergence} used a technique to show convergence of expected primal-dual gap for SPDHG of~\cite{chambolle2018stochastic}.
This rate is for ergodic sequence averaging $x_k$ and the full dual variable $\bar{y}_k$.
We can use this technique for our analysis. However, there remains another technical challenge as full dual variable is not computed in PURE-CD.
Thus, averaging $\bar{y}_k$ is not feasible in our case.

In addition to~\Cref{asmp: asmp1}, in this section we will assume separability of $h$, to be able to do an efficient averaging with the dual iterate.

Due to the asymmetric nature of Algorithm~\ref{alg:stripd}, there are fundamental difficulties for proving a rate with averaging $y_{k+1}$.
On this front, we propose a new type of analysis for the dual variable.
To start with, we define the following iterate which has the same cost to compute as $y_{k+1}$ each iteration.
Let $\breve{y}_1 = y_1 = \bar{y}_1$,
\begin{equation}
\begin{aligned}\label{eq: breve_def}
\breve{y}_{k+1}^j &= \bar{y}_{k+1}^j, &&\forall j \in J(i_{k+1}),\\
\breve{y}_{k+1}^j &= \breve{y}_k^j, &&\forall j\not\in J(i_{k+1}).
\end{aligned}
\end{equation}
We note that $\breve{y}_k$ is $\mathcal{F}_k$-measurable and more useful properties of $\breve{y}_k$ for analysis are discussed in~\Cref{lem: breve} in the appendix.

Due to the definition of $\breve{y}_k$, it is now feasible to compute and average this iterate.
We can show the convergence of expected primal-dual gap by averaging $\breve{y}_k$ and $x_k$.
We remark that we use some coarse inequalities to give simple constants for~\Cref{th:erg} and~\Cref{cor: erg} in this section, which results in suboptimal dependence with respect to dimension $n$.
In~\Cref{sec: proofs}, we give these theorems with their original, tighter bounds and we show how we transform the tighter bounds into the constants we give in this section.
\begin{theorem}\label{th:erg}
Let Assumption~\ref{asmp: asmp1} hold, $\theta, \tau, \sigma$ are chosen as in~\eqref{eq: theta_choice},~\eqref{eq: ss_choice}, and $h$ be separable.
Let $x^{av}_K = \frac{1}{K} \sum_{k=1}^K x_{k}$ and $y^{av}_K = \frac{1}{K} \sum_{k=1}^K \breve{y}_{k}$, where $\breve{y}_k$ is defined in~\eqref{eq: breve_def}. Then it holds that for any bounded set $\mathcal{C}=\mathcal{C}_x\times\mathcal{C}_y\subset\mathcal{Z}$
\begin{align}
\mathbb{E} \left[G_\mathcal{C}(x^{av}_K, y^{av}_K) \right] \leq \frac{C_{g}}{\underline pK},\notag
\end{align}
where $C_g = \sum_{i=1}^4 C_{g, i}$, $C_{\tau, \tilde V} =  \min_i C(\tau)_i\tau_i$, \\
$C_{g, 1} = \sup_{z\in \mathcal{C}} \big\{2\underline p \| x_0 - x\|^2_{\tau^{-1}P^{-1}} + 2\underline p\| y_0 - y \|^2_{\sigma^{-1}\pi^{-1}} \big\} 
+4\sqrt{\Delta_0\underline p^{-1}}\|A\| \sup_{y\in \mathcal{C}_y}\|y\|_{\tau P} \\+ \sqrt{\Delta_0(\underline p^{-1} - 2\underline p^{-2}C_{\tau,\tilde V}^{-1})}\|A\|\sup_{x\in\mathcal{C}_x}\|x\|_{\sigma\pi}$,\\
$\sum_{i=2}^4 C_{g, i} = \Delta_0\left( 5+9\underline p^{-1} + C_{\tau,\tilde V}^{-1}\left(1+10\underline p^{-1} + 8\underline p^{-2}\right) \right)  \\
+(1-\underline p)(f(x_0) + g(x_0) - f(x_\star) - g(x_\star)) + h^\ast(y_0) - h^\ast(y_\star)  \\
+\underline p\| Ax_\star \|_{\sigma\pi^{-1}} + \|A^\top y_\star \|^2_{\tau P}$.
\end{theorem}
\begin{remark}
When implementing averaging of $x_k$, and $\breve{y}_k$, one should use a technique similar to~\cite{dang2015stochastic}.
The main idea is to only update the averaged vector at the coordinates where an update occurred.
For this, one needs to remember for each coordinate, the last time it is updated, and update the averaged vector using this information, when a coordinate is selected.
\end{remark}
The result in~\Cref{th:erg} would give a rate for primal-dual gap when $\mathcal{C}=\mathcal{Z}$. 
However, in general such a rate is not desirable as taking a supremum over $\mathcal{Z}$ might result in an infinite bound.
This rate would be meaningful when both primal and dual domains are bounded in which case one would take the supremum in $C_{g,1}$ over the bounded domains. 

Alternatively, in the following theorem, we show that for two important special cases, we can extend this result to show guarantees without bounded domains.
Namely, we show the same rate for the case when $h(\cdot) = \delta_{\{b\}}(\cdot)$, $b\in\mathbb{R}^m$ to cover linearly constrained problems.
Moreover, we show the result for the case when $h$ is Lipschitz continuous.

\begin{theorem}\label{cor: erg} 
Let Assumption~\ref{asmp: asmp1} hold. We use the same parameters $\theta,\tau,\sigma$ and the definitions for $x_K^{av}$ and $y_K^{av}$ as Theorem~\ref{th:erg}. 
We consider two cases separately: \\
$\triangleright$ If $h(\cdot) = \delta_{\{ b \}}(\cdot)$, we obtain
\begin{align}
&\mathbb{E} [ f(x_K^{av}) + g(x_K^{av}) - f(x_\star) - g(x_\star)] \leq \frac{C_o}{\underline pK}. \notag\\
&\mathbb{E}[ \|Ax_K^{av}-b \|] \leq \frac{C_f}{\underline pK}.\notag
\end{align}
$\triangleright$ If $h$ is $L_h$-Lipschitz continuous, we obtain
\begin{multline}
\mathbb{E}[ f(x_K^{av}) + g(x_K^{av}) + h(Ax_K^{av}) \\
- f(x_\star) - g(x_\star) -h(Ax_\star) ] \leq \frac{C_l}{\underline pK},\notag
\end{multline}
where 
$C_f = 2c_1\|y_\star - y_0 \|_{\sigma^{-1}\pi^{-1}} + 2\sqrt{c_1C_s} + 2\sqrt{2}c_1$, \\
$C_o = C_s + \| y_\star \|_{\sigma^{-1}\pi^{-1}} C_f + 2c_1\underline p^{-1}V(z_0-z_\star)$,\\
$C_l = C_s + c_1 \| x_\star-x_0\|^2_{\tau^{-1}P^{-1}} + 4 c_1 L_h^2$,\\
$C_s =  C_{g, 2}+ C_{g, 5}+C_{g,6}$, with $c_1 = 2\underline p + 2$, $C_{g, 2}$ as defined in Theorem~\ref{th:erg} and $C_{g,5}, C_{g, 6}$ are defined in the proof in~\eqref{eq: cg3_def},~\eqref{eq: cg4_def}.
\end{theorem}


\section{Related works}\label{sec: lit}
In the deterministic setting, many primal-dual methods are proposed~\cite{chambolle2011first,vu2013splitting,condat2013primal,tan2020accelerated,latafat2019new}.
The standard results in these papers include linear convergence when $g, h^\ast$ are strongly convex, with step sizes selected by using strong convexity constants.
In addition, these papers also show sublinear $\mathcal{O}(1/k)$ rate with general convexity, which is known to be optimal~\cite{nesterov2005smooth}.
Moreover, linear rates under metric subregularity is shown for deterministic methods in~\cite{liang2016convergence,latafat2019new}.

Randomized coordinate descent is proposed in~\cite{nesterov2012efficiency} and improved by a large body of subsequent papers~\cite{richtarik2014iteration,fercoq2015accelerated}.
Primal randomized coordinate descent requires full separability on the nonsmooth parts of the objective function.
Nonsmooth and nonseparable functions are handled by primal-dual coordinate descent methods~\cite{fercoq2019coordinate}.

One of the first primal-dual coordinate descent (PDCD) methods is SPDC, which is proposed in~\cite{zhang2017stochastic}, that solves a special case of problem~\eqref{eq: prob_temp} with $f=0$.
SPDC has linear convergence when $g, h^\ast$ are strongly convex and the step sizes are selected according to strong convexity constants.
In the general convex case, SPDC has perturbation-based analysis, which needs to set an $\epsilon$, requires knowing $\| x_\star\|^2$, and shows $\epsilon$-based iteration complexity results, and not anytime convergence rates.
Almost sure convergence of the iterates of SPDC is not proven in the general convex case.
Moreover, the step sizes of SPDC are scalar and they depend on the maximum block norm of $A$.
It is shown in~\cite{zhang2017stochastic} that in the specific cases when $g(x)=\| x \|^2$ or $g(x)=\| x\|_1 + \| x \|^2$, one can use a special implementation for efficiency with sparse data.

\citet{tan2020accelerated} proposed a new method similar to SPDC with the same type of guarantees as~\cite{zhang2017stochastic}.
Due to similar analysis techniques, this method inherits the abovementioned drawbacks of SPDC.
For this method,~\citet{tan2020accelerated} showed a new implementation technique for sparse data, that can be used with any separable $g(x)$.

For solving the specific case of empirical risk minimization problems, stochastic dual coordinate ascent (SDCA) is proposed in~\cite{shalev2013stochastic,shalev2014accelerated}.
SDCA uses strong convexity constant to set step sizes and attain linear convergence.
A limitation of SDCA is to require strong convexity in the primal, to ensure smoothness of the dual objective, which is essential in the design of the method.

Another early PDCD method is by~\cite{dang2014randomized} where the authors focused on showing sublinear convergence rates.
The authors showed guarantees for a relaxed version of expected primal-dual gap function in~\eqref{eq: def_pdgap}.

Building on~\cite{dang2014randomized}, block-coordinate variants of alternating direction method of multipliers (ADMM) are proposed in~\cite{gao2019randomized,xu2018accelerated}.
These papers focus on linearly constrained problems and show ergodic sublinear convergence rates.
Moreover,~\cite{xu2018accelerated} showed that under strong convexity assumption and special decomposition of the blocks, the method achieves linear convergence.
This linear convergence result, similar to~\cite{zhang2017stochastic} requires knowing the strong convexity constants to set the algorithmic parameters.
Moreover, these results generally set step sizes depending on global Lipschitz constants and norm of $A$ rather than the norm of the blocks of $A$.

\begin{figure*}[ht]
\begin{center}
\includegraphics[width=0.6\columnwidth]{rcv1_tex}
\includegraphics[width=0.6\columnwidth]{w8a_tex}
\includegraphics[width=0.6\columnwidth]{cov_tex}
\caption{Lasso: Left: rcv1, $n=20,242, m=47,236$, density $=0.16 \%$, $\lambda=10$; Middle: w8a, $n=49,749, m=300$, density $=3.9\%$, $\lambda=10^{-1}$; Right: covtype, $n=581,012$, $m=54$, density $=22.1\%$, $\lambda=10$.}
\label{fig:1}
\end{center}
\vskip -0.2in
\end{figure*}
\begin{figure*}[ht]
\begin{center}
\includegraphics[width=0.6\columnwidth]{sec_tex}
\includegraphics[width=0.6\columnwidth]{a9a_tex}
\includegraphics[width=0.6\columnwidth]{mnist_tex}
\caption{Ridge regression: Left: sector, $n=6,412, m=55,197$, density $=0.3 \%$, $\lambda=0.1$; Middle: a9a, $n=32,561, m=123$, density $=11.3\%$, $\lambda=0.1$; Right: mnist, $n=60,000$, $m=780$, density $=19.2\%$, $\lambda=1$.}
\label{fig:2}
\end{center}
\vskip -0.2in
\end{figure*}



PDCD variants are also proposed in~\cite{combettes2015stochastic,combettes2019stochastic,pesquet2015class} and analyzed under the general setting of monotone operators. These methods use global constants of the problem such as global Lipschitz constant of smooth part and $\|A\|$, rather than blockwise constants, resulting in worse practical performance, as illustrated in the experiments of~\cite{chambolle2018stochastic}.

Another early PDCD variant to solve problem~\eqref{eq: prob_temp} in its full generality, where $f, g, h$ are all nonseparable, is by~\citet{fercoq2019coordinate}.
This method uses coordinatewise Lipschitz constants of the smooth part and it is designed to exploit sparsity of $A$.
This method has almost sure convergence guarantees as well as linear convergence when $g, h^\ast$ are strongly convex.
As opposed to most results in this nature, it is not required to know strong convexity constants to set the step sizes.
In the general convex case, the method has $\mathcal{O}(1/\sqrt{k})$ rate for a randomly selected iterate.
As argued in Section~\ref{sec: as}, main limitation of~\cite{fercoq2019coordinate} is that small step sizes are required when matrix $A$ is dense.
Moreover, the results in this paper are restricted to uniform probability law for selecting coordinates.

One of the most related works to ours, and a building block of PURE-CD is TriPD-BC from~\cite{latafat2019new}.
The authors showed almost sure convergence of the iterates and linear convergence under metric subregularity, by using global Lipschitz constants of $f$ for the step sizes. This paper did not have any sublinear convergence rates in the general convex case.
Similar to~\cite{fercoq2019coordinate}, TriPD-BC is designed for sparse setting and a naive implementation in the dense setting requires the same per iteration cost as the deterministic algorithm.
An efficient implementation is by duplication of dual variables, which as explained in~\cite{fercoq2019coordinate} results in small step sizes.

Another building block of PURE-CD is SPDHG by~\cite{chambolle2018stochastic}, to solve~\eqref{eq: prob_temp} when $f=0$.
Linear convergence result of SPDHG by~\cite{chambolle2018stochastic} is similar to~\cite{zhang2017stochastic} and requires setting step sizes with strong convexity constants.
In the general convex case and partially strongly convex case,~\cite{chambolle2018stochastic} proved optimal sublinear rates.
Recently,~\cite{alacaoglu2019convergence} analyzed SPDHG and proved additional theoretical results.
In particular, this work showed almost sure convergence of the iterates of SPDHG and linear convergence under metric subregularity.
Even though it is fast in the dense setting, the main limitation of SPDHG, as discussed in Section~\ref{sec: as} is that it needs to update all the dual coordinates, resulting in high per iterations costs in the sparse setting.

Similar algorithms are proposed in~\cite{luke2018block,alacaoglu2017smooth} where the authors focused on the linearly constrained problems and proved sublinear rates.





\section{Numerical experiments}\label{sec: numexp}
\subsection{Effect of sparsity}\label{sec:sparse}
As explained in Section~\ref{sec: as}, and Remark~\ref{rm: rm1}, PURE-CD brings together the benefits of different methods that are designed for dense and sparse cases.
We will now compare the empirical performance of PURE-CD with Vu-Condat-CD from~\cite{fercoq2019coordinate} which has desirable properties with sparse data and SPDHG from~\cite{chambolle2018stochastic} which has desirable properties with dense data.

We select uniform sampling, $p_i = 1/n$, so~\eqref{eq: ss_choice} simplifies to
\begin{align}\label{eq: ss_rule_exp}
\tau_i < \frac{1}{\sum_{j=1}^m \theta_j \sigma_j A_{j, i}^2}.
\end{align}
We provide a step size policy inspired by the step size rules chosen in~\cite{chambolle2018stochastic} and~\cite{fercoq2019coordinate}.
We use the following step sizes, for $\gamma < 1$
\begin{align}
\sigma_j = \frac{1}{\theta_j \max_{i'}\|A_{i'}\|}, ~~~~ \tau_i = \frac{\gamma \max_{i'} \|A_{i'} \|}{\| A_i \|^2}.\notag
\end{align}
We note that in contrast to~\cite{chambolle2018stochastic}, step sizes are both diagonal. In our case, it is important to utilize diagonal step sizes for both primal and dual variables since we perform coordinate-wise updates for both primal and dual variables and the step sizes need to be set appropriately to obtain good practical performance.
For SPDHG and Vu-Condat-CD, we use step sizes suggested in the papers.

In the edge cases (one nonzero element per row or fully dense), it is easy to see that our step size policy reduces to the suggested step sizes of~\cite{chambolle2018stochastic} and~\cite{fercoq2019coordinate}.

For experiments, we used the generic coordinate descent solver, implemented in Cython, by~\citet{fercoq2019generic}, which includes an implementation of Vu-Condat-CD with duplication and we implemented SPDHG and PURE-CD.
We solve Lasso and ridge regression, where we let
$g(x) = \lambda \|x\|_1$, $h(Ax) = \frac{1}{2} \| Ax-b\|^2$, $f=0$, and $g(x) = \frac{\lambda}{2} \|x\|^2$, $h(Ax) = \frac{1}{2} \| Ax-b\|^2$, $f=0$, respectively, in our template~\eqref{eq: prob_temp}. Then, we apply all the methods to the dual problems of these, to access data by rows.

\begin{figure*}[ht]
\begin{center}
\includegraphics[width=1.8\columnwidth]{a9a_ridge_ml.pdf}
\end{center}
\vskip -0.2in
\end{figure*}
\begin{figure*}[ht]
\begin{center}
\includegraphics[width=1.8\columnwidth]{sec_ridge_ml.pdf}
\caption{top: a9a, $n=32,561$, $m=123$, bottom: sector, $n=6,412$, $m=55,197$.}
\label{fig:3}
\end{center}
\vskip -0.2in
\end{figure*}

We use datasets from LIBSVM with different sparsity levels~\cite{CC01a}. The properties of each data matrix are given in the caption of the corresponding figures.
For preprocessing, we removed all-zero rows and all-zero columns of $A$ and we performed row normalization.
The results are compiled in~\Cref{fig:1,fig:2}.

We observe the behavior predicted by theory. With sparse data such as rcv1, where density level is $0.16\%$, SPDHG makes very little progress in the given time window. The reason is that the per iteration cost of SPDHG in this case is updating $47,236$ dual variables, whereas for PURE-CD and Vu-Condat-CD, the cost is updating $75$ dual variables.
We note that PURE-CD is faster than Vu-Condat-CD due to better step sizes.
On the other hand, with moderate sparsity, SPDHG and Vu-Condat-CD is comparable, whereas PURE-CD exhibits the best performance.
For denser data, SPDHG and PURE-CD exhibit similar behavior where Vu-Condat-CD is slower than both due to smaller step sizes.

\subsection{Comparison with specialized methods}
In this section, we compare the practical performance of PURE-CD with state-of-the-art algorithms that are designed for strongly convex-strongly concave problems.
Due to space constraints, we defer some of the plots and more details about experiments to the appendix.
We focus on the problem
$
\min_x \frac{1}{n} \sum_{i=1}^n h_i(A_i x) + \frac{\lambda}{2} \| x \|^2,
$
where $h_i(x) = (x - b_i)^2$.
Each $h_i$ is smooth with Lipschitz constants $L_i = 2$ and the second component is strongly convex.
This is equivalent to strong convexity in both primal and dual problems.

In this case, the algorithms SDCA~\cite{shalev2013stochastic}, ProxSVRG~\cite{xiao2014proximal}, Accelerated SVRG~\cite{zhou2018simple}, SPDC~\cite{zhang2017stochastic} are all designed to use the strong convexity to obtain linear convergence.
These algorithms use the strong convexity constant $\lambda$ for setting the algorithmic parameters.
Moreover, as all the abovementioned algorithms have special implementations to exploit sparsity, we make the comparison with respect to number of passes of the data, rather than time.
The results are compiled for two datasets in~\Cref{fig:3} and more datasets are included in~\Cref{app: exp}.

$\bullet$ PURE-CD-$\lambda$: This variant uses the non-agnostic step sizes, using $\lambda$, which still satisfy the theoretical requirement~\eqref{eq: ss_rule_exp}.
\begin{align}
\sigma_j = \frac{n}{\theta_j \sqrt{n\lambda} \max_{i'}\|A_{i'}\|}, ~~~~ \tau_i = \frac{\gamma\sqrt{n\lambda} \max_{i'} \|A_{i'} \|}{n\| A_i \|^2}.\notag
\end{align}
$\bullet$ PURE-CD: This variant is with the standard agnostic step sizes. We note that the step sizes are scaled by $n$ since the problem is scaled by $1/n$, compared to~\Cref{sec:sparse}.
\begin{align}
\sigma_j = \frac{n}{\theta_j \max_{i'}\|A_{i'}\|}, ~~~~ \tau_i = \frac{\gamma \max_{i'} \|A_{i'} \|}{n \| A_i \|^2}.\notag
\end{align}
We observe that PURE-CD has a consistent linear convergence behavior as predicted by theory.
In most of the datasets (see~\Cref{app: exp}), it has the fastest convergence behavior.
However, in some datasets, as $\lambda$ gets smaller, we observed that the linear rate of PURE-CD slowed down, which motivated us to try PURE-CD-$\lambda$, which incorporates the knowledge of $\lambda$ as the other methods.
It seems to show favorable behavior when PURE-CD slows down.

The takeaway message is that PURE-CD, which is designed for a general problem, adapts to strong convexity well with agnostic step sizes in most cases.
However, in some cases, it does not perform as good as the algorithms which are designed to exploit strong convexity.
In those cases however one can choose separating step sizes of PURE-CD accordingly, and use PURE-CD-$\lambda$ to get better performance.


\section*{Acknowledgements}
This project has received funding from the European Research Council (ERC) under the European Union's Horizon $2020$ research and innovation programme (grant agreement no $725594$ - time-data) and the Swiss National Science Foundation (SNSF) under grant number $200021\_178865 / 1$.

