# Alternating Direction Method of Multipliers (ADMM) for Frictional Contact

Vincent Acary, Yoshihiro Kanno, Nicolas Molina

► To cite this version:

Vincent Acary, Yoshihiro Kanno, Nicolas Molina. Alternating Direction Method of Multipliers (ADMM) for Frictional Contact. 2019. hal-05222639

HAL Id: hal-05222639

https://hal.science/hal-05222639v1

Preprint submitted on 25 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Alternating Direction Method of Multipliers (ADMM)
## for Frictional Contact

Vincent Acary *, Yoshihiro Kanno † , Nicolas Molina ‡

Abstract

**Abstract:** In this report, we review the quadratic programming (QP) over second-order cones formulation of the discrete frictional contact problem that arises in space and time discretized mechanical systems with unilateral contact and three-dimensional Coulomb’s friction. Thanks to this formulation, various numerical methods emerge naturally for solving the problem [1]. Here we propose the application of the Alternating Direction Method of Multipliers (ADMM) with various improvements proposed in the literature. This numerical technique is compared over a large set of test examples using performance profiles.

**Keywords:** Multibody systems, nonsmooth mechanics, unilateral constraints, Coulomb's friction, impact, numerical methods, ADMM.

<footnote>*INRIA, 655, Avenue de l'Europe, 38334 Saint Ismier Cedex, France. E-mail: vincent.acary@inria.fr.
†Laboratory for Future Interdisciplinary Research of Science and Technology, Institute of Innovative Research, Tokyo Institute of Technology, Nagatsuta 4259, Yokohama 226-8503, Japan. E-mail: kanno.y.af@m.titech.ac.jp.
‡Department of Mechanical and Metallurgical Engineering, Escuela de Ingenieria, Pontificia Universidad Catolica de Chile, Vicuna Mackenna 4860, Santiago, Chile. E-mail: namolina@uc.cl.</footnote>
# Contents

<table><tr><td>1</td><td>Introduction</td><td>3</td></tr><tr><td>1.1</td><td>Problem statement. 3D Discrete frictional contact problems</td><td>3</td></tr><tr><td>1.2</td><td>Outline of the article</td><td>5</td></tr><tr><td>2</td><td>Generalities on the ADMM method</td><td>5</td></tr><tr><td>2.1</td><td>Fundamentals</td><td>5</td></tr><tr><td>2.2</td><td>Algorithm</td><td>9</td></tr><tr><td>3</td><td>ADMM formulation of problem (P<sub>c</sub>) and (R<sub>c</sub>)</td><td>10</td></tr><tr><td>4</td><td>Numerical framework for the comparisons</td><td>11</td></tr><tr><td>4.1</td><td>The FCLib collection of examples.</td><td>11</td></tr><tr><td>4.2</td><td>Measuring errors</td><td>13</td></tr><tr><td>4.3</td><td>Comparison with performance profiles</td><td>13</td></tr><tr><td>4.4</td><td>Implementation details.</td><td>14</td></tr><tr><td>4.5</td><td>Nomenclature for the algorithms</td><td>14</td></tr><tr><td>5</td><td>Numerical comparisons for the convex problems</td><td>15</td></tr><tr><td>5.1</td><td>Convergence rate</td><td>15</td></tr><tr><td>5.2</td><td>Comparison of ADMM variants for the convex problems</td><td>16</td></tr><tr><td>6</td><td>Solving the original nonconvex problems (P) and (R).</td><td>20</td></tr><tr><td>6.1</td><td>Initial guess</td><td>20</td></tr><tr><td>6.2</td><td>Internal update</td><td>20</td></tr><tr><td>6.3</td><td>External update</td><td>21</td></tr><tr><td>6.4</td><td>SICONOS error</td><td>21</td></tr><tr><td>6.5</td><td>Numerical experience for original nonsmooth and nonconvex problems.</td><td>22</td></tr><tr><td>7</td><td>Numerical experiments.</td><td>25</td></tr><tr><td>7.1</td><td>Standard first order techniques.</td><td>25</td></tr><tr><td>7.2</td><td>Convergence rate</td><td>25</td></tr><tr><td>7.3</td><td>Comparisons of variants of ADMM solvers for the global problem (P)</td><td>25</td></tr><tr><td>7.4</td><td>Comparisons of variants of ADMM solvers</td><td>32</td></tr><tr><td>7.5</td><td>Comparisons of variants of ADMM solvers with NSGS</td><td>37</td></tr><tr><td>A</td><td>Projection onto dual second-order cone</td><td>42</td></tr><tr><td>B</td><td>Convex formulation without initial affine constraints</td><td>42</td></tr><tr><td>C</td><td>Optimality conditions of problem (6)</td><td>43</td></tr></table>
# 1 Introduction

Mechanical problems with unilateral contact and Coulomb friction are ubiquitous in a wide range of industrial applications ranging from multi-body systems such as linkages, robotic applications, or mechanics of continuum media with interfaces to granular and divided materials. Nevertheless, after the pioneering work of [2], [3], [4–6], [7, 8], [9], [10], [11], [12], there are still active research on this subject in the computational mechanics and applied mathematics communities. The difficulty lies mainly in the fact that the discrete frictional contact problems can be stated as a variational inequality that is non-monotone, and cannot be reformulated as a gentle convex optimization problem. It is well known that the frictionless case yields convex optimization problem, but non-associated character of Coulomb friction destroys this property. Furthermore, industrial applications lead to the definition of models where the unilateral constraints are redundant and then imply hyper-staticity. This results in very slow convergence rate or even loss of convergence and robustness issues. Recently, a thorough comparison of solvers has been made in [13] and the main conclusions are as follows:

* the standard first-methods (projected Jacobi, projected Gauss Seidel, Extra-gradient methods to cite a few), that do not use the information on the Hessian of the constraints, are generally robust even with large coefficients of friction and a rank-deficiency of the constraints Jacobian matrix, but suffer from very slow convergence rates. They are able to identify quickly a rough solution (around $10^{-2}$ in relative error) and the contact status for a large percentile of constraints (around 95%), but it is difficult to improve the accuracy of the solution in a efficient way.

* the second-order methods (semi-smooth Newton of interior points techniques) are, on the contrary, not robust due the redundancy of constraints, but improves a lot the convergence rate when they converge.

In this article, we propose an alternative method for solving the contact friction problem based the Accelerated Direction methods of multipliers (ADMM). The ultimate goal is to design a numerical method for large-scale problems with redundancy, that is robust and that improves the convergence rate of the standard first order methods. The ADMM is an invaluable element of the modern optimization toolbox. The ADMM decomposes complex optimization problems into sequences of simpler subproblems, often solvable in closed form; its simplicity, flexibility, and broad applicability, make ADMM a state-of-the-art solver in machine learning, signal processing, and many other areas [14].

## 1.1 Problem statement. 3D Discrete frictional contact problems

In this section, following the work in [15, 16] we formulate an abstract, algebraic finite-dimensional frictional contact problem. This problem is modeled as a complementarity problem over second-cones, and the properties of the latter are discussed. Let $n_c \in \mathbb{N}$ be the number of contact points and $n \in \mathbb{N}$ the number of degrees of freedom of a discrete mechanical system, after space and time-discretization.

**Formulation in generalized coordinates** The problem data are: a positive definite symmetric matrix $M \in \mathbb{R}^{n \times n}$, a vector $f \in \mathbb{R}^n$, a matrix $H \in \mathbb{R}^{n \times m}$ with $m = 3n_c$, a vector $w \in \mathbb{R}^m$ and a vector of coefficients of friction $\mu \in \mathbb{R}^{n_c}$. The unknowns are two vectors $v \in \mathbb{R}^n$, a velocity-like
vector and $r \in \mathbb{R}^m$, a contact reaction or impulse, solution to

