###### Abstract

Many applications using large datasets require efficient methods for
minimizing a proximable convex function subject to satisfying a set of linear
constraints within a specified tolerance. For this task, we present a proximal
projection (PP) algorithm, which is an instance of Douglas-Rachford splitting
that directly uses projections onto the set of constraints. Formal guarantees
are presented to prove convergence of PP estimates to optimizers. Unlike many
methods that obtain feasibility asymptotically, each PP iterate is feasible.
Numerically, we show PP either matches or outperforms alternatives (e.g.
linearized Bregman, primal dual hybrid gradient, proximal augmented
Lagrangian, proximal gradient) on problems in basis pursuit, stable matrix
completion, stable principal component pursuit, and the computation of earth
mover’s distances.

Key words: |  proximal, projection, convex optimization, Douglas-Rachford splitting, basis pursuit, compressed sensing, matrix completion, principal component pursuit, earth mover’s distance, Wasserstein distance
---|---

##  1 Introduction

With the profound increase in the use of extremely high-dimensional data, an
ongoing challenge is to create efficient tools for processing this data. In
several important applications, this processing takes the form of solving a
convex optimization problem. Practical value is, thus, derived from efficient
and easy-to-code methods for such tasks, which is the focus of this work.
Specifically, for a convex function111We set
$\overline{\mathbb{R}}\triangleq\mathbb{R}\cup\\{+\infty\\}$.
$f\colon\mathbb{R}^{n}\rightarrow\overline{\mathbb{R}}$, a matrix
$A\in\mathbb{R}^{m\times n}$, a vector $b\in\mathbb{R}^{m}$ and a scalar
$\varepsilon\geq 0$, this work considers the problem

$\min_{x\in\mathbb{R}^{n}}f(x)\quad\mbox{s.t.}\quad\|Ax-b\|\leq\varepsilon.$
(P)

We let $\mathcal{C}\triangleq\\{x:\|Ax-b\|\leq\varepsilon\\}$ denote the
constraint set so that (P) equates to minimizing $f$ over $\mathcal{C}$. As
many applications introduce a tolerance $\varepsilon>0$ on linear constraints
to ensure stability with respect data being corrupted by noise, herein we
refer to (P) as a stable linearly constrained optimization problem.

The problem (P) is well-studied. When $\varepsilon=0$, several algorithms can
be directly applied (e.g. primal dual hybrid gradient [1], conditional
gradient [2]). When $\varepsilon>0$, an auxiliary variable can be introduced
to decompose $\mathcal{C}$ into a linear constraint and a ball constraint; for
example, the alternating direction method of multipliers [3, 4] can be readily
applied to that formulation. In some works, the constraint is moved into the
objective as a quadratic penalty [5, 6]; for an appropriate penalty weight,
this soft-penalty variation shares the same minimizers as (P). The
superiorization methodology [7, 8] may also approximate solutions to (P) by
interweaving projection steps onto $\mathcal{C}$ (or sets whose intersection
forms $\mathcal{C}$) with subgradient steps. Other approaches [9, 10] use
various forms of smoothing or added regularization to solve problems that
approximate (P).

Much prior work aims to solve (P) by either approximating (P) or obtaining
feasibility asymptotically. In contrast, we solve (P) and maintain feasibility
at each step of our iterative algorithm. This is done at comparable per-
iteration cost to existing methods by interweaving projections onto
$\mathcal{C}$ and proximal operations (defined below).

#### Contribution.

Our main result is Algorithm 2 and its convergence to solutions of (P), which
is possible due to a novel formula for the projection onto the constraint set
$\mathcal{C}$. Our numerical examples show favorable performance of PP against
alternatives in a varied collection of practical problems, i.e. basis pursuit,
stable principal component pursuit, computation of earth mover’s distances,
and stable matrix completion.

#### Notation.

Here, $\|\cdot\|$ is the Euclidean norm, $\|\cdot\|_{F}$ is the Frobenius
norm, $\|\cdot\|_{1}$ is the $1$-norm, and $\|\cdot\|_{\star}$ is the nuclear
norm. The relative interior domain of $f$ is $\mbox{ri}(\mbox{dom}(f))$. For
an integer $n\in\mathbb{N}$, we set $[n]\triangleq\\{1,\ldots,n\\}$.

Algorithm 1 Proximal Projection (PP) for (P) with $\varepsilon=0$

1: PP($f$, $A$, $b$):

2: initialize iterate $z\in\mathbb{R}^{n}$ and parameter $\alpha>0$

3: while stopping criteria not met

4: $x\leftarrow z-A^{\top}(AA^{\top})^{-1}(Az-b)$

5: $z\leftarrow z+\mbox{prox}_{\alpha f}\left(2x-z\right)-x$

6: return $x$

##  2 Main Results

In this section, we define our proposed proximal projection (PP) algorithm,
and analytically show it generates a solution to (P). Throughout, we make use
of the following conditions:

  1. (C1)

the function $f\colon\mathbb{R}^{n}\rightarrow\overline{\mathbb{R}}$ is
closed, convex, and proper;

  2. (C2)

either the matrix $A$ has full row-rank or $\varepsilon>0$;

  3. (C3)

$b$ is in the image of $A$;

  4. (C4)

either $f$ is coercive or $\mathcal{C}$ is bounded;

  5. (C5)

there is $x\in\mbox{ri}(\mbox{dom}(f))$ such that $Ax=b$ if $\varepsilon=0$
and $\|Ax-b\|<\varepsilon$ if $\varepsilon>0$.

To minimize the function $f$, we make use of a proximal operator. For a
parameter $\alpha>0$, this is defined by

$\mbox{prox}_{\alpha
f}(x)\triangleq\operatorname*{\mbox{argmin}}_{u\in\mathbb{R}^{n}}\alpha
f(u)+\frac{1}{2}\|u-x\|^{2},$ (1)

and (C1) ensures it uniquely exists [11, 12, 13]. In many instances, explicit
formulas exist for the proximal (e.g. see [11, Chapter 6]). Proximals also
generalize projections. Specifically, if $\mathcal{C}$ is a nonempty, closed
and convex set and $\delta_{\mathcal{C}}$ is the indicator function taking
value $0$ in $\mathcal{C}$ and $\infty$ elsewhere, then

$\mbox{prox}_{\alpha\delta_{\mathcal{C}}}(x)=\mathcal{P}_{\mathcal{C}}(x)\triangleq\operatorname*{\mbox{argmin}}_{u\in\mathcal{C}}\frac{1}{2}\|u-x\|^{2}.$
(2)

That is, the proximal for $\delta_{\mathcal{C}}$ is precisely the projection
$\mathcal{P}_{\mathcal{C}}$ onto $\mathcal{C}$. This leads to our next
hypotheses (C2) and (C3), which are used to obtain our projection formula in
the following lemma. (See Appendix A for a proof.)

Algorithm 2 Proximal Projection (PP) for (P)

1: PP($f$, $A$, $b$, $\varepsilon$):

2: initialize iterate $z\in\mathbb{R}^{n}$ and parameter $\alpha>0$

3: while stopping criteria not met

4: if $\|Az-b\|\leq\varepsilon$

5: $x\leftarrow z$

6: else

7:
$\tau\leftarrow\mbox{solution}(1=\tau\|(AA^{\top}+\varepsilon\tau\mbox{I})^{-1}(Az-b)\|)$

8: $x\leftarrow z-A^{\top}(AA^{\top}+\varepsilon\tau\mbox{I})^{-1}(Az-b)$

9: $z\leftarrow z+\mbox{prox}_{\alpha f}\left(2x-z\right)-x$

10: return $x$

######  Proposition 1 (Projection Formula).

If conditions (C2) and (C3) hold and
$\mathcal{C}\triangleq\\{x:\|Ax-b\|\leq\varepsilon\\}$, then

$\mathcal{P}_{\mathcal{C}}(x)=\begin{cases}\begin{array}[]{cl}x&\mbox{if}\
\|Ax-b\|\leq\varepsilon,\\\
x-A^{\top}(AA^{\top}+\varepsilon\tau_{x}\mbox{I})^{-1}(Ax-b)&\mbox{otherwise},\end{array}\end{cases}$
(3)

where, if $\|Ax-b\|>\varepsilon$, the scalar $\tau_{x}$ is the unique positive
solution to

$1=\tau\|(AA^{\top}+\varepsilon\tau\mbox{I})^{-1}(Ax-b)\|.$ (4)

The final conditions (C4) and (C5) ensure a solution exists and total duality
holds, a key condition required by many operator splitting methods to
establish convergence [13]. To describe our method, let
$z^{1}\in\mathbb{R}^{n}$ and $\alpha>0$. We construct sequences $\\{x^{k}\\}$
and $\\{z^{k}\\}$ with the update at each index $k$ given by the formulas

$\displaystyle x^{k}$
$\displaystyle=\begin{cases}\begin{array}[]{cl}z^{k}&\mbox{if}\
\|Az^{k}-b\|\leq\varepsilon,\\\
z^{k}-A^{\top}(AA^{\top}+\varepsilon\tau_{z^{k}}\mbox{I})^{-1}(Az^{k}-b)&\mbox{otherwise,}\end{array}\end{cases}$
(5a) $\displaystyle z^{k+1}$ $\displaystyle=z^{k}+\mbox{prox}_{\alpha
f}(2x^{k}-z^{k})-x^{k},$ (5b)

with $\tau_{z^{k}}$ defined as in Proposition 1. The iteration in (5) is a
special case of a more general scheme known as Douglas-Rachford splitting
(DRS) [14, 15], which has many uses (e.g. finding the zero of a sum of
monotone operators [14, 16], feasibility problems [17], combinatorial
optimization [18]). Making use of prior DRS results, the following theorem
justifies use of Algorithm 1 (the $\varepsilon=0$ case) and Algorithm 2 (the
$\varepsilon\geq 0$ case).

######  Theorem 1 (Convergence of PP).

If conditions (C1) to (C5) hold, then the sequences $\\{x^{k}\\}$ and
$\\{z^{k}\\}$ generated by (5) converge, with $\\{x^{k}\\}$ converging to a
solution of (P). Moreover, $\|Ax^{k}-b\|\leq\varepsilon$ for all $k$.

A proof of Theorem 1 is provided in Appendix A. The rest of this section
considers per-iteration costs of PP.

######  Remark 1 (Computation of $\tau$).

Explicit formulas can sometimes be derived for $\tau_{x}$ in Line 7 of
Algorithm 2 (e.g. see Subsections 3.2 and 3.4). Otherwise, $\tau_{x}$ may be
computed via a 1D solver (e.g. bisection method). Note
$0\leq\tau_{x}\leq\sigma_{max}(A)^{2}/(\|Ax-b\|-\varepsilon)$, with
$\sigma_{max}(A)$ the largest singular value of $A$ (see Lemma A.1).

The update formulas for PP and standard alternatives involve a proximal
operation and matrix multiplications. The costs associated with multiplication
by $A$ are $\mathcal{O}(mn^{2})$ and for $A^{\top}$ they are
$\mathcal{O}(m^{2}n)$. Unlike most schemes we compare to, PP includes
multiplication by
$(AA^{\top}+\varepsilon\tau_{x}\mbox{I})^{-1}\in\mathbb{R}^{m\times m}$. In
the worst case, this adds $\mathcal{O}(m^{3})$ cost. Computing the matrix
$(AA^{\top}+\varepsilon\tau_{x}\mbox{I})^{-1}$ can also add
$\mathcal{O}(m^{3})$ cost. However, with certain structures of $A$, both of
these costs can be reduced to $\mathcal{O}(m)$, yielding little impact. The
numerical examples below investigate whether the combination of per-iteration
cost and convergence rate of PP is more efficient than alternatives.

######  Remark 2 (Inversion with $\varepsilon=0$).

In the case that $\varepsilon=0$, there is a one-time computation of
$(AA^{\top})^{-1}$ that has cost $\mathcal{O}(m^{3})$. Ammortizing this over
hundreds of iterations yields negligible per-iteration cost. Moreover, in some
applications the same matrix $A$ is repeatedly used with new measurement data
$b$, in which case $(AA^{\top})^{-1}$ may be computed in an offline setting.

######  Remark 3 (SVD Inverse).

If $USV^{\top}$ is the singular value decomposition (SVD) of $A$, with
$\Sigma=\mbox{diag}(\sigma_{i})$, then

$A^{\top}(AA^{\top}+\varepsilon\tau\mbox{I})^{-1}=V\,\mbox{diag}\left(\frac{\sigma_{i}}{\sigma_{i}^{2}+\varepsilon\tau}\right)U^{\top}.$
(6)

Thus, if the SVD of $A$ is available, then the matrix inverse in Line 8 of
Algorithm 2 is readily available too.

######  Remark 4 (Inversion of Tridiagonal Matrices).

When the matrix $A$ has tridiagonal structure, Thomas’ algorithm [19] can be
used to multiply $(AA^{\top})^{-1}$ with $\mathcal{O}(m)$ cost rather than the
$\mathcal{O}(m^{3})$ cost via Gaussian elimination.

##  3 Numerical Examples

