An Optimal ADMM for Unilateral Obstacle Problems

Abstract

We propose a new alternating direction method of multipliers (ADMM) with an optimal parameter for the unilateral obstacle problem. We first use the five-point difference scheme to discretize the problem. Then, we present an augmented Lagrangian by introducing an auxiliary unknown, and an ADMM is applied to the corresponding saddle-point problem. Through eliminating the primal and auxiliary unknowns, a pure dual algorithm is then used. The convergence of the proposed method is analyzed, and a simple strategy is presented for selecting the optimal parameter, with the largest and smallest eigenvalues of the iterative matrix. Several numerical experiments confirm the theoretical findings of this study.

1. Introduction

As we know, many variant versions methods of the alternating direction method of multipliers (ADMM) have been successfully applied to different models such as constrained convex optimization problems [1,2,3,4,5], obstacle problems [6,7], contact problems [8,9,10], nonlinear complementarity problems [11,12], etc. Since the penalty parameter plays an important role for this method, the choice of the parameter is crucial for the ADMM. However, the optimal parameter is not easily obtainable in practical applications. Therefore, a good strategy for selecting the optimal parameter is very important for the ADMM [10,13].

Recently, some optimal parameter selections of the ADMM have been considered for convex optimization problems with linear equality constraints [14,15,16] and inequality constraints [14,17,18]. In addition, the automatic parameter selection for the ADMM was designed for the unilateral contact problem [8]. This method derives a pure dual form of the ADMM to eliminate the auxiliary and primal unknowns. This method obtains the optimal parameter using the smallest and largest eigenvalues of the iterative matrix. Unfortunately, the optimal parameter is not easily obtainable in practice because of the complex computing and the high computation cost.

The unilateral obstacle problem has many applications, such as lubrication phenomena, fluid flow in porous media and option pricing, etc. [19,20]. In this paper, we propose a new ADMM [7,21] with an optimal parameter for the unilateral obstacle problem [22,23,24,25,26,27,28,29,30,31], which is based on the extreme eigenvalues of the matrix using the finite difference method (FDM) [6,22,26,32]. The main advantage to this method is the simplicity of selecting the optimal parameter [8,33,34,35,36,37]. Moreover, the method saves computational time for large-scale problems because the iterative matrix is always kept constant.

The rest of this paper is structured as follows. We present some basic results in Section 2 for the unilateral obstacle problem. In particular, we developed a new ADMM and show the convergence of the method in Section 3, and we present an optimal parameter selection scheme in detail. In Section 4, some numerical results are presented to show the performance of the method. Some conclusions and perspectives are presented in Section 5.

2. Unilateral Obstacle Problems and ADMM Algorithm in Finite Dimension

For simplicity, we consider a square domain $\mathrm{\Omega }$ in ${\mathbb{R}}^2$, and $\mathrm{\Gamma }=\partial \mathrm{\Omega }$ is its boundary. For given functions $f\left(x\right)$ and $g\left(x\right)$ and an obstacle $\psi \left(x\right)$, the unilateral obstacle problem can be described as follows: find $u\left(x\right)$ such that

