Proximal point algorithms have found numerous applications in the field of convex optimization, and their accelerated forms have also been proposed. However, the most commonly used accelerated proximal point algorithm was first introduced in 1967, and recent studies on accelerating proximal point algorithms are relatively scarce. In this paper, we propose high-resolution ODEs for the proximal point operators for both closed proper convex functions and maximally monotone operators, and present a Lyapunov function framework to demonstrate that the trajectories of our high-resolution ODEs exhibit accelerated behavior. Subsequently, by symplectically discretizing our high-resolution ODEs, we obtain new proximal point algorithms known as symplectic proximal point algorithms. By decomposing the continuous-time Lyapunov function into its elementary components, we demonstrate that symplectic proximal point algorithms possess $O(1/k^2)$ convergence rates.

§ INTRODUCTION

The burgeoning development of high-dimensional statistics and deep learning has led to an exponential surge in the demand for sparse solutions. Nevertheless, when confronted with sparse optimization problems, the absence of a smooth gradient often renders direct application of gradient descent methods untenable. In response to this challenge, numerous optimization algorithms that rely on solving optimization sub-problems have garnered significant attention from researchers. Nonetheless, the complexity of these methods is prohibitively high for use with high-dimensional problems. Therefore, this article seeks to elucidate how to accelerate these methods from the vantage point of ordinary differential equations.

In this paper, we consider the following zero-point problem
\begin{equation}
\label{eq:main}
0\in A(x),
\end{equation}

where $A$ is a maximally monotone operators.

In [21], Rockafellar introduced the use of proximal point algorithms to solve the zeros of maximally monotone operators. Over time, proximal point algorithms have found a wide range of applications in optimization problems, such as nonsmooth optimization problems[21], augmented Lagrangian methods(ALM)[19], Douglas-Rachford algorithms(DR) and alternating direction method of multipliers(ADMM)[11]. For a more comprehensive introduction to theory and application of proximal point algorithms, we refer to[10, 18]. The recursive rule of proximal point algorithms can be described as the following fixed point iteration method
\begin{equation}
\end{equation}

Here $c>0, J_{cA}=(I+cA)^{-1}$ is the resolvent of $A$ or proximal point operator of $A$. It has been proved that the worst-case of convergence rate of proximal point algorithms is $O(1/k)$ . Owing to Nesterov's accelerated gradient methods, which exhibit improved convergence rates without augmenting computational complexity per iteration, the accelerated form of proximal point algorithms has garnered significant interest in both theory and practice. Using the Halpern's iteration[13], the following recursive rule is the accelerated form of proximal point algorithms.

\begin{equation}
\label{eq:acceleratedproximal1}
\begin{split}
\end{split}
\end{equation}

In comparison to $O(1/k)$ convergence rate of proximal point algorithms, the accelerated form of proximal point algorithms yields $O(1/k^2)$ convergence rate.

However, studies of accelerated proximal point algorithms from a theoretical point of view are infrequently. In [14], Kim interpreted worst-cast convergence rate of general proximal point operator as optimal value of a optimization problem, called performance estimation problem(PEP). Although PEP framework provides a unify framework to analysis convergence rats of general proximal point algorithms and has been successfully applied to analyzing (<ref>), PEP shows limitations to obtain new accelerated form of proximal point algorithms. After Su established the link between Nesterov accelerated gradient method and ordinary differential equations(ODEs) in [25], the study of optimization algorithms from ODEs perspective began to increased. In [4, 3], Attouch innovatively considered Yosida approximation of $A$ as gradient of proximal point algorithms and proposed his ODEs for proximal point operators basic on ODEs in [1, 25]. After proper discretization, Attouch obtain new accelerated proximal point algorithms. However the index of accelerated proximal point algorithms proposed by Attouch is the same order as $k^2$, which limits its applications in accelerating Douglas-Rachford algorithms, ADMM, etc.

In our work, we utilize the ODEs in [26], which is originated from the high-resolution ODEs proposed in [23], to propose ODEs for proximal point operators. The reason that we use high-resolution ODEs techinique is that the high-resolution ODEs in [23] prove refining convergence results of Nesterov's accelerated gradient method and has many applications such as analyzing FISTA[15] and algorithms for solving saddle point problems[8]. Also by investigating the structure of our ODEs, we utlize implicit symplectic methods to discretize our ODEs and obtain new proximal point algorithms.

§.§ Structure of this paper

* In section <ref>, we establish a connection between proximal point algorithms and ordinary differential equations for both closed proper convex functions and maximally monotone operators. Furthermore, we propose corresponding Lyapunov functions to analyze the convergence rate of proximal point algorithms for both cases. These results offer valuable insights into extending the high-resolution ODEs framework to accelerate proximal point algorithms.

* In section <ref>, we generalize the high-resolution ODEs in [26], and consider the following ODEs
\begin{equation}
\begin{split}
Z&=b_t\X+c_t\nabla f(X)+X,\\
\Z&=-a_t\nabla f(X).
\end{split}
\label{eq:ode1}
\end{equation}