We provide four numerical examples. In each setting, conditions (C1) to (C5)
hold and, when applicable, a simplification of Algorithm 2 and Proposition 1
are presented. Shown methods may be accelerated (including PP) to yield better
results than shown; for simplicity of comparison, we restrict attention to
unaccelerated variants. Code222See python source code on Github:
github.com/TypalAcademy/proximal-projection-algorithm. was run on a Macbook
with an Apple M1 Pro chip and 16 GB RAM.

###  3.1 Basis Pursuit

In the field of compressed sensing, the aim is to recover a sparse signal
$x^{\star}\in\mathbb{R}^{n}$ via a collection of linear measurements
$b\in\mathbb{R}^{m}$ (see the survey [20]). If $m<n$ and the matrix
$A\in\mathbb{R}^{m\times n}$ defining the measurements satisfies certain
conditions (e.g. restricted isometry [21, 22]) the signal $x^{\star}$ is often
the solution to the problem

$\min_{x\in\mathbb{R}^{n}}\|x\|_{1}\quad\mbox{s.t.}\quad Ax=b.$ (BP)

Here we set the matrix $A$ to have i.i.d. Gaussian entries, $m=500$, and
$n=2000$. Elements of $x^{\star}$ are independently nonzero with probability
$p=0.05$ and the nonzero values are i.i.d. Gaussian. Algorithm 1 is used with
$f(x)=\|x\|_{1}$, for which

$\mbox{prox}_{\alpha\|\cdot\|_{1}}(x)=\mbox{shrink}(x,\
\alpha)\triangleq\mbox{sgn}(x)\odot\max\\{|x|-\alpha,0\\},$ (7)

where $\odot$ is the element-wise product and sgn is the sign function with
value $1$ for positive input, $-1$ for negative input, and $0$ otherwise.
Further setup details are described in Appendix B.1.

We compare PP against similar first-order methods: Linearized Bregman (LB)
[23, 24, 25], Primal Dual Hybrid Gradient (PDHG) [1], and the Linearized
Method of Multipliers333This is also referred to as an instance of the
proximal augmented Lagrangian method. (LMM) [13]. Figure 1a shows PP is
feasible at each iteration (to machine precision) while the benchmark methods
reduce iteratively reduce constraint violation. The other plots in Figure 1b
and 1c show PP linearly converges to machine precision within hundreds of
iterations while the other methods exhibit sublinear convergence. The mean
time for 10 trials of PP, LB, PDHG, and LMM to compute 2,000 iterations were,
respectively, 8.06s, 7.40s, 20.60s, and 22.85s. Specifically, PP has
comparable per-iteration cost to LB and less per-iteration cost than PDHG and
LMM. Thus, in this example, PP is fastest (converging in finite steps) and the
only method to achieve feasibility.

(a) Violation $\|Ax^{k}-b\|$ (b) Objective $\|x^{k}\|_{1}$ (c) Residual
$\|x^{k+1}-x^{k}\|$

Figure 1:  Basis pursuit (BP) plots, using median for 10 samples of $A$ and
$x^{\star}$. Each PP iterate is feasible, as shown in (a). Note LB and PDHG
update other variables for several steps before $\|x^{k+1}-x^{k}\|\neq 0$.

(a) Original Image $M$

(b) Low Rank $L$

(c) Sparse $S$

Figure 2: Example output from PP (Algorithm 3) with $M$ as a video consisting
of 250 grayscale 960$\times$540 images (each image forming a column of $M$).
The video shows a group of people entering from the right and walking to the
left. The matrix $M$ is approximately represented as the sum of a low rank
matrix $L$ (background) and a sparse matrix $S$ (people walking).

###  3.2 Stable Principal Component Pursuit

Many applications with high dimensional data (e.g. image alignment [26],
hyperspectral image restoration [27], scene triangulation [28]) consider the
problem of recovering a low-rank matrix (i.e. the principal components) from a
high-dimensional matrix that has sparse errors and entry-wise noise. For a
matrix $M\in\mathbb{R}^{n_{1}\times n_{2}}$, the task we consider is finding a
corresponding low-rank matrix $L\in\mathbb{R}^{n_{1}\times n_{2}}$ and a
sparse matrix $S\in\mathbb{R}^{n_{1}\times n_{2}}$ such that

$M=L+S+N,$ (8)

where $N$ is a noise term, for which we assume $\|N\|_{F}\leq\varepsilon$.
Following [5] who proved the effectivenss of this model, we estimate $L$ and
$S$ via

$\min_{L,S}\
\|L\|_{\star}+\lambda\|S\|_{1}\quad\mbox{s.t.}\quad\|L+S-M\|_{F}\leq\varepsilon,$
(SPCP)

where $\lambda=1/\sqrt{m}$ (as chosen in [29]). Letting $U\Sigma V^{\top}$ be
the SVD for $X$, the proximal operator for the nuclear norm
$\|\cdot\|_{\star}$ is the singular value threshold (see [30]):

$\mbox{svt}(X,\alpha)\triangleq\mbox{prox}_{\alpha\|\cdot\|_{\star}}(W)=U\,\mbox{diag}(\max\\{\sigma_{i}-\alpha,\
0\\})\,V^{\top}.$ (9)

Rather than compute an SVD, the threshold $\mbox{svt}(\cdot,\alpha)$ may also
be computed more quickly without an SVD by using a polar decomposition and
projection [10]. In the framing of (P), here the matrix $A$ takes the form
$A=[\mbox{I}\ \mbox{I}]$, where we set $X=[X_{L};X_{S}]$. This simple
structure enables the projection formula in Proposition 1 to admit a simple,
explicit expression, as outlined by the following lemma. (See Appendix B.2 for
a proof.)

######  Lemma 3.1 (SPCP Projection).

If
$\mathcal{C}\triangleq\\{X=[X_{L};X_{S}]:\|X_{L}+X_{S}-M\|_{F}\leq\varepsilon\\}$,
then

$\mathcal{P}_{\mathcal{C}}(Z)=\left[\begin{array}[]{c}Z_{L}-\mu_{Z}(Z_{L}+Z_{S}-M)\\\
Z_{S}-\mu_{Z}(Z_{L}+Z_{S}-M)\end{array}\right],$ (10)

where444We adopt the convention of using $-\varepsilon/0=-\infty$.

$\mu_{Z}\triangleq\max\left\\{0,\
\frac{\|Z_{L}+Z_{S}-M\|_{F}-\varepsilon}{2\|Z_{L}+Z_{S}-M\|_{F}}\right\\}.$
(11)

Using the shrink in (7) and singular value threshold in (9) along with the
projection formula in Lemma 3.1, a special case of Algorithm 2 for the problem
(SPCP) is given by Algorithm 3 below.

######  Remark 5 (Similarity of PP to Proximal Gradient).

When $\varepsilon=0$ and the constraint in (SPCP) is moved into the objective
as a quadratic penalty, the sequence of estimates generated by proximal
gradient for this softly constrained version are identical (for a particular
stepsize) to the sequence $\\{Z^{k}\\}$ generated by PP. The difference is PP
defines solution estimates as the projection of $Z^{k}$ onto the set of all
$X$ satisfying $X_{L}+X_{S}=M$.

Algorithm 3 Proximal Projection for Stable Principal Component Pursuit (PP-
SPCP)

1: PP-SPCP($M$, $\lambda$, $\varepsilon$):

2: initialize parameter $\alpha>0$, matrices $Z_{L}\leftarrow M$,
$Z_{S}\leftarrow 0$, $X_{L}\leftarrow M$, $X_{S}\leftarrow 0$

3: while stopping criteria not met

4: $Z_{L}\leftarrow Z_{L}+\mbox{svt}\left(2X_{L}-Z_{L},\alpha\right)-X_{L}$
$\vartriangleleft$ Use (9) or method in [10]

5: $Z_{S}\leftarrow Z_{S}+\mbox{shrink}\left(2X_{S}-Z_{S},\
\alpha\lambda\right)-X_{S}$ $\vartriangleleft$ Use (7)

6:
$\mu\leftarrow\max\left\\{\frac{\|X_{L}+X_{S}-M\|_{F}-\varepsilon}{2\|X_{L}+X_{S}-M\|_{F}},\
0\right\\}$ $\vartriangleleft$ Treat $-\varepsilon/0$ as $-\infty$

7: $X_{L}\leftarrow Z_{L}-\mu(Z_{L}+Z_{S}-M)$

8: $X_{S}\leftarrow Z_{S}-\mu(Z_{L}+Z_{S}-M)$

9: return $X_{L}$ and $X_{S}$

######  Remark 6 (Proximal Gradient and SPCP).

Prior work (e.g. [5]) uses proximal gradient (PG) for (SPCP), with the
constraint moved into the objective via a quadratic penalty. The weight of
this penalty is important. If it is too large, then the constraint is
satisfied in (SPCP) and the objective for PG is suboptimal. If it is too
small, then the constraint is not satisfied. For a particular weight,
dependent on $M$ and $\varepsilon$, PG is well-known to solve (SPCP). However,
we note PG is is generally not apt when it is unclear how to pick the penalty
weight.

We compare PP to a variant of alternating splitting augmented Lagrangian
method (VASALM) [31] and a “partially smoothed” proximal gradient555In [9],
FISTA accleration was used; however, as noted above, unaccelerated variants
are compared in this work. (PSPG) on a variation (SPCPμ) of (SPCP) wherein the
nuclear norm is smoothed (via a tunable parameter $\mu$); this smoothing was
proposed in [9]. Note VASALM, PSPG, and PP have essentially identical per-
iteration costs (dominated by $\mbox{svt}(\cdot)$). Here $M$ encodes a
sequence of images in a video derived from the PNNL Parking Lot 1 dataset
[32], with each image vectorized to form a column of $M$ (see Figure 2).
Figure 3a shows PP and PSPG666Constraint violation for PSPG (beyond floating
point precision) occurs due to inexactly solving a 1D root-finding problem.
are feasible at each iteration while VASALM approaches feasibility
asymptotically. Due to this feasibility, objective values for PP and PSPG are
bounded from below by the optimal value. On the other hand, VASALM oscillates
around the optimal value. Figures 3b and 3c shows update residuals and
objective values of PP and VASALM have comparable convergence speed, with
slight favor toward PP. Meanwhile, PSPG is notably slower. Overall, PP
performs at least as well as VASALM and PSPG.

(a) Violation
$\max\\{\|X_{L}^{k}+X_{S}^{k}-M\|_{F}-\varepsilon,0\\}/\varepsilon$ (b)
Objective $\|L^{k}\|_{\star}+\lambda\|S^{k}\|_{1}$ (c) Residual
$\|X^{k+1}-X^{k}\|_{F}/\|M\|_{F}$

Figure 3: Convergence comparison for PP, VASALM and PSPG for (SPCP). Each PP
estimate is feasible to machine precision and PSPG estimates are feasible to
accuracy of 1D root-finding solver.

Algorithm 4 Proximal Projection for Computation of Earth Mover’s Distance (PP-
EMD)

1: PP-EMD($\rho^{0}$, $\rho^{1}$):

2: initialize matrix $z\in\mathbb{R}^{2(n-1)\times n}$ and parameter
$\alpha>0$

3: while stopping criteria not met

4: if $\|\mbox{div}(z)+\rho^{1}-\rho^{0}\|_{F}\leq\varepsilon$

5: $m\leftarrow z$

6: else

7: $\tau\leftarrow\mbox{solution}\left(1=\tau\|Y_{\tau}\|_{F}\right)$
$\vartriangleleft$ Use root-finding algorithm and $Y_{\tau}$ in Lemma 3.2

8: $m\leftarrow z-K^{\top}U\,[Y_{\tau};\,Y_{\tau}^{\top}]\,U^{\top}$
$\vartriangleleft$ Use $U$ and $Y_{\tau}$ as in Lemma 3.2

9: $z\leftarrow z+\mbox{shrink}(2m-z,\alpha)-m$ $\vartriangleleft$ Use (7)

10: return $m$

###  3.3 Earth Mover’s Distance

The earth mover’s distance (EMD) is a key metric that is widely used in
several fields (e.g. image processing [33], statistics [34], image retrieval
[35], seismology [36, 37], machine learning [38]). This numerical example
shows PP can be an efficient method for estimation of EMD, outperforming some
existing alternatives.

Consider two distributions $\rho^{0}$ and $\rho^{1}$, which are here
represented as matrices in $\mathbb{R}^{n\times n}$. Let $m^{1}$ and $m^{2}$
together denote components of a 2D flux. The divergence operator is here
defined for a grid size $h>0$ by

$\mbox{div}(m)\triangleq\frac{1}{h}\left(m^{1}_{i,j}-m_{i-1,j}^{1}+m_{i,j}^{2}-m_{i,j-1}^{2}\right).$
(12)

Following similarly to [39], we estimate the earth mover’s distance777PP
immdiately extends to other distances (e.g. L2). To keep experiments concise,
we restrict scope to L1. (i.e. Wasserstein-1) between $\rho^{0}$ and
$\rho^{1}$ as the optimal objective value for the problem

$\min_{m}\|m\|_{1}\quad\mbox{s.t.}\quad\|\mbox{div}(m)+\rho^{1}-\rho^{0}\|_{F}\leq\varepsilon,$
(EMDε)

