The accelerated proximal point method

Problem setup

Consider the monotone inclusion

\[\text{find } x \in \calH \text{ such that } 0 \in G(x),\]

where \(G:\calH\rightrightarrows\calH\) is maximally monotone.

For initial points \(x^0,y^0,y^{-1}\in\calH\), step size \(\gamma\in\reals_{++}\), and \(\lambda_k = k/(k+2)\), the accelerated proximal point method [Kim21] is

\[\begin{split}(\forall k \in \naturals)\quad \left[ \begin{aligned} x^{k+1} &= J_{\gamma G}(y^k), \\ y^{k+1} &= x^{k+1} + \lambda_k(x^{k+1} - x^k) - \lambda_k(x^k - y^{k-1}). \end{aligned} \right.\end{split}\]

In this example, we use AutoLyap to compute the smallest certificate constant \(c_K\) such that

\[\|x^{K+1} - y^K\|^2 \le c_K \|y^0 - y^\star\|^2,\]

where \(K \in \mathbb{N}\) is the iteration budget (horizon), and \(y^\star \in \zer G\).

Model the problem in AutoLyap and search for the smallest c_K

• \(G\) is modeled by MaximallyMonotone.

• The monotone inclusion is represented by InclusionProblem.

• The update rule is represented by AcceleratedProximalPoint.

• Initial-horizon parameters \(\|y^0-y^\star\|^2\) are obtained with IterationDependent.get_parameters_state_component_distance_to_solution.

• Final-horizon parameters \(\|x^{K+1}-y^K\|^2\) are obtained with IterationDependent.get_parameters_state_component_cross_iteration_difference.

• The certificate constant is computed with IterationDependent.search_lyapunov.

Run the iteration-dependent analysis

This example uses the MOSEK Fusion backend (backend="mosek_fusion"). Install the optional MOSEK dependency first:

pip install "autolyap[mosek]"

from autolyap import IterationDependent, SolverOptions
from autolyap.algorithms import AcceleratedProximalPoint
from autolyap.problemclass import InclusionProblem, MaximallyMonotone

gamma = 1.0
K = 10

problem = InclusionProblem([MaximallyMonotone()])
algorithm = AcceleratedProximalPoint(gamma=gamma, type="operator")
solver_options = SolverOptions(backend="mosek_fusion")

# V(0) = ||y^0 - y^*||^2. In APP state x^k = (x^k, y^k, y^{k-1}),
# y^k is state component ell=2 (1-indexed).
Q_0 = IterationDependent.get_parameters_state_component_distance_to_solution(
    algorithm,
    k=0,
    ell=2,
)

# V(K) = ||x^{K+1} - y^K||^2 = ||x_1^{K+1} - x_2^K||^2
# with 1-indexed state components: ell=1 and ell_prime=2.
Q_K = IterationDependent.get_parameters_state_component_cross_iteration_difference(
    algorithm,
    k=K,
    ell=1,
    ell_prime=2,
)

result = IterationDependent.search_lyapunov(
    problem,
    algorithm,
    K,
    Q_0,
    Q_K,
    solver_options=solver_options,
)
if result["status"] != "feasible":
    raise RuntimeError("No feasible chained Lyapunov certificate for APP.")

c_K_autolyap = result["c_K"]
c_K_kim = 1.0 / (K + 1.0) ** 2

print(f"APP c_K (AutoLyap): {c_K_autolyap:.6e}")
print(f"APP c_K (Kim):      {c_K_kim:.6e}")

The computed value APP c_K (AutoLyap) matches (up to solver numerical tolerances) the theoretical horizon-K expression from [Kim21, Theorem 4.1], i.e.,

\[\|x^{K+1} - y^K\|^2 \le \frac{1}{(K+1)^2}\|y^0 - y^\star\|^2.\]

Sweeping over \(K \in \llbracket 1, 100\rrbracket\) gives the plot below, with Kim’s theoretical constant from [Kim21, Theorem 4.1] \(c_K^{\mathrm{Kim}} = \frac{1}{(K+1)^2}\) as a line and AutoLyap certificates as blue dots.

References

[Kim21] (1,2,3)

Donghwan Kim. Accelerated proximal point method for maximally monotone operators. Mathematical Programming, 190(1-2):57–87, March 2021. doi:10.1007/s10107-021-01643-0.