$$ \left\{ \begin{array}{ll} M\mathbf{v} + \mathbf{f} = H^\top \mathbf{r} & \mathbf{u} := H\mathbf{v} + \mathbf{w} \\ K_{e,\mu}^* \ni \tilde{\mathbf{u}} \perp \mathbf{r} \in K_{e,\mu} & \tilde{\mathbf{u}} := \mathbf{u} + \Phi(\mathbf{u}) \end{array} \right. \quad (P), \quad (1) $$

where the set $K_{e,\mu} \subseteq \mathbb{R}^m$ is the Cartesian product of Coulomb's friction second-order cone at each contact, and $K_{e,\mu}^*$ its dual cone, i.e.,

$$ K_{e,\mu} = \prod_{\alpha=1 \dots n_c} K_{e,\mu}^{\alpha} = \prod_{\alpha=1 \dots n_c} \{(x_1, x_2) \in \mathbb{R} \times \mathbb{R}^2 \mid \|x_2\| \le \mu x_1\}, \quad (2) $$

$$ K_{e,\mu}^* = \prod_{\alpha=1 \dots n_c} K_{e,\mu}^{\alpha*} = \prod_{\alpha=1 \dots n_c} \{(x_1, x_2) \in \mathbb{R} \times \mathbb{R}^2 \mid \|x_2\| \le \frac{1}{\mu} x_1\}. \quad (3) $$

The variable $\mathbf{u}$ is the local relative velocity at contact, and $\tilde{\mathbf{u}}$ is a modified velocity that enables to state the complementarity condition given in (1). The function $\Phi: \mathbb{R}^m \to \mathbb{R}^m$ is a nonsmooth function defined as

$$ \Phi(\mathbf{u}) = [\mu^{\alpha} \|E_t \mathbf{u}^{\alpha}\|_{\hat{e}_N, \alpha} = 1 \dots n_c]^{\top} \quad (4) $$

where $\hat{e}_N = [1, 0, 0]^{\top}$ is the normal component vector and $E_T = \text{diag}(0, 1, 1) \in \mathbb{R}^{3 \times 3}$ is a linear transformation that maps to the tangential part of $\tilde{\mathbf{u}}$.

This problem (P) in (1) is a quite generic problem that can obtained after a space-discretization of quasi-static problems of solids mechanics with unilateral contact (Signorini condition) and dry friction (Coulomb friction), or time and space discretization of collection of rigid or flexible bodies. The first equation in (1) models the linearized and discretized equilibrium or equations of motion. The complementarity condition over second order cones with the modified velocity $\tilde{\mathbf{u}}$ is a condensed formulation of the whole contact and nonassociated friction behavior at the interface. A more detailed description how we formulate the problem can be found in [17] and [18]. t

The problem (P) is a nonlinear and nonsmooth complementarity problem over a Cartesian product of second-order cones in $\mathbb{R}^3$. If the coefficient of friction $\mu$ vanishes, the problem is frictionless and reduces to a gentle convex Quadratic problem (QP) after straightforward manipulations.

**Reduced formulation in local coordinates at contact**

Assuming that $M$ is non singular, we can substitute the variable $\mathbf{v}$ to obtain

$$ \left\{ \begin{array}{ll} \tilde{\mathbf{u}} = W\mathbf{r} + \mathbf{q} + \Phi(\mathbf{u}) & \mathbf{u} := H\mathbf{v} + \mathbf{w} \\ K_{e,\mu}^* \ni \tilde{\mathbf{u}} \perp \mathbf{r} \in K_{e,\mu} & \tilde{\mathbf{u}} := \mathbf{u} + \Phi(\mathbf{u}) \end{array} \right. \quad (R) \quad (5) $$

where $W = H^\top M^{-1}H$ and $\mathbf{q} = \mathbf{w} + H^\top M^{-1}\mathbf{f}$. This problem is smaller since we reduce the number of unknowns but a part of the structure of the problem has been lost.

**Convex relaxation and optimisation**

As we said earlier, the problem (P) and (R) are non-linear and nonsmooth complementarity problems, or if we prefer variational inequalities, over a product of second-order cones. Since the operator involved in the definitions are not the gradient of a nonlinear functions, these problems are not the first order optimality conditions of “natural optimization” problems. Although is it possible to formulate direct nonconvex and nonsmooth optimization problems from (P) or (R), we prefer to introduce a simple convex problem by relaxation, and to propose to solve the original problems by a succession of convex problems.

Convex relaxation is common for the frictional contact problems, there are mainly three approaches:

* Formulations with two coupled convex problems, splitting the normal and the tangential directions
* Formulations that are based on the Tresca friction with a given friction threshold that does not depends on the normal reaction; we end up with an optimisation problem on semi-cylinders

* Fixing the correction of the modified velocity $\Phi(u)$ for a while, we obtain a convex optimization problem over second order cones.

Let us detail the third approach. Since the problem (P) is nonsmooth and nonconvex, the use of an associated optimization problem is interesting from the numerical point of view if we want to improve the robustness and the stability of the numerical methods. In [18, 19], a parametric convex optimization formulation is presented for the dynamical Coulomb friction problem. The problem (P) reads as

$$ \underset{\mathbf{v},\tilde{\mathbf{u}}}{\text{Minimize}} \quad \frac{1}{2} \mathbf{v}^{\top} M \mathbf{v} + \mathbf{f}^{\top} \mathbf{v} \quad (\text{P}_c) \qquad (6a) $$

$$ \text{subject to} \quad \tilde{\mathbf{u}} = H\mathbf{v} + \mathbf{w} + \phi(\mathbf{s}) \qquad (6b) $$

$$ \tilde{\mathbf{u}} \in K_{e,\mu}^* \qquad (6c) $$

where $\mathbf{s} \in \mathbb{R}^m$ is a parameter vector that plays the role of the norm of the sliding velocity and $\phi(\mathbf{s}) = [\mu^\alpha \mathbf{s}, \alpha \in [1...n_c]^\top$. Similarly, the problem (R) reads as

$$ \underset{\mathbf{r}}{\text{Minimize}} \quad \frac{1}{2} \mathbf{r}^{\top} W \mathbf{r} + \tilde{\mathbf{q}}(\mathbf{s})^{\top} \mathbf{r} \qquad (\text{R}_c) \qquad (7a) $$

$$ \text{subject to} \quad \tilde{\mathbf{r}} \in K_{e,\mu} \qquad (7b) $$

with $\tilde{\mathbf{q}}(\mathbf{s}) := \mathbf{q} + \phi(\mathbf{s})$. Both problems are convex problems and if they are feasible, it exists some solutions [19].

## 1.2 Outline of the article

In this article, we will propose accelerated first-order methods for solving (P) and (R). Firstly, in Section 2 and 3, we will design a new method solving ($P_c$) and ($R_c$) based the standard ADMM solver. Then, in Section we will recall standard first-order methods (Projected Gauss-Seidel and Projected Gradient direction methods) for the original problems (P) and (R). Then, we will try to propose accelerations of these methods In Section 6, we will also study several ways to update the $\mathbf{s}$ variables to solve the original problems.

# 2 Generalities on the ADMM method

## 2.1 Fundamentals

In this section we explain the alternating direction method of multipliers (ADMM) for convex optimization following [14]. Let $f: \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ and $g: \mathbb{R}^m \to \mathbb{R} \cup \{+\infty\}$ be a closed proper convex functions. Consider the following convex optimization problem in variables $\mathbf{x} \in \mathbb{R}^n$ and $\mathbf{y} \in \mathbb{R}^m$:

$$ \underset{\mathbf{x}}{\text{Minimize}} \quad f(\mathbf{x}) + g(\mathbf{y}) \qquad (8a) $$

$$ \text{subject to} \quad A\mathbf{x} + B\mathbf{y} = \mathbf{c} \qquad (8b) $$

where $A \in \mathbb{R}^{l \times n}$ and $B \in \mathbb{R}^{l \times m}$ are constraint constant matrices, and $\mathbf{c} \in \mathbb{R}^l$ is a constant vector. The augmented Lagrangian of problem (8) is defined as

$$ L_{\rho}(\mathbf{x}; \mathbf{y}; \mathbf{w}) = f(\mathbf{x}) + g(\mathbf{y}) + \mathbf{w}^{\top}(A\mathbf{x} + B\mathbf{y} - \mathbf{c}) + \frac{\rho}{2}\|A\mathbf{x} + B\mathbf{y} - \mathbf{c}\|^2, \qquad (9) $$
where $\rho > 0$ is the penalty parameter, and $\mathbf{w} \in \mathbb{R}^l$ is the Lagrange multiplier (or the dual variable). The ADMM consists of the following iterations

$$ \mathbf{x}^{k+1} := \arg \min_{\mathbf{x}} L_{\rho}(\mathbf{x}; \mathbf{y}^k; \mathbf{w}^k) \qquad (10) $$

$$ \mathbf{y}^{k+1} := \arg \min_{\mathbf{y}} L_{\rho}(\mathbf{x}^{k+1}; \mathbf{y}; \mathbf{w}^k) \qquad (11) $$

$$ \mathbf{w}^{k+1} := \mathbf{w}^k + \rho(A\mathbf{x}^{k+1} + B\mathbf{y}^{k+1} - \mathbf{c}) \qquad (12) $$

where the dual variable is updated via dual ascent method. To compare with, the standard method of multipliers (or the augmented Lagrangian method) updates the variables as follows:

$$ (\mathbf{x}^{k+1}, \mathbf{y}^{k+1}) := \arg \min_{\mathbf{x}, \mathbf{y}} L_{\rho}(\mathbf{x}; \mathbf{y}; \mathbf{w}^k) \qquad (13a) $$

$$ \mathbf{w}^{k+1} := \mathbf{w}^k + \rho(A\mathbf{x}^{k+1} + B\mathbf{y}^{k+1} - \mathbf{c}) \qquad (13b) $$

Namely, the method of multipliers updates the primal variables $\mathbf{x}$ and $\mathbf{y}$ simultaneously. In contrast, the ADMM updates $\mathbf{x}$ and $\mathbf{y}$ sequentially.

**Scaled form** The ADMM is often described in a slightly different form by defining a scaled multiplier $\mathbf{z}$ by $\mathbf{z} = \mathbf{w}/\rho$. Then the augmented Lagrangian in (9) is reduced to

$$ L_{\rho}(\mathbf{x}; \mathbf{y}; \mathbf{z}) = f(\mathbf{x}) + g(\mathbf{y}) + \frac{\rho}{2} \|A\mathbf{x} + B\mathbf{y} - \mathbf{c} + \mathbf{z}\|^2 - \frac{\rho}{2} \|\mathbf{z}\|^2. \qquad (14) $$

Using $L_{\rho}$ in (14), we can express ADMM as follows:

$$ \mathbf{x}^{k+1} := \arg \min_{\mathbf{x}} L_{\rho}(\mathbf{x}; \mathbf{y}^k; \mathbf{z}^k), \qquad (15) $$

$$ \mathbf{y}^{k+1} := \arg \min_{\mathbf{y}} L_{\rho}(\mathbf{x}^{k+1}; \mathbf{y}; \mathbf{z}^k), \qquad (16) $$

$$ \mathbf{z}^{k+1} := \mathbf{z}^k + A\mathbf{x}^{k+1} + B\mathbf{y}^{k+1} - \mathbf{c} \qquad (17) $$

The expression in (15), (16), and (17) is called the ADMM in the scaled form, and $\mathbf{z}$ is called the scaled dual variable.

**Residuals** The necessary and sufficient optimality conditions for the ADMM problem (8) are primal and dual feasibility, therefore, one common way to measure how well the iterates satisfy the KKT conditions is to define the primal and dual residuals:

$$ \begin{aligned} \mathbf{r}^{k+1} &:= A\mathbf{x}^{k+1} + B\mathbf{y}^{k+1} - \mathbf{c} \\ \mathbf{s}^{k+1} &:= \rho A^T B (\mathbf{y}^{k+1} - \mathbf{y}^k) \end{aligned} \qquad (18) $$

**Stopping criterion** The residuals of the optimality conditions can be related to an approximated bound on the objective suboptimality of the current point

$$ \begin{aligned} \|\mathbf{r}^k\| &\le \epsilon^{\text{pri}} \\ \|\mathbf{s}^k\| &\le \epsilon^{\text{dual}} \end{aligned} \qquad (19) $$

where $\epsilon^{\text{pri}} > 0$ and $\epsilon^{\text{dual}} > 0$ are feasibility tolerances for the primal and dual residuals. These tolerances can be chosen using an absolute and relative criterion, such as

$$ \begin{aligned} \epsilon^{\text{pri}} &= \sqrt{l}\epsilon^{\text{abs}} + \epsilon^{\text{rel}} \max\{\|A\mathbf{x}^{k+1}\|, \|B\mathbf{y}^{k+1}\|, \|\mathbf{c}\|\} \\ \epsilon^{\text{dual}} &= \sqrt{n}\epsilon^{\text{abs}} + \epsilon^{\text{rel}} \|A^T \rho \mathbf{z}^{k+1}\| \end{aligned} \qquad (20) $$

where $\epsilon^{\text{abs}} > 0$ and $\epsilon^{\text{real}} > 0$ is an absolute and relative tolerance, respectively. A reasonable value for the absolute and relative tolerance might be $10^{-6}$ and $10^{-3}$, respectively. This values typically depend on the scale of the data size.
**Acceleration and relaxation** A relaxed version of ADMM was proposed in [20]. This method applies Nesterov's type overrelaxation scheme to the updates of $y^k$ and $z^k$ which takes into account past iterations:

$$ \alpha_{k+1} := \frac{1}{2} \left( 1 + \sqrt{1 + 4\alpha_k^2} \right) \qquad (21a) $$

$$ \hat{y}^{k+1} := y^{k+1} + \frac{\alpha_k - 1}{\alpha_{k+1}} (y^{k+1} - y^k) \qquad (21b) $$

$$ \hat{z}^{k+1} := z^{k+1} + \frac{\alpha_k - 1}{\alpha_{k+1}} (z^{k+1} - \frac{\rho^{k-1}}{\rho^k} z^k). \qquad (21c) $$

It is shown that the method converges in the primal and dual residuals with the rate of $O(1/k^2)$.

For weakly convex problems, we must enforce stability using a restart rule, which often makes the relaxed algorithm more efficient. It is also proposed in [20]. The restart rule relies on a combined residual, which measures both the primal and dual error simultaneously:

$$ e_k := \rho^k \|z^{k+1} - \frac{\rho^{k-1}}{\rho^k} \hat{z}^k\|^2 + \rho^k \|B(y^{k+1} - \hat{y}^k)\|^2. \qquad (22) $$

If $e_k < \eta e_{k-1}$ the relaxed ADMM (21) is performed, otherwise, the following restart rule is applied

$$ \alpha_{k+1} := 1 \qquad (23a) $$

$$ \hat{y}^{k+1} := y^{k+1} \qquad (23b) $$

$$ \hat{z}^{k+1} := z^{k+1} \qquad (23c) $$

$$ e_k \leftarrow \frac{e_{k-1}}{\eta} \qquad (23d) $$

**Choosing the initial penalty parameter** In this section we present two main results that have been proposed for QPs in the initial guess of the parameter and we define another two ones:

* In [21], the variable $\rho^*$ that is shown to maximizes the rate of convergence is

$$ \rho^* := \left( \sqrt{\lambda_{min}(HM^{-1}H^{\top})\lambda_{max}(HM^{-1}H^{\top})} \right)^{-1} \qquad (24) $$

If the constraint matrix $H$ is not full rank, the smallest non-zero eigenvalue is taken, which still works as a heuristic to reduce the rate of convergence. It is worth noting that the matrix $HM^{-1}H^{\top}$ is called *Delassus Matrix* and is previously denoted by $W$.

* In [22], The variable $\rho^*$ that maximizes the rate of convergence is

$$ \rho^* := \sqrt{\lambda_{min}(M)\lambda_{max}(M)} \qquad (25) $$

In this case, the issue associated to the rank-deficiency of $H$ is avoided.

* In this report, we propose another possible value, that avoids the computation of eigenvalues, where the penalty parameter could be seen as the factor that balance the matrix involved in the first term of (38)

$$ \rho^* := \frac{\|M\|_1}{\|H\|_1} \qquad (26) $$

* Finally, a less expensive and standard method is to set

$$ \rho^* := 1 \qquad (27) $$
**Self-adaptive penalty parameter** The ADMM technique is surprisingly robust to poorly selected algorithm parameters: under mild conditions, the method is guaranteed to convergence for all positive values of $\rho$ [14, 23, 24]. However, a key ingredient in the efficiency of ADMM with respect to time is a smart choice of a sequence of parameters $\{\rho_k\}$. A substantial work has been done in the literature mainly motivated by the rate of convergence. This section presents the most popular approach for choosing the sequence $\{\rho_k\}$ for QPs, which is closely related to the formulation in (33). In the following, we present optimal penalty parameter selections for QP.

He et al. [25] argue that adaptively choosing the penalty parameter to balance the residuals is a reasonable heuristic for minimising the distance from convergence: increasing $\rho$ strengthens the penalty term, yielding smaller primal residuals but larger dual ones; conversely, decreasing $\rho$ leads to larger primal and smaller dual residuals. Therefore, this heuristic is implemented so that $\rho$ keeps both residuals of similar magnitude

$$ \rho^{k+1} := \begin{cases} \tau \rho^k & \text{if } \|\mathbf{r}^k\| > \varphi \|\mathbf{s}^k\| \\ \tau^{-1} \rho^k & \text{if } \|\mathbf{s}^k\| > \varphi \|\mathbf{r}^k\| \\ \rho^k & \text{otherwise,} \end{cases} \qquad (28) $$

where $\tau > 1$ and $\varphi > 1$. This method appears to becoming quite popular in the literature (see [26] for references).

The behavior of ADMM under scaling is addressed by Wohlberg [26]. For problems involving physical quantities, for example, scaling corresponds to choices of the units in which the functional value and solution are expressed. The heuristic implemented by He et al. is appropriately scaled to maintain the invariance of the algorithm with respect to scaling as follows:

$$ \rho^{k+1} := \begin{cases} \tau \rho^k & \text{if } \|\mathbf{r}_{rel}^k\| > \xi \varphi \|\mathbf{s}_{rel}^k\| \\ \tau^{-1} \rho^k & \text{if } \|\mathbf{s}_{rel}^k\| > \xi \varphi \|\mathbf{r}_{rel}^k\| \\ \rho^k & \text{otherwise,} \end{cases} \qquad (29) $$

where $\xi \in \mathbb{R}^+$ and the relative primal and dual residuals are defined by

$$ \begin{aligned} \mathbf{r}_{rel}^{k+1} &:= \frac{A\mathbf{x}^{k+1} + B\mathbf{y}^{k+1} - \mathbf{c}}{\max\{\|A\mathbf{x}^{k+1}\|, \|B\mathbf{y}^{k+1}\|, \|\mathbf{c}\|\}} \\ \mathbf{s}_{rel}^{k+1} &:= \frac{\rho A^T B (\mathbf{y}^{k+1} - \mathbf{y}^k)}{\|A^T \mathbf{w}^{k+1}\|} = \frac{A^T B (\mathbf{y}^{k+1} - \mathbf{y}^k)}{\|A^T \mathbf{\zeta}^{k+1}\|} \end{aligned} \qquad (30) $$

This approach avoids the need for explicit compensation for problem scaling when the formulation is modified. The stopping criteria in (19) is invariant to problem scaling when $\epsilon^{abs} = 0$.

The fixed multiplier $\tau$ is a potential weakness of the penalty update policies: if $\tau$ is small, then a large number of iterations may be required to reach appropriate $\rho$ value if $\rho^0$ is poorly chosen; on the other hand, if $\tau$ is large, the corrections to $\rho$ may be too large when $\rho$ is close to the optimal value. A straightforward solution is to adapt $\tau$ at each iteration

$$ \tau^{k+1} := \begin{cases} \sqrt{\xi^{-1} \frac{\|\mathbf{r}_{rel}^k\|}{\|\mathbf{s}_{rel}^k\|}} & \text{if } 1 \le \sqrt{\xi^{-1} \frac{\|\mathbf{r}_{rel}^k\|}{\|\mathbf{s}_{rel}^k\|}} < \tau_{max} \\ \sqrt{\xi \frac{\|\mathbf{s}_{rel}^k\|}{\|\mathbf{r}_{rel}^k\|}} & \text{if } \tau_{max}^{-1} < \sqrt{\xi^{-1} \frac{\|\mathbf{r}_{rel}^k\|}{\|\mathbf{s}_{rel}^k\|}} < 1 \\ \tau_{max} & \text{otherwise} \end{cases} \qquad (31) $$

where $\tau_{max}$ provides a bound on $\tau$, i.e., the convergence results in [25] still hold for this extension.

Xu et al. [27] propose to automate and speed up ADMM by using stepsize selection rules adapted from the gradient descent literature, namely the Barzilai-Borwein 'spectral' method, based in the dual of the ADMM in (8). The results of these techniques on our set of problems do not show an improvement of the convergence of the algorithm. The results will not be reported in the article.
## 2.2 Algorithm

The relaxed ADMM involves a parameter $\eta \in (0, 1)$. It is desirable to restart the method as infrequently as possible, it is recommended a value of $\eta$ close to 1 [20]. In our case $\eta = 0.999$ was used.

<caption>**Algorithm 1 ADMM.** Alternative Direction Method of Multipliers</caption>

Require: $\mathbf{y}^0 = \hat{\mathbf{y}}^0$, $\mathbf{z}^0 = \hat{\mathbf{z}}^0$, $\alpha_0 = 1$, $\eta \approx 1$, $\varphi = 10$, $\tau = 2$, $f = 1$
Require: with\_restart, with\_acceleration, with\_adaptive\_rho Boolean variables
1: Choose the initial value of $\rho^*$ with on the rules.
2: **for** $k = 0, 1, 2, \dots$ **do**
3: &nbsp;&nbsp;&nbsp;&nbsp; $\mathbf{x}^{k+1} \leftarrow \arg \min_{\mathbf{x}} L_{\rho_k}(\mathbf{x}; \hat{\mathbf{y}}^k; \hat{\mathbf{z}}^k)$
4: &nbsp;&nbsp;&nbsp;&nbsp; $\mathbf{y}^{k+1} \leftarrow \arg \min_{\mathbf{y}} L_{\rho_k}(\mathbf{x}^{k+1}; \mathbf{y}; \hat{\mathbf{z}}^k)$
5: &nbsp;&nbsp;&nbsp;&nbsp; $\mathbf{z}^{k+1} \leftarrow \frac{\rho^{k-1}}{\rho^k}(\mathbf{z}^k + A\mathbf{x}^{k+1} + B\mathbf{y}^{k+1} - \mathbf{c})$
6: &nbsp;&nbsp;&nbsp;&nbsp; Stop criterion in (19)
7: &nbsp;&nbsp;&nbsp;&nbsp; **if** with\_acceleration **then**
8: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; /* Relaxation defined in (21) */
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\alpha_{k+1} := \frac{1}{2}(1 + \sqrt{1 + 4\alpha_k^2})$
10: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\hat{\mathbf{y}}^{k+1} \leftarrow \mathbf{y}^{k+1} + \frac{\alpha_{k-1}}{\alpha_{k+1}}(\mathbf{y}^{k+1} - \mathbf{y}^k)$
11: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\hat{\mathbf{z}}^{k+1} \leftarrow \mathbf{z}^{k+1} + \frac{\alpha_{k-1}}{\alpha_{k+1}}(\mathbf{z}^{k+1} - \frac{\rho^{k-1}}{\rho^k}\mathbf{z}^k)$
12: &nbsp;&nbsp;&nbsp;&nbsp; **end if**
13: &nbsp;&nbsp;&nbsp;&nbsp; **if** with\_restart **then**
14: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; /* Restart defined in (23) */
15: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $e_k \leftarrow \rho^k \|\mathbf{z}^{k+1} - \frac{\rho^{k-1}}{\rho^k}\hat{\mathbf{z}}^k\|^2 + \rho^k \|B(\mathbf{y}^{k+1} - \hat{\mathbf{y}}^k)\|^2$
16: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **if** $e_k < \eta e_{k-1}$ **then**
17: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\alpha_{k+1} \leftarrow 1$; $\hat{\mathbf{y}}^{k+1} \leftarrow \mathbf{y}^{k+1}$; $\hat{\mathbf{z}}^{k+1} \leftarrow \mathbf{z}^{k+1}$; $e_k \leftarrow \frac{e_{k-1}}{\eta}$;
18: &nbsp;&nbsp;&nbsp;&nbsp; **end if**
19: &nbsp;&nbsp;&nbsp;&nbsp; **end if**
20: &nbsp;&nbsp;&nbsp;&nbsp; **if** with\_adaptive\_rho **then**
21: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **if** mod$(k, f) = 1$ **then**
22: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Update $\rho_{k+1}$ with one of the updating rules.
23: &nbsp;&nbsp;&nbsp;&nbsp; **else**
24: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\rho^{k+1} \leftarrow \rho^k$
25: &nbsp;&nbsp;&nbsp;&nbsp; **end if**
26: &nbsp;&nbsp;&nbsp;&nbsp; **end if**
27: **end for**
### 3 ADMM formulation of problem ($P_c$) and ($R_c$)

**ADMM formulation of problem ($P_c$)** Let $\delta_{K_{e,\mu}^*} : \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ denote the indicator function of $K_{e,\mu}^*$, i.e.,

$$ \delta_{K_{e,\mu}^*}(\boldsymbol{x}) = \begin{cases} 0 & \text{if } \boldsymbol{x} \in K_{e,\mu}^* \\ +\infty & \text{otherwise.} \end{cases} \qquad (32) $$

Then problem ($P_c$) in (6) is equivalently rewritten as follows:

$$ \underset{\boldsymbol{v}, \tilde{\boldsymbol{u}}}{\text{Minimize}} \quad \frac{1}{2} \boldsymbol{v}^{\top} M \boldsymbol{v} + \boldsymbol{f}^{\top} \boldsymbol{v} + \delta_{K_{e,\mu}^*}(\tilde{\boldsymbol{u}}) \qquad (33a) $$

$$ \text{subject to} \quad \tilde{\boldsymbol{u}} = H \boldsymbol{v} + \boldsymbol{w} + \phi(\boldsymbol{s}). \qquad (33b) $$

Therefore, in the notation of (8), we get $f(\boldsymbol{v}) = \frac{1}{2} \boldsymbol{v}^{\top} M \boldsymbol{v} + \boldsymbol{f}^{\top} \boldsymbol{v}$, $g(\tilde{\boldsymbol{u}}) = \delta_{K_{e,\mu}^*}(\tilde{\boldsymbol{u}})$, $A = H$, $B = -I$ and $\boldsymbol{c} = -\boldsymbol{w} - \phi(\boldsymbol{s})$. Problem (33) is a minimization problem of a convex function under linear equality constraints. Its augmented Lagrangian (in the scaled form) is given by

$$ L_{\rho}(\boldsymbol{v}; \tilde{\boldsymbol{u}}; \boldsymbol{\zeta}) = \frac{1}{2} \boldsymbol{v}^{\top} M \boldsymbol{v} + \boldsymbol{f}^{\top} \boldsymbol{v} + \delta_{K_{e,\mu}^{*}}(\tilde{\boldsymbol{u}}) + \frac{\rho}{2} \|H\boldsymbol{v} + \boldsymbol{w} + \phi(\boldsymbol{s}) + \boldsymbol{\zeta} - \tilde{\boldsymbol{u}}\|^2 - \frac{\rho}{2} \|\boldsymbol{\zeta}\|^2, \quad (34) $$

where $\rho$ is the penalty parameter, and $\boldsymbol{\zeta} \in \mathbb{R}^m$ is the scaled dual variable (the scaled Lagrange multiplier).

The ADMM for solving problem (33) consists of iterating the updates (see (15), (16), and (17))

$$ \boldsymbol{v}^{k+1} := \arg \min_{\boldsymbol{v}} L_{\rho}(\boldsymbol{v}; \tilde{\boldsymbol{u}}^k; \boldsymbol{\zeta}^k) \qquad (35) $$

$$ \tilde{\boldsymbol{u}}^{k+1} := \arg \min_{\tilde{\boldsymbol{u}}} L_{\rho}(\boldsymbol{v}^{k+1}; \tilde{\boldsymbol{u}}; \boldsymbol{\zeta}^k) \qquad (36) $$

$$ \boldsymbol{\zeta}^{k+1} := \boldsymbol{\zeta}^k + H \boldsymbol{v}^{k+1} + \boldsymbol{w} + \phi(\boldsymbol{s}) - \tilde{\boldsymbol{u}}^{k+1} \qquad (37) $$

The first step in (35) corresponds to unconstrained minimization of a convex quadratic function. It follows from the stationarity condition of this convex quadratic function that $\boldsymbol{v}^{k+1}$ is the solution to the following system of linear equations:

$$ [M + \rho H^{\top} H] \boldsymbol{v}^{k+1} = -\boldsymbol{f} + \rho H^{\top} (\tilde{\boldsymbol{u}}^k - \boldsymbol{w} - \phi(\boldsymbol{s}) - \boldsymbol{\zeta}^k) \qquad (38) $$

The second step in (36) can be computed independently for each $\alpha = 1...n_c$ as

$$ \tilde{\boldsymbol{u}}^{k+1} := \arg \min_{\tilde{\boldsymbol{u}}} \delta_{K_{e,\mu}^*}(\tilde{\boldsymbol{u}}) + \frac{\rho}{2} \|H\boldsymbol{v}^{k+1} + \boldsymbol{w} + \phi(\boldsymbol{s}) + \boldsymbol{\zeta}^k - \tilde{\boldsymbol{u}}\|^2 \qquad (39) $$

which is identified as the proximity operator or the Moreau-Yosida regularization of $\delta_{K_{e,\mu}^*}(\tilde{\boldsymbol{u}})$ and this can be written explicitly by using the projection onto the second-order cone as

$$ \tilde{\boldsymbol{u}}^{k+1} := \Pi_{K_{e,\mu}^*}(H\boldsymbol{v}^{k+1} + \boldsymbol{w} + \phi(\boldsymbol{s}) + \boldsymbol{\zeta}^k) \qquad (40) $$

moreover, it is recognized in (40) the fixed-point iteration solution of the VI formulation of (1), described by Acary in [1], which is related to the zeroes of the natural map, i.e.,

$$ F_{\rho}^{nat}(\tilde{\boldsymbol{u}}) = \tilde{\boldsymbol{u}} - \Pi_{K_{e,\mu}^*}(\tilde{\boldsymbol{u}} - \rho \boldsymbol{r}) = 0 \qquad (41) $$

where $\boldsymbol{r}$ is the contact reaction or impulse (see Appendix B for the relationship between $\boldsymbol{r}$ and $\boldsymbol{\zeta}$). This formulation helps us to notice the fixed-point iteration for the VI formulation is embedded into the ADMM algorithm.

Finally, it is worthy to note that for (38) we carry out a Cholesky factorization and (40) can be solved explicitly by using the formulae given in the next section.
**ADMM formulation of problem (R_c)** The problem (R_c) can be equivalently written in the form

$$ \text{Minimize} \quad \frac{1}{2} \mathbf{r}^\top W \mathbf{r} + \tilde{q}(\mathbf{s})^\top \mathbf{r} + \delta_K(z) \qquad (42a) $$

$$ \text{subject to} \quad z = r \qquad (42b) $$

Its augmented Lagrangian (in the scaled form) is formulated as

$$ L_{\rho}(\mathbf{r}; z; \zeta) = \frac{1}{2} W \mathbf{r} + \tilde{q}(\mathbf{s})^{\top} \mathbf{r} + \delta_{K}(z) + \frac{\rho}{2} \|\mathbf{r} + \zeta - z\|^{2} - \frac{\rho}{2} \|\zeta\|^{2}, \qquad (43) $$

where $\rho$ is the penalty parameter, and $\zeta \in \mathbb{R}^{d \times n}$ is the scaled dual variables (the scaled Lagrange multipliers).

The ADMM solving problem (33) consists of iterating the updates (see (??), (??), and (??))

$$ \mathbf{r}^{k+1} := \arg \min_{\mathbf{r}} L_{\rho}(\mathbf{r}; z^k; \zeta^k), \qquad (44) $$

$$ z^{k+1} := \arg \min_{z} L_{\rho}(\mathbf{r}^{k+1}; z; \zeta^{k}), \qquad (45) $$

$$ \zeta^{k+1} := \zeta^k + \mathbf{r}^{k+1} - z^{k+1}, \quad i = 1, \dots, n. \qquad (46) $$

The first step in (56) corresponds to unconstrained minimization of a convex quadratic function. It follows from the stationarity condition of this convex quadratic function that $\mathbf{r}^{k+1}$ is the solution to the following system of linear equations:

$$ [W + \rho I] \mathbf{r} = -\tilde{q}(\mathbf{s}) - \rho(\zeta^k - z^k). \qquad (47) $$

The second step in (57) can be computed independently for each $i = 1, \dots, n$ as

$$ z^{k+1} := \arg \min_{z} \delta_{K}(z) + \frac{\rho}{2} \|\mathbf{r}^{k+1} + \zeta^{k} - z\|^{2}. \qquad (48) $$

(64) can be rewritten by using the projection onto the second-order cone as

$$ z^{k+1} := \Pi_{K}(\mathbf{r}^{k+1} + \zeta^{k}) \qquad (49) $$

# 4 Numerical framework for the comparisons

Since the convergence is only proved for convex problems, and the convergence rate is problem dependent, the comparison of the variants of the ADMM algorithm is made on a large set of problems and we use perfomance profiles to draw conclusions. In this section, we explain the set of problems that we use, how we measure errors between different algorithms and how the performance profiles are computed.

## 4.1 The FCLib collection of examples.

The FCLib (Frictional Contact Library) is a open-source collection of problems stored in an hdf5 format described in [28] and publicly available at https://github.com/FrictionalContactLibrary/fclib-library. The library which contains problems in the standard from (P) or in its reduced form (R), is used in this paper in version v2.0. The test sets is illustrated in Figure 1.
<figure><img src="image_2.png" /><figcaption>Figure 1: Illustrations of the FClib test problems</figcaption></figure>
**Algorithm 2** Fast ADMM for problem with restart (33) condensed form.

Require: $y^0 = \tilde{y}^0$, $z^0 = \tilde{z}^0$, $\eta^0 = \tilde{\eta}^0$, $\zeta^0 = \tilde{\zeta}^0$, $\tau_0 = 1$, and $\rho > 0$.

1: **for** $k = 0, 1, 2, \dots$ **do**

2: $\mathbf{r}^{k+1}$ solves $[W + \rho I]\mathbf{r} = -\tilde{q}(\mathbf{s}) + \rho(\hat{\zeta}^k - \hat{z}^k)$

3: $\mathbf{z}^{k+1} := \Pi_K(\mathbf{r}^{k+1} + \hat{\zeta}^k)$

4: $\hat{\zeta}^{k+1} := \hat{\zeta}^k + \mathbf{r}^{k+1} - \mathbf{z}^{k+1}$

5:

6: $e_k := \|\zeta^{k+1} - \hat{\zeta}^k\|^2 + \|\mathbf{z}^{k+1} - \hat{\mathbf{z}}^k\|^2$

7: **if** $e_k < \eta e_{k-1}$ **then**

8: $\tau_{k+1} := \frac{1}{2}(1 + \sqrt{1 + 4\tau_k^2})$

9: $\hat{\mathbf{z}}^{k+1} := \mathbf{z}^{k+1} + \frac{\tau_k - 1}{\tau_{k+1}}(\mathbf{z}^{k+1} - \mathbf{z}^k)$

10: $\hat{\zeta}_i^{k+1} := \zeta_i^{k+1} + \frac{\tau_k - 1}{\tau_{k+1}}(\zeta_i^{k+1} - \zeta_i^k)$

11: **else**

12: $\tau_{k+1} := 1$

13: $\hat{\mathbf{z}}^{k+1} = \mathbf{z}^k$

14: $\hat{\zeta}_i^{k+1} := \zeta_i^k$

15: $e_k \leftarrow e_{k-1}/\eta$

16: **end if**

17: **end for**

## 4.2 Measuring errors

The stopping criteria of the variants of the ADMM is based on (19). Nevertheless, to get a uniform measure of the error for all the numerical methods we compare, we use the following measure

$$ \text{Error}(v, r) = \frac{\|Mv + f - H^\top r\|}{\max(\|Mv\|, \|H^\top r\|, \|f\|)} + \frac{\|r - \Pi_{K_\mu}(r - (Hv + w + \Phi(Hv + w)))\|}{\max(\|r\|, \|Hv + w\|)} \quad (50) $$

## 4.3 Comparison with performance profiles

In this section, following [1] we explain the concept of performance profiles that was introduced in [29] for benchmarking optimization solvers on a large set of problems. For a set $P$ of $n_p$ problems, and a set $S$ of $n_s$ solvers, we define a performance criterion for a solver $s$, a problem $p$ and a required precision $tol$ by

$$ t_{p,s} = \text{computing time required for } s \text{ to solve } p \text{ at precision } tol \quad (51) $$

A performance ratio over all the solvers is defined by

$$ r_{p,s} = \frac{t_{p,s}}{\min\{t_{p,s}, s \in S\}} \ge 1 \quad (52) $$

For $\tau > 1$, we define a distribution function $\rho_s$ for the performance ratio for a solvers as

$$ \rho_s(\tau) = \frac{1}{n_p} \mathrm{card}\{p \in P, r_{p,s} \le \tau\} \le 1 \quad (53) $$

This distribution computes the number of problems $p$ that are solved with a performance ratio below a given threshold $\tau$. In other words, $\rho_s(\tau)$ represents the probability that the solver $s$ has a performance ratio not larger than a factor $\tau$ of the best solver. It is worth noting that $\rho_s(1)$
represents the probability that the solver *s* beats the other solvers, and $\rho_s(\tau)$ characterizes the robustness of the method for large values of $\tau$. To summarize: the higher $\rho_s$ is, the better the method is. In the sequel, the term performance profile denotes a graph of the functions $\rho_s(\tau)$, $\tau > 1$. The computational time is used to measure performance in the algorithms.

## 4.4 Implementation details.

## 4.5 Nomenclature for the algorithms

For the sake of readability, we introduce a nomenclature for the variants of the algorithms in Table ?? shows the notation used in the results:

<table><caption>Table 1: Nomenclature for the algorithms.</caption><thead><tr><th>Name</th><th>Algorithm</th><th>Additional informations</th></tr></thead><tbody><tr><td>ADMM-CST</td><td>Algorithm 1</td><td>Constant ρ with ρ* = 1</td></tr><tr><td>ADMM-CST-NO</td><td>Algorithm 1</td><td>Constant ρ with ρ* = 1<br>No acceleration and restart</td></tr><tr><td>ADMM-CST-NORM-INF</td><td>Algorithm 1</td><td>Constant ρ with ρ* given by norm inf</td></tr><tr><td>ADMM-BR</td><td>Algorithm 1</td><td>Residual balancing for ρ<sub>k</sub> with ρ* = 1</td></tr><tr><td>ADMM-BR-NO</td><td>Algorithm 1</td><td>Residual balancing for ρ<sub>k</sub> with ρ* = 1<br>No acceleration and restart</td></tr><tr><td>ADMM-BR-SCALED</td><td>Algorithm 1</td><td>Scaled Problem. Residual balancing for ρ<sub>k</sub> with ρ* = 1</td></tr><tr><td>ADMM-SBR</td><td>Algorithm 1</td><td>Scaled residual balancing for ρ<sub>k</sub> with ρ* = 1</td></tr></tbody></table>
# 5 Numerical comparisons for the convex problems

## 5.1 Convergence rate

To start with, we compare the convergence rate of variants of algorithms for ADMM in the convex case. The example Spheres-i1000-3196-8516, chosen in the FCLib collection has 3196 contacts with $\mu = 0.7$ and 12000 degrees of freedom and corresponds to a pile of 2000 spheres into a box. Similar results are obtained for other examples. The user tolerance is set $10^{-12}$. The results

<figure><img src="image_3.png" /><figcaption>Figure 2: Convergence rate of variants of ADMM for the convex problem ($P_c$).</figcaption></figure>

are depicted in Figure 2. We observe that the acceleration technique improves a lot the rate of convergence as expected. Furthermore, the update of the penalty parameter $\rho$ with the residual balancing method speeds up the convergence. We can also observe that the residual balancing renders the convergence non-monotone. When the residual balancing technique is used with the acceleration, a sensitive part of the iterations are restarted due to the increase of the residual.

**Influence of the initial choice for $\rho^*$** Let us now consider the effect of the initial choice of $\rho$, denoted by $\rho^*$ on the same example. In Figure 3a, the convergence rate is depicted with various choices of $\rho^*$ without acceleration and restart. This choice of $\rho^*$ has a strong influence on the convergence rate:

* quite surprisingly, the best result is obtained with $\rho^* = 1.0$; this may be explained by the fact
that the problem is well-scaled in the sense that the 1-norm and the $\infty$-norm of the matrices $H$ and $M$ are of order 1;

* the choice based on the computation of the eigenvalues of $M$ as in (25) that gives $\rho^* \approx 0.44$ and the norm of the matrices $\rho^* = 0.17$ for the equation (26) gives also good results

* choosing relatively large or small values for $\rho^*$ entails the rate of convergence;

This results are in accordance with the literature [21] where an optimal choice of $\rho^*$ must be found.

In 3b, we also observe that this choice an influence on the long term behaviour of the algorithm when the acceleration and restart technique is used. This result are more difficult to explain since there is an interplay between the choice of $\rho^*$ and the acceleration strategy.

## 5.2 Comparison of ADMM variants for the convex problems
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3: Influence of the initial choice of $\rho_*$ for ADMM without adaptive $\rho$ strategy</figcaption></figure>
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>Figure 4: Performance profile for Spheres example.</figcaption></figure>
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>Figure 5: Performance profile for Spheresimm example.</figcaption></figure>
# 6 Solving the original nonconvex problems (P) and (R).

To solve the friction problem, the parameter $s$ in problem (6) must be updated. One obvious way is that, once we solve problem (6) with fixed $s$ by ADMM, then update $s$ by using the obtained solution, and repeat this procedure. This is what we called *external update*. Another possibility is to update $s$ at each iteration of the ADMM. This is what we called *internal update*. Intuitively, the latter saves the total computational cost, but stability of the algorithm in this is not clear, there is no proof of existence or uniqueness in comparison to the first one [16].

## 6.1 Initial guess

Another issue is the initial guess of $s$ to improve the rate of convergence. One simple method that performs well in our tests (results not shown) is $s^0 = 0 \in \mathbb{R}^m$.

## 6.2 Internal update

For the next development, we consider

$$ \Phi(\mathbf{u}) = [\mu^{\alpha} \|E_t \mathbf{u}^{\alpha} \|_{\hat{e}_n}, \alpha = 1 \dots n_c]^T = [\mu^{\alpha} \|E_t (H\mathbf{v} + \mathbf{w})^{\alpha} \|_{\hat{e}_t}, \alpha = 1 \dots n_c]^T \quad (54) $$

therefore $\Phi(\mathbf{u})$ could be seen as $\Phi(\mathbf{v})$. Then, the problem (34) is a minimization problem of a convex function under non-linear equality constraints. Its augmented Lagrangian (in the scaled form) is formulated as

$$ L_{\rho}(\mathbf{v}; \tilde{\mathbf{u}}; \zeta) = \frac{1}{2} \mathbf{v}^{\top} M \mathbf{v} + \mathbf{f}^{\top} \mathbf{v} + \delta_{K_{e,\mu}^{*}}(\tilde{\mathbf{u}}) + \frac{\rho}{2} \|H\mathbf{v} + \mathbf{w} + \Phi(\mathbf{v}) + \zeta - \tilde{\mathbf{u}}\|^2 - \frac{\rho}{2} \|\zeta\|^2 \quad (55) $$

where $\rho$ is the penalty parameter, and $\zeta \in \mathbb{R}^m$ is the scaled dual variable (the scaled Lagrange multiplier). The ADMM solving problem (33) consists of iterating the updates

$$ \mathbf{v}^{k+1} := \arg \min_{\mathbf{v}} L_{\rho}(\mathbf{v}; \tilde{\mathbf{u}}^k; \zeta^k) \quad (56) $$

$$ \tilde{\mathbf{u}}^{k+1} := \arg \min_{\tilde{\mathbf{u}}} L_{\rho}(\mathbf{v}^{k+1}; \tilde{\mathbf{u}}; \zeta^k) \quad (57) $$

$$ \zeta^{k+1} := \zeta^k + H\mathbf{v}^{k+1} - \tilde{\mathbf{u}}^{k+1} + \Phi(\mathbf{v}^{k+1}) \quad (58) $$

The first step in (56) corresponds to unconstrained minimization of a convex function. It follows from the stationarity condition of this convex function that $\mathbf{v}^{k+1}$ is the solution to the following system of non-linear equations:

$$ [M + \rho (H^{\top} + \partial_{\mathbf{v}} \Phi(\mathbf{v}^{k+1})) H] \mathbf{v}^{k+1} = -\mathbf{f} + \rho (H^{\top} + \partial_{\mathbf{v}} \Phi(\mathbf{v}^{k+1})) (\tilde{\mathbf{u}}^k - \mathbf{w} - \Phi(\mathbf{v}^{k+1}) - \zeta^k) \quad (59) $$

this non-linearity is difficult to handle, however, if the previous iteration is considered in the $\Phi(\mathbf{v})$ term, we get the following system of linear equations:

$$ [M + \rho (H^{\top} + \partial_{\mathbf{v}} \Phi(\mathbf{v}^k)) H] \mathbf{v}^{k+1} = -\mathbf{f} + \rho (H^{\top} + \partial_{\mathbf{v}} \Phi(\mathbf{v}^k)) (\tilde{\mathbf{u}}^k - \mathbf{w} - \Phi(\mathbf{v}^k) - \zeta^k) \quad (60) $$

which follows the same structure of (38) and allows us to perform a LU or Cholesky factorization. If $\|E_t \mathbf{u}^{\alpha}\| \neq 0$

$$ \partial_{\mathbf{v}} \Phi(\mathbf{v}) = \left[ \mu^{\alpha} \hat{\mathbf{e}}_{n}^{\top} \frac{E_t (H \mathbf{v} + \mathbf{w})^{\alpha}}{\|E_t (H \mathbf{v} + \mathbf{w})^{\alpha}\|} E_t^{\top} H^{\alpha^{\top}}, \alpha = 1 \dots n_c \right]^{\mathrm{T}} = \Psi \quad (61) $$
Therefore

$$[M + \rho (H^{\top} + \Psi) H] \mathbf{v}^{k+1} = -\mathbf{f} + \rho (H^{\top} + \Psi) (\tilde{\mathbf{u}}^k - \mathbf{w} - \Phi(\mathbf{v}^k) - \mathbf{\zeta}^k) \quad (62)$$

Note that if in (56) it is evaluated $\Phi(\mathbf{v})$ in the previous iteration before compute the stationary condition, it results in

$$[M + \rho H^{\top} H] \mathbf{v}^{k+1} = -\mathbf{f} + \rho H^{\top} (\tilde{\mathbf{u}}^k - \mathbf{w} - \Phi(\mathbf{v}^k) - \mathbf{\zeta}^k) \quad (63)$$

which is a pseudo-update of (56) in comparison with (62), but it aims to be faster without the computation of $\Psi$, therefore, (63) will be used.

The second step in (57) can be computed independently for each $\alpha = 1...n_c$ as

$$\tilde{\mathbf{u}}^{k+1} := \arg \min_{\tilde{\mathbf{u}}} \delta_{K_{e,\mu}^*} (\tilde{\mathbf{u}}) + \frac{\rho}{2} \|H\mathbf{v}^{k+1} + \mathbf{w} + \Phi(\mathbf{v}^k) + \mathbf{\zeta}^k - \tilde{\mathbf{u}}\|^2 \quad (64)$$

(64) can be rewritten by using the projection onto the second-order cone as

$$\tilde{\mathbf{u}}^{k+1} := \Pi_{K_{e,\mu}^*} (H\mathbf{v}^{k+1} + \mathbf{w} + \mathbf{\zeta}^k + \Phi(\mathbf{v}^{k+1})) \quad (65)$$

## 6.3 External update

In this case, the external error measurement is given for each $\alpha = 1...n_c$ by

$$\left[ \frac{\Phi(\mathbf{v}^{k+1}) - \Phi(\mathbf{v}^k)}{\Phi(\mathbf{v}^k)} \right]^\alpha \le \epsilon^{ext} \quad (66)$$

where $\epsilon^{ext}$ is the external tolerance, with a value $10^{-3}$.

## 6.4 SICONOS error

The primal residual measures the error in the constraint in (1), but also the error could be measured from the dynamic equation and the natural map (41) as

$$\epsilon^{SIC} = \sqrt{\|M\mathbf{v} + \mathbf{f} - H^{\top}\mathbf{r}\|^2 + \|\tilde{\mathbf{u}} - \Pi_{K_{e,\mu}^*}(\tilde{\mathbf{u}} - \rho\mathbf{r})\|^2} \quad (67)$$
## 6.5 Numerical experience for original nonsmooth and nonconvex problems.

In most of the test examples the behavior depicted in Figure 5-7 is observed. This is analyzed with the best solver of Section 7 and the test example *Box Stacks-i1000-352-6*.

In Figure 5 is observed a convergence after 7 external iterations. Note that the first value is not shown because it is zero.

<figure><img src="image_12.png" /><figcaption>Figure 6: External update with vp-RR-He (Di Cairano).</figcaption></figure>
In Figure 6 is observed a convergence after 43 internal iterations. From iteration 25 the same value holds, but the stop criterion (19) is not satisfied.

<figure><img src="image_13.png" /><figcaption>Figure 7: Internal update with vp-RR-He (Di Cairano).</figcaption></figure>
In Figure 7 is depicted the same internal update, but the stop criterion (19) is changed for the external one (66). It is observed that $||\Phi(v)||$ holds from iteration 25, but $\Phi(v)^{\alpha}$ is not the same from iteration to iteration.

<figure><img src="image_14.png" /><figcaption>Figure 8: Internal update with vp-RR-He (Di Cairano) and stop criterion in (66).</figcaption></figure>

In Table 2 are shown the values of the final $||\Phi(s)||$ with its respective first and last contact vector $\Phi(s)^{\alpha}$:

<table><caption>Table 2: Comparison of s-update</caption><thead><tr><th></th><th>External</th><th>Internal</th><th>Internal-s</th></tr></thead><tbody><tr><td>||Φ(s)||</td><td>7.10e-02</td><td>7.10e-02</td><td>7.10e-02</td></tr><tr><td>Φ(s)<sup>α=1</sup></td><td>[8.31e-06 0.00e+00 0.00e+00]</td><td>[6.49e-06 0.00e+00 0.00e+00]</td><td>[7.45e-06 0.00e+00 0.00e+00]</td></tr><tr><td>Φ(s)<sup>α=n<sub>c</sub></sup></td><td>[1.02e-05 0.00e+00 0.00e+00]</td><td>[1.19e-05 0.00e+00 0.00e+00]</td><td>[7.85e-06 0.00e+00 0.00e+00]</td></tr></tbody></table>

Both external and internal update converges to the same value of $||\Phi(s)||$, but this is an apparent result of uniqueness when $\Phi(s)^{\alpha}$ is compared in each case. Therefore, the internal update does not guarantee the same result for $\Phi(s)$. However, all of them have low values of SICONOS error and residuals, i.e., they solve the problem (1). This opens the question if these approaches converge at some level to a different solutions of problem (1).
# 7 Numerical experiments.

## 7.1 Standard first order techniques.

## 7.2 Convergence rate

<figure><img src="image_15.png" /><figcaption>Figure 9: Convergence rate of variants of ADMM for the problem (P) for $\mu = 0.7$.</figcaption></figure>

## 7.3 Comparisons of variants of ADMM solvers for the global problem (P)
<figure><img src="image_16.png" /><figcaption>Figure 10: Convergence rate of variants of ADMM for the problem (P) for $\mu = 0.3$.</figcaption></figure>
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /><figcaption>Figure 11: Performance profile for Spheres example.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure 12: Performance profile for Spheres1mm example.</figcaption></figure>
<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /><figcaption>Figure 13: Performance profile for KaplasTower example.</figcaption></figure>
<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /><figcaption>Figure 14: Performance profile for Chute example.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Figure 15: Comparison of ADMM-$\star$ algorithms for the global problem (P).</figcaption></figure>
reduced coefficient of friction

7.4 Comparisons of variants of ADMM solvers
<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /><figcaption>Figure 16: Performance profile for Spheres example with $\mu = 0.3$.</figcaption></figure>
<figure><img src="image_31.png" /><figcaption>Figure 17: Comparison of ADMM-$\star$ algorithms for the reduced problem (R).</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>Figure 17: Comparison of ADMM-$\star$ algorithms.</figcaption></figure>

ADM
ADMM-SCALE
<figure><img src="image_33.png" /><figcaption>Figure 17: Comparison of ADMM-⋆ algorithms.</figcaption></figure>
## 7.5 Comparisons of variants of ADMM solvers with NSGS

<figure><img src="image_34.png" /><figcaption>Figure 18: Comparison of ADMM-$\star$ algorithms with NSGS</figcaption></figure>
<figure><img src="image_35.png" /></figure>

<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /></figure>

<figure><img src="image_38.png" /></figure>

<figure><img src="image_39.png" /></figure>

<figure><img src="image_40.png" /></figure>

<figure><img src="image_41.png" /><figcaption>Figure 18: Comparison of ADMM-$\star$ algorithms.</figcaption></figure>
<figure><img src="image_42.png" /></figure>

<figure><img src="image_43.png" /></figure>

<figure><img src="image_44.png" /></figure>

<figure><img src="image_45.png" /></figure>

<figure><img src="image_46.png" /></figure>

<figure><img src="image_47.png" /></figure>

<figure><img src="image_48.png" /></figure>

<figure><img src="image_49.png" /><figcaption>Figure 18: Comparison of ADMM-$\star$ algorithms.</figcaption></figure>
References

[1] Vincent Acary, Maurice Brémond, and Olivier Huber. On solving contact problems with Coulomb friction: formulations and numerical comparisons. Research Report RR-9118, INRIA, November 2017.

[2] P.D. Panagiotopoulos. A nonlinear programming approach to the unilateral contact-, and friction-boundary value problem in the theory of elasticity. *Ingenieur-Archiv*, 44(6):421-432, 1975.

[3] J. Nečas, J. Jarušek, and J. Haslinger. On the solution of the variational inequality to the Signorini problem with small friction. *Bollettino U.M.I.*, 5(17-B):796-811, 1980.

[4] J. Haslinger. Approximation of the signorini problem with friction, obeying the coulomb law. *Mathematical Methods in the Applied Sciences*, 5:422-437, 1983.

[5] J. Haslinger. Least square method for solving contact problems with friction obeying coulomb's law. *Applications of mathematics*, 29(3):212-224, 1984.

[6] J. Haslinger and P. D. Panagiotopoulos. The reciprocal variational approach to the signorini problem with friction. approximation results. *Proceedings of the Royal Society of Edinburgh: Section A Mathematics*, 98:365-383, 1984.

[7] G. Del Piero and F. Maceri, editors. *Unilateral Problems in Structural Analysis*, volume 288 of CISM courses and lectures, Ravello, Italy, September 22-24 1983. Springer.

[8] G. Del Piero and F. Maceri, editors. *Unilateral Problems in Structural Analysis – II.*, volume 304 of CISM courses and lectures, Prescudin, Italy, June 17–20 1985. Springer.

[9] Michael G. Katona. A simple contact–friction interface element with applications to buried culverts. *International Journal for Numerical and Analytical Methods in Geomechanics*, 7(3):371–384, 1983.

[10] A. Chaudhary and K.J. Bathe. A solution method for static and dynamic analysis of three-dimensional contact problems with friction. *Computers & Structures*, 24(6):855-873, 1986.

[11] M. Jean and J.J. Moreau. Dynamics in the presence of unilateral contacts and dry friction: a numerical approach. In G. Del Pietro and F. Maceri, editors, *Unilateral problems in structural analysis. II*, pages 151–196. CISM 304, Spinger Verlag, 1987.

[12] E.N. Mitsopoulou and I.N. Doudoumis. On the solution of the unilateral contact frictional problem for general static loading conditions. *Computers & Structures*, 30(5):1111-1126, 1988.

[13] Vincent Acary, Maurice Brémond, and Olivier Huber. On solving contact problems with Coulomb friction: formulations and numerical comparisons. In Leine, R, V. Acary, and O. Bruls, editors, *Advanced Topics in Nonsmooth Dynamics - Transactions of the European Network for Nonsmooth Dynamics*, pages 375-457. Springer International Publishing, June 2018.

[14] Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, Jonathan Eckstein, et al. Distributed optimization and statistical learning via the alternating direction method of multipliers. *Foundations and Trends in Machine Learning*, 3(1):1-122, 2011.
[15] Vincent Acary and Florent Cadoux. *Applications of an Existence Result for the Coulomb Friction Problem*, pages 45–66. Springer Berlin Heidelberg, Berlin, Heidelberg, 2013.

[16] Vincent Acary, Florent Cadoux, Claude Lemaréchal, and Jérôme Malick. A formulation of the linear discrete coulomb friction problem via convex optimization. *ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik*, 91(2):155–175, 2011.

[17] V. Acary and B. Brogliato. *Numerical methods for nonsmooth dynamical systems. Applications in mechanics and electronics*. Lecture Notes in Applied and Computational Mechanics 35. Berlin: Springer. xxi, 525 p. , 2008.

[18] V. Acary and F. Cadoux. *Recent Advances in Contact Mechanics, Stavroulakis, Georgios E. (Ed.)*, volume 56 of *Lecture Notes in Applied and Computational Mechanics*, chapter Applications of an existence result for the Coulomb friction problem. Springer Verlag, 2013.

[19] V. Acary, F. Cadoux, C. Lemaréchal, and J. Malick. A formulation of the linear discrete coulomb friction problem via convex optimization. *ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik*, 91(2):155–175, 2011.

[20] Tom Goldstein, Brendan O'Donoghue, Simon Setzer, and Richard Baraniuk. Fast alternating direction optimization methods. *SIAM Journal on Imaging Sciences*, 7(3):1588–1623, 2014.

[21] Euhanna Ghadimi, André Teixeira, Iman Shames, and Mikael Johansson. Optimal parameter selection for the alternating direction method of multipliers (admm): quadratic problems. *IEEE Transactions on Automatic Control*, 60(3):644–658, 2015.

[22] Arvind U Raghunathan and Stefano Di Cairano. Alternating direction method of multipliers for strictly convex quadratic programs: Optimal parameter selection. In *American Control Conference (ACC)*, 2014, pages 4324–4329. IEEE, 2014.

[23] Wei Deng and Wotao Yin. On the global and linear convergence of the generalized alternating direction method of multipliers. *Journal of Scientific Computing*, 66(3):889–916, 2016.

[24] Jonathan Eckstein and W Yao. Augmented lagrangian and alternating direction methods for convex optimization: A tutorial and some illustrative computational results. *RUTCOR Research Reports*, 32:3, 2012.

[25] B. S. He, H. Yang, and S. L. Wang. Alternating direction method with self-adaptive penalty parameters for monotone variational inequalities. *Journal of Optimization Theory and Applications*, 106(2):337–356, Aug 2000.

[26] Brendt Wohlberg. Admm penalty parameter selection by residual balancing. *arXiv preprint arXiv:1704.06209*, 2017.

[27] Zheng Xu, Mário AT Figueiredo, and Tom Goldstein. Adaptive admm with spectral penalty parameter selection. *arXiv preprint arXiv:1605.07246*, 2016.

[28] V. Acary, M. Brémond, T. Koziara, and F. Pérignon. FCLIB: a collection of discrete 3D Frictional Contact problems. Technical Report RT-0444, INRIA, February 2014.

[29] Elizabeth D. Dolan and Jorge J. Moré. Benchmarking optimization software with performance profiles. *Mathematical Programming*, 91(2):201–213, Jan 2002.

[30] Shunsuke Hayashi, Nobuo Yamashita, and Masao Fukushima. A combined smoothing and regularization method for monotone second-order cone complementarity problems. *SIAM Journal on Optimization*, 15(2):593–615, 2005.
[31] Masao Fukushima, Zhi-Quan Luo, and Paul Tseng. Smoothing functions for second-order-cone complementarity problems. *SIAM Journal on optimization*, 12(2):436–460, 2002.

## A Projection onto dual second-order cone

For vector $x = (x_1, x_2) \in \mathbb{R} \times \mathbb{R}^{n-1}$, its spectral factorization with respect to $K_{e,\mu}^*$ by means of Jordan algebra (in a general SOCC-function) is defined by [30]

$$x = \lambda_1 u^1 + \lambda_2 u^2 \qquad (68)$$

Here, $\lambda_1, \lambda_2 \in \mathbb{R}$ are the spectral values given by

$$\lambda_i = x_1 + (-1)^i \mu^{(-1)^{i+1}} \|x_2\| \qquad (69)$$

and $u^1, u^2 \in \mathbb{R}^n$ are the spectral vectors given by

$$u^i = \begin{cases} \frac{\mu^2}{1+\mu^2} \left[ (-1)^i \frac{\mu^{2(i-2)}}{\mu} x_2 / \|x_2\| \right] & \text{if } x_2 \neq \mathbf{0} \\ \frac{\mu^2}{1+\mu^2} \left[ (-1)^i \frac{\mu^{2(i-2)}}{\mu} \omega \right] & \text{if } x_2 = \mathbf{0} \end{cases} \qquad (70)$$

with $\omega \in \mathbb{R}^{n-1}$ satisfying $\|\omega\| = 1$.

For $x \in \mathbb{R}^n$, let $\Pi_{K_{e,\mu}^*}(x) \in \mathbb{R}^n$ denote the projection of $x$ onto $K_{e,\mu}^*$, i.e.,

$$\Pi_{K_{e,\mu}^*}(x) = \arg \min \{\|x' - x\| \mid x' \in K_{e,\mu}^*\}. \qquad (71)$$

This can be computed explicitly as [31]

$$\Pi_{K_{e,\mu}^*}(x) = \max\{0, \lambda_1\}u^1 + \max\{0, \lambda_2\}u^2 \qquad (72)$$

Therefore the projection of $x$ onto $K_{e,\mu}^*$ could be written as follows

$$\Pi_{K_{e,\mu}^*}(x) = \begin{cases} 0 & \text{if } -x \in K_{e,\mu} \rightarrow \lambda_i \le 0 \\ x & \text{if } x \in K_{e,\mu}^* \rightarrow \lambda_i \ge 0 \\ \frac{\mu^2 \left(x_1 + \frac{1}{\mu} \|x_2\|\right)}{1 + \mu^2} \left[ \frac{1}{\mu} x_2 / \|x_2\| \right] & \text{if } -x \notin K_{e,\mu} \land x \notin K_{e,\mu}^* \rightarrow \lambda_1 < 0 \land \lambda_2 > 0, \end{cases} \qquad (73)$$

## B Convex formulation without initial affine constraints

The problem reads

$$\text{Minimize} \quad \frac{1}{2} r^\top W r + r^\top (q + s) \qquad (74a)$$

$$\text{subject to} \quad r \in K_{e,\mu} \qquad (74b)$$

This problem can be rewritten in ADMM form introducing a lack variable

$$\text{Minimize} \quad \frac{1}{2} r^\top W r + r^\top (q + s) \qquad (75a)$$

$$\text{subject to} \quad r = p \qquad (75b)$$

$$p \in K_{e,\mu} \qquad (75c)$$
This problem is a minimization problem of a convex function under linear equality constraints. Its augmented Lagrangian (in the scaled form) is formulated as

$$L_{\rho}(\mathbf{r}; \mathbf{p}; \zeta) = \frac{1}{2} \mathbf{r}^{\top} W \mathbf{r} + \mathbf{r}^{\top} (\mathbf{q} + \mathbf{s}) + \delta_{K_{e,\mu}}(\mathbf{p}) + \frac{\rho}{2} \|\mathbf{r} - \mathbf{p} + \zeta\|^2 - \frac{\rho}{2} \|\zeta\|^2 \quad (76)$$

where $\rho$ is the penalty parameter, and $\zeta \in \mathbb{R}^m$ is the scaled dual variable (the scaled Lagrange multiplier). The ADMM solving problem consists of iterating the updates

$$\mathbf{r}^{k+1} := \arg \min_{\mathbf{r}} L_{\rho}(\mathbf{r}; \mathbf{p}^k; \zeta^k), \quad (77)$$

$$\mathbf{p}^{k+1} := \arg \min_{\mathbf{p}} L_{\rho}(\mathbf{r}^{k+1}; \mathbf{p}; \zeta^k), \quad (78)$$

$$\zeta^{k+1} := \zeta^k + (\mathbf{r}^{k+1} - \mathbf{p}^{k+1}). \quad (79)$$

The first step corresponds to unconstrained minimization of a convex quadratic function. It follows from the stationarity condition of this convex quadratic function that $\mathbf{v}^{k+1}$ is the solution to the following system of linear equations:

$$[W + \rho I] \mathbf{r}^{k+1} = -(\mathbf{q} + \mathbf{s}) + \rho(\mathbf{p}^k - \zeta^k). \quad (80)$$

The second step can be computed independently for each $\alpha = 1...n_c$ as

$$\mathbf{p}^{k+1} := \arg \min_{\mathbf{p}} \delta_{K_{e,\mu}}(\mathbf{p}) + \frac{\rho}{2} \|\mathbf{r}^{k+1} - \mathbf{p} + \zeta^k\|^2. \quad (81)$$

it can be rewritten by using the projection onto the second-order cone as

$$\mathbf{p}^{k+1} := \Pi_{K_{e,\mu}}(\mathbf{r}^{k+1} + \zeta^k) \quad (82)$$

It is worthy to note that for (74) we carry out a Cholesky or LU factorization and (76) can be solved explicitly. All this procedure is simpler than the one applied in Section 3.

## C Optimality conditions of problem (6)

The Lagrangian reads

$$L(v; \tilde{u}; \lambda) = \frac{1}{2} v^{\top} M v + f^{\top} v + \delta_{K_{e,\mu}^{*}}(\tilde{u}) + \lambda(Hv + w + \Phi(s) - \tilde{u}) \quad (83)$$

where $\lambda$ is the unscaled Lagrange multiplier, i.e., $\rho\zeta$. Its gradient (or subgradient) is

$$\partial_v L(v; \tilde{u}; \lambda) = M v + f + H^{\top} \lambda \quad (84)$$

$$\partial_u L(v; \tilde{u}; \lambda) = \partial_{\tilde{u}} \delta_{K_{e,\mu}^*}(\tilde{u}) - \lambda \quad (85)$$

$$\partial_{\lambda} L(v; \tilde{u}; \lambda) = H v + w + \Phi(s) - \tilde{u} \quad (86)$$

Therefore, the optimality conditions are

$$M v + f = -H^{\top} \lambda \quad (87)$$

$$\lambda \in \partial_{\tilde{u}} \delta_{K_{e,\mu}^*}(\tilde{u}) \quad (88)$$

$$\tilde{u} = H v + w + \Phi(s) \quad (89)$$

But, by definition (see [1])

$$-\mathbf{r} \in \partial_{\tilde{u}} \delta_{K_{e,\mu}^*}(\tilde{u}) \quad (90)$$

Finally

$$\mathbf{r} = -\lambda = -\rho\zeta \quad (91)$$