where $\varepsilon\geq 0$ is small. (See the note in Remark 7 below about
$\varepsilon$.) Neumann boundary conditions are enforced by taking
$m_{0,j}^{1}=m_{i,0}^{2}=m_{n,j}^{2}=m_{i,n}^{2}=0$ for all
$i,j=1,2,\ldots,n$. To write the divergence in matrix form, let $K$ be the
backward differencing operator, i.e.

$K\triangleq\frac{1}{h}\left[\begin{array}[]{rrrrr}1&\\\ -1&1&\\\
&\ddots&\ddots\\\ &&-1&1\\\ &&&-1\\\
\end{array}\right]\in\mathbb{R}^{n\times(n-1)}.$ (13)

By eliminating flux entries that are zero due to boundary conditions, we may
set $m=[m^{1};(m^{2})^{\top}]\in\mathbb{R}^{2(n-1)\times n}$. Thus,
application of the divergence and its transpose can be expresed via888The
notation with left multiplication by $A$ and $A^{\top}$ is somewhat “abusive”
as right multiplications and transposes are included.

$Am=\mbox{div}(m)=Km^{1}+m^{2}K^{\top}\in\mathbb{R}^{n\times
n}\quad\mbox{and}\quad A^{\top}b=\left[\begin{array}[]{c}K^{\top}b\\\
(bK)^{\top}\end{array}\right]=\left[\begin{array}[]{c}K^{\top}b\\\
K^{\top}b^{\top}\end{array}\right]\in\mathbb{R}^{2(n-1)\times n}.$ (14)

With this notation established, the projection formula in Proposition 1 can be
rewritten for this setting as below (see Appendix B.3 for a proof).
Importantly, we note the divergence operator here does not have full row rank;
the fact $\varepsilon>0$ is what ensures the needed matrix inversion is
possible.

######  Lemma 3.2 (EMD Projection).

If
$\mathcal{C}\triangleq\\{m:\|\mbox{div}(m)+\rho^{1}-\rho^{0}\|_{F}\leq\varepsilon\\}$
and $U\Sigma V^{\top}$ is the SVD for $K$, then

$\mathcal{P}_{\mathcal{C}}(z)=\begin{cases}\begin{array}[]{cl}z&\mbox{if}\
\|\mbox{div}(z)+\rho^{1}-\rho^{0}\|_{F}\leq\varepsilon\\\
z-K^{\top}U\left[\begin{array}[]{c}Y_{\tau}\\\
Y_{\tau}^{\top}\end{array}\right]U^{\top}&\mbox{otherwise,}\end{array}\end{cases}$
(15)

where $\tau>0$ satisfies $1=\tau\|(Y_{\tau})\|_{F}$ and $Y_{\tau}$ is the
matrix defined element-wise via

$(Y_{\tau})_{ij}\triangleq\frac{[U^{\top}(Kz+zK^{\top}+\rho^{1}-\rho^{0})U]_{ij}}{\sigma_{i}^{2}+\sigma_{j}^{2}+\varepsilon\tau},\quad\mbox{for
all}\ i,j\in[n].$ (16)

(a) Standing Cat $\rho^{0}$

(b) Crouched Cat $\rho^{1}$

(c) PP Flux Estimate $m^{k}$

Figure 4: The black portion of the cats gives the supports $\rho^{0}$ and
$\rho^{1}$. Arrows in (c) show the 2D flux $m^{k}$.

######  Remark 7.

Setting $\varepsilon=0$ gives the true EMD. For this reason, here we refer to
violation as $\|\mbox{div}(m^{k})+\rho^{1}-\rho^{0}\|_{F}$. However, we
emphasize using PP with small $\varepsilon>0$ yields lower violation than
using PDHG with $\varepsilon=0$. Although counterintuitive, this occurs since
each iterate of PP is feasible while PDHG obtains feasibility asymptotically.

We compare PP to PDHG [39] and G-Prox PDHG [40]. Here, $\varepsilon=10^{-10}$,
and $\rho^{0}$ and $\rho^{1}$ are shown by cat images in Figure 4. (See
Appendix B.3 for details.) The updates for PP and G-Prox PDHG are closely
related, which is reflected by significant overlap of their plots in Figure 5.
Figure 5a shows PP has violation no more than $\varepsilon$ for each
iteration. Although the same $\varepsilon$ is used for projections onto
divergence free vector fields in G-Prox PDHG, its violation grows as that
scheme accumulates projection errors. Hence PP may better handle constraints
than G-Prox PDHG. The times to run PP, G-Prox PDHG, and PDHG for 20,000
iterations were, respectively 107.20s, 43.71s, and 107.38s (averaged over 10
trials). The per-iteration cost of PP roughly equals that of G-Prox PDHG and
2.5X that of PDHG. Despite these higher per-iteration costs, PP and G-Prox
PDHG are more efficient than PDHG here as PDHG requires order(s) of magnitude
more iterations to converge.

(a) Violation $\|\mbox{div}(m^{k})+\rho^{1}-\rho^{0}\|_{F}$ (b) Objective
$\|m^{k}\|_{1}$ (c) Residual $\|m^{k+1}-m^{k}\|_{F}$

Figure 5: Convergence comparison for PP, PDHG, and G-Prox PDHG for computation
of earth mover’s distance between cat images. PP and G-Prox PDHG have
significantly overlapping plots in (b) and (c).

###  3.4 Stable Matrix Completion

A key task in many applications consists of filling in missing entries of a
partially observed matrix. Forms of this problem arise in machine learning,
system identification, recommendations systems, localization in IoT networks,
image restoration, and more (e.g. see [41, 42, 43] and the references
therein). Although infinitely many matrices are consistent with partial
observations (causing recovery to be ill-posed), a common approach is to find
the matrix with minimal nuclear norm that is consistent with the observations.
Indeed, the seminal work [44] shows many low-rank matrices
$M\in\mathbb{R}^{n_{1}\times n_{2}}$ can be exactly recovered by solving the
convex problem

$\min_{X\in\mathbb{R}^{n_{1}\times n_{2}}}\|X\|_{\star}\quad\mbox{s.t.}\quad
X_{ij}=M_{ij}\ \mbox{for all}\ (i,j)\in\Omega,$ (MC)

where $\Omega\subseteq[n_{1}]\times[n_{2}]$ is the set of observed indices.
For settings with observations corrupted by Gaussian noise, we consider the
stable matrix completion variant999Note (SMC) and (SPCP) can both be expressed
as special cases of a common, more general formulation. proposed by [6]:

$\min_{X\in\mathbb{R}^{n_{1}\times
n_{2}}}\|X\|_{\star}\quad\mbox{s.t.}\quad\|\mathcal{P}_{\Omega}(X-M)\|_{F}\leq\varepsilon.$
(SMC)

where we let $\mathcal{P}_{\Omega}$ denote the projection onto the the span of
matrices vanishing outside of $\Omega$ and $\mathcal{P}_{\Omega^{\perp}}$ be
the projection onto the complement of this space, i.e.

$\mathcal{P}_{\Omega}(X)_{ij}\triangleq\begin{cases}\begin{array}[]{cl}X_{ij}&\mbox{if}\
(i,j)\in\Omega\\\
0&\mbox{otherwise}\end{array}\end{cases}\quad\mbox{and}\quad\quad\mathcal{P}_{\Omega^{\perp}}(X)_{ij}\triangleq\begin{cases}\begin{array}[]{cl}X_{ij}&\mbox{if}\
(i,j)\notin\Omega\\\ 0&\mbox{otherwise.}\end{array}\end{cases}$ (17)

The projection formula from Proposition 1 simplifies in this setting as
follows (see Appendix B.4 for a proof).

######  Lemma 3.3 (SMC Projection).

If
$\mathcal{C}\triangleq\\{X:\|\mathcal{P}_{\Omega}(X-M)\|_{F}\leq\varepsilon\\}$,
then

$\mathcal{P}_{\mathcal{C}}(Z)=\begin{cases}\begin{array}[]{cl}Z&\mbox{if}\
\|\mathcal{P}_{\Omega}(Z-M)\|_{F}\leq\varepsilon\\\
\mathcal{P}_{\Omega}\left(\frac{[\|\mathcal{P}_{\Omega}(Z-M)\|_{F}-\varepsilon]M+\varepsilon
Z}{\|\mathcal{P}_{\Omega}(Z-M)\|_{F}}\right)+\mathcal{P}_{\Omega^{\perp}}(Z)&\mbox{otherwise,}\end{array}\end{cases}$
(18)

With this notation and projection formula established, the special case of PP
for (SMC) is Algorithm 5 below.

Algorithm 5 Proximal Projection for Stable Matrix Completion (PP-SMC)

1: PP-SMC($\mathcal{P}_{\Omega}(M),\ \varepsilon$):

2: initialize parameter $\alpha>0$, matrices
$Z\leftarrow\mathcal{P}_{\Omega}(M)$, $X\leftarrow\mathcal{P}_{\Omega}(M)$
$\vartriangleleft$ Initialize via observed entries

3: while stopping criteria not met

4: $Z\leftarrow Z+\mbox{svt}(2X-Z,\ \alpha)-X$ $\vartriangleleft$ Use (9)

5: if $\|\mathcal{P}_{\Omega}(Z-M)\|_{F}\leq\varepsilon$ $\vartriangleleft$
Use (17)

6: $X\leftarrow Z$

7: else

8:
$X\leftarrow\mathcal{P}_{\Omega}\left(\frac{[\|\mathcal{P}_{\Omega}(Z-M)\|_{F}-\varepsilon]M+\varepsilon
Z}{\|\mathcal{P}_{\Omega}(Z-M)\|_{F}}\right)+\mathcal{P}_{\Omega^{\perp}}(Z)$
$\vartriangleleft$ Use (17)

9: return $X$

This numerical example consists of solving (SMC) for various Gaussian matrices
$M$. Letting $n=1,000$, we generated $n\times n$ matrices $M$ of rank $r$ by
sampling two $n\times r$ factors $M_{L}$ and $M_{R}$ independently, each
having i.i.i.d Gaussain entries and then setting $M=M_{L}M_{R}^{\top}$. The
set of observed entries $\Omega$ is uniformly sampled among sets of
cardinality $s$. The observations of $M$ are corrupted by noise $N$, which is
also Gaussian, and $\varepsilon$ is set to be the Frobenius norm of the noise
in the observations. Results are shown using various choices of $s$ and $r$,
and we note the degrees of freedom $d_{r}$ of an $n\times n$ matrix of rank
$r$ is $d_{r}=r\cdot(2n-r)$.

(a) Violation
$\dfrac{\max\\{\|\mathcal{P}_{\Omega}(X^{k}-M)\|_{F}-\varepsilon,0\\}}{\varepsilon}$
(b)  Objective $\|X^{k}\|_{\star}$ (c) Residual
$\dfrac{\|X^{k+1}-X^{k}\|_{F}}{\|M\|_{F}}$

Figure 6: Comparison of PP, VASALM, and SPG for solving (SMC) with
$\mbox{rank}(M)=r=10$ and $|\Omega|=s=5d_{r}$. Plots show medians from 10
trials with distinct random seeds.

(a) Violation
$\dfrac{\max\\{\|\mathcal{P}_{\Omega}(X^{k}-M)\|_{F}-\varepsilon,0\\}}{\varepsilon}$
(b)  Objective $\|X^{k}\|_{\star}$ (c) Residual
$\dfrac{\|X^{k+1}-X^{k}\|_{F}}{\|M\|_{F}}$

Figure 7: Comparison of PP, VASALM, and SPG for solving (SMC) with
$\mbox{rank}(M)=r=50$ and $|\Omega|=s=4d_{r}$. Plots show medians from 10
trials with distinct random seeds.

In similar fashion to the example for (SPCP), here we compare PP to (VASALM)
[31] and a “smoothed” proximal gradient (SPG) [9] on a variation (SMCμ) of
(SMC) wherein the nuclear norm is smoothed (via a tunable parameter $\mu$).
Note PP, VASALM, and SPG all have essentially identical per-iteration costs,
which are dominated by the evaluation of $\mbox{svt}(\cdot)$. The parameters
for each algorithm were chosen to give good performance for the case where
$r=50$ and $s/d_{r}=4$. To demonstrate the flexibility of these algorithms,
these parameters were held fixed across all trials and values of $r$ and
$s/d_{r}$.

Numerical results are plotted in Figures 6 and 7 and also shared in Table 1.
As shown in Figures 6a and 7a, PP and SPG are feasible at each iteration while
VASALM iteratively reduces its constraint violation. Figures 6b and 7b show
the objective $\|X^{k}\|_{\star}$ of PP and VASALM converge quickly, and that
SPG converged to a suboptimal limit in each case. Figures 6c and 7c show PP
converges faster than VASALM and both of these converge faster than SPG.
Similar results are reported in Table 1. The following statements hold for all
three configurations shown in the table. PP converges in fewer steps than
VASALM and SPG (i.e. about 55% as many steps as VASALM and 40% as many as
SPG). The objective values $\|X^{k}\|_{\star}$ of PP and VASALM agree to three
digits of accuracy. The output of VASALM has a violation exceeding machine
precision. SPG performs worse than PP with respect to every metric. Overall,
PP outperforms VASALM and SPG in these examples.