$$\left\{\begin{array}{c}-\Delta u\left(x\right)\ge f\left(x\right),x\mathrm{in}\mathrm{\Omega },\hfill \\ u\left(x\right)\ge \psi \left(x\right),x\mathrm{in}\mathrm{\Omega },\hfill \\ \left(u\right(x)-\psi (x\left)\right)(-\Delta u(x)-f(x\left)\right)=0,x\mathrm{in}\mathrm{\Omega },\hfill \\ u\left(x\right)=g\left(x\right),x\mathrm{on}\mathrm{\Gamma },\hfill \end{array}\right.$$

(1)

where $\Delta$ is the Laplace operator. Then, obstacle problem (1) admits a unique solution u [13].

To approximate problem (1) using the FDM, we divide $\mathrm{\Omega }$ by a regular mesh, and then the mesh consists of boundary points on $\partial G$ and $n^2$ interior points in G. Then, we obtain a linear complementarity problem (LCP):

$$\left\{\begin{array}{c}-{\Delta }_h{u}_h\left(x\right)\ge f\left(x\right),x\mathrm{in}G,\hfill \\ u_h\left(x\right)\ge \psi \left(x\right),x\mathrm{in}G,\hfill \\ (u_h\left(x\right)-\psi \left(x\right))(-{\Delta }_h{u}_h\left(x\right)-f\left(x\right))=0,x\mathrm{in}G,\hfill \\ u_h\left(x\right)=g\left(x\right),x\mathrm{on}\partial G,\hfill \end{array}\right.$$

where $u_h\left(x\right)$ is an approximate solution to $u\left(x\right)$ at the grid points x, and ${\Delta }_h$ is a difference operator that approximates $\Delta$.

In particular, we use the following five-point difference scheme in $\mathrm{\Omega }$ to discrete the first inequality (1):

$$4u_{ij}-u_{i-1,j}-u_{i+1,j}-u_{i,j-1}-u_{i,j+1}\ge h^2{f}_{i,j},i,j=1,2,\cdots ,n,$$

(2)

where h is the mesh width. We use (2) and substitute the known values of $u_{ij}$ and $f_{ij}$ on the boundary grid and domain, grid respectively. Tthe LCP reads as follows:

$$\left\{\begin{array}{c}u\ge \psi ,\hfill \\ Au+b\ge 0,\hfill \\ {(u-\psi )}^T(Au+b)=0,\hfill \end{array}\right.$$

(3)

where u is the vector N ($N=n^2$) of unknown interior mesh values, b is the given data vector, and A is a ${\mathbb{R}}^{N\times N}$ matrix expressed as

$$A=\left[\begin{array}{ccccc}B+2\mathbb{I}& -\mathbb{I}& & & \\ -\mathbb{I}& B+2\mathbb{I}& -\mathbb{I}& & \\ & \ddots & \ddots & \ddots & \\ & & -\mathbb{I}& B+2\mathbb{I}& -\mathbb{I}\\ & & & -\mathbb{I}& B+2\mathbb{I}\end{array}\right]$$

where B = tridiag$(-1,2,-1)\in {\mathbb{R}}^{n\times n}$, and $\mathbb{I}\in {\mathbb{R}}^{n\times n}$ is the identity matrix. It is clear that A is diagonally dominant. Moreover, the matrix A has the following special features:

(a)

A is symmetric, positive-definite.

(b)

A has n diagonal blocks in form $B+2\mathbb{I}$ and subdiagonal blocks $-\mathbb{I}$.

To design efficient algorithms using these features for LCP (3), some iterative algorithms have been presented [6,22,26]. Moreover, LCP (3) is equivalent to the following constrained minimization problem: Find $u\ge \psi$ such that

$$F\left(u\right)\le F\left(v\right),\forall v\ge \psi ,$$

(4)

where $F\left(v\right)=\frac{1}{2}{v}^{T}Av+v^{T}b$. Consequently, we can solve optimization problem (4) and obtain the solution of LCP (3).

To achieve the solution of problem (4) using the ADMM algorithm [7,21], we introduce an indicator function I as follows:

$$I\left(q\right)=\left\{\begin{array}{ccc}0\hfill & \hfill & \mathrm{if}q\ge \psi ,\hfill \\ +\infty \hfill & \hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

We introduce an auxiliary unknown p and consider the following constraint minimization problem:

$$\left\{\begin{array}{c}F\left(u\right)+I\left(p\right)\le F\left(v\right)+I\left(q\right),\forall v,q\in {\mathbb{R}}^N,\hfill \\ \mathrm{subject} \mathrm{to}p-u=0.\hfill \end{array}\right.$$

(5)

Then, (4) is equivalent to problem (5), and we have the following preliminary result [6,21].

 Lemma 1. 

$u\in {\mathbb{R}}^N$ is the solution of (4) if and only if there exists $u\ge \psi$ such that $\{u,p\}$ is the solution of (5).

To solve problem (5), we present the saddle-point problem:

$$L_{\rho }(v,q,\lambda )\ge L_{\rho }(u,p,\lambda )\ge L_{\rho }(u,p,\mu ),\forall \{v,q,\mu \}\in {\mathbb{R}}^N\times {\mathbb{R}}^N\times {\mathbb{R}}^N,$$

(6)

with the augmented Lagrangian $L_{\rho }$ defined by

$$L_{\rho }(u,p,\lambda )=F\left(u\right)+I\left(p\right)+〈\lambda ,p-u〉+\frac{\rho }{2}{\parallel p-u\parallel }^2,$$

(7)

and the penalty parameter $\rho >0$. Let $〈·,·〉$ denote the inner product, and we obtain the following results for problems (5) and (6) [13,21].

 Lemma 2. 

Let $\{u,p,\lambda \}\in {\mathbb{R}}^N\times {\mathbb{R}}^N\times {\mathbb{R}}_{+}^N$ be the solution of (6). Then, u is the solution of (4), and $p=u$; moreover, $\{u,\lambda \}$ is the unique solution of the following problem:

$$\left\{\begin{array}{}\hfill \mathrm{(8)}& Au+b=\lambda ,\hfill \hfill \mathrm{(9)}& \langle \lambda ,u\rangle =0.\hfill \end{array}\right.$$

 Proof. 

From the two inequalities (6) and (7), it follows that

$$〈\lambda -\mu ,u-p〉\le 0,\forall \mu \in {\mathbb{R}}_{+}^N,$$

which leads to

$$p=u.$$

From the first inequality (6), we have

$$L_{\rho }(v,q,\lambda )\ge L_{\rho }(u,p,\lambda )=E\left(u\right),\forall (v,q)\in {\mathbb{R}}^N\times {\mathbb{R}}^N.$$

In taking $q=v$, it follows that

$$E\left(u\right)\le L_{\rho }(v,v,\lambda )=E\left(v\right),\forall v\in {\mathbb{R}}^N.$$

Hence, u is the solution of (4); in taking $q=p$, it follows that

$$Au=f+\lambda -\rho (u-p).$$

Since $p=u$, we obtain (8). We take $v=u$ in the first inequality (6) and have

$$〈\lambda ,q-p〉\ge 0,\forall q\in {\mathbb{R}}^N,$$

It follows from $p=u$ that

$$〈\lambda ,q-u〉\ge 0,\forall q\in {\mathbb{R}}^N.$$

Then, we take $q=0$ and $q=2u$, and we have

$$〈\lambda ,u〉\le 0,〈\lambda ,u〉\ge 0,$$

which proves (9).    □

For given $u^0,{\lambda }^0\in {\mathbb{R}}_{+}^N$, the ADMM algorithm can now be described as the following three-step pattern:

$$\left\{\begin{array}{}\hfill \mathrm{(10)}& u^{k+1}=arg{\displaystyle \underset{u\in {\mathbb{R}}^N}{min}}{L}_{\rho }\left(u,p^k,{\lambda }^k\right),\hfill \hfill \mathrm{(11)}& p^{k+1}=arg{\displaystyle \underset{p\in {\mathbb{R}}_{+}^N}{min}}{L}_{\rho }\left(u^{k+1},p,{\lambda }^k\right),\hfill \hfill \mathrm{(12)}& {\lambda }^{k+1}={\lambda }^k+\rho \left(p^{k+1}-u^{k+1}\right),\hfill \end{array}\right.$$

until some stopping criterion is satisfied. Minimizing (10) results in an equilibrium problem that admits a unique solution. The minimization of (11) is solved explicitly. Let ${\left(x\right)}^{+}=max(0,x)(x\in {\mathbb{R}}^N)$. Then, we describe Algorithm 1 as follows [8].

Algorithm 1 ADMM

Step 0. ${\lambda }^0\in {\mathbb{R}}^N$, $p^0\in {\mathbb{R}}^N$ and $\rho >0$ are given; set $k:=0$. Step 1. Find $u^{k+1}\in {\mathbb{R}}^N$ by solving $$(A+\rho \mathbb{I})u^{k+1}=\rho p^k+{\lambda }^k-b.$$ (13) Step 2. Find $p^{k+1}\in {\mathbb{R}}^N$ by solving $$p^{k+1}=u^{k+1}-\frac{1}{\rho }[{\lambda }^k-{({\lambda }^k-\rho (u^{k+1}-\psi ))}^{+}].$$ (14) Step 3. Find ${\lambda }^{k+1}\in {\mathbb{R}}^N$ by solving $${\lambda }^{k+1}={\lambda }^k-\rho (u^{k+1}-p^{k+1}).$$ (15) If the stopping criterion is satisfied then stop; else, go to Step 1 with $k:=k+1$.

In letting ${\lambda }_{-}^k=min(0,{\lambda }_{+}^{k-1}-\rho (u^k-\psi ))$ and ${\left(x\right)}^{-}=min(0,x)(x\in {\mathbb{R}}^N)$, Algorithm 1 without $p^k$ and $p^{k+1}$ can be described as Algorithm 2 [9].

Algorithm 2 ADMM without auxiliary unknown

Step 0. ${\lambda }^0\in {\mathbb{R}}^N$, $p^0\in {\mathbb{R}}^N$ and $\rho >0$ are given; set $k:=0$. Step 1. Find $u^{k+1}\in {\mathbb{R}}^N$ by solving $$(A+\rho \mathbb{I})u^{k+1}=({\lambda }_{-}^k-{\lambda }_{+}^k)-b.$$ (16) Step 2. Find ${\lambda }_{+}^{k+1},{\lambda }_{-}^{k+1}\in {\mathbb{R}}^N$ by solving $${\lambda }_{+}^{k+1}={({\lambda }_{+}^k-\rho (u^{k+1}-\psi ))}^{+},$$ (17) $${\lambda }_{-}^{k+1}={({\lambda }_{+}^k-\rho (u^{k+1}-\psi ))}^{-},$$ (18) and go to Step 1 with $k:=k+1$. To obtain the pure dual algorithm for the unilateral obstacle problem, we from (16) to obtain $$u^{k+1}=A_{\rho }^{-1}({\lambda }_{-}^k-{\lambda }_{+}^k-b).$$ (19) where $$A_{\rho }=A+\rho \mathbb{I}.$$ (20)

We eliminate the vector $u^{k+1}$ in (17) and (18) using (19). We then obtain Algorithm 3.

Algorithm 3 Dual ADMM

Step 0. ${\lambda }^0\in {\mathbb{R}}^N$ and $\rho >0$ are given; set $a=-A_{\rho }^{-1}b$ and $k:=0$. Step 1. Compute ${\lambda }_{+}^{k+1}$ and ${\lambda }_{-}^{k+1}$ as $${\lambda }_{+}^{k+1}={((\mathbb{I}-\rho A_{\rho }^{-1}){\lambda }_{+}^k+\rho A_{\rho }^{-1}{\lambda }_{-}^k+\rho a)}^{+},$$ (21) $${\lambda }_{-}^{k+1}={((\mathbb{I}-\rho A_{\rho }^{-1}){\lambda }_{+}^k+\rho A_{\rho }^{-1}{\lambda }_{-}^k+\rho a)}^{-},$$ (22) and repeat Step 1 with $k:=k+1$.

3. Proof of Convergence

In this section, we prove the convergence of Algorithm 3. From $min(0,x)=-max(0,-x)$, it follows that iteration (22) can be rewritten as

$${\lambda }_{-}^{k+1}=-{(-(\mathbb{I}-\rho A_{\rho }^{-1}){\lambda }_{+}^k-\rho A_{\rho }^{-1}{\lambda }_{-}^k-\rho a)}^{+}.$$

(23)

Since $\parallel x^{+}\parallel \le \parallel x\parallel$ for any $x\in {\mathbb{R}}^N$, from (21) and (23), we obtain

$$\parallel {\lambda }_{+}^{k+1}\parallel =\parallel ((\mathbb{I}-\rho A_{\rho }^{-1}){\lambda }_{+}^k+\rho A_{\rho }^{-1}{\lambda }_{-}^k\parallel +\parallel \rho a\parallel ,$$

(24)

$$\parallel {\lambda }_{-}^{k+1}\parallel =\parallel -(-(\mathbb{I}-\rho A_{\rho }^{-1}){\lambda }_{+}^k-\rho A_{\rho }^{-1}{\lambda }_{-}^k\parallel +\parallel \rho a\parallel .$$

(25)

By adding (24) and (25), we have

$$\begin{array}{cc}\hfill \parallel {\lambda }_{+}^{k+1}\parallel +\parallel {\lambda }_{-}^{k+1}\parallel \le & \parallel ((\mathbb{I}-\rho A_{\rho }^{-1}){\lambda }_{+}^k+\rho A_{\rho }^{-1}{\lambda }_{-}^k\parallel \hfill \\ \hfill & +\parallel -(-(\mathbb{I}-\rho A_{\rho }^{-1}){\lambda }_{+}^k-\rho A_{\rho }^{-1}{\lambda }_{-}^k\parallel +2\parallel \rho a\parallel .\hfill \end{array}$$

(26)

Let us set

$${\lambda }_{\pm }^k=\left[\begin{array}{c}{\lambda }_{+}^k\\ {\lambda }_{-}^k\end{array}\right]\mathrm{and}M=\left[\begin{array}{cc}\mathbb{I}-\rho A_{\rho }^{-1}& \rho A_{\rho }^{-1}\\ -(\mathbb{I}-\rho A_{\rho }^{-1})& -\rho A_{\rho }^{-1}\end{array}\right]$$

and define the normal $\parallel {\lambda }_{\pm }^k\parallel =\parallel {\lambda }_{+}^k\parallel +\parallel {\lambda }_{-}^k\parallel$. Then, from (26), we obtain

$$\parallel {\lambda }_{\pm }^{k+1}\parallel =\parallel M{\lambda }_{\pm }^k\parallel +2\parallel \rho a\parallel .$$

Therefore, the spectral radius of M determines the convergence of Algorithm 3. To prove the convergence, the following lemma will be needed [8].

 Lemma 3. 

M has the same non-zero eigenvalues as $\mathbb{I}-2\rho A_{\rho }^{-1}$.

 Proof. 

Let $M_1=\mathbb{I}-\rho A_{\rho }^{-1}$ and $M_2=\rho A_{\rho }^{-1}$. We have

$$M=\left[\begin{array}{cc}{M}_1& M_2\\ -M_1& -M_2.\end{array}\right]$$

We need to prove that M has the same eigenvalues as $M_1-M_2=\mathbb{I}-2\rho A_{\rho }^{-1}$. Suppose that $\alpha$ is an eigenvalue of M and $\left[\begin{array}{c}{\xi }_1\\ {\xi }_2\end{array}\right]$ (${\xi }_1,{\xi }_2\in {\mathbb{R}}^N$) is its eigenvector. It follows that

$$M_1{\xi }_1+M_2{\xi }_2=\alpha {\xi }_1,$$

(27)

$$-M_1{\xi }_1-M_2{\xi }_2=\alpha {\xi }_2,$$

(28)

Adding (27) and (28), we obtain $\alpha ({\xi }_1+{\xi }_2)=0$. For non-zero $\alpha$, we then have ${\xi }_1=-{\xi }_2$, and it follows from (27) and (28) that

$$(M_1-M_2){\xi }_1=\alpha {\xi }_1,$$

(29)

$$(M_1-M_2){\xi }_2=\alpha {\xi }_2,$$

(30)

Then, we prove that $\alpha$ is an eigenvalue of $M_1-M_2$. From the relations above, it is also the eigenvalues of M.    □

From Lemma 3 above, we have the convergence analysis as follows.

 Theorem 1. 

ADMM Algorithm 3 converges when $\rho >0$. Let ${\lambda }_M$ and ${\lambda }_m$ be the largest and smallest eigenvalues of A, respectively. Then, the optimal penalty parameter is given by

$${\rho }^{*}=\sqrt{{\lambda }_M{\lambda }_m}.$$

(31)

 Proof. 

To study the eigenvalues of $\mathbb{I}-2\rho A_{\rho }^{-1}$, we consider the following sequence:

$${\mu }^{k+1}=(\mathbb{I}-2\rho A_{\rho }^{-1}){\mu }^k,$$

(32)

Since

$$A_{\rho }^{-1}={(A+\rho \mathbb{I})}^{-1}={(\mathbb{I}+\rho A^{-1})}^{-1}{A}^{-1},$$

(33)

from (32) and (33), we have

$${\mu }^{k+1}=(\mathbb{I}-2\rho A^{-1}{(\mathbb{I}+\rho A^{-1})}^{-1}){\mu }^k,$$

(34)

From Lemma 1 and the convergence of the sequence $\left\{{\mu }^k\right\}$, we can obtain the convergence of the sequences $\{{\lambda }_{+}^k,{\lambda }_{-}^k\}$ and $\left\{u^k\right\}$. To obtain the eigenvalues of the iterative matrix $G=\mathbb{I}-2\rho A^{-1}{(\mathbb{I}+\rho A^{-1})}^{-1}$ from the eigenvalues of $A^{-1}$, let ${\lambda }_i\left(A^{-1}\right)(i=1,\cdots ,N)$ be the eigenvalues of $A^{-1}$ sorted in ascending order. Then, the eigenvalues of G satisfy

$${\lambda }_i\left(G\right)=1-\frac{2\rho {\lambda }_i\left(A^{-1}\right)}{1+\rho {\lambda }_i\left(A^{-1}\right)}=\frac{1-\rho {\lambda }_i\left(A^{-1}\right)}{1+\rho {\lambda }_i\left(A^{-1}\right)}.$$

(35)

Because all eigenvalues ${\lambda }_i\left(A^{-1}\right)$ and penalty parameter $\rho$ are positive, from (35), we obtain

$$-1<{\lambda }_i\left(G\right)<1,i=1,2,\cdots ,N.$$

(36)

Hence, we proved the convergence of Algorithm 3. Moreover, from (35), it follows that the optimal choice for $\rho$ is

$$\frac{1-\rho {\lambda }_m\left(A^{-1}\right)}{1+\rho {\lambda }_m\left(A^{-1}\right)}=-\frac{1-\rho {\lambda }_M\left(A^{-1}\right)}{1+\rho {\lambda }_M\left(A^{-1}\right)},$$

(37)

so that the optimal parameter is given by

$${\rho }^{*}=\frac{1}{\sqrt{{\lambda }_M\left(A^{-1}\right){\lambda }_m\left(A^{-1}\right)}}=\sqrt{{\lambda }_M\left(A\right){\lambda }_m\left(A\right)}.$$

(38)

Then, we complete the proof.    □

4. Optimal Penalty Parameter Approximation

In this section, we introduce Kronecker Products to compute the largest and smallest eigenvalues of $A^{-1}$. For the sake of simplicity, we use the discretization with a regular grid in the square. As the tridiagonal matrix B = tridiag$(-1,2,-1)$ plays a central role in our method, here, we consider its eigenvalues and eigenvectors [33].

 Lemma 4. 

The eigenvalues of B are given as ${\lambda }_i=2(1-cos\frac{i\pi }{n+1})$, $i=1,2,\cdots ,n$. The corresponding eigenvectors can be given by $z_i\left(k\right)=\sqrt{\frac{2}{n+1}}(1-cos\frac{i\pi }{n+1})sin\left(\frac{ik\pi }{n+1}\right)$, where $i,k=1,2,\cdots ,n$. Let the diagonal matrix $\mathrm{\Lambda }=diag({\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_n)$, and Z be the orthogonal matrix with columns as eigenvectors. Then, $B=Z\mathrm{\Lambda }{Z}^T$.

It follows that the greatest eigenvalue of B is

$${\lambda }_M={\lambda }_n\approx 4,$$

and the smallest eigenvalue of B is

$${\lambda }_m={\lambda }_1\approx \frac{{\pi }^2}{{(n+1)}^2}.$$

As the matrix A can be rewritten in terms of Kronecker products,

$$\begin{array}{cc}\hfill A& =(B\otimes \mathbb{I}+\mathbb{I}\otimes B)\hfill \\ \hfill & =\left[\begin{array}{ccccc}2\mathbb{I}& -\mathbb{I}& & & \\ -\mathbb{I}& 2\mathbb{I}& -\mathbb{I}& & \\ & \ddots & \ddots & \ddots & \\ & & -\mathbb{I}& 2\mathbb{I}& -\mathbb{I}\\ & & & -\mathbb{I}& 2\mathbb{I}\end{array}\right]+\left[\begin{array}{ccccc}B& & & & \\ & B& & & \\ & & \ddots & & \\ & & & B& \\ & & & & B\end{array}\right],\hfill \end{array}$$

we refer to the following lemma to compute the eigenvalue of the matrix A [33].

 Lemma 5. 

Let $B=Z\mathrm{\Lambda }{Z}^T$ be the eigendecomposition of B, where $\mathrm{\Lambda }=diag({\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_n)$, and $Z=[z_1,z_2,\cdots ,z_n]$ is the orthogonal matrix with columns as eigenvectors. Then, the eigendecomposition of $A=(B\otimes \mathbb{I}+\mathbb{I}\otimes B)$ is

$$B\otimes \mathbb{I}+\mathbb{I}\otimes B=(Z\otimes Z)·(\mathrm{\Lambda }\otimes \mathbb{I}+\mathbb{I}\otimes \mathrm{\Lambda }){(Z\otimes Z)}^T.$$

$\mathrm{\Lambda }\otimes \mathbb{I}+\mathbb{I}\otimes \mathrm{\Lambda }$ is a diagonal matrix, and its $(in+j)th$ diagonal component is obtained using the $(i,j)th$ eigenvalue of B${\lambda }_{i,j}={\lambda }_i+{\lambda }_j$. $Z\otimes Z$ is an orthogonal matrix, and its $(in+j)th$ columns are obtained using the corresponding eigenvalue $z_i\otimes z_j$.

It follows from Lemma 5 that the largest eigenvalue of A is

$${\lambda }_M={\lambda }_n+{\lambda }_n=4(1-cos\frac{n\pi }{n+1})\approx 8,$$

and the smallest eigenvalue of A is

$${\lambda }_m={\lambda }_1+{\lambda }_1=4(1-cos\frac{\pi }{n+1})=4(1-(1-2{sin}^2\frac{\pi }{2(n+1)}))\approx \frac{2{\pi }^2}{{(n+1)}^2}.$$

Therefore, from the eigenvalues of A, we obtain the optimal penalty parameter approximation as follows:

$${\rho }^{*}=\sqrt{{\lambda }_m\left(A\right){\lambda }_M\left(A\right)}=4\sqrt{(1-cos\frac{n\pi }{n+1})(1-cos\frac{\pi }{n+1})}\approx \frac{4\pi }{n+1}.$$

(39)

5. Results and Discussion

In this section, we present some preliminary numerical results to assess the performance of the ADMM (Algorithm 1). Note that the semi-smooth Newton method (SSNM) is superlinear convergent and equivalent to the primal-dual active set algorithm. This approach is widely used for unilateral obstacle problems [10,24,25]. We set the set of indices $N_h=\{1,2,\cdots ,N\}$ and describe Algorithm 4 as follows.

Algorithm 4 SSNM

Step 0. Given ${\lambda }^0\in {\mathbb{R}}^N$, $u^0\in {\mathbb{R}}^N$, and $\rho >0$, set $k:=0$. Step 1. Set $A_k=\{i,{\lambda }_i^k-\rho {\displaystyle \sum _{j=1}^N}{A}_{ij}{u}_j^k>0\}$. Step 2. Find $u^{k+1}\in {\mathbb{R}}^N$ and ${\lambda }^{k+1}\in {\mathbb{R}}^N$ such that $$A{u}^{k+1}+b={\lambda }^{k+1},$$ (40) where ${\lambda }_i^{k+1}={\lambda }_i^k-\rho {\displaystyle \sum _{j=1}^N}{A}_{ij}{u}_j^{k+1}$ for $i\in A_k$, and ${\lambda }_i^{k+1}=0$ for $i\in N_h\setminus A_k$. Step 3. If $A_{k+1}=A_k$, then stop; else, return to Step 1 with $k:=k+1$.

It is well known that the active set $A_k$ changes during the iterative process of the SNNM, and the matrix of linear system (40) varies at each iteration. Different from the SSNM, the iterative matrix of our proposed ADMM is constant in the iteration, saving time on matrix computation. This advantage is particularly noticeable for large systems where N is enough large.

For the SSNM algorithm, the convergence speed also depends on the penalty parameter $\rho$, and the change is very moderate when $\rho$ is large enough. For the examples that we tested, the number of iterations did not change anymore for $\rho \ge {10}^4$. In all tests, we took parameters $\rho ={10}^4$ for the SSNM and ${\rho }^{*}$, using (39) for the ADMM. All computations were performed with Matlab R2016b on a Notebook personal computer Lenovo with a 13th Gen Intel(R) Core(TM) i7-13700H CPU 2.40 GHz, and the stopping criterion was $\parallel u^{k+1}-u^k{\parallel }_2\le {10}^{-5}{\parallel u^{k+1}\parallel }_2$. Here, we present several numerical examples to show the efficiency of the ADMM and compare the results with those of the SSNM.

5.1. Example 1

We first considered unilateral obstacle problem (1) on a unit square $\mathrm{\Omega }=(0,1)\times (0,1)$ with $f\left(x\right)=-50$ and $\psi \left(x\right)=-0.5$. We choose the parameter $\rho =\frac{4\pi }{81}$ and applied the proposed algorithm to this problem with step size $h=\frac{1}{80}$. We show the numerical solution $u_h$ and the contact zone $u=\psi$ in Figure 1 and Figure 2, respectively. Our numerical results appear to be consistent with those in [22]. Let Iterations and CPU(s) represent the number of iterations and amount of CPU time in seconds, respectively. We compare the number of iterations and amount of CPU time with that of the SSNM in

Table 1. Performances of ADMM with ${\rho }^{*}$ and SSNM with $\rho ={10}^4$ (example 1).

h	ADMM (Algorithm 1)		SSNM (Algorithm 4)
	Iterations	CPU (s)	Iterations	CPU (s)
$\frac{1}{10}$	20	0.0242	7	0.0245
$\frac{1}{20}$	16	0.0260	8	0.0331
$\frac{1}{40}$	22	0.0497	9	0.0723
$\frac{1}{80}$	30	0.3891	9	0.7741
$\frac{1}{160}$	55	3.6584	9	8.9213

. One can see that the ADMM requires less CPU time than the SSN [8] because of the constant matrix for the ADMM and the matrix factorizations for the SSNM in every iteration.

To further investigate the convergence of the proposed ADMM, we choose different values of $\rho$, and Figure 3 shows the number of iterations for different step sizes h. The curves in Figure 3 show that the penalty parameter $\rho$ can heavily affect the convergence rate of the ADMM. We can easily use (39) to compute all optimal parameters ${\rho }^{*}\approx \frac{4\pi }{11},\frac{4\pi }{21},\frac{4\pi }{41},\frac{4\pi }{81},\frac{4\pi }{161}$ for $h=\frac{1}{10},\frac{1}{20},\frac{1}{40},\frac{1}{80},\frac{1}{160}$ ($N={(\frac{1}{h}-1)}^2$), respectively. These numerical results also demonstrate the accuracy and effectiveness of the proposed ADMM.

5.2. Example 2

Let $\mathrm{\Omega }=(-1.5,1.5)\times (-1.5,1.5)$. We computed obstacle problem (1) with $f\left(x\right)=-2$ and $\psi \left(x\right)=0$. For this example, the analytical solution is

$$u\left(x\right)=\left\{\begin{array}{cc}\frac{r^2-1}{2}-lnr\hfill & \hfill \mathrm{if}r\ge 1,\\ 0\hfill & \hfill \mathrm{if}r<1,\end{array}\right.$$

where $x=(x_1,x_2)\in {\mathbb{R}}^2$, and $r=\sqrt{x_1^2+x_2^2}$ [24,25,28,29,31].

In this example, we considered the same parameters as in Example 1. Figure 4 and Figure 5, respectively, plot the numerical solution $u_h$ and the contact zone with $h=\frac{3}{80}$. We observe that our numerical solution is in good agreement with the analytical solution.

Table 2. Performances of ADMM with ${\rho }^{*}$ and SSNM with $\rho ={10}^4$ (example 2).

h	ADMM (Algorithm 1)		SSNM (Algorithm 4)
	Iterations	CPU (s)	Iterations	CPU (s)
$\frac{3}{10}$	23	0.0248	6	0.0311
$\frac{3}{20}$	18	0.0270	8	0.0333
$\frac{3}{40}$	23	0.0485	11	0.0869
$\frac{3}{80}$	38	0.3730	13	1.0514
$\frac{3}{160}$	63	3.4886	14	12.5551

presents the number of iterations and amount of CPU time of both algorithms, with different step sizes h. Figure 6 shows the convergence speed versus the parameter $\rho$ with different step sizes h. We also observe that the proposed ADMM also requires less CPU time than the SSNM, and the optimal parameters of the numerical results are in good agreement with the ${\rho }^{*}$ obtained using (39).

5.3. Example 3

Let $\mathrm{\Omega }=(-1,1)\times (-1,1)$. We computed obstacle problem (1) with $\psi \left(x\right)=0$ and

$$f\left(x\right)=\left\{\begin{array}{cc}-8(2r^2-r_0^2)\hfill & \hfill \mathrm{if}r>r_0,\\ -8r_0^2(1-r^2+r_0^2)\hfill & \hfill \mathrm{otherwise}.\end{array}\right.$$

For any $0<r_0<1$, the analytical solution of this example is

$$u\left(x\right)={(max\{r^2-r_0^2,0\})}^2.$$

(41)

As in [23,27], we took $r_0=0.7$ for the computation. Figure 7 and Figure 8, respectively, plot the numerical solution $u_h$ and the contact zone with $h=\frac{1}{40}$.

Table 3. Performances of ADMM with ${\rho }^{*}$ and SSNM with $\rho ={10}^4$ (example 4).

h	ADMM (Algorithm 1)		SSNM (Algorithm 4)
	Iterations	CPU (s)	Iterations	CPU (s)
$\frac{1}{5}$	27	0.0365	6	0.0386
$\frac{1}{10}$	32	0.0426	8	0.0431
$\frac{1}{20}$	53	0.0857	11	0.0976
$\frac{1}{40}$	79	0.6882	17	1.3378
$\frac{1}{80}$	116	4.710	29	24.5476

displays the number of iterations and the CPU time. Figure 9 shows the convergence speed versus the parameter $\rho$ and step size h. It can be seen that the proposed ADMM is very effective.

5.4. Example 4

Finally, we considered the unilateral obstacle problem (1) in $\mathrm{\Omega }=(-1.5,1.5)\times (-1.5,1.5)$ with $f=0$ and

$$\psi \left(x\right)=\left\{\begin{array}{cc}\sqrt{1-r^2}\hfill & \hfill \mathrm{if}r\le 1,\\ -1\hfill & \hfill \mathrm{if}r>1.\end{array}\right.$$

The exact solution for this example is

$$u\left(x\right)=\left\{\begin{array}{cc}\sqrt{1-r^2}\hfill & \hfill \mathrm{if}r\le r_0,\\ -\frac{r_0^{2}ln(r/2)}{\sqrt{1-r_0^2}}\hfill & \hfill \mathrm{otherwise},\end{array}\right.$$

where $r_0=0.6979651482$ [24,25,28,29,31].

Figure 10 and Figure 11 plot the numerical solution $u_h$ and the contact zone, respectively.

Table 4. Performances of ADMM with ${\rho }^{*}$ and SSNM with $\rho ={10}^4$ (example 4).

h	ADMM (Algorithm 1)		SSNM (Algorithm 4)
	Iterations	CPU (s)	Iterations	CPU (s)
$\frac{3}{10}$	21	0.0285	5	0.0292
$\frac{3}{20}$	26	0.0326	5	0.0336
$\frac{3}{40}$	43	0.0691	8	0.0768
$\frac{3}{80}$	75	0.6283	12	1.0372
$\frac{3}{160}$	133	5.4149	21	18.8991

displays the number of iterations and CPU time. Figure 12 shows the impact of parameters $\rho$ in the convergence of the ADMM. It can be seen that the numerical results are also in good agreement with our theoretical analysis.

6. Conclusions

In this paper, we propose an ADMM with an optimal parameter for unilateral obstacle problems. To solve the LCP arising from the discretization of free boundary problems with the FDM, we use an optimal penalty parameter selection to improve the performance of the ADMM. The numerical results demonstrate that the proposed ADMM is feasible and efficient for the numerical solution. In the future, we will apply this method to bilateral obstacle problems.