Also, we extend the Lyapunov function in [my previous work, hasn't been published yet] to prove the convergence rate of (<ref>). Next, we take advantage of symplectic method to discretize (<ref>) and obtain a new accelerated proximal point, called symplectic proximal point algorithms. Inspired by the work in [21], we applied the symplectic proximal point algorithms to Lagrangian duality and obtain symplectic augmented Lagrangian methods.

* In section <ref>, we first build up the high-resolution ODEs for (<ref>), which inspire us to consider the following ODEs
\begin{equation}\label{eq:ode2}
\begin{split}
\Z&=-\frac{C}{r}\tilde{A}(X).
\end{split}
\end{equation}

Above $\tilde{A}(X)=X-J_A(X)$ is the Yosida-approximation of $A$. Basic on the convergence result discussed in section <ref>, we proposed the corresponding Lyapunov function, which shows the accelerated phenomenon of (<ref>). Also, by utilizing Lyapunov function, we shows symplectic discretization of (<ref>) can achieve new accelerated proximal point algorithms, also called symplectic proximal point algorithms.

* In section <ref>, we apply symplectic proximal point algorithms to several famous instances of proximal point algorithms, and provide discussion on numerical performance of symplectic proximal point algorithms.

§ CONNECTION BETWEEN PROXIMAL POINT ALGORITHMS AND ORDINARY DIFFERENTIAL EQUATIONS

§.§ Proximal Point Algorithms for Closed Proper Convex Functions

First, we study the convergence rate of proximal point algorithms when $A$ is the subgradient of a closed proper convex function $f$. When $A=\partial f$, the proximal point algorithms is equivalent to recursively solving the following optimization subproblems
\begin{equation}
\end{equation}

In practice, the parameter $c$ will change with iteration. Therefore, we consider the following iteration rule
\begin{equation}\label{eq:proximalf}
\end{equation}

If $f$ is continuous differentiable, the first-order equivalent characterization of (<ref>) is
\begin{equation}
\frac{x_{k+1}-x_k}{c_k}=-\nabla f(x_{k+1}).
\label{eq:proximalffirst}
\end{equation}

(<ref>) can be seen as apply implicit Euler method to standard gradient flow $\X=-\nabla f(X)$. It is has been shown in [24] that the convergence rate of standard gradient flow is $O(1/t)$. However, The effect of the parameter $\{c_k\}$ cannot be reflected if we use $\X=-\nabla f(X)$ to analyze (<ref>). Such problem can be overcome by considering

\begin{equation}
\label{eq:gradientflow}
\X=-c_t\nabla f(X).
\end{equation}

To analyze trajectories of (<ref>), we propose the following Lyapunov function
\begin{equation}
\E(t)=\int_{0}^{t}c_udu[f(X(t))-f(x^*)]+\frac{1}{2}\norm{X(t)-x^*}^2.
\end{equation}

Here $\int_0^t c_udu[f(X(t))-f(x^*)]$ is potential term and $\dfrac{1}{2}\norm{X(t)-x^*}^2$ plays a mix. With $\E(t)$, the convergence rate of (<ref>) is given as followed.

Let $X(t)$ be the solution to (<ref>) with initial $X(0)=x_0$, $x^*\in\argmin_x f(x)$. If $f$ is a continuous differentiable closed proper convex function, then
\begin{align*}
\inf_{0\leqslant v\leqslant t}\norm{\nabla f(X(v))}^2&\leqslant\frac{\norm{x_0-x^*}^2}{2\int_0^t\int_0^vc_ududv}.
\end{align*}

The derivative of $\E(t)$ is given as followed.

\begin{align*}
\EE(t)&=c_t[f(X(t))-f(x^*)]+\int_{0}^{t}c_udu\inner{\nabla f(X(t)),\X(t)}+\inner{\X(t),X(t)-x^*}\\
&\leqslant c_t\inner{\nabla f(X(t)),X(t)-x^*}-c_t\inner{\nabla f(X(t)),X(t)-x^*}-\int_{0}^{t}c_udu\norm{\nabla f(X(t))}^2\\
&=-\int_{0}^{t}c_udu\norm{\nabla f(X(t))}^2.
\end{align*}

The above inequality shows that $\E(t)$ is non-increasing, which implies
\[
\int_0^tc_udu[f(X(t))-f(x^*)]\leqslant\E(t)\leqslant\E(0)=\frac{1}{2}\norm{x_0-x^*}^2.
\]

Diving $\int_0^t c_udu$ on both side of above inequality, we obtain the convergence rate of $f(X(t))-f(x^*)$.

Next. integrate $\EE(v)$ from $0$ to $t$, we have
\begin{align*}
\E(t)-\E(0)&=\int_0^t\EE(v)dv\\
&\leqslant-\int_0^t\int_{0}^{v}c_udu\norm{\nabla f(X(v))}^2dv\\
&\leqslant-\inf_{0\leqslant v\leqslant t}\norm{\nabla f(X(v)}^2\int_0^t\int_0^vc_ududv.
\end{align*}

\[
\inf_{0\leqslant v\leqslant t}\norm{\nabla f(X(v))}^2\leqslant\frac{\E(0)-\E(t)}{\int_0^t\int_0^vc_ududv}\leqslant\frac{\norm{x_0-x^*}^2}{2\int_0^t\int_0^vc_ududv}.
\]

Now, we can apply Lyapunov function technique from the continuous case to the discrete case.

Let $f$ be a closed proper convex function, $\{x_k\}$ be the sequence generated by (<ref>), $x^*\in\argmin_x f(x)$. Then we have
\begin{align*}
\min_{0\leqslant j\leqslant k}\norm{x_{j+1}-x_j}^2&\leqslant\frac{\norm{x_0-x^*}^2}{2\sum_{j=0}^k\sum_{i=0}^j\frac{c_i}{c_j}}.
\end{align*}

Define the following Lyapunov function.
\[
\E(k)=\sum_{i=0}^{k-1}c_i[f(x_k)-f(x^*)]+\frac{1}{2}\norm{x_k-x^*}^2.
\]

Difference of $\{\E(k)\}$ is given as followed.
\begin{align*}
\E(k+1)-\E(k)=&c_k[f(x_{k+1})-f(x^*)]+\sum_{i=0}^{k-1}c_i[f(x_{k+1})-f(x_k)]\\
\end{align*}

Since $x_{k}-x_{k+1}\in c_k\partial f(x_{k+1})$, we have
\[
\]

\[
\]

The upper bound of $\E(k+1)-\E(k)$ can be estimated as followed.
\begin{align*}
\E(k+1)-\E(k)\leqslant&\inner{x_k-x_{k+1},x_{k+1}-x^*}+\frac{1}{2}\norm{x_{k+1}-x^*}^2-\frac{1}{2}\norm{x_k-x^*}^2\\
\end{align*}

\[
\sum_{i=0}^{k-1}c_i[f(x_k)-f(x^*)]\leqslant\E(k)\leqslant\E(0)=\frac{1}{2}\norm{x_0-x^*}^2,
\]

Diving $\sum_{i=0}^{k-1}c_i$ on both sides of above inequality, we have
\[
\]

Summing $\E(j+1)-\E(j) $ from $0$ to $k$, we have

\begin{align*}
\E(k+1)-\E(0)&=\sum_{j=0}^k\E(j+1)-\E(j)\\
&\leqslant-\min_{0\leqslant j\leqslant k}\norm{x_{j+1}-x_j}^2\sum_{j=0}^k\sum_{i=0}^j\frac{c_i}{c_j}.
\end{align*}

\[
\min_{0\leqslant j\leqslant k}\norm{x_{j+1}-x_j}^2\leqslant\frac{\E(0)-\E(k+1)}{\sum_{j=0}^k\sum_{i=0}^j\frac{c_i}{c_j}}\leqslant\frac{\norm{x_0-x^*}^2}{2\sum_{j=0}^k\sum_{i=0}^j\frac{c_i}{c_j}}.
\]

§.§ Proximal Point Algorithms for Maximally Monotone Operators

In this section, we build up ODEs framework for
\begin{equation}
\end{equation}

For convenience, we only consider
\begin{equation}
\label{eq:proximalmain}
\end{equation}
That is because (<ref>) can be seem as applying (<ref>) to maximally monotone operator $cA$. Analysis in section <ref> seem to suggest that we consider the ODEs $\X=-A(X)$ with assumption that $A$ is single-value. However, this approach is not feasible due to the reasons described below. Unlike proximal point algorithms for closed proper convex functions, the convergence rate of (<ref>) that researchers concern is $\norm{A(x)}^2$, which is generally not differentiable. If we include $\norm{A(x)}^2$ as potential term in Lyapunov function for $\X=-A(X)$, the Lyapunov function will be generally not differentiable. Inspired by [4, 3], we consider the following ODEs.

\begin{equation}
\X=J_A(X)-X=-\tilde{A}(X).
\label{eq:proximalfirst}
\end{equation}

Above $\tilde{A}$ is the Yosida approximation of $A$. $\tilde{A}$ possesses some properties that is useful in the following text.

Let $\tilde{A}$ be the Yosida approximation of maximally monotone operator $A$.

* $\tilde{A}$ is 1-Lipschitz continuous;

* $\tilde{A}$ is 1-cocoercive, i. e.
\[
\inner{\tilde{A}(x)-\tilde{A}(y), x-y}\geqslant\norm{\tilde{A}(x)-\tilde{A}(y)}^2\quad\forall x,\ y.
\]

The proof of Proposition <ref> can be found in the proof of Corollary 23.11 in [5]. The corresponding Lyapunov function is given as followed.
\begin{equation}
\label{eq:lyapunovfirst}
\E(t)=t\|\X(t)\|^2+\frac{1}{2}\norm{X(t)-x^*}^2
\end{equation}

Here $t\|\X(t)\|^2=t\|\tilde{A}(X)\|^2$ is the potential term, and $\dfrac{1}{2}\norm{X(t)-x^*}^2$ plays a mix. Since $\tilde{A}$ is a single-value 1-Lipschitz continuous operator, the generalized derivative of $\tilde{A}$ exists. Thus, we can analyze the convergence rate of (<ref>) by presuming that $\tilde{A}$ is differentiable.

Let $X(t)$ be the solution to (<ref>) with initial $X(0)=x_0$, $x^*$ be zero of $A$. If $\tilde{A}$ is differentiable, then
\[
\|\X(t)\|^2=\norm{\tilde{A}(X(t))}^2\leqslant\frac{\norm{x_0-x^*}^2}{2t}.
\]

The derivative of $\E(t)$ is estimated as followed.

\begin{align*}
\EE(t)&=-2t\inner{\X(t),\nabla\tilde{A}(X(t))\X(t)}+\norm{\tilde{A}(X(t))}^2-\inner{\tilde{A}(X(t)),X(t)-x^*}\\
&\leqslant -2t\inner{\X(t),\nabla\tilde{A}(X(t))\X(t)}.
\end{align*}

Above inequality is held due to Proposition <ref>. Since $\tilde{A}$ is a monotone operator, $\nabla\tilde{A}(X)$ is semi-positive definite. Then we have $\EE(t)\leqslant 0$,
\[
\]

Diving $t$ on both sides of above inequality, we achieve the required inequality.

Also, the discrete Lyapunov function for (<ref>) can be proposed by using the Lyapunov function (<ref>) as a reference.

Let $\{x_k\}$ be the sequence generated by (<ref>). Then we have

\begin{align*}
\norm{\tilde{A}(x_k)}^2=\norm{x_{k+1}-x_k}^2&\leqslant\frac{\norm{x_0-x^*}^2}{k}.
\end{align*}

Lyapunov function

\begin{equation}
\E(k)=k\norm{x_{k+1}-x_k}^2+\norm{x_k-x^*}^2.
\end{equation}

\begin{align*}
\E(k+1)-\E(k)=&(k+1)\left(\norm{x_{k+2}-x_{k+1}}^2-\norm{x_{k+1}-x_k}^2\right)+\norm{x_{k+1}-x_{k}}^2\\
\leqslant&(k+1)\inner{x_{k+2}-x_{k+1}-(x_{k+1}-x_k),x_{k+2}-x_{k+1}}-2\inner{x_k-x_{k+1},x_{k+1}-x^*}.
\end{align*}

Since $x_k-x_{k+1}\in A(x_{k+1}), x_{k+1}-x_{k+2}\in A(x_{k+2})$, we have $\E(k+1)-\E(k)\leqslant 0$. Thus
\[
\]

Diving $k$ on both sides of above inequality, we obtain required inequality.

§ ACCELERATED PROXIMAL POINT ALGORITHMS FOR CLOSED PROPER CONVEX FUNCTIONS

§.§ Generalized High-resolution ODEs

In section <ref>, we study convergence rate of proximal point algorithms for closed proper convex functions from perspective of gradient flow. Compared to standard gradient flow, the following ODEs proposed in [26] admits faster convergence rate.

\begin{equation}
\begin{split}
\Z&=-\frac{t}{2}\nabla f(X),\\
Z&=\frac{t}{2}\X+\frac{t}{2}\sqrt{s}\nabla f(X)+X,\\
\end{split}
\label{eq:acceleratedmirrorode2}
\end{equation}

In [my previous work], (<ref>) was applied to analyze accelerated mirror descent proposed in [16] and higher-order mirror descent. However, because the coefficient of $\nabla f(X)$ in (<ref>) is an arbitrary positive function $c_t$, directly deriving accelerated proximal point algorithms by using (<ref>) does not meet actual needs. To overcome these problems, we propose the following generalized high-resolution ODEs.
\begin{equation}
\begin{split}
\Z&=-a_t\nabla f(X),\\
Z&=b_t\X+c_t\nabla f(X)+X,\\
\end{split}
\label{eq:acceleratedodeg}
\end{equation}

Above $a_t, b_t, c_t$ are positive on $(0,+\infty)$, $b_0=c_0=0$. Also, we propose the following generalized Lyapunov function.

\begin{equation}
\E(t)=A_t[f(X(t))-f(x^*)]+\frac{1}{2}\norm{Z(t)-x^*}^2.
\label{eq:lyapunovg}
\end{equation}

Here $A_t$ is positive on $(0,+\infty), A_0=0$, $x^*\in\argmin_x f(x)$. Above $A_t[f(X(t))-f(x^*)]$ is the potential terms, $ A_t $ indicates the convergence rate of function value of trajectories to (<ref>), and $\dfrac{1}{2}\norm{Z(t)-x^*}^2$ is the mixed term. It is crucial to study under what conditions $\E(t)$ becomes a Lyapunov function, so we need to study the derivative of $\E(t)$.
\begin{align*}
\EE(t)=&\overset{\cdot}{A_t}[f(X(t))-f(x^*)]+A_t\inner{\nabla f(X(t)),\X(t)}+\inner{\Z(t),Z(t)-x^*}\\
=&\overset{\cdot}{A_t}[f(X(t))-f(x^*)]+A_t\inner{\nabla f(X(t)),\X(t)}+\inner{b_t\X(t)+X_t-x^*+c_t\nabla f(X(t)),-a_t\nabla f(X(t))}.
\end{align*}

Due to convexity of $f$,  we have
\[
f(X(t))-f(x^*)\leqslant\inner{\nabla f(X(t)),X(t)-x^*}.
\]

\begin{align*}
\EE(t)\leqslant&(A_t-a_tb_t)\inner{\nabla f(X(t)),\X(t)}+(\overset{\cdot}{A_t}-a_t)\inner{\nabla f(X(t)),X(t)-x^*}-a_tc_t\norm{\nabla f(X_t)}^2.
\end{align*}

In order to ensure that $ \EE(t)\leqslant 0 $, we impose the following conditions:

\begin{equation}
A_t=a_tb_t,\qquad \overset{\cdot}{A_t}\leqslant a_t.
\label{eq:lyapunonvcondition}
\end{equation}

So far, we can draw the following conclusion.

Let $ f $ be a differentiable closed proper convex function, $ X(t), Z(t) $ be the solution to (<ref>).. If $ A_t=a_tb_t, \overset{\cdot}{A_t}\leqslant a_t, A_0=b_0=c_0=0 $, then we have

\begin{align*}
\inf_{0\leqslant v\leqslant t}\norm{\nabla f(X(v))}^2_2&\leqslant\frac{\norm{x_0-x^*}^2}{2\int_0^t a_vc_vdv}.
\end{align*}

Since $A_t=a_tb_t, \overset{\cdot}{A_t}\leqslant a_t$, we have
\[
\EE(t)\leqslant-a_tc_t\norm{\nabla f(X(t))}^2\leqslant 0.
\]

Then we have
\[
A_t[f(X(t))-f(x^*)]\leqslant\E(t)\leqslant \E(0)=\frac{1}{2}\norm{x_0-x^*}^2.
\]

Integrate $\EE(v)$ from $0$ to $t$, we have
\[
\E(t)-\E(0)=\int_0^t\EE(v)dv\leqslant-\int_0^t a_vc_v\norm{\nabla f(X(v))}^2dv\leqslant-\int_0^ta_vc_vdv\inf_{0\leqslant v\leqslant t}\norm{\nabla f(X(v))}^2.
\]

Then we have
\[
\inf_{0\leqslant v\leqslant t}\norm{\nabla f(X(v))}^2\leqslant\frac{\E(0)-\E(t)}{\int_0^t a_vc_vdv}\leqslant\frac{\norm{x_0-x^*}^2}{2\int_0^t a_vc_vdv}.
\]

Now we discuss the selection of parameters $ a_t, b_t, c_t. $
[Polynomial convergence rate]
Our first goal is to propose a high-resolution ODEs with $ p- $th order convergence rates, where $ p\geqslant 2 $. Let $ A_t=t^p $, $ a_t $ can be chosen as the derivative of $ A_t $, i. e. $ pt^{p-1}. $ Since $ A_t=a_tb_t $, $ b_t=\dfrac{t}{p}. $ The corresponding convergence rate is
\[
\begin{split}
\inf_{0\leqslant v\leqslant t}\norm{\nabla f(X(v))}^2_2\leqslant\frac{\norm{x_0-x^*}^2}{2p\int_0^t v^{p-1}c_vdv}.
\end{split}
\]

[Exponential convergence rate]
Next we want to obtain a high-resolution ODE with exponential convergence rate. Let $ A_t=e^{\lambda t}-1 $, $ a_t $ should be the derivative of $ A_t $, i. e. $ \lambda e^{\lambda t} $. Owing to $ A_t=a_tb_t $, $ b_t=\lambda^{-1}-\lambda^{-1}e^{-\lambda t} $. According to theorem <ref>, the corresponding convergence rate is given
\begin{align*}
f(X_t)-f(x^*)&\leqslant\frac{\norm{x_0-x^*}^2}{2e^{\lambda t}-2},\\
\inf_{0\leqslant v\leqslant t}\norm{\nabla f(X(v))}_2^2&\leqslant\frac{\norm{x_0-x^*}^2}{2\lambda\int_0^t e^{\lambda v}c_vdv}.
\end{align*}

§.§ Discretization  via Symplectic methods

Among the many implicit numerical methods to ODEs, our rationale for adopting Symplectic method is based on the following observation. The low-resolution version of (<ref>) is

\begin{equation}
\label{eq:mirrorlow}
\begin{split}
\Z&=-\frac{t}{2}\nabla f(X),\\
\frac{t}{2}\X+X&=Z.
\end{split}
\end{equation}

We will show that (<ref>) can be translated into a Hamitonian system.
By multiplying $ t $ on both side of second equality of (<ref>), we have
\[
\begin{split}
\Z&=-\frac{t}{2}\nabla f(X),\\
\left(\frac{t^2}{2}X\right)'&=tZ.
\end{split}
\]

Let $ Y=\dfrac{t^2}{2}X $, then we have
\[
\begin{split}
\Z&=-\frac{t}{2}\nabla f\left(\frac{2Y}{t^2}\right)=-\nabla_YH(Y,Z),\\
\overset{\cdot}{Y}&=tZ=\nabla_ZH(Y,Z).
\end{split}
\]

The above Hamiltonian function $ H(Y,Z)=\dfrac{t^3}{4}f\left(\dfrac{2Y}{t^2}\right)+\dfrac{t}{2}Z^2 $. Inspired by this observation and the fact that symplectic method is powerful in numerically solving Hamiltonian systems, we take advantage of implicit symplectic method instead of Euler method and propose the following recursive framework

\begin{equation}
\begin{split}
z_k&=(b_k+1)x_{k+1}+c_k\nabla f(x_{k+1})-b_kx_k,\\
z_{k+1}&=z_k-a_k\nabla f(x_{k+1}),\\
\end{split}
\label{eq:symplecticproximal1}
\end{equation}

Above parameters $ a_k, b_k, c_k $ are positive. Since the resolvent of $\nabla f$ is the proximal operator, (<ref>) is equivalent to the following algorithm.

Symplectic proximal point algorithms, SPPA

Initialize $ x_0$ $, z_0=x_0 $
Termination condinitions aren't hold
$ y_{k+1}=\dfrac{1}{b_k+1}z_k+\dfrac{b_k}{b_k+1}x_k $
$ x_{k+1}=\mathrm{prox}_{\frac{c_k}{b_k+1}f}(y_{k+1}) $
$ z_{k+1}=z_k+\dfrac{a_k}{c_k}(b_k+1)(x_{k+1}-y_{k+1})$

Output $ x_{k+1}. $

Compared to proximal point algorithms, Algorithm <ref> yields an additional computational cost of 3 scalar products and 2 vector additions per iterations, which means Algorithm <ref> does not increase computational complexity per iteration. To investigate convergence results of Algorithm <ref>, we propose the corresponding discrete-time Lyapunov function given by
\begin{equation}
\E(k)=A_k[f(x_k)-f(x^*)]+\frac{1}{2}\norm{z_k-z^*}^2.
\end{equation}

Suggested by (<ref>), we assume that

\begin{equation}\label{eq:discretelyapunovkey}
A_k= a_kb_k,\qquad A_{k+1}-A_k\leqslant a_k.
\end{equation}

Now we can show that the convergence result of Algorithm <ref> is quite analogous to the continuous time case. First, we need the following lemma.

If $ A_k= a_kb_k, A_{k+1}-A_k\leqslant a_k $, then

\[
b_kA_{k+1}\leqslant (b_k+1)A_k.
\]

Since $ A_{k+1}-A_k\leqslant a_k $, we have
\[
(A_{k+1}-A_k)b_k\leqslant a_kb_k=A_k.
\]

Reorder the inequality, we have
\[
\]

Now we can prove the desired convergence result.

Let $f$ be a closed proper convex function, $\{x_k\}, \{z_k\}$ be the sequences generated by
\begin{align*}
\end{align*}
If the parameters $\{A_k\}, \{a_k\}, \{b_k\}, \{c_k\}$ satisfy
\[
A_k=a_kb_k,\quad A_{k+1}-A_k\leqslant a_k,\quad a_k\leqslant c_k,\quad A_0=b_0=0,
\]

then we have

\begin{equation}
\begin{gathered}
f(x_k)-f(x^*)\leqslant\frac{\norm{x_0-x^*}^2}{2A_k},\quad\forall k\geqslant 1,\\
\min_{0\leqslant j\leqslant k}\norm{\frac{(b_k+1)x_{k+1}-b_kx_k-x_k}{c_k}}^2\leqslant\frac{\norm{x_0-x^*}^2}{\sum_{j=0}^ka_j^2}.
\end{gathered}
\end{equation}

Recall the Lyapunov function
\[
\E(k)=A_k[f(x_k)-f(x^*)]+\frac{1}{2}\norm{z_k-x^*}^2.
\]

Let $\tilde{\nabla}f(x_{k+1})=\dfrac{b_kx_k+z_k-(b_k+1)x_{k+1}}{c_k}$. By first order characterization of proximal operators, we have $\tilde{\nabla}f(x_{k+1})\in\partial f(x_{k+1})$.

Step 1: Separate the difference $ \E(k+1)-\E(k) $ into three parts.

\begin{align*}
\E(k+1)-\E(k)=&\underbrace{A_k [f(x_{k+1})-f(x_k)]}_{\text{I}}+\underbrace{(A_{k+1}-A_k)[f(x_{k+1})-f(x^*)]}_{\text{II}}\\
\end{align*}

Step 3: Reckon the upper bound of three parts. Since $\tilde{\nabla}f(x_{k+1})\in\partial f(x_{k+1})$, we have
\begin{align*}
\text{I}&\leqslant \underbrace{A_k\inner{\tilde{\nabla} f(x_{k+1}),x_{k+1}-x_k}}_{\text{I}_1},\\
\text{II}&\leqslant \underbrace{(A_{k+1}-A_k)\inner{\tilde{\nabla} f(x_{k+1}),x_{k+1}-x^*}}_{\text{II}_1}.
\end{align*}

For III, we have

\begin{align*}
\text{III}&=\inner{z_{k+1}-z_k,z_k-x^*}+\frac{1}{2}\norm{z_{k+1}-z_k}^2\\
\end{align*}

Step 2: Estimate the upper bound of $ \E(k+1)-\E(k) $. From previous arguments, we have
\begin{align*}
\E(k+1)-\E(k)\leqslant&\underbrace{A_k\inner{\tilde{\nabla} f(x_{k+1}),x_{k+1}-x_k}}_{\text{I}_1}+\underbrace{(A_{k+1}-A_k)\inner{\tilde{\nabla}f(x_{k+1}),x_{k+1}-x^*}}_{\text{II}_1}\\
&+\underbrace{\inner{z_{k+1}-z_k,z_k-x^*}}_{\text{III}_1}+\underbrace{\frac{a_k^2}{2}\norm{\tilde{\nabla} f(x_k)}^2}_{\text{III}_2}.
\end{align*}

First, consider  $ \text{I}_1+\text{II}_1 $. By Lemma <ref>, we have
\begin{align*}
\text{I}_1+\text{II}_2\leqslant & A_k\inner{\tilde{\nabla} f(x_{k+1}),x_{k+1}-x_k}+\frac{A_k}{b_k}\inner{\tilde{\nabla}f(x_{k+1}),x_{k+1}-x^*}\\
=&A_k\inner{\tilde{\nabla} f(x_{k+1}),\frac{b_k+1}{b_k}x_{k+1}-x_k-\frac{1}{b_k}x^*}.
\end{align*}

Since $ (b_k+1)x_{k+1}=b_kx_k+z_k-c_k\tilde{\nabla} f(x_{k+1}) $, we have
\begin{align*}
\text{I}_1+\text{II}_1=&\frac{A_k}{b_k}\inner{\tilde{\nabla} f(x_{k+1}),z_k-x^*-c_k\tilde{\nabla} f(x_{k+1})}\\
=&a_k\inner{\tilde{\nabla} f(x_{k+1}),z_k-x^*-c_k\tilde{\nabla} f(x_{k+1})}\\
=&-\inner{z_{k+1}-z_k,z_k-x^*}-a_kc_k\norm{\tilde{\nabla} f(x_{k+1})}_2^2.
\end{align*}

\begin{align*}
\E(k+1)-\E(k)\leqslant\left(\frac{a_k^2}{2}-a_kc_k\right)\norm{\tilde{\nabla} f(x_{k+1})}^2_2.
\end{align*}

Step 3: Deduce the convergence rates. Since $c_k\geqslant a_k$, we have
\[
\E(k+1)-\E(k)\leqslant-\frac{a_k^2}{2}\norm{\tilde{\nabla}f(x_{k+1})}^2\leqslant 0.
\]

Then we have
\[
\]

Diving $A_k$ on both sides of the above inequality, we have
\[
\]

Summing $\E(j+1)-\E(j)$ from $0$ to $k$, we have
\begin{align*}
\E(k+1)-\E(0)&=\sum_{j=0}^k\E(j+1)-\E(j)\leqslant-\sum_{j=0}^k\frac{a_k^2}{2}\norm{\tilde{\nabla}f(x_{k+1})}^2\leqslant-\sum_{j=0}^k\frac{a_k^2}{2}\min_{0\leqslant j\leqslant k}\norm{\tilde{\nabla}f(x_{k+1})}^2.
\end{align*}

Then we have
\[
\min_{0\leqslant j\leqslant k}\norm{\tilde{\nabla}f(x_{k+1})}^2\leqslant\frac{2\E(0)}{\sum_{j=0}^ka_j^2}=\frac{\norm{x_0-x^*}^2}{\sum_{j=0}^ka_j^2}.
\]

Next, we discuss some options for the parameters.

First, we require that the index of proximal point operator in Algorithm <ref> is constant. Let $ A_k=\dfrac{c}{2}k(k+1), c>0 $, according to (<ref>), we have $ a_k=c(k+1), b_k=\dfrac{k}{2} $. If $ c_k=c(k+2) $, then $ c_k\geqslant a_k $ is held, and the index of proximal point operator is $c$. The resulting convergence rates is

\begin{equation}
\begin{gathered}
\min_{0\leqslant j\leqslant k}\norm{\frac{(j+2)x_{j+1}-jx_j-2z_j}{2(j+1)}}^2_2\leqslant\frac{6\norm{x_0-x^*}^2}{c^2(k+1)(k+2)(2k+3)}.
\end{gathered}
\end{equation}

Next, we try to obtain $ p- $th order convergence rates accelerated proximal point algorithm. Let $ A_k=k^{(p)}=k(k+1)\cdots(k+p-1) $ rather than $ A_k=k^{(p)} $ with integer $ p\geqslant 2 $. Followed by (<ref>), $ a_k=p(k+1)^{(p-1)} $, $ b_k=\dfrac{k}{p} $. Finally, we let $ c_k=a_k. $ The corresponding convergence rates is

\begin{equation}
\begin{gathered}
\min_{0\leqslant j\leqslant k}\norm{\frac{(j+p)x_{j+1}-jx_j-pz_j}{p^2(j+1)^{(p-1)}}}^2\leqslant\frac{\norm{x_0-x^*}^2}{p^2\sum_{j=0}^{k}[(j+1)^{(p-1)}]^2}.
\end{gathered}
\end{equation}

In the end, we derive exponential convergence rates accelerated proximal point algorithm. Let $ A_k=\rho^k-1 $, due to (<ref>), we have $ a_k=\rho^{k}(\rho-1) $, $ b_k=\dfrac{\rho^k-1}{\rho^k(\rho-1)}. $ The parameter $ c_k $ is chosen to be $ a_k $. The corresponding convergence rates is

\begin{equation}
\begin{gathered}
\min_{0\leqslant j\leqslant k}\norm{\frac{(\rho^{k+1}-1) x_{j+1}-(\rho^k-1)x_j-z_j}{\rho^{2j}(\rho-1)^2}}^2\leqslant\frac{(\rho^2-1)\norm{x_0-x^*}^2}{(\rho-1)^2(\rho^{2k+2}-1)}.
\end{gathered}
\end{equation}

In the end of this section, we discuss significance of the conditions $A_{k+1}-A_k\leqslant a_k$, $A_k=a_kb_k$. Let $b_k\equiv 0, a_k\equiv c_k$, then Algorithm <ref> becomes the proximal point algorithms (<ref>). Because $A_k=a_kb_k, b_k\equiv 0$, then $A_k\equiv 0$ is held, which contradicts to $A_{k+1}-A_k=a_k=c_k\ne 0$. In conclusion, the conditions $A_{k+1}-A_k\leqslant a_k, A_k=a_kb_k$ must be held so that our symplectic proximal point algorithms converge faster than proximal point algorithms.

§.§ Application: Symplectic Augmented Lagrangian Multipliers Methods

In this section, we will develop a new augmented Lagrangian multipliers methods via symplectic proximal point algorithms. For a given objective function $f:\mathbb{R}^n\to\overline{\mathbb{R}}$, let $\varphi:\mathbb{R}^n\times\mathbb{R}^m\to\overline{\mathbb{R}}$ be perturbation of $f$, i. e. $\varphi(x,0)=f(x), \forall x\in\mathbb{R}^n.$ We assume that the perturbation function $\varphi$ is a closed proper convex function. The Lagrangian function associated with $\varphi$ is
\begin{equation}
L(x,\lambda)=\inf_u \varphi(x,u)-\inner{\lambda,u}.
\end{equation}

For primal problem
\[
\min_x f(x),
\]

the Lagrangian duality is given as followed
\begin{equation}
\label{eq:lagrangianduality}\tag{LP}
\max_\lambda\min_x L(x,\lambda).
\end{equation}

The Lagrangian duality theory we need below is stated as followed.

Assume $\varphi$ is a closed proper convex function. If $(x^*,\lambda^*)$ is the saddle point of $L$, then
\[
\min_x f(x)=f(x^*)=L(x^*,\lambda^*).
\]

The proof of Theorem <ref> can be found in section 2.5 in [6]. If we consider $\min_x L(x,\lambda)$ as an objective function respect to $\lambda$, (<ref>) becomes an optimization problem, then we can apply Algorithm <ref> to (<ref>), and obtain the following equation.
\begin{equation}
\begin{split}
\tilde{\lambda}_{k+1}&=\frac{1}{b_k+1}z_k+\frac{b_k}{b_k+1}\lambda_k,\\
\lambda_{k+1}&=\mathop{\arg\max}_\lambda\min_x L(x,\lambda)-\frac{b_k+1}{2c_k}\norm{\lambda-\tilde{\lambda}_k}^2,\\
\end{split}
\label{eq:symplecticalm1}
\end{equation}

Due to Theorem 34.3 in [20], the iteration formula of $\lambda_{k+1}$ can be simplified as

\begin{equation}
\begin{split}
\lambda_{k+1}&=\lambda_k-\frac{c_k}{b_k+1}u_{k+1}.
\end{split}
\label{eq:symplecticalm2}
\end{equation}

In practice, above $\norm{u}^2$ can be replaced by $\norm{u}_H^2=u\tr Hu$, where $H$ is a positive-definite matrix. To meet practical requirements, we transform (<ref>) and (<ref>) into the following Algorithm <ref>.

Symplectic augmented Lagrangian multipliers methods, S-ALM

Initialize $ \lambda_0$ $, z_0=\lambda_0 $
Termination condinitions aren't hold
$ \tilde{\lambda}_{k+1}=\dfrac{1}{b_k+1}z_k+\dfrac{b_k}{b_k+1}\lambda_k $
$ (x_{k+1},u_{k+1})=\argmin_{x,u}\varphi(x,u)-\inner{\tilde{\lambda}_{k+1},u}+\dfrac{c_k}{2(b_k+1)}\norm{u}_H^2 $
$ z_{k+1}=z_k+\dfrac{a_k}{c_k}(b_k+1)H(\lambda_{k+1}-\tilde{\lambda}_{k+1}). $

Output $ x_{k+1}, \lambda_{k+1}. $

We consider the convex programming with equality constraints
\begin{equation}
\label{eq:equalityconstrants}
\begin{gathered}
\min_x f(x),\\
s.\ t. \quad Ax=b,
\end{gathered}
\end{equation}

where $ f:\mathbb{R}^n\to\overline{\mathbb{R}} $ is a closed proper convex function, $ A\in\mathbb{R}^{m\times n} $ is a matrix and $ \mathrm{ri}(\mathrm{dom}f)\cap\mathrm{im}(A)\ne\emptyset $. The corresponding perturbation function is given by
\[
\varphi(x,u)=f(x)+I(Ax-b=u).
\]

Here, $I$ is the indicator function. It is easily to check that $\varphi$ is closed and proper. According to Algorithm <ref>, the corresponding recursive rule is

\begin{equation}\label{eq:salmequalityconstraints}
\begin{split}
\tilde{\lambda}_{k+1}&=\frac{1}{b_k+1}z_k+\frac{b_k}{b_k+1}\lambda_k,\\
x_{k+1}&=\argmin_x f(x)-\inner{\tilde{\lambda}_{k+1},Ax-b}+\frac{c_k}{2(b_k+1)}\norm{Ax-b}^2_H,\\
\lambda_{k+1}&=\tilde{\lambda}_{k+1}-\frac{c_k}{b_k+1}H(Ax_{k+1}-b),\\
\end{split}
\end{equation}

Next, we discuss convergence rates of symplectic augmented Lagrangian method. Since Algorithm <ref> is derived from applying Algorithm <ref> to Lagrangian duality, the convergence rates of Algorithm <ref> we required is that
\[
L(x^*,\lambda^*)-\min_x L(x,\lambda_k).
\]

Here $(x^*,\lambda^*)$ is saddle point of concave-convex function $L$. Moreover, $\min_x L(x,\lambda_k)$ can be represented by $x_k, u_k$.

Let $\{x_k\}, \{u_k\}, \{\lambda_k\}, \{\tilde{\lambda}_k\}$ be the sequences generated by Algorithm <ref>. Then
\[
\min_x L(x,\lambda_k)=L(x_k,\lambda_k)=\varphi(x_k,u_k)-\inner{\lambda_k,u_k},\quad\forall k\geqslant 1.
\]

By the first-order characterization of $(x_{k+1},u_{k+1})$, we have
\[
\]

Thus $(x_{k+1},u_{k+1})\in\argmin_{x,u}\varphi(x,u)-\inner{\lambda_{k+1},u}$. By the definition of Lagrangian function, we have
\[
\min_x L(x,\lambda_{k+1})=L(x_{k+1},\lambda_{k+1})=\varphi(x_{k+1},u_{k+1})-\inner{\lambda_{k+1},u_{k+1}},\quad\forall k\geqslant 0.
\]

With Lemma <ref>, we can prove the convergence results of Algorithm <ref>.

Let $\{x_k\}, \{u_k\}, \{\lambda_k\}, \{\tilde{\lambda}_k\}$ be the sequences generated by
\begin{align*}
\tilde{\lambda}_{k+1}&=\dfrac{1}{b_k+1}z_k+\dfrac{b_k}{b_k+1}\lambda_k \\
(x_{k+1},u_{k+1})&=\argmin_{x,u}\varphi(x,u)-\inner{\tilde{\lambda}_{k+1},u}+\dfrac{c_k}{2(b_k+1)}\norm{u}_H^2 \\
\lambda_{k+1}&=\tilde{\lambda}_{k+1}-\dfrac{c_k}{b_k+1}Hu_{k+1}\\
\end{align*}

Let $(x^*,\lambda^*)$ be the saddle point of Lagrangian function $L$. If $\{A_k\}, \{a_k\}, \{b_k\}, \{c_k\}$ be the sequences that satisfy
\[
A_{k+1}-A_k=a_k,\quad A_k=a_kb_k,\quad A_0=b_0=0, \quad c_k\geqslant a_k,
\]

then we have
\begin{align*}
\min_{0\leqslant j\leqslant k}\norm{u_{j+1}}^2_H&\leqslant\frac{\norm{x_0-x^*}^2_{H^{-1}}}{\sum_{j=0}^k a_j^2}.
\end{align*}

Consider the following Lyapunov function
\[
\E(k)=A_k[L(x^*,\lambda^*)-L(x_k,\lambda_k)]+\frac{1}{2}\norm{z_k-x^*}_{H^{-1}}^2.
\]

Step 1: Diving the difference $\E(k+1)-\E(k)$ into 3 parts.
\begin{align*}
\E(k+1)-\E(k)=&\underbrace{A_k[L(x_k,\lambda_k)-L(x_{k+1},\lambda_{k+1})]}_{\text{I}}+\underbrace{(A_{k+1}-A_k)[L(x^*,\lambda^*)-L(x_{k+1},\lambda_{k+1})]}_{\text{II}}\\
\end{align*}

If $k\geqslant 1$, by Lemma <ref>, we have
\begin{align*}
\text{I}=A_k[\varphi(x_k,u_k)-\inner{\lambda_k,u_k}-\varphi(x_{k+1},u_{k+1})+\inner{\lambda_{k+1},u_{k+1}}].
\end{align*}

By the proof of Lemma <ref>, we have
\[
\binom{0}{\lambda_k}\in\partial\varphi(x_k,u_k).
\]

\begin{align*}
\text{I}&\leqslant A_k(\inner{\lambda_k,u_k-u_{k+1}}-\inner{\lambda_k,u_k}+\inner{\lambda_{k+1},u_{k+1}})\\
\end{align*}

If $k=0$, $0=\text{I}\leqslant A_k\inner{\lambda_{k+1}-\lambda_k.u_{k+1}}=0$ is also held.

By Theorem <ref>, we have $\binom{0}{\lambda^*}\in\partial\varphi(x^*,0)$. Owing to the facts that $\binom{0}{\lambda^*}\in\partial\varphi(x^*,0)$ and $A_{k+1}-A_k\leqslant a_k$, the upper bound of II can be estimated as followed.

\begin{align*}
\text{II}&\leqslant a_k[L(x^*,\lambda^*)-L(x_{k+1},\lambda_{k+1})]\\
&\leqslant a_k(-\inner{\lambda^*,u_{k+1}}+\inner{\lambda_{k+1},u_{k+1}})\\
\end{align*}

For III, we have
\[
\text{III}=-\underbrace{a_k\inner{z_k-\lambda^*,u_{k+1}}}_{\text{III}_1}+\underbrace{\frac{a_k^2}{2}\norm{u_{k+1}}_H^2}_{\text{III}_2}.
\]

Step 2: Reckon the upper bound of $\E(k+1)-\E(k)$. Combining the upper bound of I, II, III, we have
\begin{align*}
\E(k+1)-\E(k)\leqslant&\underbrace{A_k\inner{\lambda_{k+1}-\lambda_k,u_{k+1}}}_{\text{I}_1}+\underbrace{a_k\inner{\lambda_{k+1}-\lambda^*,u_{k+1}}}_{\text{II}_1}\\
\end{align*}

Since $A_k=a_kb_k$, we have
\begin{align*}
\text{I}_1+\text{II}_1&\leqslant a_k\inner{(b_k+1)\lambda_{k+1}-b_k\lambda_k-\lambda^*,u_{k+1}}\\
\end{align*}

\[
\E(k+1)-\E(k)\leqslant\left(-a_kc_k+\frac{a_k^2}{2}\right)\norm{u_{k+1}}_H^2\leqslant-\frac{a_k^2}{2}\norm{u_{k+1}}_H^2.
\]

Step 3: Deduce the convergence rates. Since $\E(k+1)-\E(k)\leqslant 0$, we have
\[
\]

Diving $A_k$ on both sides of above inequality, we have
\[
\]

Since $\E(k+1)-\E(k)\leqslant-\dfrac{a_k^2}{2}\norm{u_{k+1}}_H^2$, we have
\begin{align*}
\E(k+1)-\E(0)&=\sum_{j=0}^k\E(j+1)-\E(j)\leqslant\sum_{j=0}^k-\frac{a_j^2}{2}\norm{u_{j+1}}^2_H\leqslant-\min_{0\leqslant j\leqslant k}\norm{u_{j+1}}^2_H\sum_{j=0}^k\frac{a_j^2}{2}.
\end{align*}

Then we have
\[
\min_{0\leqslant j\leqslant k}\norm{u_{j+1}}^2_H\leqslant\frac{\norm{x_0-x^*}^2_{H^{-1}}}{\sum_{j=0}^ka_j^2}.
\]

Let $\{x_k\}, \{z_k\}, \{\lambda_k\}, \{\tilde{\lambda}_k\} $ be the sequences generated by (<ref>). Let $(x^*,\lambda^*)$ be the saddle point of Lagrangian function $L$. If $\{A_k\}, \{a_k\}, \{b_k\}, \{c_k\}$ be the sequences that satisfy
\[
A_{k+1}-A_k\leqslant a_k,\quad A_k=a_kb_k,\quad A_0=b_0=0, \quad c_k\geqslant a_k,
\]

then we have
\begin{align*}
\min_{0\leqslant j\leqslant k}\norm{Ax_{j+1}-b}^2_H&\leqslant\frac{\norm{\lambda_0-\lambda^*}^2_{H^{-1}}}{\sum_{j=0}^k a_j^2}.
\end{align*}

§ ACCELERATED PROXIMAL POINT ALGORITHMS FOR MAXIMALLY MONOTONE OPERATORS

§.§ High-resolution ODEs Framework

Since it is not possible to directly transform the Nesterov accelerated gradient method into an accelerated proximal point algorithm, a direct application of (<ref>) to proximal point operator for maximally monotone operator is not feasible. In order to build up high-resolution ODEs for proximal point operator, first we need to study continuous-limit of (<ref>). For convinice, we also assume that $c=1$ in (<ref>). Combining the two equations of (<ref>), we have
\begin{equation}
\label{eq:temp1}
\end{equation}

By introducing the assumptions that $\tilde{A}$ is differentiable, $y_k\approx X(k)$ for some smooth curve $X(t)$ defined for $t\geqslant 0$, we obtain the following Taylor expansions.
\begin{align*}
\tilde{A}(y_{k-1})&=\tilde{A}(X(k))-\nabla\tilde{A}(X(k))+O(1).
\end{align*}

Substitute $y_{k+1}, y_{k-1}, \tilde{A}(y_{k-1})$ in (<ref>) by above equations, and omit infinitesimal terms, we obtain the following ODEs
\begin{equation}
\XX+\frac{2}{t}\X+\frac{2}{t}\tilde{A}(X)+2\nabla\tilde{A}(X)\X=0.
\label{eq:acceleratedode1}
\end{equation}

If $X(0)=x_0$, it is easily to verify that (<ref>) is equivalent to
\begin{equation}
\label{eq:acceleratedode2}
\end{equation}

The convergence results of trajectories of (<ref>) are shown below.

Let $X(t)$ be the solution to (<ref>). If $\tilde{A}$ is differentiable, then we have
\begin{align*}
\norm{\tilde{A}(X(t))}^2\leqslant\frac{\norm{x_0-x^*}^2}{t^2+2t},\\
\inner{\tilde{A}(X(t)),X(t)-x^*}\leqslant\frac{\norm{x_0-x^*}^2}{2t}.
\end{align*}

Consider the following Lyapunov function
\[
\E(t)=t^2\norm{\tilde{A}(X(t))}^2+t\inner{\tilde{A}(X(t)),X(t)-x_0}.
\]

The derivative of $\E(t)$ is
\begin{align*}
\EE(t)=&2t^2\inner{\nabla\tilde{A}(X(t))\X(t),\tilde{A}(X(t))}+2t\norm{\tilde{A}(X(t))}^2\\
\end{align*}

In view of the fact that $\tilde{A}$ is monotone, $\nabla\tilde{A}$ is semi-positive definite. Then we have
\[
\EE(t)\leqslant 0, \quad \E(t)\leqslant\E(0)=0.
\]

Because of Proposition <ref>, we have
\begin{align*}
\E(t)=&t^2\norm{\tilde{A}(X(t))}^2+t\inner{\tilde{A}(X(t)),X(t)-x^*}+t\inner{\tilde{A}(X(t)),x^*-x_0}\\
\geqslant& \frac{t^2}{2}\norm{\tilde{A}(X(t))}^2+t\inner{\tilde{A}(X(t)),X(t)-x^*}-\frac{1}{2}\norm{x_0-x^*}^2.
\end{align*}

Then we have
\[
\frac{t^2}{2}\norm{\tilde{A}(X(t))}^2+t\inner{\tilde{A}(X(t)),X(t)-x^*}\leqslant\frac{1}{2}\norm{x_0-x^*}^2.
\]

By Proposition <ref>, we have
\begin{align*}
\left(\frac{t^2}{2}+t\right)\norm{\tilde{A}(X(t))}^2\leqslant\frac{1}{2}\norm{x_0-x^*}^2
\end{align*}

Diving $\dfrac{t^2}{2}+t$ on both sides of above inequality, we obtain
\[
\norm{\tilde{A}(X(t))}^2\leqslant\frac{\norm{x_0-x^*}^2}{t^2+2t}.
\]

Also, we have
\[
\]

Diving $t$ on both sides of above inequality, we obtain the desired inequality.

For applying discrete-time Lyapunov function to (<ref>), we refer to section 12.2 in [22]. It is worthwhile to mention that the discussion of $\inner{\tilde{A}(X(t)),X(t)-x^*}$, which is omitted in most previous work, plays crucial role in the following analysis. First of all, in the proof of Theorem <ref>, we prove that
\[	\frac{t^2}{2}\norm{\tilde{A}(X(t))}^2+t\inner{\tilde{A}(X(t)),X(t)-x^*}\leqslant\frac{1}{2}\norm{x_0-x^*}^2,
\]

which provides us valuable clues for designing Lyapunov function in the following text. Also, the role of $\inner{\tilde{A}(X(t)),X(t)-x^*}$ in (<ref>) is similar to the role of $f(X(t))-f(x^*)$ in (<ref>). For example, if $\tilde{A}$ is the gradient of a closed proper convex function $f$, then
\[
\]

which means that an upper bound of $\inner{\tilde{A}(X(t)),X(t)-x^*}$ can also be an upper bound of $f(X(t))-f(x^*)$.

Now, we discuss how to modify (<ref>) to fit $\tilde{A}$ without lossing faster convergence rates. It is worth noting that (<ref>) can be rewritten as
\begin{align*}
\Z&=0,\\
\end{align*}

Above observation suggests us to consider the following ODEs
\begin{equation}\label{eq:acceleratedproximalode}
\begin{split}
\Z&=-\frac{C}{r}\tilde{A}(X),\\
\end{split}
\end{equation}

Here $C>0, r>0, a_t$ is positive on $(0,+\infty], a_0=0$. The corresponding Lyapunov function is given by

\begin{equation}
\begin{split}
\E(t)=&\frac{Ct(ra_t-Ct)}{2}\norm{\tilde{A}(X(t))}^2+Crt\inner{\tilde{A}(X(t)),X(t)-x^*}\\
\end{split}
\label{eq:acceleratedproximallyapunov}
\end{equation}

The sum of first two term is the potential term, and the sum of last two term plays a mix. In the proof of Theorem <ref>, we receive
\[
\frac{t^2}{2}\norm{\tilde{A}(X(t))}^2+t\inner{\tilde{A}(X(t)),X(t)-x^*}\leqslant\frac{1}{2}\norm{x_0-x^*}^2,
\]

which inspires us that the potential term of Lyapunov function should be linear combination of $\norm{A(X(t))}^2$ and $\inner{\tilde{A}(X(t)),X(t)-x^*}$. Mixed term are chosen to ensure that $\EE(t)$ is consist of linear combination of three non-negetive terms. Before describing convergence results of (<ref>), we need to estimate $\EE(t)$.

Let $X(t), Z(t)$ be the solutions to (<ref>). If $\tilde{A}$ is differentiable, then we have
\begin{align*}
\EE(t)\leqslant&-Ct^2\inner{\nabla\tilde{A}(X(t))\X(t),\X(t)}-C(r^2-r)\inner{\tilde{A}(X(t)),X(t)-x^*}\\
\end{align*}

Directly calculate the derivative of $\E(t)$, then we have

\begin{align*}
\EE(t)=&Ct(ra_t-Ct)\inner{\nabla\tilde{A}(X(t))\X(t),\tilde{A}(X(t))}+\frac{C(ra_t-Ct)+Ct(r\overset{\cdot}{a}_t-C)}{2}\norm{\tilde{A}(X(t))}^2\\
\end{align*}

To ensure that $\EE(t)\leqslant 0$, the function $a_t$ should satisfy
\[
Cr^2a_t-\frac{C(ra_t-Ct)+Ct(r\overset{\cdot}{a}_t-C)}{2}\geqslant 0.
\]

Determining all $a_t$ that satisfy above condition is difficult. For convenience, we focus on the special case $a_t=1+\dfrac{t}{r}$, which plays central role in the following section.

Let $X(t), Z(t)$ be the solutions to

\begin{equation}
\begin{split}
\Z&=-\frac{C}{r}\tilde{A}(X),\\
\end{split}
\end{equation}

If $\tilde{A}$ is differentiable, $0<C\leqslant 1, r\geqslant 2$, then we have
\begin{align*}
\norm{\tilde{A}(X(t))}^2&\leqslant\frac{(r^3-r^2)\norm{x_0-x^*}^2}{Ct(3r+t-Ct)},\\
\inner{\tilde{A}(X(t)),X(t)-x^*}&\leqslant\frac{(r^2-r)\norm{x_0-x^*}^2}{2Ct}.
\end{align*}

By Lemma <ref>, we have

\begin{align*}
\EE(t)\leqslant&-Ct^2\inner{\nabla\tilde{A}(X(t))\X(t),\X(t)}-C(r^2-r)\inner{\tilde{A}(X(t)),X(t)-x^*}\\
\end{align*}

Owing to $r\geqslant 2, C>0$, we have $Cr^2-\dfrac{Cr}{2}-C(r+1-C)t\geqslant 0$. Thus  $\EE(t)\leqslant 0$ by monotonicity of $\tilde{A}$. Using Proposition <ref>, we have
\[
\frac{1}{2}Ct(3r+t-Ct)\norm{\tilde{A}(X(t))}^2\leqslant\E(t)\leqslant\E(0)=\frac{r^3-r^2}{2}\norm{x_0-x^*}^2.
\]

Diving $\dfrac{1}{2}Ct(3r+t-Ct)$ on both sides of the above inequality, we have
\[
\norm{\tilde{A}(X(t))}^2\leqslant\frac{(r^3-r^2)\norm{x_0-x^*}^2}{Ct(3r+t-Ct)}.
\]

Also, we have
\begin{align*}
Crt\inner{\tilde{A}(X(t)),X(t)-x^*}\leqslant\E(t)\leqslant \E(0)=\frac{r^3-r^2}{2}\norm{x_0-x^*}^2.
\end{align*}

Diving $Crt$ on both sides of the above inequality, we obtain the desired inequality.

§.§ Symplectic Discretization

Similar to section <ref>, we also apply implicit symplectic method to (<ref>), and obtain

\begin{equation}
\label{eq:temp3}
\begin{split}
\end{split}
\end{equation}

To simplify (<ref>), the following property of $\tilde{A}$ is needed.

Let $\tilde{A}$ be Yosida approximation of maximally monotone operator $A$. The resolvent of $\tilde{A}$ is given by
\[
\]

By Proposition <ref>, (<ref>) is equivalent to
\begin{equation}
\begin{split}
\tilde{x}_{k+1}&=\frac{r}{k+r}z_k+\frac{k}{k+r}x_k,\\
\end{split}\label{eq:symplecticproximal}
\end{equation}

If $a_k=1+\dfrac{k}{r}$, which makes (<ref>) become directly symplectic discretization of ODEs in Theorem <ref>, then the index of proximal operator in (<ref>) is constant 2. Next, we deduce convergence rates of (<ref>) with $a_k=1+\dfrac{k}{r}$.

Let $\{x_k\}, \{z_k\}$ be the sequences generated by
\begin{equation}
\label{eq:symplecticproximalconstant}
\begin{split}
\end{split}
\end{equation}

If $r\geqslant 2, 0<C\leqslant 1$, then we have
\begin{align*}
\norm{\tilde{A}(x_k)}^2&\leqslant\frac{(r^3-r^2)\norm{x_0-x^*}^2}{Ck[k+3r-C(k+1)]},\\
\inner{\tilde{A}(x_k),x_k-x^*}&\leqslant\frac{(r^2-r)\norm{x_0-x^*}^2}{2Ck}.
\end{align*}

Consider the following Lyapunov function:

\begin{align*}
\E(k)=&\frac{Ck[k+r-C(k+1)]}{2}\norm{\tilde{A}(x_k)}^2+Crk\inner{\tilde{A}(x_k),x_k-x^*}\\
\end{align*}

Step 1: Diving the Difference of $\{\E(k)\}$ into 4 parts:
\begin{align*}
\E(k+1)-\E(k)=&\underbrace{\frac{C(k+1)[k+1+r-C(k+2)]}{2}\norm{\tilde{A}(x_{k+1})}^2-\frac{Ck[k+r-C(k+1)]}2\norm{\tilde{A}(x_k)}^2}_{\text{I}}\\
\end{align*}

By the inequality
\[
\frac{1}{2}\norm{a}^2-\frac{1}{2}\norm{b}^2\leqslant\inner{a,a-b},\quad\forall a,b\in\mathbb{R}^n,
\]

we have the following estimation:
\begin{align*}
\text{I}=&\frac{Ck[k+r-C(k+1)]}{2}\left(\norm{\tilde{A}(x_{k+1})}^2-\norm{\tilde{A}(x_k)}^2\right)+\frac{C[2k(1-C)+r+1-2C]}{2}\norm{\tilde{A}(x_{k+1})}^2\\
\leqslant&\underbrace{Ck[k+r-C(k+1)]\inner{\tilde{A}(x_{k+1})-\tilde{A}(x_k),\tilde{A}(x_{k+1})}}_{\text{I}_1}+\underbrace{\frac{C[2k(1-C)+r+1-2C]}{2}\norm{\tilde{A}(x_{k+1})}^2}_{\text{I}_2},\\
\text{III}\leqslant&\inner{Ck[\tilde{A}(x_{k+1})-\tilde{A}(x_k)]+2C\tilde{A}(x_{k+1}),C(k+1)\tilde{A}(x_{k+1})-rz_k+rx^*}\\
\end{align*}

Also, by
\[
\frac{1}{2}\norm{a}^2-\frac{1}{2}\norm{b}^2=\inner{a-b,b}+\frac{1}{2}\norm{a-b}^2,
\]

we have
\begin{align*}
\text{IV}=&-\underbrace{C(r^2-2r)\inner{\tilde{A}(x_{k+1}),z_k-x^*}}_{\text{IV}_1}+\underbrace{\frac{C^2(r-2)}{2}\norm{\tilde{A}(x_{k+1})}^2}_{\text{IV}_2}.
\end{align*}

For II, we have
\begin{align*}
\text{II}=&Crk\inner{\tilde{A}(x_{k+1}),x_{k+1}-x^*}-Crk\inner{\tilde{A}(x_k),x_k-x^*}+Cr\inner{\tilde{A}(x_{k+1}),x_{k+1}-x^*}\\
\end{align*}

Step 2: Estimate the upper bound of $\E(k+1)-\E(k)$.
First, consider $\text{I}_1+\text{II}_1+\text{III}_1$.

\begin{align*}
\text{I}_1+\text{II}_1+\text{III}_1=&Ck\inner{\tilde{A}(x_{k+1})-\tilde{A}(x_k),r(x_{k+1}-z_k)+(k+r)\tilde{A}(x_{k+1})}\\
\end{align*}

\[
\text{I}_1+\text{II}_1+\text{II}_2+\text{III}_1=-C(k^2+rk)\inner{\tilde{A}(x_{k+1})-\tilde{A}(x_k),x_{k+1}-x_k}.
\]

Now, consider $\text{II}_3$. By (<ref>), we have
\[
\text{II}_3=Cr\inner{\tilde{A}(x_{k+1}),r(z_k-x_{k+1})-(k+r)\tilde{A}(x_{k+1})}.
\]

\begin{align*}
\text{II}_3+\text{III}_2+\text{IV}_1=&Cr^2\inner{\tilde{A}(x_{k+1}),z_k-x_{k+1}}-Cr(k+r)\norm{\tilde{A}(x_{k+1})}^2-Cr^2\inner{\tilde{A}(x_{k+1}),z_k-x^*}\\
\end{align*}

In conclusion,
\begin{align*}
\E(k+1)-\E(k)\leqslant&-C(k^2+rk)\inner{\tilde{A}(x_{k+1})-\tilde{A}(x_k),x_{k+1}-x_k}-C(r^2-r)\inner{\tilde{A}(x_{k+1}),x_{k+1}-x^*}\\
\end{align*}

Step 3: Deduce the convergence rates. If $r\geqslant 2, 0<C\leqslant 1$, then all parts of $\E(k)$ are non-negetive, and $Crk\geqslant (C+C^2)k, Cr^2\geqslant\dfrac{C(r+1)+C^2r}{2}$, which implies
\[
\E(k+1)-\E(k)\leqslant 0.
\]

By Proposition <ref>, we have
\begin{align*}
\frac{Ck[k+3r-C(k+1)]}{2}\norm{\tilde{A}(x_{k})}^2\leqslant\E(k)\leqslant\E(0)=\frac{r^3-r^2}{2}\norm{x_0-x^*}^2.
\end{align*}

Diving $\dfrac{1}{2}Ck[k+3r-C(k+1)]$ on both sides of above inequality, we have
\[
\norm{\tilde{A}(x_k)}^2\leqslant\frac{(r^3-r^2)\norm{x_0-x^*}^2}{Ck[k+r-C(k+1)]}.
\]

Also we have

\[Crk\inner{\tilde{A}(x_k),x_k-x^*}\leqslant\E(k)\leqslant\E(0)=\frac{r^3-r^2}{2}\norm{x_0-x^*}^2.\]

Diving $Crk$ on both sides of the above inequality, then we receive the desired inequality.

It is worth noting that the fact $\E(k+1)-\E(k)\leqslant 0$ in the proof of Theorem <ref> relies only on that $\tilde{A}$ is monotone, which prompts us to consider replacing $\tilde{A}$ in (<ref>) by $A$ when applying symplectic discretization techinique. The corresponding discretization form is summerized as followed.

Symplectic proximal point algorithm, SPPA

Initialize $x_0, z_0=x_0$
Termination conditions aren't hold

Output $x_{k+1}$.

The convergence rates of Algorithm <ref> is given as followed.

Let $\{x_k\}, \{z_k\}$ be the sequences generated by
\begin{align*}
\tilde{x}_{k+1}&=\dfrac{r}{k+r}z_k+\frac{k}{k+r}x_k,\\
\end{align*}

If $0<C\leqslant 1, r\geqslant 2$, then we have
\begin{align*}
\tilde{x}_{k+1}-x_{k+1}&\in A(x_{k+1}),\\
\norm{\tilde{x}_k-x_{k}}^2&\leqslant\frac{(r^3-r^2)\norm{x_0-x^*}^2}{Ck[k+3r-C(k+1)]},\quad\forall k\geqslant 1,\\
\inner{\tilde{x}_k-x_k,x_k-x^*}&\leqslant\frac{(r^3-r^2)\norm{x_0-x^*}^2}{2Ck},\quad\forall k\geqslant 1.
\end{align*}

The proof of Corollary <ref> is just to replace $\tilde{A}(x_k)$ in the Lyapunov function in the proof of Theorem <ref> with $\tilde{x}_{k+1}-x_{k+1}$, and the remaining part of the proof is similar to the proof of Theorem <ref>.

In the end, we demonstrate that the constrain of parameter $C$ is necessary. If $C=r$, then Algorithm <ref> is equivalent to proximal point algorithms (<ref>).

§ APPLICATIONS OF SYMPLECTIC PROXIMAL POINT ALGORITHMS AND ITS NUMERICAL EXPERIMENTS

§.§ Symplectic Augmented Lagrangian Multipliers Methods

In section <ref>, we have derived symplectic augmented Lagrangian multipliers methods via symplectic proximal point algorithms. In this section, we test numerical performance of symplectic augmented Lagrangian multipliers methods on LASSO problems
\[
\min_x F(x)=\frac{1}{2}\norm{Ax-b}^2+\mu\norm{x}_1.
\]

To apply (<ref>), we first transform above unconstrained optimization problems into the following constrained optimization problems

\begin{equation}
\label{eq:lasso}
\begin{gathered}
\min_{x,y} \frac{1}{2}\norm{Ax-b}^2+\mu\norm{y}_1,\\
s.\ t. \quad x=y.
\end{gathered}
\end{equation}

Applying (<ref>) with $H=I$, we obtain the following recursive rule.
\begin{equation}
\label{eq:lassosymplectic1}
\begin{split}
\tilde{\lambda}_{k+1}&=\frac{1}{b_k+1}z_k+\frac{b_k}{b_k+1}\lambda_k,\\
\lambda_{k+1}&=\tilde{\lambda}_{k+1}-\frac{c_k}{b_k+1}(x_{k+1}-y_{k+1}),\\
\end{split}
\end{equation}

(<ref>) can be transformed into
\begin{equation}
\begin{split}
\tilde{\lambda}_{k+1}&=\frac{1}{b_k+1}z_k+\frac{b_k}{b_k+1}\lambda_k,\\
\lambda_{k+1}&=\tilde{\lambda}_{k+1}-\frac{c_k}{b_k+1}(x_{k+1}-y_{k+1}),\\
\end{split}
\label{eq:lassosymplectic2}
\end{equation}

Above $y_{k+1}(x)=\mathrm{soft}\left(\dfrac{b_k+1}{c_k}\tilde{\lambda}_{k+1}+x,\dfrac{\mu(b_k+1)}{c_k}\right)$, here $\mathrm{soft}(x,\mu)=\argmin_y \dfrac{1}{2}\norm{y-x}^2+\mu\norm{y}_1$ is the soft-threshold operator. In our experiment, the matrix $A$ and the initial point $\lambda_0$ are randomly generated, also $A\in\mathbb{R}^{100\times 200}$. The penalty parameter $\rho$ is seted to be 10. As a comparison, we include augmented Lagrangian multipliers methods in our experiment. The optimization subproblem in (<ref>) is solved by using FISTA. We run 1000 iterations of symplectic augmented Lagrangian multipliers methods with parameters $a_k=\dfrac{k+1}{4}, b_k=\dfrac{k}{2}, c_k=\dfrac{k+2}{2}$.

Evolution of $ F(x_k)-F(x^*) $ and $\norm{x_k-y_k}$ versus iterations on LASSO problems , where the optimal solution $x^*$ is generated by using FISTA.

As we can see from Figure <ref>, symplectic augmented Lagrangian multipliers method numerically perform better than augmented Lagrangian multipliers method. Also, symplectic augmented Lagrangian multipliers method exhibits oscillations when iteration grows large. To address oscillation phenomenon, we restart symplectic augmented Lagrangian multipliers method every 50 iterations. Figure <ref> shows that with restarting, oscillation phenomenon is alleviated and convergence results become better as iteration grows large. Moreover, as iteration grows large, symplectic augmented Lagrangian multipliers methods with restarting exhibits linear convergence rate.

§.§ Symplectic Douglas-Rachford Algorithms

Douglas-Rachford algorithm proposed by Lions and Mercier in 1979 is an unconditionally convergent algorithms for solving the following zero-point problems.
\begin{equation}\label{eq:drmain}
0\in A(x)+B(x),
\end{equation}

where $A, B$ are maximally monotone operators. The standard Douglas-Rachford algotithms is
\begin{equation}\label{eq:dr}
x_{k+1}=T_{\rho A,\rho B}(x_k).
\end{equation}

Here $T_{\rho A,\rho B}$ is the Douglas-Rachford operator defined as $T_{\rho A,\rho B}=\dfrac{1}{2}I+\dfrac{1}{2}(2J_{\rho A}-I)\circ(2J_{\rho B}-I)$. [11] proved that $T_{\rho A,\rho B}$ is resolvent of a maximally monotone operator. Thus we can replace resolvent of $A$ in Algorithm <ref> by $T_{\rho A,\rho B}$, and obtain the following symplectic Douglas-Rachford algorithms.

\begin{equation}\label{eq:symplecticdr}
\begin{split}
\tilde{x}_{k+1}&=\dfrac{r}{k+r}z_k+\dfrac{k}{k+r}x_k,\\
u_{k+1}&=J_{\rho B}(\tilde{x}_k),\\
v_{k+1}&=J_{\rho A}(2u_{k+1}-\tilde{x}_k),\\
\end{split}
\end{equation}

§.§ Symplectic Alternating Direction Method of Multipliers

In this section, we consider the following convex optimization problems
\begin{equation}\label{eq:admmmain}
\begin{gathered}
\min_{x,y} f(x)+g(y),\\
s.\ t. \quad Ax+By=c.
\end{gathered}
\end{equation}

Alternating direction method of multipliers(ADMM) may be the most populer algorithms for solving (<ref>). Theoretically, ADMM is equivalent to applying (<ref>) to Lagrangian duality of (<ref>). The Lagrangian duality of (<ref>) is given by
\[
\min_\lambda f^*(-A\tr \lambda)+g^*(-B\tr\lambda)+\inner{c,\lambda}.
\]

The symplectic Douglas-Rachford algorithm (<ref>) with $A=\partial(f^*\circ -A\tr)+c, B=\partial(g^*\circ-B\tr)$ is given as followed.

\begin{align*}
\tilde{u}_{k+1}&=\frac{r}{k+r}z_k+\frac{k}{k+r}u_k,\\
\lambda_{k+1} &=\argmin_\lambda f^*(-A\tr\lambda)+\inner{c,\lambda}+\frac{1}{2\rho}\norm{\lambda-\tilde{u}_{k+1}}^2, \\
v_{k+1} &=\argmin_v g^*(-B\tr v)+\frac{1}{2\rho}\norm{v-2\lambda_{k+1}+\tilde{u}_{k+1}}^2,\\
\end{align*}

By utilizing Fenchel-Rockafellar duality theory, we obtain the following recursive rule.

\begin{equation}
\label{eq:symplecticadmm}
\begin{split}
\tilde{u}_{k+1}&=\frac{r}{k+r}z_k+\frac{k}{k+r}u_k,\\
x_{k+1}&=\argmin_x f(x)+\inner{\tilde{u}_{k+1},Ax-c}+\frac{\rho}{2}\norm{Ax-c}^2,\\
\lambda_{k+1}&=\tilde{u}_{k+1}+\rho(Ax_{k+1}-c),\\
y_{k+1}&=\argmin_y g(y)+\inner{2\lambda_{k+1}-\tilde{u}_{k+1},By}+\frac{\rho}{2}\norm{By}^2,\\
v_{k+1}&=2\lambda_{k+1}-\tilde{u}_{k+1}+\rho By_{k+1},\\
\end{split}
\end{equation}

There are two ways to simplify (<ref>). By representing $\lambda_{k+1}, v_{k+1}$ by $x_{k+1}, y_{k+1}$ in (<ref>), we obtain the following recursive rule.
\begin{equation}
\label{eq:symplecticadmm1}
\begin{split}
\tilde{u}_{k+1}&=\frac{r}{k+r}z_k+\frac{k}{k+r}u_k,\\
x_{k+1}&=\argmin_x f(x)+\inner{\tilde{u}_{k+1},Ax-c}+\frac{\rho}{2}\norm{Ax-c}^2,\\
y_{k+1}&=\argmin_y g(y)+\inner{\tilde{u}_{k+1},By}+\frac{\rho}{2}\norm{2(Ax_{k+1}-c)+By}^2,\\
\end{split}
\end{equation}

Also by representing $u_{k+1}, \tilde{u}_{k+1}, v_{k+1}$ by $x_{k+1}, y_{k+1}$ and reordering the recursive rule, we obtain the following iteration rule.

\begin{equation}
\label{eq:symplecticadmm2}
\begin{split}
x_{k+1}&=\argmin_x f(x)+\inner{\lambda_k,Ax+By_k-c}+\frac{\rho}{2}\norm{Ax+By_k-c}^2\\
\tilde{\lambda}_{k+1}&=\frac{r}{k+r}z_{k+1}+\frac{k}{k+r}\lambda_k\\
y_{k+1}&=\argmin_y g(y)+\inner{\tilde{\lambda}_{k+1},Ax_{k+1}+By-c}+\frac{\rho}{2}\norm{\frac{k}{k+r}(Ax_{k+1}-c)+By}^2,\\
\lambda_{k+1}&=\tilde{\lambda}_{k+1}+\frac{k\rho}{k+r}(Ax_{k+1}-c)+\rho By_{k+1}.
\end{split}
\end{equation}

Both (<ref>) and (<ref>) can be seem as instances of symplectic ADMM(S-ADMM). In the following numerical experiments, we test numerical preformance of (<ref>) rather than (<ref>) because (<ref>) is easier to impletement than (<ref>). By Corollary <ref>, we have the following convergence results of (<ref>).

Let $\{x_k\}, \{y_k\}, \{z_k\}, \{u_k\}, \{\tilde{u}_k\}$ be the sequences generated by (<ref>). If $0<C\leqslant 1, r\geqslant 2$, then we have
\[
\norm{Ax_k+By_k-c}^2\leqslant\frac{(r^3-r^2)\norm{x_0-x^*}^2}{\rho^2Ck[k+3r-C(k+1)]},\quad\forall k\geqslant 1.
\]

Since $u_k-\tilde{u}_k=\rho(Ax_k+By_k-c)$, then we have
\[
\rho^2\norm{Ax_k+By_k-c}^2\leqslant\frac{(r^3-r^2)\norm{x_0-x^*}^2}{Ck[k+3r-C(k+1)]},\quad\forall k\geqslant 1.
\]

Diving $\rho^2$ on above inequality, we obtain desired inequality.

We apply (<ref>) to the problem
\begin{equation}
\begin{gathered}
\min_{x,y} \norm{x}_1+I(Ay=b),\\
s. t. \quad x=y,
\end{gathered}
\end{equation}
which is associated with basis pursuit problem
\begin{equation}
\begin{gathered}
\min_x \norm{x}_1,\\
s. t. \quad Ax=b.
\end{gathered}
\end{equation}

Let $f(x)=\norm{x}_1, g(y)=I(Ay=b), A=I, B=-I, c=0$, the symplectic ADMM schemes is given as followed.
\begin{align*}
\tilde{u}_{k+1}&=\frac{r}{k+r}z_k+\frac{k}{k+r}u_k,\\
\end{align*}

In our experiment, the matrix $A$ and the vector $b$ are randomly generated, also $A\in\mathbb{R}^{100\times 200}, b\in\mathbb{R}^{100\times 1}$. The ADMM method is comparable in this experiment. After running 5000 iterations of symplectic ADMM with parameters $\rho=10, r=2, C=1$, we plot $\norm{x_k-y_k}$ in the following Figure <ref>.

Evolution of $ \norm{x_k-y_k} $ versus iterations on basis pursuit problems.

Figure <ref> shows that, compared to ADMM, numerical preformance of symplectic ADMM is better and oscillation phenomenon of symplectic ADMM is lighter. In addition, if we restart symplectic ADMM every 50 iterations, the algorithms exhibits faster convergence results as iteration grows large, and linear convergence rate can be observed.

We consider fused LASSO problem
\[
\min_x F(x)=\frac{1}{2}\norm{Ax-b}^2+\mu_1\norm{x}_1+\mu_2\sum_{i=1}^{d-1}|x_{i+1}-x_i|.
\]

Here $x\in\mathbb{R}^d, x_i$ is $i-$th entry of $x$. Fused Lasso problem can be simplified as

\begin{equation}
\begin{gathered}
\min_x F(x)=\frac{1}{2}\norm{Ax-b}^2+\mu_1\norm{x}_1+\mu_2\norm{Dx}_1,
\end{gathered}
\end{equation}

where $D$ is given as
\[
\vdots&\ddots&\ddots&\ddots&\ddots&\vdots\\
\vdots&\ddots&0&1&-1&0\\
\end{pmatrix}.
\]

To apply (<ref>), we transform the fused LASSO problem into
\begin{equation}
\begin{gathered}
\min_{x^1,x^1,x^1,y} \frac{1}{2}\norm{Ax^1-b}^2+\mu_1\norm{x^2}_1+\mu_2\norm{Dx^3}_1\\
s. t. \begin{pmatrix}
\end{pmatrix}\begin{pmatrix}
\end{pmatrix}-\begin{pmatrix}
\end{pmatrix}y=0.
\end{gathered}
\end{equation}

By letting $x=(x^1,x^1,x^1), f(x)=\dfrac{1}{2}\norm{Ax^1-b}^2+\mu_1\norm{x^1}_1+\mu_2\norm{Dx^3}_1, g(y)=0$, $A=\begin{pmatrix}
\end{pmatrix}, B=-\begin{pmatrix}
\end{pmatrix},$ $ c=0$, we have the following symplectic ADMM schemes.

\begin{align*}
\tilde{u}^1_{k+1}&=\frac{r}{k+r}z^1_k+\frac{k}{k+r}u^1_k,\\
\tilde{u}^2_{k+1}&=\frac{r}{k+r}z^2_k+\frac{k}{k+r}u^2_k,\\
\tilde{u}^3_{k+1}&=\frac{r}{k+r}z^3_k+\frac{k}{k+r}u^3_k,\\
x^1_{k+1}&=(A\tr A+\rho I)^{-1}(A\tr b+y_k-\tilde{u}^1_{k+1}),\\
y_{k+1}&=(\rho D\tr D+2\rho I)^{-1}(\tilde{u}_{k+1}^1+\tilde{u}^2_{k+1}+D\tr\tilde{u}^3_{k+1}+\rho x^1_{k+1}+\rho x^2_{k+1}+\rho D\tr x^3_{k+1}),\\
\end{align*}

In our experiment, the matrix $A$, vector $b$ and the initial point $u^1_0, u^2_0, u^3_0$ are randomly generated, also $A\in\mathbb{R}^{100\times 200}, b\in\mathbb{R}^{100\times 1}$. The parameters $\mu_1, \mu_2$ are choosen to be $5, 10$. Also, ADMM is a comparator in our experiment. We run 20000 iterations of symplectic augmented Lagrangian multipliers methods with $\rho, C, r$ are chosen to be $1, 1, 2$ and plot our numerical results in Figure <ref>.

Evolution of $F(y_k)-F^*$ and primal residual $\norm{x^1_k-y_k}^2+\norm{x^2_k-y_k}^2+\norm{x^3-Dy_k}^2$, where $F^*=F(y')$ and $y'$ is generated by recursively iterate ADMM until primal residual less than $10^{-12}$.

As we can see from Figure <ref>, at the beginning of the iteration, symplectic ADMM did not perform as well as ADMM. But with increasing number of iterations, symplectic ADMM converge faster than ADMM. Compared to symplectic ADMM, symplectic ADMM with restarting can slightly improve numerical performance.

§.§ Symplectic Primal-dual Hybrid Gradient Methods

The primal-dual hybrid gradient methods, also known as the Chambolle-Pock methods, is a powerful first-order methods to solve the following convex-concave saddle point problem
\begin{equation}
\min_x\max_y F(x,y)=f(x)+\inner{x,Ay}-g(y).
\label{eq:pdhgmain}
\end{equation}

Here $f, g$ are two closed proper convex functions. In [7], they shows that primal-dual hybrid gradient methods are a preconditioned proximal point methods for the maximally monotone operator $\partial F$. In addition, the equivalence between primal-dual hybrid gradient methods and Douglas-Rachford algorithms is been shown in [17]. These two previous works inspire us to apply Algorithm <ref> or symplectic Douglas-Rachford algorithms (<ref>) to (<ref>) as below.

\begin{equation}
\label{eq:symplecticpdhg}
\begin{split}
\tilde{u}_{k+1} &=(\tilde{x}_{k+1},\tilde{y}_{k+1})=\frac{r}{k+r}z_k+\frac{k}{k+r}u_k \\
x_{k+1}&=\argmin_x f(x)+\inner{x,A\tilde{y}_{k+1}}+\frac{1}{2\mu_1}\norm{x-\tilde{x}_{k+1}},\\
y_{k+1}&=\argmin_y g(y)-\inner{2(x_{k+1}-\tilde{x}_{k+1}),Ay}+\frac{1}{2\mu_2}\norm{y-\tilde{y}_{k+1}},\\
\end{split}
\end{equation}

Above (<ref>) is the symplectic primal-dual hybrid gradient method. By Theorem <ref>, the convergence results of (<ref>) is described as followed.
Let $\{u_k\}, \{\tilde{u}_k\}, \{x_k\}, \{y_k\}, \{z_k\}$ be the sequences generated by (<ref>). If $ 0<C\leqslant 1, r\geqslant 2, \mu_1\mu_2\norm{A}^2< 1 $, then we have
\[
\norm{u_k-\tilde{u}_k}^2\leqslant\frac{(r^3-r^2)\norm{x_0-x^*}^2}{Ck[k+3r-C(k+1)]},\quad\forall k\geqslant 1.
\]

The following min-max problems
\begin{equation}
\label{eq:nashequilibrium}
\min_{x\in\Delta_m}\max_{y\in\Delta_n}\inner{a,x}+\inner{x,Ay}-\inner{b,y}
\end{equation}

is what we need to test next. Here $\Delta_m=\{x\in\mathbb{R}^m|\sum\limits_{i=1}^mx_i=1, x\geqslant 0\}$ is the $m-1$dimensional unit simplex. (<ref>) can be converted into
\begin{equation}
\min_{x}\max_{y} I_{\Delta_m}(x)+\inner{a,x}+\inner{x,Ay}-\inner{b,y}-I_{\Delta_n}(y).
\end{equation}

Let $f(x)=I_{\Delta_m}(x)+\inner{a,x}, g(y)=\inner{b,y}+I_{\Delta_n}(y)$, the corresponding symplectic primal-dual hybrid gradient methods is
\begin{align*}
\tilde{u}_{k+1} &=(\tilde{x}_{k+1},\tilde{y}_{k+1})=\frac{r}{k+r}z_k+\frac{k}{k+r}u_k, \\
x_{k+1}&=\argmin_{x\in I_{\Delta_m}} \inner{a,x}+\inner{x,A\tilde{y}_{k+1}}+\frac{1}{2\mu_1}\norm{x-\tilde{x}_{k+1}}^2\\
y_{k+1}&=\argmin_y \inner{b,y}-\inner{2x_{k+1}-\tilde{x}_{k+1},Ay}+\frac{1}{2\mu_2}\norm{y-\tilde{y}_{k+1}}^2\\
\end{align*}

Before the numerical experiment, we need to mention that the projection onto unit simplex can be solved exactly in polynomial time. In our experiment, the matrix $A$, the vector b and the initial point $u_0$ are randomly generated. The parameters $\mu_1, \mu_2, C, r$ are chosen to be $0.99/\norm{A}, 0.99/\norm{A}, 1, 2$. In comparision, primal-dual hybrid gradient methods is also tested in our experiment. Since (<ref>) can be altered into the following constrained optimization
\[
\min_{x\in\Delta_m} \inner{a,x}+\max_{i=1,\cdots,n}A_ix-b,
\]

where $A_i$ is $i-$row of $A$, and $x_{k+1}\in\Delta_m$, $\max_y F(x_k,y)-F(x^*,y^*)$ can measure the speed of convergence, where $(x^*,y^*)$ is the saddle point of $F$. Also, we let $\overline{A}(u_k)=u_k-\tilde{u}_k$ in symplectic primal-dual hybrid gradient methods, $\overline{A}(u_k)=u_k-u_{k-1}$ in primal-dual hybrid gradient methods, both of which are measurements of convergence. After running 10000 iterations, the quantative results are shown in Figure <ref>.

Evolution of $\norm{\tilde{A}(u_k)}^2$ and $\max_y F(x_k,y)-F(x^*,y^*)$ versus iterations on (<ref>), where $(x^*,y^*)$ is the $(x_k,y_k)$ generated by recursive iterate symplectic primal-dual hybrid gradient methods until $\norm{\overline{A}(u_k)}^2\leqslant 10^{-12}$.

Figure <ref> demonstrates that symplectic primal-dual hybrid gradient method are faster than primal-dual hybrid methods. Unlike previous numerical experiments, symplectic primal-dual hybrid gradient methods with restarting does not speed up the convergence, but is still faster than primal-dual hybrid gradient methods.

§ CONCLUSION AND DISCUSSION

Throughout this paper, we obtain high-resolution ODEs for accelerated proximal point algorithms respect to both closed proper convex functions and maximally monotone operators by modifying high-resolution ODEs for accelerated mirror descent. Also by investigating convergence results of proximal point algorithms and existing accelerated proximal point algorithms, we proposed Lyapunov functions framework in both continuous-time and discrete-time. By using Lyapunov functions, we proved that the high-resolution ODEs for proximal point algorithms and its corresponding symplectic discretization scheme admit accelerated phenomenon.

There are also many theoretical issues about accelerated proximal point algorihthms worth investigating. To obtain accelerated schemes, the parameter $r$ in (<ref>) and (<ref>) must be not less than $2$. Also the accelerated phenomenon exists for Nesterov's accelerated gradient methods when $r\geqslant 2$. Thus, it is interesting to search the effect of $2$ in accelerated schemes. In [2, 9], they discussed convergence rates of accelerated gradient methods and FISTA with $0<r<2$. The convergence results of (<ref>) and (<ref>) with $0<r<2$ are worthy of discussion.

For practice, there are some valuable problems. Theorem <ref> and Corollary <ref> tells us that the upper bound of $\inner{\tilde{A}(x_k),x_k-x^*}$ and $\inner{\tilde{x}_k-x_k,x_k-x^*}$, which are interpreted as objective function in finding zeros of $A$ in section <ref>, are monotone decreasing respect to $C$. This fact suggests a conjecture that the symplectic proximal point algorithms converge faster as the parameter $C$ grows larger. Numerical experiments give the basis for this conjecture.

Evolution of $F(y_k)-F^*$ and $\norm{x_k-y_k}$ versus iterations on LASSO problems.

Figure <ref> illustrates that when $C=1, r=2 $, the symplectic ADMM does not preform as well as ADMM, which can be seem as symplectic ADMM with $C=r=2$. That is because the convergence rates of ADMM is $O\left(\dfrac{1}{k^2}\right)$ when it is applied to LASSO problem. In contrast, when $C=16, r=2$, the symplectic ADMM converge significantly faster than ADMM. But in section <ref>, we present two numerical experiments that symplectic ADMM with $C=1, r=2$ converge faster than ADMM. Also in section <ref>, the symplectic primal-dual hybrid gradient methods converge faster than primal-dual hybrid gradient methods. A problem arise natually.
For what kind of operators does symplectic proximal point algorithms converge faster than proximal point algorithms.
Also, the following problem is interesting in practice.
If proximal point algorithms converge faster than symplectic proximal point algorithm, what is the maximum value of $C$ respect to specific maximally monotone operator $A$ that gurantee $O\left(\dfrac{1}{k^2}\right)$ convergence rates of symplectic proximal point algorithms?

Determine maximum value of $C$ respect to a great many of maximally monotone operators may cost too much time and effort. Develop an algorithms which produces as large $C$ that maintain accelerated nature of Algorithm <ref> as possible is vital in practice.
[1] F. Alvarez, H. Attouch, J. Bolte and P. Redont, A second-order gradient-like dissipative dynamical system with hessian-driven damping. Application to optimization and mechanics, Journal de mathématiques pures et appliquées, 81 (2002), pp. 747-779.
[2]	H. Attouch, Z. Chbani and H. Riahi, Rate of convergence of the Nesterov accelerated gradient method in the subcritical case $\alpha\leqslant 3$, ESAIM: Control, Optimisation and Calculus of Variations, 25 (2019), pp. 2.
[3]	H. attouch and S. C. László, Newton-like inertial dynamics and proximal algorithms governed by maximally monotone operators, SIAM Journal on Optimization, 30 (2020), pp. 3252-3283.
[4]H. attouch and J. Peypouquet, Convergence of inertial dynamics and proximal algorithms governed by maximally monotone operators, Mathematical Programming, 174 (2019), pp. 391-432.
[5]H. H. Bauschke and P. L. Combettes, Correction to: Convex Analysis and Monotone Operator Theory in Hilbert Spaces, Springer, Cham, 2017.
[6] J. F. Bonnans and A. Shapiro, Perturbation analysis of optimization problems, Springer, New York, 2013.
[7]A. Chambolle and T. Pock, On the ergodic convergence rates of a first-order primal–dual algorithm, Mathematical Programming, 159 (2016), pp. 253-287.
[8]	T. Chavdarova, M. I. Jordan and M. Zampetakis, Last-iterate convergence of saddle point optimizers via high-resolution differential equations, arXiv preprint arXiv:2112.13826, 2021.
[9]	S. Chen, B. Shi and Y.-X. Yuan, On Underdamped Nesterov's Acceleration, arXiv preprint arXiv:2304.14642, 2023.
[10]	P. L. Combettes, Monotone operator theory in convex optimization, Mathematical Programming, 170 (2018), pp. 177-206.
[11]	J. Eckstein and D. P. Bertsekas, On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators, Mathematical programming, 55 (1992), pp. 293-318.
[12]	K. Feng and M. Qin, Symplectic Geometric Algorithms for Hamiltonian Systems, Springer, Berlin, 2010.

[13]	B. Halpern, Fixed Points of Nonexpanding Maps. (1967). Bulletin of the American Mathematical Society,  73(6), pp. 957-961.
[14]	D. KIM, Accelerated proximal point method for maximally monotone operators, Mathematical Programming, 190 (2021), pp. 57-87.
[15]	B. Li, B. Shi and Y.-X. Yuan, Proximal Subgradient Norm Minimization of ISTA and FISTA, arXiv preprint arXiv:2211.01610, 2022.
[16]Y. Nesterov, Smooth Minimization of Non-smooth Functions, Mathematical programming, 103 (2005), pp. 127-152.
[17]D. O'Connor and L. Vandenberghe, On the equivalence of the primal-dual hybrid gradient method and Douglas–Rachford splitting, Mathematical Programming, 179 (2020), pp. 85-108.
[18]	N. Parikh and S. Boyd, Proximal algorithms, Foundations and trends® in Optimization, 1 (2013), pp. 127-239.
[19]	R. ROCKAFELLAR, Monotone operators and augmented Lagrangian methods in nonlinear programming, Nonlinear Programming 3, 1978, pp. 1-25.
[20]	R. T. Rockafellar, Convex Analysis, Princeton University Press, Princeton, 1997.
[21]	R. T. Rockafellar, Monotone Operators and the Proximal Point Algorithm, SIAM Journal on Control and Optimization, 14 (1976), pp. 877-898.
[22]E. Ryu and W. Yin, Large-scale convex optimization via monotone operators, Cambridge University Press, Cambridge, 2022.
[23]	B. Shi, S. S. Du, M. I. Jordan and W. J. Su, Understanding the acceleration phenomenon via high-resolution differential equations, Mathematical Programming, 195 (2021), pp. 79-148.
[24]	B. Shi, S. S. Du, W. Su and M. I. Jordan, Acceleration via symplectic discretization of high-resolution differential equations, Advances in Neural Information Processing Systems, 32 (2019).
[25]	W. Su, S. Boyd and E. Candes, A Differential Equation for Modeling Nesterov’s Accelerated Gradient Method: Theory and Insights, Journal of Machine Learning Research, 17 (2016) , pp. 1-43.
[26]Y.-X. Yuan and Y. Zhang, Analysis Accelerated Mirror Descent via High-resolution ODEs, arXiv preprint arXiv:2308.03242,2023.