Unknown Matrix $M$ | Result
---|---
rank $(r)$ | $s/d_{r}$ | $s/n^{2}$ | PP | VASALM | SPG
# | Viol | $\|X^{k}\|_{\star}$ | # | Viol | $\|X^{k}\|_{\star}$ | # | Viol | $\|X^{k}\|_{\star}$
10 | 5 | 0.0995 | 105 | 1.80e-16 | 9.49e+03 | 185 | 5.45e-04 | 9.49e+03 | 247 | 3.24e-16 | 1.09e+04
50 | 4 | 0.3900 | 61 | 3.28e-16 | 4.82e+04 | 110 | 1.35e-03 | 4.82e+04 | 145 | 9.29e-16 | 4.96e+04
100 | 3 | 0.5700 | 54 | 1.16e-15 | 9.56e+04 | 96 | 1.30e-03 | 9.56e+04 | 133 | 6.02e-16 | 9.69e+04
Table 1: Results for stable matrix completion example with
$M\in\mathbb{R}^{n\times n}$ and $n=1,000$. Each algorithm executed until
$\|X^{k+1}-X^{k}\|_{F}/\|M\|_{F}\leq 10^{-5}$. The number of iterations $k$ to
meet this stopping condition is denoted by #. The relative constraint
violation
$\max\\{\|\mathcal{P}_{\Omega}(X^{k}-M)\|_{F}-\varepsilon,0\\}/\varepsilon$ is
denoted by “Viol.” Note an $n\times n$ matrix of rank $r$ depends upon
$r(2n-r)$ degrees of freedom $d_{r}$. Results shown are averages for 10 random
seeds.

##  4 Discussion

This work proposes the proximal projection (PP) algorithm for minimizing
proximable functions subject to satisfication of linear constraints within a
specified tolerance. It has a simple and easy-to-code formulation. The primary
novelty herein is showing how to project onto the constraint set. This can be
efficiently computed by solving a 1D root-finding problem; in some important
instances, explicit formulas are available. This projection is interwoven with
proximal operations via Douglas-Rachford splitting to obtain PP. Our formal
analysis shows PP converges, with any choice of step size, to an optimal
solution. Moreover, unlike many algorithms that only achieve feasibility
asymptotically, the output of PP is feasible for any stopping iteration.

Numerically, we show PP performs favorably against alternatives in a few
important applications: basis pursuit, stable principal component pursuit,
earth mover’s distances, and matrix completion. In each example provided, PP
is numerically verified to be feasible and to converge to an optimal solution.
Moreover, PP performs at least as well as each alternative method in these
examples with respect to every metric considered.

The presented projection formula provided can be included in other algorithms
(e.g. projected gradient when the objective is smooth); we leave investigation
of this to future work. Acceleration techniques [45, 46, 47] can be used on
Douglas-Rachford splitting, which therefore apply to PP; future work could
examine the performance of these. Lastly, extensions may also involve PP’s
incorporation as an optimization layer in machine learning applications (e.g.
the realm of learning-to-optimize and predict-then-optimize [48, 49, 50, 51,
52]).

## Acknowledgements

We thank Samy Wu Fung for his helpful discussions and feedback on early drafts
of this manuscript.

## References

  * [1] Antonin Chambolle and Thomas Pock.  A first-order primal-dual algorithm for convex problems with applications to imaging.  Journal of mathematical imaging and vision, 40(1):120–145, 2011.
  * [2] Marguerite Frank, Philip Wolfe, et al.  An algorithm for quadratic programming.  Naval research logistics quarterly, 3(1-2):95–110, 1956.
  * [3] Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, Jonathan Eckstein, et al.  Distributed optimization and statistical learning via the alternating direction method of multipliers.  Foundations and Trends® in Machine learning, 3(1):1–122, 2011.
  * [4] Wei Deng and Wotao Yin.  On the global and linear convergence of the generalized alternating direction method of multipliers.  Journal of Scientific Computing, 66:889–916, 2016.
  * [5] Zihan Zhou, Xiaodong Li, John Wright, Emmanuel Candes, and Yi Ma.  Stable Principal Component Pursuit.  In 2010 IEEE international symposium on information theory, pages 1518–1522. IEEE, 2010.
  * [6] Emmanuel J Candes and Yaniv Plan.  Matrix completion with noise.  Proceedings of the IEEE, 98(6):925–936, 2010.
  * [7] Yair Censor, Ran Davidi, Gabor T Herman, Reinhard W Schulte, and Luba Tetruashvili.  Projected subgradient minimization versus superiorization.  Journal of Optimization Theory and Applications, 160:730–747, 2014.
  * [8] Perturbation resilience and superiorization of iterative algorithms.  Inverse problems, 26(6):065008, 2010.
  * [9] Necdet Serhat Aybat, Donald Goldfarb, and Shiqian Ma.  Efficient algorithms for robust and stable principal component pursuit problems.  Computational Optimization and Applications, 58:1–29, 2014.
  * [10] Jian-Feng Cai and Stanley Osher.  Fast Singular Value Thresholding Without Singular Value Decomposition.  Methods and Applications of Analysis, 20(4):335–352, 2013.
  * [11] Amir Beck.  First-Order Methods in Optimization.  SIAM, 2017.
  * [12] Heinz H Bauschke and Patrick L Combettes.  Convex Analysis and Monotone Operator Theory in Hilbert Spaces.  Springer, 2017.
  * [13] Ernest Ryu and Wotao Yin.  Large-Scale Convex Optimization: Algorithm Designs via Monotone Operators.  Cambridge University Press, 2022.
  * [14] Jonathan Eckstein and Dimitri P Bertsekas.  On the Douglas—Rachford Splitting Method and the Proximal Point Algorithm for Maximal Monotone Operators.  Mathematical programming, 55:293–318, 1992.
  * [15] Pierre-Louis Lions and Bertrand Mercier.  Splitting Algorithms for the Sum of Two Nonlinear Operators.  SIAM Journal on Numerical Analysis, 16(6):964–979, 1979.
  * [16] Jonathan Eckstein.  Splitting Methods for Monotone Operators with Applications to Parallel Optimization.  PhD thesis, Massachusetts Institute of Technology, 1989.
  * [17] Scott B Lindstrom and Brailey Sims.  Survey: Sixty Years of Douglas–Rachford.  Journal of the Australian Mathematical Society, 110(3):333–370, 2021.
  * [18] Francisco J Aragón Artacho, Jonathan M Borwein, and Matthew K Tam.  Recent results on Douglas–Rachford methods for combinatorial optimization problems.  Journal of Optimization Theory and Applications, 163:1–30, 2014.
  * [19] Llewellyn Hilleth Thomas.  Elliptic problems in linear difference equations over a network.  Watson Sci. Comput. Lab. Rept., Columbia University, New York, 1:71, 1949.
  * [20] Saad Qaisar, Rana Muhammad Bilal, Wafa Iqbal, Muqaddas Naureen, and Sungyoung Lee.  Compressive sensing: From theory to applications, a survey.  Journal of Communications and networks, 15(5):443–456, 2013.
  * [21] Emmanuel J Candes and Terence Tao.  Decoding by linear programming.  IEEE transactions on information theory, 51(12):4203–4215, 2005.
  * [22] Emmanuel J Candes and Terence Tao.  Near-optimal signal recovery from random projections: Universal encoding strategies?  IEEE transactions on information theory, 52(12):5406–5425, 2006.
  * [23] Jian-Feng Cai, Stanley Osher, and Zuowei Shen.  Linearized Bregman Iterations for Compressed Sensing.  Mathematics of computation, 78(267):1515–1536, 2009.
  * [24] Wotao Yin.  Analysis and generalizations of the linearized bregman method.  SIAM Journal on Imaging Sciences, 3(4):856–877, 2010.
  * [25] Stanley Osher, Yu Mao, Bin Dong, and Wotao Yin.  Fast linearized Bregman iteration for compressive sensing and sparse denoising.  2010\.
  * [26] Yigang Peng, Arvind Ganesh, John Wright, Wenli Xu, and Yi Ma.  RASL: Robust alignment by sparse and low-rank decomposition for linearly correlated images.  IEEE transactions on pattern analysis and machine intelligence, 34(11):2233–2246, 2012.
  * [27] Hongyan Zhang, Wei He, Liangpei Zhang, Huanfeng Shen, and Qiangqiang Yuan.  Hyperspectral image restoration using low-rank matrix recovery.  IEEE transactions on geoscience and remote sensing, 52(8):4729–4743, 2013.
  * [28] Zhengdong Zhang, Arvind Ganesh, Xiao Liang, and Yi Ma.  Tilt: Transform invariant low-rank textures.  International journal of computer vision, 99:1–24, 2012.
  * [29] Emmanuel J Candès, Xiaodong Li, Yi Ma, and John Wright.  Robust Principal Component Analysis?  Journal of the ACM (JACM), 58(3):1–37, 2011.
  * [30] Jian-Feng Cai, Emmanuel J Candès, and Zuowei Shen.  A singular value thresholding algorithm for matrix completion.  SIAM Journal on optimization, 20(4):1956–1982, 2010.
  * [31] Min Tao and Xiaoming Yuan.  Recovering low-rank and sparse components of matrices from incomplete and noisy observations.  SIAM Journal on Optimization, 21(1):57–81, 2011.
  * [32] Pacific Northwest National Lab.  PNNL Parking Lot 1 Dataset, 2012.  www.crcv.ucf.edu/research/data-sets/pnnl-parking-lot.
  * [33] Sylvain Boltz, Frank Nielsen, and Stefano Soatto.  Earth mover distance on superpixels.  In 2010 IEEE International Conference on Image Processing, pages 4597–4600. IEEE, 2010.
  * [34] Victor M Panaretos and Yoav Zemel.  Statistical aspects of Wasserstein distances.  Annual review of statistics and its application, 6(1):405–431, 2019.
  * [35] Yossi Rubner, Carlo Tomasi, and Leonidas J Guibas.  The earth mover’s distance as a metric for image retrieval.  International journal of computer vision, 40:99–121, 2000.
  * [36] Application of the Wasserstein metric to seismic signals.  Communications in Mathematical Sciences, 2014.
  * [37] L. Métivier, R. Brossier, Q. Mérigot, E. Oudet, and J. Virieux.  Measuring the misfit between seismograms using an optimal transport distance: application to full waveform inversion.  Geophysical Journal International, 205(1):345–377, 2016.
  * [38] Ilya Tolstikhin, Olivier Bousquet, Sylvain Gelly, and Bernhard Schoelkopf.  Wasserstein auto-encoders.  arXiv preprint arXiv:1711.01558, 2017.
  * [39] Wuchen Li, Ernest K Ryu, Stanley Osher, Wotao Yin, and Wilfrid Gangbo.  A parallel method for earth mover’s distance.  Journal of Scientific Computing, 75(1):182–197, 2018.
  * [40] Matt Jacobs, Flavien Léger, Wuchen Li, and Stanley Osher.  Solving large-scale optimization problems with a convergence rate independent of grid size.  SIAM Journal on Numerical Analysis, 57(3):1100–1123, 2019.
  * [41] Andy Ramlatchan, Mengyun Yang, Quan Liu, Min Li, Jianxin Wang, and Yaohang Li.  A survey of matrix completion methods for recommendation systems.  Big Data Mining and Analytics, 1(4):308–323, 2018.
  * [42] Mark A Davenport and Justin Romberg.  An overview of low-rank matrix recovery from incomplete observations.  IEEE Journal of Selected Topics in Signal Processing, 10(4):608–622, 2016.
  * [43] Luong Trung Nguyen, Junhan Kim, and Byonghyo Shim.  Low-Rank Matrix Completion: A Contemporary Survey.  IEEE Access, 7:94215–94237, 2019.
  * [44] Emmanuel Candes and Benjamin Recht.  Exact matrix completion via convex optimization.  Communications of the ACM, 55(6):111–119, 2012.
  * [45] Andreas Themelis and Panagiotis Patrinos.  SuperMann: a superlinearly convergent algorithm for finding fixed points of nonexpansive operators.  IEEE Transactions on Automatic Control, 64(12):4875–4890, 2019.
  * [46] Junzi Zhang, Brendan O’Donoghue, and Stephen Boyd.  Globally convergent type-I Anderson acceleration for nonsmooth fixed-point iterations.  SIAM Journal on Optimization, 30(4):3170–3197, 2020.
  * [47] Jisun Park and Ernest K Ryu.  Exact optimal accelerated complexity for fixed-point iterations.  In International Conference on Machine Learning, pages 17420–17457. PMLR, 2022.
  * [48] Tianlong Chen, Xiaohan Chen, Wuyang Chen, Howard Heaton, Jialin Liu, Zhangyang Wang, and Wotao Yin.  Learning to optimize: A Primer and a Benchmark.  arXiv preprint arXiv:2103.12828, 2021.
  * [49] Brandon Amos and J Zico Kolter.  Optnet: Differentiable optimization as a layer in neural networks.  In International conference on machine learning, pages 136–145. PMLR, 2017.
  * [50] Nir Shlezinger, Jay Whang, Yonina C Eldar, and Alexandros G Dimakis.  Model-based deep learning.  Proceedings of the IEEE, 111(5):465–499, 2023.
  * [51] Daniel McKenzie, Samy Wu Fung, and Howard Heaton.  Differentiating through integer linear programs with quadratic regularization and davis-yin splitting.  Transactions on Machine Learning Research, 2024.
  * [52] Adam N Elmachtoub and Paul Grigas.  Smart “predict, then optimize”.  Management Science, 68(1):9–26, 2022.
  * [53] Dimitri Bertsekas, Angelia Nedic, and Asuman Ozdaglar.  Convex Analysis and Optimization, volume 1.  Athena Scientific, 2003.
  * [54] Yu Nesterov.  Smooth minimization of non-smooth functions.  Mathematical programming, 103:127–152, 2005.
  * [55] Zhouchen Lin, Arvind Ganesh, John Wright, Leqin Wu, Minming Chen, and Yi Ma.  Fast convex optimization algorithms for exact recovery of a corrupted low-rank matrix.  Coordinated Science Laboratory Report no. UILU-ENG-09-2214, DC-246, 2009.

##  Appendix A Proofs of Main Results

Each of our results from Section 2 is proven below. For ease of reference,
results are restated before their proof. We begin with an auxiliary lemma used
for the proof of Lemma 1.

######  Lemma A.1 (Unique Positive $\tau_{x}$).

If (C2) holds and $x\in\mathbb{R}^{n}$ is such that $\|Ax\|>\varepsilon$, then
there is unique $\tau_{x}\in\mathbb{R}$ satisfying

$0<\tau_{x}\leq\frac{\sigma_{max}(A)^{2}}{\|Ax\|-\varepsilon}$ (19)

for which

$\tau_{x}\|(AA^{\top}+\varepsilon\tau_{x})^{-1}Ax\|=1.$ (20)

###### Proof.

We proceed in the following manner. First, we obtain an alternative expression
for the left hand side of (20) in terms of the SVD of $A$, which is continuous
in $\tau$ (Step 1). For the $\varepsilon=0$ case, we then show the desired
range holds for $\tau_{x}$ (Step 2). For the $\varepsilon>0$ case, we show the
left hand side of (20) goes to zero as $\tau\rightarrow 0^{+}$ and that it is
at least one for $\tau$ sufficently large. This enable application of the
intermediate value theorem to assert existence of the desired $\tau_{x}$ (Step
3). Uniqueness of $\tau_{x}$ is then established (Step 4).

Step 1. Let $U\Sigma V^{\top}$ be the singular value decomposition of $A$,
with $\Sigma=\mbox{diag}(\sigma_{i})$. Then

$AA^{\top}+\varepsilon\tau\mbox{I}=U\Sigma
V^{\top}V\Sigma^{\top}U^{\top}+\varepsilon\tau\mbox{I}=U\,\mbox{diag}(\sigma_{i}^{2}+\varepsilon\tau)U^{\top},\quad\mbox{for
all}\ \tau>0,$ (21)

where the final equality holds since $U$ and $V$ are orthogonal. Whenever
$\tau>0$, (C2) ensures $\sigma_{i}^{2}+\varepsilon\tau>0$ for each index
$i\in[m]$, and so this matrix is invertible. In particular, right mutiplying
this inverse by $\tau A$ yields

$(AA^{\top}+\varepsilon\tau\mbox{I})^{-1}\tau
A=U\,\mbox{diag}\left(\frac{\tau}{\sigma_{i}^{2}+\varepsilon\tau}\right)U^{\top}A=U\,\mbox{diag}\left(\frac{\sigma_{i}\tau}{\sigma_{i}^{2}+\varepsilon\tau}\right)V^{\top},\quad\mbox{for
all}\ \tau>0.$ (22)

Together (22) and the orthogonality of $U$ imply

$\|\tau(AA^{\top}+\varepsilon\tau\mbox{I})^{-1}Ax\|=\left\|U\,\mbox{diag}\left(\frac{\sigma_{i}\tau}{\sigma_{i}^{2}+\varepsilon\tau}\right)V^{\top}x\right\|=\left\|\mbox{diag}\left(\frac{\sigma_{i}\tau}{\sigma_{i}^{2}+\varepsilon\tau}\right)V^{\top}x\right\|.$
(23)

Note each fraction $\sigma_{i}\tau/(\sigma_{i}^{2}+\varepsilon\tau)$ is
continuous on $(0,\infty)$, and the entire expression, being a composition of
continuous functions, is also continuous.

Step 2. If $\varepsilon=0$, then $A$ is full-rank by (C2), which implies
$AA^{\top}$ is invertible. Additionally, this also yields
$(AA^{\top})^{-1}Ax\neq 0$ since $Ax\neq 0$. Thus, letting $\sigma_{+}$ denote
the maximum singular value of $A$,

$0<\tau_{x}=\frac{1}{\|(AA^{\top})^{-1}Ax\|}=\frac{1}{\left\|U\,\mbox{diag}\left(\sigma_{i}^{-2}\right)U^{\top}Ax\right\|}\leq\frac{1}{\sigma_{+}^{-2}\|UU^{\top}Ax\|}=\frac{\sigma_{+}^{2}}{\|Ax\|}.$
(24)

Step 3. For the remainder of the proof, we assume $\varepsilon>0$. Set

$\tau_{+}\triangleq\frac{\sigma_{+}^{2}}{\|Ax\|-\varepsilon}$ (25)

so that

$\frac{\tau_{+}}{\sigma_{i}^{2}+\varepsilon\tau_{+}}\geq\frac{\tau_{+}}{\sigma_{+}^{2}+\varepsilon\tau_{+}}=\frac{\sigma_{+}^{2}}{(\|Ax\|-\varepsilon)\sigma_{+}^{2}+\varepsilon\sigma_{+}^{2}}=\frac{1}{\|Ax\|},\quad\mbox{for
all indices $i$.}$ (26)

This implies

$\|\tau_{+}(AA^{\top}+\varepsilon\tau_{+}\mbox{I})^{-1}Ax\|=\left\|\mbox{diag}\left(\frac{\sigma_{i}\tau_{+}}{\sigma_{i}^{2}+\varepsilon\tau_{+}}\right)V^{\top}x\right\|\geq\frac{\|\mbox{diag}(\sigma_{i})V^{\top}x\|}{\|Ax\|}=\frac{\|U\,\mbox{diag}(\sigma_{i})V^{\top}x\|}{\|Ax\|}=1.$
(27)

On the other hand,

$\lim_{\tau\rightarrow
0^{+}}\|\tau_{+}(AA^{\top}+\varepsilon\tau_{+}\mbox{I})^{-1}Ax\|=\lim_{\tau\rightarrow
0^{+}}\left\|\mbox{diag}\left(\frac{\sigma_{i}\tau}{\sigma_{i}^{2}+\varepsilon\tau}\right)V^{\top}x\right\|=0,$
(28)

where we note, if $\sigma_{i}=0$ for any index $i$, then
$\sigma_{i}\tau/(\sigma_{i}^{2}+\varepsilon\tau)=0$ for all $\tau>0$. Thus, by
the intermediate value theorem, there is $\tau_{x}\in(0,\tau_{+}]$ such that
(20) holds, from which (19) then follows.

Step 4. All that remains is to verify $\tau_{x}$ is unique. To do this, it
suffices to show the square of the norm expression is strictly increasing in
$\tau$. Note

$\frac{\mbox{d}}{\mbox{d}\tau}\left[\frac{\sigma_{i}\tau}{\sigma_{i}^{2}+\varepsilon\tau}\right]=\frac{\sigma_{i}^{3}}{(\sigma_{i}^{2}+\varepsilon\tau)^{2}}\geq
0,\quad\mbox{for all $i$ and}\ \tau>0.$ (29)

Using this fact, differentiating reveals

$\displaystyle\frac{\mbox{d}}{\mbox{d}\tau}\left\|\mbox{diag}\left(\frac{\sigma_{i}\tau}{\sigma_{i}^{2}+\varepsilon\tau}\right)V^{\top}x\right\|^{2}$
$\displaystyle=\frac{\mbox{d}}{\mbox{d}\tau}\sum_{i=1}^{m}\left(\frac{\sigma_{i}\tau}{\sigma_{i}^{2}+\varepsilon\tau}\cdot(V^{\top}x)_{i}\right)^{2}$
(30a)
$\displaystyle=2\cdot\sum_{i=1}^{m}\frac{\sigma_{i}\tau}{\sigma_{i}^{2}+\varepsilon\tau}\cdot\frac{\sigma_{i}^{3}}{(\sigma_{i}^{2}+\varepsilon\tau)^{2}}\cdot(V^{\top}x)_{i}^{2}$
(30b)
$\displaystyle\geq\underbrace{\frac{2\tau}{(\sigma_{+}+\varepsilon\tau)^{3}}}_{>0}\cdot\underbrace{\sum_{i=1}^{m}\sigma_{i}^{4}(V^{\top}x)_{i}^{2}}_{\geq
0},\quad\mbox{for all}\ \tau>0.$ (30c)

Since

$0<\varepsilon<\|Ax\|^{2}=\|\mbox{diag}(\sigma_{i})V^{\top}x\|^{2}=\sum_{i=1}^{m}\sigma_{i}^{2}(V^{\top}x)_{i}^{2},$
(31)

there is an index $\ell$ for which $\sigma_{\ell}>0$ and
$(V^{\top}x)_{\ell}\neq 0$, and so at least one term in the sum in (30c) is
positive, making the entire derivative expression in (30) positive, as
desired. ∎

The next lemma’s proof draws heavily from [11, Lemma 6.68].

Proposition 1 (Projection Formula).  If conditions (C2) and (C3) hold and
$\mathcal{C}\triangleq\\{x:\|Ax-b\|\leq\varepsilon\\}$, then

$\mathcal{P}_{\mathcal{C}}(x)=\begin{cases}\begin{array}[]{cl}x&\mbox{if}\
\|Ax-b\|\leq\varepsilon,\\\
x-A^{\top}(AA^{\top}+\varepsilon\tau_{x}\mbox{I})^{-1}(Ax-b)&\mbox{otherwise},\end{array}\end{cases}$
(32)

where, if $\|Ax-b\|>\varepsilon$, the scalar $\tau_{x}$ is the unique positive
solution to

$1=\tau\|(AA^{\top}+\varepsilon\tau\mbox{I})^{-1}(Ax-b)\|.$ (33)

###### Proof.

The set $\mathcal{C}$ is nonempty by (C3). As $\mathcal{C}$ is also closed and
convex, projections onto $\mathcal{C}$ exist and are well-defined [??]. Let
$u^{\star}$ be the projection of $x$ onto $\mathcal{C}$ so that $u^{\star}$ is
the unique solution to

$\min_{u}\left\\{\frac{1}{2}\|u-x\|^{2}\ :\ \|Au-b\|\leq\varepsilon\right\\}.$
(34)

If $\|Ax-b\|\leq\varepsilon$, then the solution to (34) is $u^{\star}=x$
since, in that case, the objective value is zero and norms are nonnegative.
Henceforth, we assume $\|Ax-b\|>\varepsilon$. Introducing
$z\in\mathbb{R}^{m}$, (34) may be rewritten as

$\min_{u,z}\left\\{\frac{1}{2}\|u-x\|^{2}\ :\ z=Au-b,\
\|z\|\leq\varepsilon\right\\}.$ (35)

The Lagrangian
$\mathcal{L}\colon\mathbb{R}^{n}\times\mathbb{R}^{m}\times\mathbb{R}_{\geq
0}\times\mathbb{R}^{m}\rightarrow\mathbb{R}$ for this problem is

$\displaystyle\mathcal{L}(u,z;\alpha,y)$
$\displaystyle=\frac{1}{2}\|u-x\|^{2}+y^{\top}(z-[Au-b])+\alpha(\|z\|-\varepsilon)$
(36a)
$\displaystyle=\left[\frac{1}{2}\|u-x\|^{2}-(A^{\top}y)^{\top}u\right]+\left[\alpha(\|z\|-\varepsilon)+y^{\top}(z+b)\right].$
(36b)

Since the Lagrangian $\mathcal{L}$ is separable with respect to $u$ and $z$,
the dual objective can be written as

$\min_{u,z}\mathcal{L}(u,z;\alpha,y)=\min_{u}\left[\frac{1}{2}\|u-x\|^{2}-(A^{\top}y)^{\top}u\right]+\min_{z}\left[\alpha(\|z\|-\varepsilon)+y^{\top}(z+b)\right].$
(37)

The minimizer of the minimization problem for $u$ is $u^{\star}=x+A^{\top}y$,
which yields

$\min_{u}\left[\frac{1}{2}\|u-x\|^{2}-(A^{\top}y)^{\top}u\right]=\frac{1}{2}\|u^{\star}-x\|^{2}-(A^{\top}y)^{\top}u^{\star}=-\frac{1}{2}\|Ay\|^{2}-(Ax)^{\top}y.$
(38)

Due to the linear term $y^{\top}z$, the minimization of $\mathcal{L}$ with
respect to $z$ is obtained when $z$ is anti-parallel to $y$, for which the
expression in the $z$ minimization of (37) becomes

$\left(\alpha-\|y\|\right)\|z\|-\alpha\varepsilon+y^{\top}b.$ (39)

Consequently, if $\|y\|>\alpha$, the expression in (39) goes to $-\infty$ as
$\|z\|\rightarrow\infty$. Otherwise, it is minimized by picking $z=0$. In
summary,

$\min_{z}\left[\alpha(\|z\|-\varepsilon)+y^{\top}(z+b)\right]=\begin{cases}\begin{array}[]{cl}-\alpha\varepsilon+y^{\top}b&\mbox{if}\
\|y\|\leq\alpha\\\ -\infty&\mbox{otherwise}.\end{array}\end{cases}$ (40)

Combining (37), (38), and (40), we obtain the dual problem

$\max_{\alpha\geq 0,\
y\in\mathbb{R}^{m}}\left\\{-\frac{1}{2}\|Ay\|^{2}-(Ax)^{\top}y-\alpha\varepsilon+y^{\top}b\
:\ \|y\|\leq\alpha\right\\}.$ (41)

Importantly, strong duality holds for the primal-dual pair of problems (34)
and (41) [53, Proposition 6.4.4]. Thus, upon finding an optimal dual solution
$y^{\star}$, the primal solution is $x^{\star}=x=A^{\top}y^{\star}$.

The dual objective is strictly decreasing as $\alpha$ increases beyond
$\|y\|$, and so the optimal choice for $\alpha$ is always $\alpha=\|y\|$.
Thus, the dual problem (41) may be simplified to

$\min_{y\in\mathbb{R}^{m}}\left\\{\frac{1}{2}\|A^{\top}y\|^{2}+(Ax)^{\top}y+\varepsilon\|y\|-y^{\top}b\right\\}.$
(42)

If the optimal dual variable $y^{\star}$ were zero, it would follow that
$\|Au^{\star}-b\|=\|A(x+A^{\top}0)-b\|=\|Ax-b\|>\varepsilon$, a contradiction
to the constraint in (34). Thus, $y^{\star}$ is nonzero and satisfies the
first-order optimality condition

$0=\frac{\mbox{d}}{\mbox{d}y}\left[\frac{1}{2}\|A^{\top}y\|^{2}+(Ax)^{\top}y+\varepsilon\|y\|-y^{\top}b\right]_{y=y^{\star}}=AA^{\top}y^{\star}+Ax+\frac{\varepsilon
y^{\star}}{\|y^{\star}\|}-b=\left(AA^{\top}+\frac{\varepsilon\mbox{I}}{\|y^{\star}\|}\right)y^{\star}+Ax-b.$
(43)

By (C2), either $A$ has full row rank or $\varepsilon>0$, and so the matrix
$AA^{\top}+\varepsilon I/\|y^{\star}\|$ is invertible. Thus,

$y^{\star}=-\left(AA^{\top}+\frac{\varepsilon\mbox{I}}{\|y^{\star}\|}\right)^{-1}(Ax-b).$
(44)

Taking norms of both sides reveals

$\|y^{\star}\|=\left\|\left(AA^{\top}+\frac{\varepsilon\mbox{I}}{\|y^{\star}\|}\right)^{-1}(Ax-b)\right\|.$
(45)

By Lemma A.1,101010Since $b\in\mbox{im}(A)$ by (C3), there is $\xi$ such that
$A\xi=b$ and the lemma can be applied with the $x-\xi$ in place of $x$. there
is a unique scalar $\tau_{x}>0$ for which

$\frac{1}{\tau_{x}}=\left\|(AA^{\top}+\varepsilon\tau_{x}\mbox{I})^{-1}(Ax-b)\right\|.$
(46)

Hence $\|y^{\star}\|=1/\tau_{x}$ and (44) becomes

$y^{\star}=-\left(AA^{\top}+\varepsilon\tau_{x}\mbox{I}\right)^{-1}(Ax-b),$
(47)

from which we conclude

$u^{\star}=x+A^{\top}y^{\star}=x-A^{\top}(AA^{\top}+\varepsilon\tau\mbox{I})^{-1}Ax,$
(48)

as desired. ∎

We now verify our convergence result, which is a special case of existing
results for Douglas Rachford splitting.

Theorem 1 (Convergence of PP).  If conditions (C1) to (C5) hold, then the
sequences $\\{x^{k}\\}$ and $\\{z^{k}\\}$ generated by (5) converge, with
$\\{x^{k}\\}$ converging to a solution of (P). Moreover,
$\|Ax^{k}-b\|\leq\varepsilon$ for all $k$.

###### Proof.

By (C5), $\mbox{ri}(\mbox{dom}(f))\cap\mbox{ri}(\mathcal{C})\neq\varnothing$
and, in particular, $\mbox{dom}(f)\cap\mathcal{C}\neq\varnothing$. This fact,
together with (C4) and $\mathcal{C}$ being closed and convex, enable [12,
Proposition 11.15] to be applied to deduce the existence of a solution to (P).
Futhermore, by [12, Proposition 6.19],
$0\in\mbox{sri}(\mathcal{C}-\mbox{dom}(f)),$111111Here $\mbox{sri}(\cdot)$
denotes the strong relative interior (see [12, Definition 6.9]). and so [12,
Proposition 27.8] may be applied to deduce

$\\{x:0\in\partial f(x)+\partial
N_{\mathcal{C}}(x)\\}=\operatorname*{\mbox{argmin}}_{x\in\mathcal{C}}f(x)=\operatorname*{\mbox{argmin}}_{x\in\mathbb{R}^{n}}\Big{\\{}f(x)\,:\,\|Ax-b\|\leq\varepsilon\Big{\\}},$
(49)

where $N_{\mathcal{C}}$ is the normal cone operator for $\mathcal{C}$. By (C2)
and (C3) and Proposition 1, the update for each $x^{k}$ is precisely the
projection of $z^{k}$ onto $\mathcal{C}$. Thus, the iteration (5) is an
instance of Douglas-Rachford splitting [12, 15]. By [12, Theorem 26.11],
$\\{z^{k}\\}$ and $\\{x^{k}\\}$ converge, and the limit $x^{\star}$ of
$\\{x^{k}\\}$ satisfies

$0\in\partial f(x^{\star})+N_{\mathcal{C}}(x^{\star}).$ (50)

By (49) and (50), we conclude $x^{\star}$ is a solution to (P), as desired.
Lastly, note $\|Ax^{k}-b\|\leq\varepsilon$ since, as noted above,
$x^{k}=\mathcal{P}_{\mathcal{C}}(z^{k})\in\mathcal{C}$ for all
$k\in\mathbb{N}$. ∎

##  Appendix B Numerical Examples Supplement

A subsection is dedicated herein to providing further details for each of the
numerical examples, particularly formulations of the algorithms to which PP is
compared and proofs for the special cases of the projection formula in
Proposition 1 to the various settings.

###  B.1 Basis Pursuit

We initialize iterates to the zero vector (e.g. $z^{1}=0$ for PP). Entries of
$A$ are drawn from $A_{ij}\sim\mathcal{N}(0,1/m)$. In each case, we attempted
to pick parameters that yield best performance while respecting conditions
needed to ensure convergence guarantees. Here 10 trials were used, with the
mean time reported in the main text and medians used for the plots in Figure
1.

Proximal Projection (PP). We applied Algorithm 1 with $\alpha=0.1$ and the
shrink operator.

Linearized Bregman (LB). Rather than directly minimize $\ell_{1}$, linearized
Bregman solves

$\min_{x}\mu\|x\|_{1}+\frac{1}{2\alpha}\|u\|^{2}\quad\mbox{s.t.}\quad Ax=b,$
(51)

which yields the same result as (BP) when $\alpha$ is sufficiently large.
Following [23], we use the iteration

$\displaystyle v^{k+1}$ $\displaystyle=v^{k}-A^{\top}(Ax^{k}-b),$ (52a)
$\displaystyle x^{k+1}$ $\displaystyle=\mbox{shrink}\left(\alpha v^{k+1},\
\alpha\mu\right).$ (52b)

By [23, Theorem 2.4], this iteration converges for
$\alpha\in(0,2/\|AA^{\top}\|)$; we used $\mu=2\|AA^{\top}\|$ and
$\alpha=2/\|AA^{\top}\|$.

Linearized Method of Multipliers (LMM). For step sizes $\alpha,\lambda>0$, the
linearized method of multipliers solves (BP) using iterates of the form

$\displaystyle x^{k+1}$ $\displaystyle=\mbox{shrink}\left(x^{k}-\alpha
A^{\top}[v^{k}+\lambda(Ax^{k}-b)],\ \alpha\right),$ (53a) $\displaystyle
v^{k+1}$ $\displaystyle=v^{k}+\lambda(Ax^{k+1}-b).$ (53b)

This converges for $\alpha\lambda\|A^{\top}A\|<1$ (e.g. see [13, Section
3.5]). We used $\lambda=100\|A^{\top}A\|$ and
$\alpha=1/(\lambda\|A^{\top}A\|)$.

Primal Dual Hybrid Gradient (PDHG). For step sizes $\alpha,\lambda>0$, the
primal dual hybrid gradient (PDHG) iteration solves (BP) using the iteration

$\displaystyle x^{k+1}$ $\displaystyle=\mbox{shrink}(x^{k}-\alpha
A^{\top}v^{k},\ \alpha),$ (54a) $\displaystyle v^{k+1}$
$\displaystyle=v^{k}+\lambda[A(2x^{k+1}-x^{k})-b].$ (54b)

Here $\alpha\lambda\|A^{\top}A\|<1$ ensures convergence (see [13, Section
3.2]). We used $\lambda=100\|A^{\top}A\|$ and
$\alpha=1/(\lambda\|A^{\top}A\|)$.

###  B.2 Stable Principal Component Pursuit

For each method, we initialize the low rank term to $M$ and the sparse term to
the zero matrix (e.g. for proximal projection
$X^{1}=(X_{L}^{1},X_{S}^{1})=(M,0)$).

#### Proximal Projection (PP).

Letting $X=(X_{L},X_{S})$, the problem (SPCP) may be rewritten as

$\min_{X}\|X_{L}\|_{\star}+\lambda\|X_{S}\|_{1}\quad\mbox{s.t.}\quad\|AX-M\|_{F}\leq\varepsilon,$
(55)

where $A=[\mbox{I}\ \mbox{I}]$. The proximal for (SPCP) can be written as

$\mbox{prox}_{\alpha
f}(X)=\left[\begin{array}[]{c}\mbox{prox}_{\alpha\|\cdot\|_{\star}}(X_{L})\\\
\mbox{prox}_{\alpha\lambda\|\cdot\|_{1}}(X_{S})\end{array}\right]=\left[\begin{array}[]{c}\mbox{svt}(X_{L},\
\alpha)\\\ \mbox{shrink}_{1}(X_{S},\ \alpha\lambda)\end{array}\right].$ (56)

We next verify the projection formula used by PP for (SPCP). For ease of
reference, the result is restated before its proof.

Lemma 3.1 (SPCP Projection).  If
$\mathcal{C}\triangleq\\{X=[X_{L};X_{S}]:\|X_{L}+X_{S}-M\|_{F}\leq\varepsilon\\}$,
then

$\mathcal{P}_{\mathcal{C}}(Z)=\left[\begin{array}[]{c}Z_{L}-\mu_{Z}(Z_{L}+Z_{S}-M)\\\
Z_{S}-\mu_{Z}(Z_{L}+Z_{S}-M)\end{array}\right],$ (57)

where121212We adopt the convention of using $-\varepsilon/0=-\infty$.

$\mu_{Z}\triangleq\max\left\\{0,\
\frac{\|Z_{L}+Z_{S}-M\|_{F}-\varepsilon}{2\|Z_{L}+Z_{S}-M\|_{F}}\right\\}.$
(58)

###### Proof.

When $Z$ is feasible, $\mathcal{P}_{\mathcal{C}}(Z)=Z$. In this case,
$\mu_{Z}=0$ and so the result holds. For the remainder of the proof, we assume
$Z$ is not feasible. By Proposition 1, when $Z$ is not feasible, $\tau_{Z}$ is
the positive solution to

$\tau\left\|(AA^{\top}+\varepsilon\tau\mbox{I})^{-1}[AZ-M]\right\|_{F}=1.$
(59)

For each $\tau\geq 0$, note

$\left(AA^{\top}+\tau\mbox{I}\right)^{-1}=\left((2+\varepsilon\tau)\mbox{I}\right)^{-1}=\frac{\mbox{I}}{2+\varepsilon\tau}.$
(60)

Thus, (59) may be rewritten as

$\frac{\tau}{2+\varepsilon\tau}\|AZ-M\|_{F}=1.$ (61)

Rearranging to isolate $\tau$ yields

$\tau_{Z}=\frac{2}{\|AZ-M\|_{F}-\varepsilon}.$ (62)

Thus, if $\|AZ-M\|_{F}>\varepsilon$, then

$\displaystyle A^{\top}(AA^{\top}+\varepsilon\tau_{X}\mbox{I})^{-1}(AZ-M)$
$\displaystyle=\frac{A^{\top}(AZ-M)}{2+\varepsilon\tau_{Z}}$ (63a)
$\displaystyle=\frac{\|AZ-M\|_{F}-\varepsilon}{2\|AZ-M\|_{F}}A^{\top}(AZ-M)$
(63b)
$\displaystyle=\frac{\|AZ-M\|_{F}-\varepsilon}{2\|AZ-M\|_{F}}\left[\begin{array}[]{c}Z_{L}+Z_{S}-M\\\
Z_{L}+Z_{S}-M\end{array}\right]$ (63e)
$\displaystyle=\frac{\|Z_{L}+Z_{S}-M\|_{F}-\varepsilon}{2\|Z_{L}+Z_{S}-M\|_{F}}\left[\begin{array}[]{c}Z_{L}+Z_{S}-M\\\
Z_{L}+Z_{S}-M\end{array}\right].$ (63h)

Hence

$\displaystyle\mathcal{P}_{\mathcal{C}}(X)$
$\displaystyle=\begin{cases}\begin{array}[]{cl}X-\frac{\|X_{L}+X_{S}-M\|_{F}-\varepsilon}{2\|X_{L}+X_{S}-M\|_{F}}\left[\begin{array}[]{c}X_{L}+X_{S}-M\\\
X_{L}+X_{S}-M\end{array}\right]&\mbox{if}\
\|X_{L}+X_{S}-M\|_{F}>\varepsilon\\\ X&\mbox{otherwise}\end{array}\end{cases}$
(64a)
$\displaystyle=X-\max\left\\{\frac{\|X_{L}+X_{S}-M\|_{F}-\varepsilon}{2\|X_{L}+X_{S}-M\|_{F}},0\right\\}\left[\begin{array}[]{c}X_{L}+X_{S}-M\\\
X_{L}+X_{S}-M\end{array}\right],$ (64d)

and the proof is complete. ∎

Substituting the proximal formula (56) and projection formula from Lemma 3.1
into (5) yields the updates

$\displaystyle Z_{L}^{k+1}$
$\displaystyle=\mbox{svt}(Z_{L}^{k}-2\mu_{k}(Z_{L}^{k}+Z_{S}^{k}-M),\ \alpha)$
(65a) $\displaystyle Z_{S}^{k+1}$
$\displaystyle=\mbox{shrink}_{1}(Z_{S}^{k}-2\mu_{k}(Z_{L}^{k}+Z_{S}^{k}-M),\
\lambda\alpha)$ (65b)

Thus, Algorithm 2 for (SPCP) simplifes to Algorithm 3. Furthermore, in the
case where $\varepsilon=0$, the sequence $\\{Z^{k}\\}$ is the same as the
sequence $\\{X^{k}\\}$ in (79) below.

#### Variant of Alternating Splitting Augmented Lagrangian Method (VASALM).

For a parameter $\alpha>0$, the augmented Lagrangian $\mathcal{L}$ for (SPCP)
is

$\mathcal{L}_{\alpha}(L,S,N;\Lambda)=\|L\|_{\star}+\lambda\|S\|_{1}+\delta_{B(0,\varepsilon)}(N)+\left<\Lambda,L+S+N-M\right>+\frac{\alpha}{2}\|L+S+N-M\|_{F}^{2}.$
(66)

Note

$\displaystyle\frac{\alpha}{2}\left\|L+S+N-M+\frac{\Lambda}{\alpha}\right\|_{F}^{2}$
$\displaystyle=\frac{1}{2}\left\|\sqrt{\alpha}(L+S+N-M)+\frac{\Lambda}{\sqrt{\alpha}}\right\|_{F}^{2}$
(67)
$\displaystyle=\frac{\alpha}{2}\|L+S+N-M\|_{F}^{2}+\left<\frac{\Lambda}{\sqrt{\alpha}},\sqrt{\alpha}(L+S+N-M)\right>+\frac{1}{2\alpha}\|\Lambda\|_{F}^{2}$
(68)
$\displaystyle=\frac{\alpha}{2}\|L+S+N-M\|_{F}^{2}+\left<\Lambda,L+S+N-M\right>+\frac{1}{2\alpha}\|\Lambda\|_{F}^{2},$
(69)

and so

$\mathcal{L}_{\alpha}(L,S,N;\Lambda)=\|L\|_{\star}+\lambda\|S\|_{1}+\delta_{B(0,\varepsilon)}(N)+\frac{\alpha}{2}\|L+S+N-M+\alpha\Lambda\|_{F}^{2}-\frac{1}{2\alpha}\|\Lambda\|_{F}^{2}.$

VASALM does proximal steps for $L$, $S$, and $N$ separately, with two dual
variable $\Lambda$ updates. Specifically, for $\alpha>0$ and $\eta>2$ it
generates a sequence of updates via

$\displaystyle N^{k+1}$
$\displaystyle=\mathcal{P}_{B(0,\varepsilon)}\left(\frac{\Lambda^{k}}{\alpha}+M-X_{L}^{k}-X_{S}^{k}\right)$
(70a) $\displaystyle\hat{\Lambda}^{k}$
$\displaystyle=\Lambda^{k}-\alpha(X_{L}^{k}+X_{S}^{k}+N^{k+1}-M)$ (70b)
$\displaystyle X_{S}^{k+1}$
$\displaystyle=\mbox{shrink}\left(X_{S}^{k}+\frac{\hat{\Lambda}^{k}}{\alpha\eta},\
\frac{\lambda}{\alpha\eta}\right)$ (70c) $\displaystyle X_{L}^{k+1}$
$\displaystyle=\mbox{svt}\left(X_{L}^{k}+\frac{\hat{\Lambda}^{k}}{\alpha\eta},\
\frac{1}{\alpha\eta}\right)$ (70d) $\displaystyle\Lambda^{k+1}$
$\displaystyle=\hat{\Lambda}^{k}+\alpha(X_{L}^{k}-X_{L}^{k+1})+\alpha(X_{S}^{k}-X_{S}^{k+1}).$
(70e)

We used $\eta=3$ and $\alpha=10^{-5}$.

#### Partially Smoothed Proximal Gradient.

Following the Nesterov smoothing technique [54], for a parameter $\mu>0$, the
work [9] considers the smoothed the version of (SPCP) given by

$\min_{L,S}\lambda\|S\|_{1}+\left(\max_{\|W\|\leq
1}\left<L,W\right>-\frac{\mu}{2}\|W\|_{F}^{2}\right)\quad\mbox{s.t.}\quad\|L+S-M\|_{F}\leq\varepsilon,$
(SPCPμ)

which approaches (SPCP) as $\mu\rightarrow 0^{+}$. The nuclear norm
approximation $1/\mu$-smooth and has gradient given by

$W_{\mu}(L)=U\,\mbox{diag}\left(\min\left\\{\frac{\sigma_{i}}{\mu},\
1\right\\}\right)V^{\top},$ (71)

where $U\Sigma V^{\top}$ is here the SVD of $L$. Proximal gradient updates
take the form

$X^{k+1}=\operatorname*{\mbox{argmin}}_{X}\lambda\|S\|_{1}+\left<W_{\mu}(X_{L}^{k}),\,X_{L}-X_{L}^{k}\right>+\frac{1}{2\mu}\|X_{L}-X_{L}^{k}\|_{F}^{2}\quad\mbox{s.t.}\quad\|X_{L}+X_{S}-M\|_{F}\leq\varepsilon.$
(72)

Following [9, Lemma 6.1], here proximal gradient update steps are explicitly
given by

$\displaystyle X_{S}^{k+1}$
$\displaystyle=\mbox{shrink}\left(M-\mu(X_{L}^{k}-W_{\mu}(X_{L}^{k})),\frac{\lambda[1+\mu\theta^{\star}]}{\theta^{\star}}\right)$
(73a) $\displaystyle X_{L}^{k+1}$
$\displaystyle=\frac{1}{1+\mu\theta^{\star}}\Big{[}\mu\theta^{\star}\left(M-X_{S}^{k+1}\right)+\big{(}X_{L}^{k}-\mu
W_{\mu}(X_{L}^{k})\big{)}\Big{]},$ (73b)

where $\theta^{\star}$ is the unique positive solution to131313Here we assume
$\varepsilon>0$.

$\varepsilon=\left\|\min\left\\{\frac{\lambda}{\theta},\
\frac{|M-X_{L}^{k}+\mu
W_{\mu}(X_{L}^{k})|}{1+\mu\theta}\right\\}\right\|_{F}.$ (74)

Since $L-\mu W_{\mu}(L)=\mbox{svt}(L)$, the proximal gradient updates can be
rewritten as

$\displaystyle X_{S}^{k+1}$
$\displaystyle=\mbox{shrink}\left(M-\mbox{svt}(X_{L}^{k},\
\mu),\frac{\lambda[1+\mu\theta^{\star}]}{\theta^{\star}}\right)$ (75a)
$\displaystyle X_{L}^{k+1}$
$\displaystyle=\frac{1}{1+\mu\theta^{\star}}\Big{[}\mu\theta^{\star}\left(M-X_{S}^{k+1}\right)+\mbox{svt}(X_{L}^{k},\
\mu)\Big{]}$ (75b)

with $\theta^{\star}$ the solution to

$\varepsilon=\left\|\min\left\\{\frac{\lambda}{\theta},\
\frac{|M-\mbox{svt}(X_{L}^{k},\ \mu)|}{1+\mu\theta}\right\\}\right\|_{F}.$
(76)

Since $M\in\mathbb{R}^{n_{1}\times n_{2}}$, we have the bound

$\theta^{\star}\leq\min\left\\{n_{1}n_{2}\lambda\varepsilon,\
\left|\frac{\|M-\mbox{svt}(X_{L}^{k},\mu)\|_{F}-\varepsilon}{\mu\varepsilon}\right|\right\\}.$
(77)

By the comment following [9, Theorem 2.1], setting
$\mu=\delta/\min\\{n_{1},n_{2}\\}$ ensures a $\delta/2$-optimal solution to
(SPCPμ) is $\delta$-optimal solution to (SPCP). We set $\delta$ to be about
$0.1$ times the optimal value for (SPCP) to ensure the objective of the limit
is within $\sim$10% of optimal. In Figure 3b, it appears that this choice of
$\mu$ leaves a visible gap between the limit of PSPG and the optimal value.
Reducing $\mu$ reduces the size of this gap, but also hinders the convergence
rate of PSPG.

#### Proximal Gradient (PG).

In the original work [5] on SPCP, the authors follow the example of [55] to
approximate (SPCP) in their numerical experiments by a soft-penalty variation

$\min_{X}\|X_{L}\|_{\star}+\lambda\|X_{S}\|_{1}+\frac{1}{2\mu}\|X_{L}+X_{S}-M\|_{F}^{2}.$
(78)

To apply proximal gradient, here we do the same. The update is

$X^{k+1}=\mbox{prox}_{\alpha
f}\left(X^{k}-\frac{\alpha}{\mu}\left[\begin{array}[]{c}X_{L}^{k}+X_{S}^{k}-M\\\
X_{L}^{k}+X_{S}^{k}-M\end{array}\right]\right),$ (79)

where the proximal is given in (56). Here the quadratic term is
$1/\mu$-smooth, and so a stepsize of $\alpha=\mu$ can be used to ensure (79)
converges to a solution to (78). In this case, the iteration simplifies to

$\displaystyle X_{L}^{k+1}$
$\displaystyle=\mbox{svt}\left(M-X_{S}^{k},\alpha\right)$ (80a) $\displaystyle
X_{S}^{k+1}$
$\displaystyle=\mbox{shrink}\left(M-X_{L}^{k},\alpha\lambda\right).$ (80b)

###  B.3 Earth Mover’s Distance

For simplicity, in the EMD example, we use $h=1.0$. The proof for the EMD
projection is given below. We emphasize, in this subsection, lowercase
variables are represented in matrix form. Moreover, although the divergence
operator (denoted by $A$) is linear, its application in this form involves
application of left and right matrix multiplications. To keep notation
concise, we will write $A$ on the leftmost side with its application
understood to be as described in the main text.

Lemma 3.2 (EMD Projection).  If
$\mathcal{C}\triangleq\\{m:\|\mbox{div}(m)+\rho^{1}-\rho^{0}\|_{F}\leq\varepsilon\\}$
and $U\Sigma V^{\top}$ is the SVD for $K$, then

$\mathcal{P}_{\mathcal{C}}(z)=\begin{cases}\begin{array}[]{cl}z&\mbox{if}\
\|\mbox{div}(z)+\rho^{1}-\rho^{0}\|_{F}\leq\varepsilon\\\
z-K^{\top}U\left[\begin{array}[]{c}Y_{\tau}\\\
Y_{\tau}^{\top}\end{array}\right]U^{\top}&\mbox{otherwise,}\end{array}\end{cases}$
(81)

where $\tau>0$ satisfies $1=\tau\|(Y_{\tau})\|_{F}$ and $Y_{\tau}$ is the
matrix defined element-wise via

$(Y_{\tau})_{ij}\triangleq\frac{[U^{\top}(Kz+zK^{\top}+\rho^{1}-\rho^{0})U]_{ij}}{\sigma_{i}^{2}+\sigma_{j}^{2}+\varepsilon\tau},\quad\mbox{for
all}\ i,j\in[n].$ (82)

###### Proof.

This proof is a corollary of Proposition 1. For $z\in\mathcal{C}$, the result
directly follows from the proposition. In the remainder of the proof, we
assume $z\notin\mathcal{C}$. The two tasks at at hand are to obtain a formula
for the term multiplied by the matrix inverse in Propostion 1 and to verify
the choice for $\tau$ matches that in Proposition 1. Set

$q\triangleq[AA^{\top}+\varepsilon\tau\mbox{I}]^{-1}(Am+\rho^{1}-\rho^{0})$
(83)

so that, by Proposition 1 and (14),

$\mathcal{P}_{\mathcal{C}}(z)=z-A^{\top}q=z-\left[\begin{array}[]{c}K^{\top}q^{1}\\\
K^{\top}(q^{2})^{\top}\end{array}\right].$ (84)

It follows that

$[AA^{\top}+\varepsilon\tau\mbox{I}]q=Am+\rho^{1}-\rho^{0}.$ (85)

By the formulas in (14),

$[AA^{\top}+\varepsilon\tau\mbox{I}]q=\left[KK^{\top}+\frac{\varepsilon\tau}{2}\mbox{I}\right]q+q\left[KK^{\top}+\frac{\varepsilon\tau}{2}\mbox{I}\right].$
(86)

For $U\Sigma V^{\top}$ the SVD of $K$, direct multiplication reveals

$KK^{\top}+\frac{\varepsilon\tau}{2}\mbox{I}=U\mbox{diag}(\sigma_{i}^{2})U^{\top}+\frac{\varepsilon\tau}{2}\mbox{I}=U\,\mbox{diag}\left(\sigma_{i}^{2}+\frac{\varepsilon\tau}{2}\right)U^{\top},$
(87)

where the final equality holds by the orthogonality of $U$. Set $Y=U^{\top}qU$
so that

$U\,\mbox{diag}\left(\sigma_{i}^{2}+\frac{\varepsilon\tau}{2}\right)YU^{\top}+UY\,\mbox{diag}\left(\sigma_{i}^{2}+\frac{\varepsilon\tau}{2}\right)U^{\top}=[AA^{\top}+\varepsilon\tau\mbox{I}]q=Am+\rho^{1}-\rho^{0}.$
(88)

Left multiplying each term by $U^{\top}$ and then right multiplying by $U$,
(88) becomes

$\mbox{diag}\left(\sigma_{i}^{2}+\frac{\varepsilon\tau}{2}\right)Y+Y\mbox{diag}\left(\sigma_{i}^{2}+\frac{\varepsilon\tau}{2}\right)=U^{\top}(Am+\rho^{1}-\rho^{0})U.$
(89)

In element-wise form, (89) may be equivalently written as

$\left(\sigma_{i}^{2}+\frac{\varepsilon\tau}{2}\right)Y_{ij}+Y_{ij}\left(\sigma_{j}^{2}+\frac{\varepsilon\tau}{2}\right)=(U^{\top}(Am+\rho^{1}-\rho^{0})U)_{ij},\quad\mbox{for
all}\ i,j\in[n].$ (90)

Thus,

$Y_{ij}=\frac{(U^{\top}(Am+\rho^{1}-\rho^{0})U)_{ij}}{\sigma_{i}^{2}+\sigma_{j}^{2}+\varepsilon\tau},$
(91)

where the division is well-defined since (C2) ensures
$\sigma_{i}^{2}+\sigma_{j}^{2}+\varepsilon\tau>0$ for all $\tau>0$. This
verifies the formula for $Y_{\tau}$ in the lemma statement. Lastly, note the
orthogonality of $U$ ensures

$\|q\|_{F}=\|UY_{\tau}U^{\top}\|_{F}=\|Y_{\tau}U^{\top}\|_{F}=\|Y_{\tau}\|_{F},$
(92)

which verifies the condition on $\tau$ may be expressed as
$1=\tau\|Y_{\tau}\|_{F}$. ∎

#### Proximal Projection (PP).

We used Algorithm 4 with $\alpha=10^{-4}$.

#### Primal Dual Hybrid Gradient (PDHG).

We use the iteration in (54) with $\lambda=5$ and $\alpha=1/(5\|A^{\top}A\|)$.

#### G-Prox PDHG.

Using the Hodge decomposition, the flux $m$ can be decomposed as
$m=u+\nabla\psi$, where $u$ is a divergence free vector field and $\nabla\psi$
is a gradient field for which

$\mbox{div}(\nabla\psi)+\rho^{1}-\rho^{0}=0.$ (93)

In this example, we used

$\nabla\psi=\mathcal{P}_{\mathcal{C}}(0).$ (94)

With $m=u+\nabla\psi$, the EMD problem can be rewritten as

$\min_{u}\|u+\nabla\psi\|_{1}\quad\mbox{s.t.}\quad\mbox{div}(u)=0.$ (95)

To solve this, for parameters $\sigma>0$ and $\tau>0$, [40] proposes the
iteration141414Here we compress (3.14), (3.15), and (3.16) from that work.

$\displaystyle p^{n+1}$
$\displaystyle=\frac{p^{n}+\sigma[u^{n}+\nabla\psi]}{\max\\{1,|p^{n}+\sigma[u^{n}+\nabla\psi|]\\}}$
(96a) $\displaystyle u^{n+1}$
$\displaystyle=u^{n}-\tau\mathcal{P}_{\nabla^{\perp}}(2p^{n+1}-p^{n}),$ (96b)

where $\nabla^{\perp}$ is the set of divergence free vector fields. To compute
this projection, we use Lemma 3.2, but with
$\mathcal{C}=\\{m:\|\mbox{div}(m)\|\leq\varepsilon\\}$. This iteration is
guaranteed to converge for $\sigma\tau<1$; we used $\tau=10^{-4}$ and
$\sigma=10^{4}$.

###  B.4 Stable Matrix Completion

We begin with notation. In vectorized form of $X$, the matrix $A$ for (SMC)
would be a diagonal matrix with $mn$ rows and $mn$ columns. The $ij$-th entry
on the diagonal would be $1$ if $(i,j)\in\Omega$ and $0$ otherwise.
Consequently, $A=A^{\top}$ and $A=\mathcal{P}_{\Omega}$. With this choice for
$A$, we may prove the projection formula.

Lemma 3.3 (SMC Projection).  If
$\mathcal{C}\triangleq\\{X:\|\mathcal{P}_{\Omega}(X-M)\|_{F}\leq\varepsilon\\}$,
then

$\mathcal{P}_{\mathcal{C}}(Z)=\begin{cases}\begin{array}[]{cl}Z&\mbox{if}\
\|\mathcal{P}_{\Omega}(Z-M)\|_{F}\leq\varepsilon\\\
\mathcal{P}_{\Omega}\left(\frac{[\|\mathcal{P}_{\Omega}(Z-M)\|_{F}-\varepsilon]M+\varepsilon
Z}{\|\mathcal{P}_{\Omega}(Z-M)\|_{F}}\right)+\mathcal{P}_{\Omega^{\perp}}(Z)&\mbox{otherwise,}\end{array}\end{cases}$
(97)

###### Proof.

Consider a matrix $X$. By Proposition 1, if $X\in\mathcal{C}$, then
$\mathcal{P}_{\mathcal{C}}(X)=X$. In what remains, we assume
$X\notin\mathcal{C}$. By Proposition 1, the projection is given by

$\mathcal{P}_{\mathcal{C}}(X)=X-A^{\top}(AA^{\top}+\varepsilon\tau_{X}\mbox{I})(AX-\mathcal{P}_{\Omega}(M)).=X-\mathcal{P}_{\Omega}\left(\left[\mathcal{P}_{\Omega}+\varepsilon\tau_{X}\mbox{I}\right]^{-1}\left[\mathcal{P}_{\Omega}(X-M)\right]\right),$
(98)

where again we note $A=A^{\top}=\mathcal{P}_{\Omega}$ and
$\mathcal{P}_{\Omega}^{2}=\mathcal{P}_{\Omega}$. Set

$Y\triangleq(\mathcal{P}_{\Omega}+\varepsilon\tau_{X}\mbox{I})^{-1}\mathcal{P}_{\Omega}(X-M).$
(99)

This implies

$\mathcal{P}_{\Omega}(Y)+\varepsilon\tau_{X}Y=\mathcal{P}_{\Omega}(X-M)\quad\iff\quad
Y_{ij}=\begin{cases}\begin{array}[]{cl}\frac{X_{ij}-M_{ij}}{1+\varepsilon\tau_{X}}&\mbox{if}\
(i,j)\in\Omega,\\\ 0&\mbox{otherwise},\end{array}\end{cases}\quad\mbox{for
all}\ (i,j)\in[m]\times[n].$ (100)

Consequently,

$\mathcal{P}_{\Omega}(X)=X-\mathcal{P}_{\Omega}(Y)=X-\mathcal{P}_{\Omega}\left(\frac{X-M}{1+\varepsilon\tau_{X}}\right)=\mathcal{P}_{\Omega^{\perp}}(X)+\mathcal{P}_{\Omega}\left(\frac{M+\varepsilon\tau_{X}X}{1+\varepsilon\tau_{X}}\right).$
(101)

By Proposition 1, $\tau_{X}$ is the unique positive scalar satisfying

$1=\tau_{X}\|Y\|_{F}=\tau_{X}\left\|\mathcal{P}_{\Omega}\left(\frac{X-M}{1+\varepsilon\tau_{X}}\right)\right\|_{F}=\frac{\tau_{X}}{1+\varepsilon\tau_{X}}\|\mathcal{P}_{\Omega}(X-M)\|_{F},$
(102)

and so rearranging reveals

$\tau_{X}=\frac{1}{\|\mathcal{P}_{\Omega}(X-M)\|_{F}-\varepsilon}.$ (103)

Plugging this choice for $\tau_{X}$ into the projection formula (101) reveals

$\mathcal{P}_{\Omega}(X)=\mathcal{P}_{\Omega^{\perp}}(Z)+\mathcal{P}_{\Omega}\left(\frac{[\|\mathcal{P}_{\Omega}(Z-M)\|_{F}-\varepsilon]M+\varepsilon
Z}{\|\mathcal{P}_{\Omega}(Z-M)\|_{F}}\right)$ (104)

as desired. ∎

Smoothed Proximal Gradient (SPG). Following [9], for a parameter $\mu>0$ we
consider the smoothed problem

$\min_{X}\left(\max_{\|W\|\leq
1}\left<X,W\right>-\frac{\mu}{2}\|W\|_{F}^{2}\right)\quad\mbox{s.t.}\quad\|\mathcal{P}_{\Omega}(X-M)\|_{F}\leq\varepsilon,$
(SMCμ)

In much the same fashion as the stable principal component pursuit problem in
Appendix B.2, here the smooth proximal gradient updates can be rewritten as

$\displaystyle
X^{k+1}=\mathcal{P}_{\Omega}\left(\frac{1}{1+\mu\theta^{\star}}\Big{[}\mu\theta^{\star}M+\mbox{svt}(X^{k},\
\mu)\Big{]}\right)+\mathcal{P}_{\Omega^{\perp}}\left(\mbox{svt}(X^{k},\
\mu)\right),$ (105)

with

$\varepsilon=\left\|\frac{\mathcal{P}_{\Omega}(M-\mbox{svt}(X^{k},\
\mu))}{1+\mu\theta^{\star}}\right\|_{F}\quad\iff\quad\theta^{\star}=\frac{\|\mathcal{P}_{\Omega}(M-\mbox{svt}(X^{k},\mu))\|_{F}-\varepsilon}{\mu\varepsilon}$
(106)

Here we use $\mu=10$. As noted before, reducing $\mu$ reduces the sub-
optimality gap of the limit of SPG, but reduces the smoothness and, thus,
step-size (which hinders convergence rate).

#### Variant of Alternating Splitting Augmented Lagrangian Method (VASALM).

Similar to (SPCP), for a parameter $\alpha>0$, the augmented Lagrangian
$\mathcal{L}$ for (SMC) is

$\mathcal{L}_{\alpha}(X,N;\Lambda)=\|X\|_{\star}+\delta_{B(0,\varepsilon)}(\mathcal{P}_{\Omega}(N))+\left<\Lambda,\mathcal{P}_{\Omega}(X+N-M)\right>+\frac{\alpha}{2}\|\mathcal{P}_{\Omega}(X+N-M)\|_{F}^{2}.$
(107)

Define the set

$\mathcal{K}\triangleq\\{Y:\|\mathcal{P}_{\Omega}(Y)\|_{F}\leq\varepsilon\\}.$
(108)

The projection $\mathcal{P}_{\mathcal{K}}$ onto $\mathcal{K}$ is simply a
Euclidean projection onto the $\varepsilon$-ball about the origin for the
submatrix $Y_{\Omega}$ and the remainder of $Y$ is unchanged by the
projection. The VASALM algorithm in this context uses the iterates

$\displaystyle N^{k+1}$
$\displaystyle=P_{\mathcal{K}}\left(\frac{\Lambda^{k}}{\alpha}+\mathcal{P}_{\Omega}(M)-X\right)$
(109a) $\displaystyle\hat{\Lambda}^{k}$
$\displaystyle=\Lambda^{k}-\alpha\left(X^{k}+N^{k}-\mathcal{P}_{\Omega}(M^{k})\right)$
(109b) $\displaystyle X^{k+1}$
$\displaystyle=\mbox{svt}\left(X^{k}+\frac{\hat{\Lambda}^{k}}{\eta\alpha},\
\frac{1}{\eta\alpha}\right)$ (109c) $\displaystyle\Lambda^{k+1}$
$\displaystyle=\hat{\Lambda}^{k}+\alpha(X^{k}-X^{k+1}),$ (109d)

where $\eta>2$. We use $\eta=3$ and $\alpha=10^{-2}$.