Vincent Acary
email: vincent.acary@inria.fr

Paul Armand
email: paul.armand@unilim.fr

Minh Hoang

Nguyen

Maksym Shpakovych
email: maksym.shpakovych@unilim.fr

Second order cone programming for frictional contact mechanics using interior point algorithm

Keywords: Second order cone optimization, Interior-point algorithm, Contact mechanics, Coulomb friction

We report experiments of an implementation of a primal-dual interior point algorithm for solving mechanical models of one-sided contact problems with Coulomb friction. The objective is to recover an optimal solution with high accuracy and as quickly as possible. These developments are part of the design of Siconos 1 , an open-source software for modeling and simulating non-smooth dynamical systems. Currently, Siconos uses mainly first order methods for the numerical solution of these systems. These methods are very robust, but suffer from a linear rate of convergence and are therefore too slow to find accurate solutions in a reasonable time. Our main objective is to apply second-order methods to speed up convergence. In this paper, we focus on solving a relaxed formulation of the initial mechanical model, corresponding to the optimality conditions of a convex quadratic minimization problem with second-order cone constraints. We will present in detail a primal-dual interior point algorithm for solving this type of problem. The main difficulty in implementing this algorithm arises from the fact that at each iteration of the algorithm, a change of variable, called a scaling, has to be performed in order to guarantee the non-singularity of the linear system to be solved, as well as to recover a symmetric system. Although this scaling strategy is very nice from a theoretical point of view, it leads to an enormous deterioration of the conditioning of the linear system when approaching the optimal solution and therefore to all the numerical difficulties that result from it. We will describe in detail the numerical algebra that we have developed in our implementation, in order to overcome these problems of numerical instability. We will also present the solution of the models resulting from the problems of rolling friction, for which the cone of constraints is no longer self-dual like the Lorentz cone.

Introduction

Contact problems with dry Coulomb friction are present in many design and validation processes in mechanical engineering. As soon as several objects are involved, the question of computing contact forces arises. Examples include multibody systems and mechanisms [START_REF] Brogliato | Nonsmooth mechanics[END_REF][START_REF] Pfeiffer | Multibody dynamics with unilateral contacts[END_REF], robotic systems and grasping problems [START_REF] Bicchi | Robotic grasping and contact: A review[END_REF][START_REF] Buss | Recursive algorithms for real-time grasping force optimization[END_REF][START_REF] Buss | Dextrous hand grasping force optimization[END_REF][START_REF] Han | Grasp analysis as linear matrix inequality problems[END_REF], deformable solid mechanics [START_REF] Johnson | Contact mechanics[END_REF][START_REF] Popov | Contact mechanics and friction[END_REF]Wriggers, 2006), and granular materials [START_REF] Cambou | Micromechanics of granular materials[END_REF][START_REF] Radjai | Discrete-element modeling of granular materials[END_REF]. A common model of one-sided contact is to consider a set of inequality constraints on the configuration parameters of the mechanical problem (positions, rotations) associated with forces which must themselves be positive. Coulomb's friction, which governs how objects slide relative to each other, imposes a constraint on the contact forces that must remain within the Lorentz cone. These laws lead us to write second order cone constraints on the contact forces. By introducing a non-linear change of variables of the contact relative velocities which defines the so-called modified velocities, the contact laws can be written as a complementarity condition on second order cones between the modified velocities and the contact forces.

1

In solid and structural dynamics, the discrete mechanical model is usually supplemented by equations of motion that relate the velocities of the system to the applied and contact forces. After discretization in space and time, the resulting system can be expressed as nonlinear Second Order Cone Complementarity Problem SOCCP [START_REF] Acary | A formulation of the linear discrete Coulomb friction problem via convex optimization[END_REF][START_REF] Acary | Numerical methods for nonsmooth dynamical systems: applications in mechanics and electronics[END_REF][START_REF] Acary | Applications of an existence result for the Coulomb friction problem[END_REF]. Section 2 details the problem formulation when the discrete equation of motions are assumed to be linear.

To solve these SOCCP problems numerically, a large number of methods have been used in the computational contact mechanics community. Many of these methods are in fact adaptations of well-known mathematical programming methods for solving variational inequalities and complementarity problems. Some of these methods have been jointly developed by optimization specialists. There are two inherent difficulties in frictional contact mechanics problems. The first difficulty is related to the non-monotone nature of the complementarity problem. Identifying an optimization problem for which the complementarity problem would be the optimality conditions is difficult, and leads to non-convex optimization problems. The second difficulty arises from the fact that the constraints are not of full rank and, in the case of rigid multi-body systems, may be of very low rank with respect to the number of constraints. In [START_REF] Acary | On solving contact problems with Coulomb friction: formulations and numerical comparisons[END_REF], a review of the main literature is made and methods are compared with performance profiles. When the second-order constraints are of full rank, second-order methods such as semi-smooth Newton methods, e.g., Alart-Curnier's method [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution method[END_REF], are generally robust and accurate. If not, as it is usual in multi-body systems made of rigid parts (robots, granular material), then the first-order methods such as projected successive over-relaxation gradient methods are robust, but slow and with a limited accuracy in practice. Second order techniques cannot solve the problem due to robustness issues. As far as we know, there is no high accuracy second order method capable of solving the friction contact problems with redundant constraints.

One of the objectives of this work is to propose a second-order method, based on an interior point method that is accurate, robust and efficient for problems where the constraints are rank deficient, but, in the first step, on a convex optimization problem. Some applications of interior point methods for contact problems have already been attempted in the literature [START_REF] Kleinert | An inexact interior point method for the large-scale simulation of granular material[END_REF][START_REF] Kučera | An interior-point algorithm for the minimization arising from 3D contact problems with friction[END_REF][START_REF] Mangoni | A primal-dual predictor-corrector interior point method for non-smooth contact dynamics[END_REF][START_REF] Melanz | A comparison of numerical methods for solving multibody dynamics problems with frictional contact modeled via differential variational inequalities[END_REF]. In the precursor work of [START_REF] Kučera | An interior-point algorithm for the minimization arising from 3D contact problems with friction[END_REF], the contact problem with Tresca friction (purely quadratic cylindrical constraints) is solved with Mehrotra's algorithm. The work in [START_REF] Kleinert | An inexact interior point method for the large-scale simulation of granular material[END_REF] is the most advanced for the case of conic constraints. The problem considered is the relaxed convex problem as in this work, which is further regularized by adding a constant diagonal term to the Jacobian matrices of the interior point algorithm. To solve large problems, linear systems are solved by iterative methods. The degeneracy of the conditioning as the iterations progress leads to the use of preconditioners and makes it costly to obtain high accuracy solutions. In [START_REF] Melanz | A comparison of numerical methods for solving multibody dynamics problems with frictional contact modeled via differential variational inequalities[END_REF], a comparison is made between an interior point method and first order methods. The interior point method on second-order cone shows a better convergence rate but without reaching higher accuracies than first order methods.

In this article, we consider a relaxed problem where the actual local velocity is complementary to the reaction force at each contact point. This yields a reformulation of the original problem as a convex second order cone optimization (SOCO) problem. A dedicated interior point method based on the primal-dual algorithm of Mehrotra with the Nesterov-Todd scaling strategy is tailored to solve this problem. In particular, we show on a large bench set of examples that a tight accuracy can be achieved at optimum without regularizing the Jacobian matrices, provided that their conditioning is controlled.

Since our main objective is to have a robust and efficient computer code capable of solving our problem with a high degree of accuracy, we have first performed experiments with existing solvers. At first, we have tried to solve the convex relaxation problem by using a differentiable nonlinear optimization as proposed by [START_REF] Benson | Solving problems with semidefinite and related constraints using interior-point methods for nonlinear programming[END_REF]. They suggest four ways of reformulating the second-order cone constraint to get round the non-differentiability of the norm at zero. Our experiments were performed by means of the modeling langage AMPL [START_REF] Gay | The AMPL modeling language: an aid to formulating and solving optimization problems[END_REF] and the optimization solver KNITRO [START_REF] Byrd | KNITRO: An integrated package for nonlinear optimization[END_REF]. Our conclusion is that, except in a few rare cases, none of these reformulations is able to solve efficiently and accurately the convex model of our problem. The second experiments were with the interior point solver SDPT3 [START_REF] Tütüncü | Solving semidefinite-quadratic-linear programs using SDPT3[END_REF]. This is an implementation of the Mehrotra predictor-corrector algorithm, which we hoped would be effective and robust in solving our problem. From our point of view, the great advantage of SDPT3 is that it is open source and we can therefore hope to modify it to solve the initial problem. Unfortunately, the results of the experiments with SDPT3 were not conclusive. In particular, the accuracy we wanted to achieve was not sufficient and too many failures were observed. But the experiments were interesting for us, because we identified that the main difficulty in implementing this algorithm is in solving the linear system at each iteration. This led us to carry out our own implementation of this algorithm, devoting our efforts in solving the linear system and in the scaling method. We also made experiments with the interior-point solver MOSEK [START_REF] Andersen | The Mosek interior point optimizer for linear programming: an implementation of the homogeneous algorithm[END_REF] whose algorithm is based on a homogeneous and self-dual model [START_REF] Andersen | On implementing a primal-dual interior-point method for conic quadratic optimization[END_REF]. The numerical results of MOSEK are better than those of SDPT3. However, in order to obtain accurate solutions, the number of iterations can become large or failures occur as soon as the stopping tolerance is reduced (below 10 -7 ). Unfortunately, MOSEK is not open source, so it has not been possible to understand the origin of these failures, nor will it be possible to modify this code to solve our initial problem. These first numerical experiments motivated us to make our own implementation of the interior point algorithm. Especially since in the conclusion of [START_REF] Tütüncü | Solving semidefinite-quadratic-linear programs using SDPT3[END_REF] the authors state that "For problems with second-order (quadratic) cone constraints, experiments indicate that there is room for improvement in SDPT3".

The outline and the contribution of the article are as follows. In Section 2, we recall the basics of mechanical models and how the convex SOCO problem is obtained. The corresponding primal-dual pair is formulated in Section 3. General results are given on the non-emptiness and compactness of the solution set under the Slater hypothesis. In Section 4, the properties of the central path are detailed. In particular, under the assumption of strict complementarity, we establish a characterization of the limit point of the central path, the analytic center of the optimal set, which to our knowledge is new in the SOCO context. This property of the algorithm is important from the mechanical point of view, since the selection of dual variables, that corresponds to the reaction forces, is completely controlled by the interior point method. The details of the numerical implementation of the algorithm are given in Section 5. Several alternative equivalent formulations of the linear system to be solved at each iteration are detailed and the comparison of their conditioning over the iterations is illustrated. Our experiments show that the choice of formulation is fundamental for the robustness of the algorithm. In Section 6, the interior point method is extended to the case of rolling friction, where the cone of constraints is no longer a Lorentz cone and is not self-dual.

Notation

Let x and y be two vectors of R n . The Euclidean scalar product of x and y is denoted x ⊤ y. The associated ℓ 2 norm is ∥x∥ = √

x ⊤ x. The perp notation x ⊥ y means that x ⊤ y = 0. Let E be a subset of R n . The interior, the relative interior and the boundary of E are respectively denoted by int(E), ri(E) and bd(E). The dual of E is denoted and defined by

E * = {y ∈ R n : ∀x ∈ E, x ⊤ y ≥ 0}. The indicator function of E is δ E (x) = 0 if x ∈ E and δ E (x) = +∞ if x / ∈ E. The domain of an extended-valued function f : R n → R ∪ {+∞} is domf = {x ∈ R n : f (x) ∈ R}.
The definitions and notations related to the Jordan algebra are described in Appendix A.

Mechanical models

Let us first introduce the original mechanical models that we are interested in. Let d be the dimension of the space of the mechanical system. Two classes of problems are considered. The frictional contact (FC) problems [START_REF] Acary | On solving contact problems with Coulomb friction: formulations and numerical comparisons[END_REF] for which d = 3 and the problems with rolling friction (RF) at contact [START_REF] Acary | Coulomb friction with rolling resistance as a cone complementarity problem[END_REF] for which d = 5. Let n ∈ N be the number of contact points and let m ∈ N be the number of degrees of freedom. The mechanical system is described by means of three vectors: the global velocity v ∈ R m , the velocity u ∈ R nd and the reaction r ∈ R nd . After discretizing the dynamical system in time and space, the general model to be considered is a conic complementarity problem of the form

M v + f = H ⊤ r, Hv + w = u, K * ∋ u + Φ(u) ⊥ r ∈ K, (1)
where the data of the problem are M ∈ R m×m a symmetric and positive-definite matrix, f ∈ R m , H ∈ R nd×m , w ∈ R nd , K is a cone and K * := {u ∈ R nd : u ⊤ r ≥ 0, for all r ∈ K} is its dual cone. The cone is of the form K = n i=1 K i , each K i being a cone whose definition depends on the model. For a FC problem, we have

K i = r i = (r i,N , r ⊤ i,T ) ⊤ ∈ R × R d-1 : µ i r i,N ≥ ∥r i,T ∥ .
This is the second order cone of Coulomb friction at the ith contact point. The scalar r i,N (resp. u i,N ) and the vector r i,T (resp. u i,T ) are the normal and tangential components of the reaction (resp. velocity) vector of the ith contact point and µ ∈ R n is the vector of friction coefficients. The dual cone of K, is given by

K * = n i=1 K * i , where K * i = u i = (u i,N , u ⊤ i,T ) ⊤ ∈ R × R d-1 : u i,N ≥ µ i ∥u i,T ∥ . The function Φ : R nd → R nd is defined by n components Φ i (u) = (µ i ∥u i,T ∥, 0 ⊤ ) ⊤ ∈ R d .
For a RF problem, the cones are of the form

K i = r i = (r i,N , r ⊤ i,T , r ⊤ i,R ) ⊤ ∈ R × R d-1 × R d-1 : µ i r i,N ≥ ∥r i,T ∥, µ R,i r i,N ≥ ∥r i,R ∥ .
The vector r i,R is the rolling friction reaction moment at contact and µ R ∈ R n is the vector of rolling friction coefficients. The dual cone of K i is

K * i = u i = (u i,N , u ⊤ i,T , u ⊤ i,R ) ⊤ ∈ R × R d-1 × R d-1 : u i,N ≥ µ i ∥u i,T ∥ + µ R,i ∥u i,R ∥ .
The vector u i,R is the relative angular velocity at contact. For this model, the components of the function Φ are defined by Φ

i (u) = (µ i ∥u i,T ∥ + µ R,i ∥u i,R ∥, 0 ⊤ ) ⊤ ∈ R d .
By observing that Φ(u) = Φ(u + Φ(u)) and by making the change of variable u ← u + Φ(u), the problem (1) is reformulated as

M v + f = H ⊤ r, Hv + w + Φ(u) = u, K * ∋ u ⊥ r ∈ K.
(2)

In [START_REF] Acary | A formulation of the linear discrete Coulomb friction problem via convex optimization[END_REF], for the case of a FC problem, an iterative solution of ( 2) is proposed as follows. The second equation is rewritten as

Hv + w + ϕ(s) = u and s i = µ i ∥u i,T ∥, i ∈ {1, . . . , n},
where ϕ : R n → R nd is defined by n components of the form (s i , 0 ⊤ ) ⊤ ∈ R d . The idea is to solve the nonlinear system as a parametric system of the form

M v + f = H ⊤ r, Hv + w + ϕ(s) = u, K * ∋ u ⊥ r ∈ K, ( 3
)
where s is periodically updated according to

s i = µ i ∥u i,T ∥, i ∈ {1, . . . , n}.
There is no guarantee of global convergence of this procedure, but the advantage of this formulation is that, for a fixed s, the parametric system corresponds to the first order optimality conditions of the following convex optimization problem: min

v,u 1 2 v ⊤ M v + f ⊤ v s.t. Hv + w + ϕ(s) = u, u ∈ K * .
This problem belongs to the class of second order cone optimization (SOCO) problems. Our study focuses on the numerical solution of this problem. Obviously, the same parametric strategy can be applied to the solution of an RF problem.

Basic properties

To simplify the models, let us remove the friction coefficients by means of a simple change of variables. Let us define the diagonal matrix P µ whose ith diagonal block is diag(1, µ i , µ i ) for the FC model and diag(1, µ i , µ i , µ R,i , µ R,i ) for the RF model. The parameter s is assumed to be fixed. Let us make the change of variables:

u ← P µ u, H ← P µ H, w ← P µ (w + ϕ(s)) and r ← P -1 µ r.

We keep the same names for the variables, but change the notation of the components. For the RF model we define

u i =   u i,0 ūi ũi   =   u i,N µ i u i,T µ R,i u i,R   and r i =   r i,0 ri ri   =   r i,N 1 µi r i,T 1 µR,i r i,R   .
The same notation is used for the FC model, except that the components ũi and ri do not appear. The system (3) becomes

M v + f = H ⊤ r, Hv + w = u, F * ∋ u ⊥ r ∈ F, (4)
where F is the product of cones. The definition of F depends on the model. For the FC model (d = 3) we have

F = n i=1 L i := {r i = (r i0 , ri ) ∈ R d : r i0 ≥ ∥r i ∥}. ( 5
)
This cone is self-dual (i.e., L * i = L i ), therefore

F * = n i=1 L * i := {u i = (u i0 , ūi ) ∈ R d : u i0 ≥ ∥ū i ∥}. (6)
For the RF model (d = 5), the cone of constraints is

F = n i=1 R i := {r i = (r i0 , ri , ri ) ∈ R d : r i0 ≥ max{∥r i ∥, ∥r i ∥}}. ( 7
)
This cone is not self-dual. The dual cone is then

F * = n i=1 R * i := {u i = (u i0 , ūi , ũi ) ∈ R d : u i0 ≥ ∥ū i ∥ + ∥ũ i ∥}. (8)
The system (4) can be viewed as the first order optimality conditions of the following second order cone optimization problem: min

v,u 1 2 v ⊤ M v + f ⊤ v s.t. Hv + w = u, u ∈ F * . (9)
Let r ∈ R nd be the dual variable associated to the linear constraint of (9). The Lagrangian function associated to problem (9) is defined by

L(v, u, r) = 1 2 v ⊤ M v + f ⊤ v + r ⊤ (u -Hv -w).
Since L is separable in v and u, convex in v and linear in u, and since F * * = F, the dual function is readily obtained and defined by

inf (v,u)∈R m ×F * L(v, u, r) = -1 2 r ⊤ W r -q ⊤ r -1 2 f ⊤ M -1 f if r ∈ F, -∞ otherwise,
where

W = HM -1 H ⊤ and q = w -HM -1 f. ( 10
)
The dual problem is then min

r 1 2 r ⊤ W r + q ⊤ r s.t. r ∈ F. ( 11
)
Let us denote by R the optimal set of this problem. It is a convex set, characterized by the first order optimality conditions:

F * ∋ W r + q ⊥ r ∈ F. (12)
Note that these conditions are equivalent to (4), thanks to the definitions (10) and the equalities

v = M -1 (H ⊤ r -f ) and u = W r + q. ( 13
)
By weak duality, the duality gap (i.e., the difference between the value of the primal and the dual function) is nonnegative. Indeed, let (v, u, r) be a primal-dual feasible solution, we have

1 2 v ⊤ M v + f ⊤ v -(-1 2 r ⊤ W r -q ⊤ r -1 2 f ⊤ M -1 f ) = u ⊤ r
≥ 0, the inequality arising from the fact that r ∈ F and u ∈ F * . The strong duality and the existence of an optimal solution for (11), we need a constraints qualification assumption, known as the Slater hypothesis.

Assumption 3.1. There exists v ∈ R m such that Hv + w ∈ int(F * ).
With respect to the reduced problem ( 11), the Slater hypothesis can be equivalently formulated as follows:

for all nonzero d ∈ F, if W d = 0 then q ⊤ d > 0.
The following result summarizes the connection between the solutions of ( 4) and ( 11). They are well known and are derived from basic properties of convex optimization. We provide a proof in Appendix B to be complete. Proposition 3.2. Consider the primal-dual pair of SOCO problems (9)-( 11), where M ∈ R m×m is symmetric and positive-definite, f ∈ R m , H ∈ R nd×m , w ∈ R nd and F is the product of second order cones of the form (5) or (7).

(i) Any solution of (4), provides a primal-dual optimal solution of (9)-( 11).

(ii) If the problem (9) is feasible, then it has a unique optimal solution (v, û). (iii) Assumption 3.1 is satisfied if and only if the optimal set of (11) is non-empty and compact.

Central path related to a frictional contact problem

In this section, we are interested in FC problems. We recall some known results about the convergence of the central path in SOCO and put them in our context. For the sake of completeness we also provide the proofs.

For FC problems, the friction contact cone is equal to the product of n Lorentz cones, i.e., F = L n , where L is the Lorentz cone of dimension d. We recall that this cone is self-dual, i.e., F * = F. There is an algebra, the so-called Euclidean Jordan algebra, associated to symmetric cones, which allows an almost direct extension of the interior-point algorithms for linear optimization to the case of SOCO (Alizadeh and Goldfarb, 2003). A summary of Euclidean Jordan algebra is given in Appendix A. Thus, the orthogonality condition in (12) becomes u • r = 0, see (Alizadeh and Goldfarb, 2003, Lemma 15). This leads to a square system of equations with the additional conical constraints. Under Assumption 3.1, (r, u) is a primal-dual optimal solution of the problem (11) if and only if

W r + q = u, r • u = 0, (r, u) ∈ L 2n . ( 14
)
It is important to keep in mind that the matrix W is not an arbitrary positive semidefinite matrix, but has a structure given by ( 10). This structure will be useful for solving the linear system done at each iteration. In interior-point methods, a perturbation is introduced in the complementarity equation, so that the system to be solved becomes

W r + q = u, u • r = 2µe, (u, r) ∈ int(L 2n ), ( 15
)
where e is the unit vector associated with the Jordan product and µ > 0 is a parameter that is driven progressively to zero along the iterations, so that in the end a solution of ( 14) is found. The parameter µ is called barrier parameter, because (15) can be interpreted as the optimality condition of the barrier problem min

r φ µ (r) := 1 2 r ⊤ W r + q ⊤ r -µ n k=1 log det r k . ( 16
)
The first order optimality condition of ( 16) is W r + q -2µr -1 = 0. By introducing the variable u = 2µr -1 , we get (15). This system (15) defines a curve, called the central path. The following result states that under Slater's hypothesis, the central path is well defined and remains bounded for bounded values of µ. The proof is given in Appendix C

Proposition 4.1. Under Assumption 3.1, for all µ > 0, the perturbed KKT system (15) has a unique solution (r(µ), u(µ)) ∈ int(L 2n ). For all μ > 0, the set {(r(µ), u(µ)) : 0 ≤ µ ≤ μ} is bounded.

From Proposition 3.2-(ii), the optimal solution of ( 9) is unique, therefore lim µ→0 u(µ) = û. For the curve r(•), it can be shown that lim µ→0 r(µ) = r, where r ∈ ri( R). Such a solution is called a maximally complementary optimal solution of (11). It can be characterized as follows. Regarding to the problem (11), the index set {1, . . . , n} of the Lorentz cones is partitioned into six sets (Mohammad-Nezhad and Terlaky, 2021):

B = {i : ∃r ∈ R, r i ∈ int(L)}, N = {i : ûi ∈ int(L)}, R = {i : ûi ∈ bd(L)\{0} and ∃r ∈ R, r i ∈ bd(L)\{0}}, T 0 = {i ∈ T : ∀r ∈ R, ûi = r i = 0}, T 1 = {i ∈ T : ûi = 0 and ∃r ∈ R, r i ∈ bd(L)\{0}}, T 2 = {i ∈ T : ûi ∈ bd(L)\{0} and ∀r ∈ R, r i = 0}.
An optimal solution r ∈ R is maximally complementary if and only if for all i ∈ B, r i ∈ int(L) and for all

i ∈ R ∪ T 1 , r i ∈ bd(L)\{0}.
The next result shows the convergence of the central path to a maximally complementary solution. A proof is given in Appendix C for self-completeness.

Proposition 4.2. Under Assumption 3.1, the central path (r(•), u(•)) converges to (r, û), where r is a maximally complementary optimal solution of (11).

In linear optimization, it is known that the central path converges to the analytic center of the optimal set [START_REF] Monteiro | On the existence and convergence of the central path for convex programming and some duality results[END_REF]. In semidefinite optimization, it is also known that this result holds under strict complementarity (De [START_REF] De Klerk | Infeasible-start semidefinite programming algorithms via selfdual embeddings[END_REF][START_REF] Halická | On the convergence of the central path in semidefinite optimization[END_REF]. But in SOCO, as far as we know, we have not found in the literature any characterization of the analytic center, even under the strict complementarity assumption. In order to complete the theory, we provide such a characterization.

Assumption 4.3. There exists r ∈ R such that û + r ∈ int(L n ).
This assumption is equivalent to T i = ∅, for i = 1, 2, 3. In that case, (B, N, R) is a partition of {1, . . . , n}. Under this assumption, we can define the analytic center of the optimal set R. If R is a singleton, then the analytic center of R is this point. Otherwise, it is the unique optimal solution of the problem

min r∈ri( R) ψ(r) := - i∈B log det r i - i∈R log r i,0 . (17)
The next proposition shows that the analytic center is well defined.

Proposition 4.4. Under Assumptions 3.1 and 4.3, if R is not reduced to a singleton, then the problem (17) has a unique optimal solution r ∈ ri( R), characterized by the following property: for all r ∈ ri( R),

i∈B r ⊤ i r-1 i + 1 2 i∈R r i,0 ri,0 ≤ |B| + 1 2 |R|. ( 18
)
Proof. Let r ∈ ri( R). For all i ∈ B, det r i > 0, and for all i ∈ R,

r i,0 = ∥r i ∥ > 0. Since ri( R) ̸ = ∅, one has domψ ∩ ri( R) ̸ = ∅.
To show that problem (17) has at least one solution, it suffices to show that the function ψ + δ ri( R) is coercive. This last property is a direct consequence of the compactness of the set R.

The uniqueness of the minimum comes from the strict convexity of ψ on ri( R). Indeed, for a conic component i ∈ B, we have ∇ 2 r (-log det r) r=ri = 2Q r -1 i , which is positive definite for all r i ∈ int(L), see Appendix A. For all conic components i ∈ R, there exists h i ∈ R d-1 , with ∥h i ∥ = 1, such that for all r ∈ ri( R), we have r i = r i,0 (1, h ⊤ i ) ⊤ . Indeed, suppose that for i ∈ R, there exist r and r ′ in ri( R), such that r i and r ′ i are not collinear. By the triangle inequality, it is easy to see that 1 2 (r i + r ′ i ) ∈ int(L), and thus i ∈ B, which would contradict i ∈ R. Finally, for all r, r ′ ∈ ri( R), if r ̸ = r ′ then there exists i ∈ B such that r i ̸ = r ′ i or there exists i ∈ R such that r i,0 ̸ = r ′ i,0 . In both cases, we have ψ(r) ̸ = ψ(r ′ ), which implies that ψ is strictly convex on ri( R). Since ψ is convex on ri( R), r is optimal if and only if -∇ψ(r) is in the normal cone to R at r, i.e., for all r ∈ ri( R), ∇ψ(r) ⊤ (r -r) ≥ 0. For r ∈ ri( R), we have

∇ψ(r) ⊤ (r -r) = -2 i∈B (r -1 i ) ⊤ (r i -ri ) - i∈R r i,0 -ri,0 ri,0 = -2 i∈B r ⊤ i r-1 i + 2|B| - i∈R r i,0 ri,0 + |R|,
from which we deduce that ( 18) is satisfied.

Let us state and prove the main result of this section.

Theorem 4.5. Under Assumptions 3.1 and 4.3, the central path r(•) converges to the analytic center of R.

Proof. Suppose that R is not reduced to a singleton, otherwise the result is a direct consequence of Proposition 4.2. Let r ∈ ri( R). As in the proof of Proposition 4.2, for all µ > 0, (44) is satisfied. By using r(µ) • u(µ) = 2µe and the definition of the partition (B, N, R), we have i∈B

r ⊤ i r -1 i (µ) + i∈N û⊤ i u -1 i (µ) + i∈R (r ⊤ i r -1 i (µ) + û⊤ i u -1 i (µ)) ≤ n. ( 19
)
For any index i, by using the Cauchy-Schwarz inequality and the fact that r i,0 ≥ ∥r i ∥, we have

r ⊤ i r -1 i (µ) = r ⊤ i Rr i (µ) det r i (µ) = r i,0 r i,0 (µ) -r⊤ i ri (µ) r 2 i,0 (µ) -∥r i (µ)∥ 2 ≥ r i,0 r i,0 (µ) -∥r i ∥∥r i (µ)∥ (r i,0 (µ) -∥r i (µ)∥)(r i,0 (µ) + ∥r i (µ)∥) ≥ r i,0 r i,0 (µ) + ∥r i (µ)∥ . ( 20
)
In the same manner, for all i ∈ {1, . . . , n} we have

û⊤ i u -1 i (µ) ≥ ûi,0 u i,0 (µ) + ∥ū i (µ)∥ . (21)
From ( 19), ( 20) and ( 21), we deduce that i∈B

r ⊤ i r -1 i (µ) + i∈N û⊤ i u -1 i (µ) + i∈R ( r i,0 r i,0 (µ) + ∥r i (µ)∥ + ûi,0 u i,0 (µ) + ∥ū i (µ)∥ ) ≤ n.
When µ tends to zero, the first sum tends to i∈B r ⊤ i r-1

i and the second one to |N|. By using the fact that ri,0 = ∥ ri ∥ and ûi,0 = ∥ ūi ∥ for i ∈ R, each term of the third sum tends to ri,0 2ri,0 + 1 2 . By taking the limit as µ → 0 and using n

= |B| + |N| + |R|, we obtain i∈B r ⊤ i r-1 i + 1 2 i∈R r i,0 ri,0 ≤ |B| + 1 2 |R|,
which shows by Proposition 4.4 that r is the optimal solution of ( 17), the analytic center of R.

We would also like to emphasise an important point concerning the characterization of the analytic center from the point of view of computational mechanics. If the matrix H is not of full row rank, the reaction forces r are not unique. The interior point algorithm allows the optimal reaction forces to be uniquely defined. In addition, the analytic center is of mechanical interest and interpretation. By choosing a solution that is furthest from the boundaries of the solution set and therefore from the boundaries of the cones, the chosen solution maximizes the sum of the normal stresses r 0 and the distance to the edges of the cones. For simple cases, this provides a better distribution of reactions for active constraints rather than particular solutions where reaction forces are concentrated on a few constraints.

Numerical solution of the friction contact problems

We have implemented the primal-dual interior point algorithm developed by Tütüncü, Toh and Todd [START_REF] Tütüncü | Solving semidefinite-quadratic-linear programs using SDPT3[END_REF] and adapted it to our context. The algorithm is an extension of the predictor-corrector algorithm of Mehrotra [START_REF] Mehrotra | On the implementation of a primal-dual interior point method[END_REF] to the solution of a SOCO problem. The differences with the algorithm implemented in SDPT3 are as follows:

• The convex quadratic objective function is taken into account directly, without reformulating the problem with a linear objective function and a new conic constraint. One consequence is that the primal and dual steps must be equal, which is not necessarily the case with a linear objective function. This is explained below.

• The linear system to solve at each iteration, resulting from the linearization of the perturbed KKT system (15), is transformed via the Nesterov-Todd scaling strategy as in SDPT3, but in a different manner.

The scaling matrix (i.e., the quadratic representation of the vector p defined by formula (25) below), which introduces bad conditioning when the iterations are close to an optimal solution, is never explicitly computed and stored in memory. Only matrix-vector products are performed. This method allows us to achieve a high degree of accuracy in the computation of the directions and thus in the solution of the quadratic convex problem. We examine in detail equivalent formulations of the linear system to justify our choice and detail all the computations performed.

• We also exploit the data structure of the problems, by using the matrices H and M instead of the reduced matrix W (see formula (10)), in order to preserve the sparsity of the linear system. This transforms the initial 2 × 2 block linear system into a 3 × 3 block system. We also show, by numerical experiments, that it is better to reduce again this system back to a 2 × 2 block system, while preserving the sparse structure, but avoiding some of the calculations with the scaling matrix.

• Finally, the experiments show that using the C data type long double for all the calculations related to the Nesterov-Todd scaling scheme, improves the robustness of the implementation, while maintaining reasonable computation times.

The main part of the algorithm is the solution of two linear systems that result from the linearization of the equation ( 15) at the current iterate (u, r) ∈ int(L 2n ). They only differ on the right-hand side and are of the form:

W -I U R d r d u = -W r -q + u -u • r -[d u a • d r a -2σµe] ( 22
)
where U = Arw(u) and R = Arw(r) are the arrow-shaped matrix associated to u and r. See Appendix A for all definitions related to the Euclidean Jordan algebra that are used in this paper. The first direction, denoted (d u a , d r a ) and called the affine scaling direction, is the solution of ( 22) without the square bracketed term in the right-hand side. It then satisfies the linear equation

u • d r a + r • d u a = -u • r.
The affine scaling direction is a Newton step on the original optimality system (14). The barrier parameter is set to µ = u ⊤ r n . The centralization parameter σ ∈ (0, 1] is fixed by comparing the current value of µ with its expected reduction obtained along the affine step. The second direction is a linear combination of the affine scaling direction and of a corrector step, to keep the iterates near the central path. It then satisfies the following linear equation u

• d r + r • d u = -u • r -d u a • d r a + 2σ u ⊤ r n e. The next iterate is set to (u + , r + ) = (u, r) + α(d u , d r )
, where α is the largest value in (0, 1] such that

(u + , r + ) ∈ (1 -τ )(u, r) + L 2n ,
for some value τ ∈ (0, 1] (typically τ = 0.995). Contrary to common practice, different primal and dual steplengths are not taken, because of the first equation in ( 22), which is linear and includes both primal and dual variables. Indeed, suppose that the current iterate is such that W r+q-u = 0 and that different steplengths α ̸ = α ′ are taken with d u ̸ = 0. We then have

W (r + αd r ) + q -(u + α ′ d u ) = (α -α ′ )d u .
If d u ̸ = 0, then at the next iteration the residual of the linear equation is no longer zero.

The algorithm will be well defined if the matrix in ( 22) is nonsingular at each iteration. Its determinant is equal to det(W + R -1 U ) det R. Although the vectors r and u are kept inside the interior of the second order cones, and thus R and U are positive definite, the matrix W +R -1 U can be singular. This is because the matrix R -1 U is not necessarily symmetric, since in general r and u do not commute. The following example is given by (Peng et al., 2002, p.143

): If U = Arw([1, 0.7, 0.7] ⊤ ), R = Arw([1, 0.8, 0.5] ⊤ ) and W = diag([0.3, 1, 0] ⊤ ), then det(W + R -1 U ) = 0.
In addition to this singularity problem, there is the symmetry problem. In the context of interior-point algorithms for linear or nonlinear optimization, where the cone of constraints is the nonnegative orthant, the matrices U and R are diagonal and so the matrix of ( 22) can be symmetrized, for example by left-multiplying the second row by -U -1 . See, e.g., [START_REF] Ghannad | Linear systems arising in interior methods for convex optimization: a symmetric formulation with bounded condition number[END_REF] for several symmetrization techniques in interior point methods. The main advantages of a symmetric system are lower factorization costs and an effective control of the inertia of the factorized matrix. In addition, very efficient codes such as MA57 [START_REF] Duff | MA57-a code for the solution of sparse symmetric definite and indefinite systems[END_REF] or MUMPS [START_REF] Amestoy | A Fully Asynchronous Multifrontal Solver Using Distributed Dynamic Scheduling[END_REF] can be used for this task.

To overcome these problems of singularity and symmetry, a change of variables, called a scaling scheme, is used to obtain a symmetric nonsingular system. The idea is to make a change of variables that leaves the Lorentz cone invariant and such that in the new space the vectors u and r commute. However, this change of variable depends on the current iterate and must be done at each iteration. Let p ∈ int K and let Q p be the associated quadratic representation (see Appendix A). Since p is in the interior of the cone K, the matrix Q p is positive definite. From (Alizadeh and Goldfarb, 2003, Theorem 9), we have

Q p (L n ) = L n and Q p (int L n ) = int L n . Let us consider the change of variables ř = Q p -1 r and û = Q p u.
The problem (11) becomes min

ř 1 2 ř⊤ Q p W Q p ř + (Q p q) ⊤ ř s.t. ř ∈ L n . (23)
The corresponding perturbed KKT system is W r + q = u and û • ř = 2µe, with (û, ř) ∈ int L 2n . The linearization of these equations leads the following linear system:

W -I Û Q p -1 ŘQ p d r d u = -W r -q + u -û • ř -[ du a • ďr a -2σµe] , (24)
where Û = Arw(û) and Ř = Arw(ř). The choice of the vector p ∈ int K is made so that û and ř commute, which implies that the matrix of the linear system (24) is nonsingular. Indeed, this matrix is nonsingular if and only if det(W + ( ŘQ p ) -1 Û Q p -1 ) ̸ = 0. Since Û and Ř are positive definite, û and ř commute, and

Q p -1 = Q -1 p , we have ( ŘQ p ) -1 Û Q p -1 = Q p -1 Ř-1/2 Û 1/2 Û 1/2 Ř-1/2 Q p -1 ,
which is symmetric and positive definite.

Several choices for the vector p are possible, see [START_REF] Hal. Alizadeh | Second-order cone programming[END_REF]. As mentioned in [START_REF] Tütüncü | Solving semidefinite-quadratic-linear programs using SDPT3[END_REF], the most efficient scaling technique for the solution of a SOCO problem, is the one using the Nesterov-Todd (NT) direction [START_REF] Nesterov | Self-scaled barriers and interior-point methods for convex programming[END_REF]:

p = Q u 1/2 (Q u 1/2 r) -1/2 -1/2 = Q r -1/2 (Q r 1/2 u) 1/2 -1/2 . ( 25
)
The main property of the NT direction is that û = ř, which implies that

Û -1 Ř = I.
The symmetrization of the system ( 24) is done by left-multiplying the last row by -Q p Û -1 , leading to a symmetric matrix of the form

W -I -I -Q p 2 .
To take advantage of the sparsity of the matrices M and H, the system ( 22) is considered in the following equivalent augmented form:

  M -H ⊤ 0 -H 0 I 0 U R     d v d r d u   =   -M v -f + H ⊤ r Hv + w -u -u • r -[d u a • d r a -2σµe]   . ( 26
)
The system (26) can be interpreted as the linearization of the perturbed KKT conditions of the problem (9). Applying the scaling scheme, the system becomes

  M -H ⊤ 0 -H 0 I 0 I Q p 2     d v d r d u   =   -M v -f + H ⊤ r Hv + w -u -r -[Q p (û -1 • ( du a • ďr a )) -2σµu -1 ]   . (27)
A reduction of ( 27) can be done by eliminating the variable d u , while keeping the sparse structure. This leads to the reduced symmetric system

M -H ⊤ -H -Q p -2 d v d r = -M v -f + H ⊤ r Hv + w + [Q p -1 (û -1 • ( du a • ďr a )) -2σµr -1 ] . ( 28
)
The big flaw of the scaling strategy is the ill-conditioning of the matrix Q p 2 when the solution pair (u, r) approaches an optimal solution. Indeed, suppose that (u * , r * ) is a primal-dual optimal solution of (11), which satisfies the strict complementarity condition. Let (u, r) be an interior point iterate near the optimal solution and let p be defined by (25). For i ∈ {1, . . . , n}, three situations can occur [START_REF] Cai | Solving second order cone programming via a reduced augmented system approach[END_REF]:

• u * i ∈ int(L) and r * i = 0, then all the eigenvalues of Q p 2 i are of order µ := u ⊤ r; • u * i = 0 and r * i ∈ int(L), then all the eigenvalues of Q p 2 i are of order 1/µ; • u * i ∈ bd(L), r * i ∈ bd(L) and (u * i , r * i ) ̸ = (0, 0), then the largest eigenvalue of Q p 2 i is of order 1/µ and the smallest is of order µ.

To overcome the difficulties due to ill-conditioning, we propose to solve the linear systems ( 27) and ( 28) under the following equivalent form:

  M -Ĥ⊤ 0 -Ĥ 0 I 0 I I     d v ďr du   =   -M v -f + H ⊤ r Ĥv + ŵ - û -ř -[û -1 • ( du a • ďr a ) -2σµû -1 ]   (29) and M -Ĥ⊤ -Ĥ -I d v ďr = -M v -f + H ⊤ r Ĥv + ŵ + [û -1 • ( du a • ďr a ) -2σµû -1 ] , (30)
where In our numerical experiments, we also consider the reduced system

Ĥ = Q p H, ŵ = Q p w, du = Q p d u and ďr = Q p -1 d r .
( ĤM -1 Ĥ⊤ + I) ďr = -ĤM -1 (f + H ⊤ r) -ŵ -[û -1 • ( du a • ďr a ) -2σµû -1 ], (31)
for which the matrix is positive definite.

Figure 1 shows the behavior of the condition number of the matrices of the six linear systems ( 26)-( 31) along the iterations of the interior-point algorithm for two examples. The first example (left figure) has a single contact point:

n = 1, m = 3, M = I, w = 0, f = [3, 3, 3, 1, -1, -3, 1, -1, -3] ⊤ , H = [D, 0, -D] where D = diag(1, 0.1, 0.1)
. Since H is of full rank, the primal-dual solution is unique, the optimal reaction and relative velocity vectors are non-zero and on the boundary of the Lorentz cone.

The second example is a "Box Stacks" problem from FCLIB, a collection of discrete three-dimensional frictional contact problems [START_REF] Acary | FCLIB: a collection of discrete 3D Frictional Contact problems[END_REF] 2 . There are n = 69 contact points and m = 450 degrees of freedom, H ∈ R 207×450 and rank(H) = 157. The optimal solution satisfies the strict complementarity condition and (|B|, |N|, |R|) = (18,5,46). The matrix of (26) at the endpoint of the minimization procedure, is nearly rank deficient, there are 15 singular values less than √ ϵ, where ϵ is the machine epsilon. These two examples show that the matrices in ( 29) and ( 30) remain the least ill-conditioned.

An advantage of the systems ( 29) and ( 30) is that the quadratic representation matrices are never explicitly built in memory for the entire computation. Only the product of these matrices times a vector needs to be performed. Indeed, for a pair of vectors (x, y) ∈ L 2 , the product of a vector y by the quadratic representation of x can be done via the formula Q x y = 2(x ⊤ y)x -(det x)R d y. Therefore, the product of a vector by a matrix Q p , where p is the , can be done by performing only three products of a quadratic representation matrix by a vector. Similarly, the computation of the inverse or the square root of a vector in the Jordan algebra, is done by using the spectral decomposition of that vector. Moreover, even if the number of cones can be large, the computational cost of a spectral decomposition per cone is very small, since the dimension of a Lorentz cone is only three. 32) is satisfied then return (v, u, r) as primal-dual solution of ( 9);

Family

η 1 ∈ (0, 1), η 2 ≥ 1, η 3 ≥ 1, τ 1 ∈ (0, 1), τ 2 ∈ (0, 1 -τ 1 ), tol > 0 1: if the stopping criterion (
2: Set µ ← u ⊤ r/n; 3: Compute (d v a , d r a , d u a
) solution of (26) (resp. ( 27)-( 31)) without the square bracketed term; 4: Find the greatest α a ∈ (0, 1] such that (u, r) 27)-( 31)); 9: Set τ ← τ 1 + α a τ 2 ; 10: Find the greatest α ∈ (0, 1] such that (u, r)

+ α a (d u a , d r a ) ∈ L 2n ; 5: Set µ a ← (u + α a d u a ) ⊤ (r + α a d r a )/n; 6: if µ > η 1 then set e ← max{1, η 2 α 2 a } else set e ← η 3 ; 7: Set σ ← min{1, (µ a /µ) e }; 8: Compute (d v , d r , d u ) solution of (26) (resp. (
+ α(d u , d r ) ∈ (1 -τ )(u, r) + L 2n ; 11: Set (v, u, r) ← (v, u, r) + α(d v , d u , d r ) and goto 1.
Algorithm 1 details one iteration of the implemented algorithm based on that of SDPT3 [START_REF] Tütüncü | Solving semidefinite-quadratic-linear programs using SDPT3[END_REF]. The stopping criterion used in GFC3D is res := max ∥Hv + w -u∥ max{∥Hv∥, ∥w∥, ∥u∥}

, ∥M v + f -H ⊤ r∥ max{∥M v∥, ∥f ∥, ∥H ⊤ r∥} , |u ⊤ r| ≤ tol, ( 32
)
The values of the parameters are fixed to η 1 = 10 -10 , η 2 = 3, η 3 = 1, τ 1 = 0.9, τ 2 = 0.09 and tol = 10 -10 . The starting point is set as follows: for all i ∈ {1, . . . , n}, u i = r i = (0.1, 0.01, 0.01) ⊤ and v = M -1 (H ⊤ r + f ). The experiments were done on 1091 problems of the FCLIB collection [START_REF] Acary | FCLIB: a collection of discrete 3D Frictional Contact problems[END_REF]. There are seven families of problems, whose dimensions are described in Table 1, where n is the number of three-dimensional cones and m is six times the number of bodies. The experiments were done on a MacBook Pro 2.3GHz equipped with a procesor Intel Core2 i9 8 coeurs and 32Gb of memory, running OSX 10.15.7.

The linear system ( 26) is solved by means of a LU factorization, the symmetric ones with a LDL ⊤ factorization with MA57 [START_REF] Duff | MA57-a code for the solution of sparse symmetric definite and indefinite systems[END_REF]. Even for the solution of the positive definite system (31) MA57 is used. Two types of failures are returned during a run:

• Failure 1: the stopping criterion (32) is not satisfied after a maximum of 100 iterations.

• Failure 2: A NaN (Not a Number) is detected during the computation of the new iterate.

Table 2 indicates the number of successes and failures when solving the 1091 problems of the FCLIB collection, with a tolerance fixed to tol = 10 -10 . Each row corresponds to a run of Algorithm 1 with the numerical solution of the indicated linear system. Figure 2 shows the corresponding performance profiles [START_REF] Dolan | Benchmarking optimization software with performance profiles[END_REF]. With the system (26) the failures are due to a nearly singular system. In these cases, either the algorithm stalls to a spurious solution (13 out of 22 cases) or the convergence becomes very slow (9 out of 22 cases). It should be noted, however, that the "no-scaling" strategy yields an optimal solution for almost 98% of the problems. The systems ( 27) and ( 28) return a large number of failures of type 2. This is mainly due to the ill-conditioning of the matrix Q p 2 when approaching an optimal solution. The reduction of the system worsens the results. Surprisingly, the worst results are obtained with the system (29). A deterioration of the residual of the second linear equation in (29) over the iterations is observed, when the matrix Q p becomes increasingly ill-conditioned. This leads to a loss of the primal feasibility of the iterates. This is mainly due to the scaling of the linear equation -Hd v + d u = u -Hv -w. To overcome this problem, the refinement procedure described in the MA57 documentation can be used to solve (29). We performed a run with a refinement tolerance fixed to the tolerance tol and a maximum of 10 refinement iterations. This results in only six type 2 failures and no more failure of type 1, but it takes longer to run than with (30) as shown in Figure 2. The best performance in terms of robustness is obtained with the system (30). No failures were detected. The average number of iterations to solve one of the 1091 problems is equal to 18, while the minimum and maximum numbers are equal to 8 and 34. The positive definite system (31) gives a good performance in terms of efficiency, but its robustness is not sufficient, even when refinement is applied. 2: Performance profiles for seven different linear system choices. For τ ≥ 0, ρ s (τ ) is the fraction of problems for which the performance of a given version of the algorithm is within a factor 2 τ of the best one.

At last, it should be mentioned that all the computations related to scaling are done in C using the floating point datatype long double. This data type is also used in the computation of the steplengths. This always results in a better accuracy, although the computing time is slightly longer. Table 3 shows the results of solving the FCLIB problems with Algorithm 1, the linear system under the form (30), with a stopping tolerance tol = 10 -11 . The comparison is between the use of the datatype long double versus double. We can see that the number of failures of type 2 is more than doubled with the type double. For these runs, even in case of a failure, the residual (left term in (32)) is small, which means that an optimal solution has been found. The column max res indicates the maximum value of the residual norm (32) of all the 1091 residuals and shows that the type long double allows to obtain a better precision. Even with the data type double, all problems are successfully solved for tol = 10 -10 . The last column of this table shows the total computational time needed to solve all the FCLIB problems with this tolerance. It can be seen that the increase in computational time is less than 10% with long double. Table 3: Performance comparison double versus long double when solving the FCLIB problems with system (30) and tol = 10 -11 .

Rolling friction contact problem

We now consider the solution of (4) in the framework of the RF model (d = 5) defined by the rolling friction cones (7) and ( 8). The main difficulty is that an elementary cone R i is no more self-dual and therefore not symmetric. There is no Jordan product such that R i is a cone of squares with respect to this product. A potential approach is to transform the cone of constraints related to problem (9) into the product of Lorentz cones by means of the following usual trick of introducing artificial variables. For real numbers, a ≥ b + c if and only if there exist two real numbers t ≥ b and t ′ ≥ c such that a = t + t ′ . By setting u i,0 = ti + ti for all i ∈ {1, . . . , n}, the primal-dual pair of problems ( 9)-( 11) can be rewritten as

min v,z 1 2 v ⊤ M v + f ⊤ v s.t. Hv + w = Jz, z ∈ L 2n , (33)
and min

r 1 2 r ⊤ W r + q ⊤ r s.t. J ⊤ r ∈ L 2n , (34)
where d+1) is a block diagonal matrix with n blocks of the form

z = (z ⊤ 1 , . . . , z ⊤ n ) ⊤ and z i = ( ti , ū⊤ i , ti , ũ⊤ i ) ⊤ ∈ R d+1 , J ∈ R nd×n(
ȷ =   1 1 I I   ∈ R d×(d+1) .
For all i ∈ {1, . . . , n}, we have ȷz i = ( ti + ti , ū⊤ i , ũ⊤ i ) ⊤ and ȷ ⊤ r i = (r i0 , r⊤ i , r i0 , r⊤ i ) ⊤ . The perturbed KKT system associated with the problem (34) is then

W r + q = Jz, z • J ⊤ r = 2µe, (z, J ⊤ r) ∈ int(L 4n ). (35)
This optimality system is associated with the barrier problem min r φ µ (r) := 1 2 r ⊤ W r + q ⊤ r -µ n i=1 log(det(r i,0 , ri ) det(r i,0 , ri )). ( 36)

As for Proposition 4.1 and by coercivity of the barrier function, under Assumption 3.1, for all µ > 0 the system (35) has a unique solution such that (z(µ), J ⊤ r(µ)) ∈ int(L 4n ). Although the optimal solution of ( 33) is not unique, because J is non-injective, it can be shown, like for Proposition 4.2, that the central path converges to a relative interior point of the primal-dual optimal set. Under the hypothesis of strict complementarity, it can also be shown that the central path r(•) converges to the analytic center of the dual optimal set. For the sake of completeness, we state the result, but without proof in order to lighten the paper.

Assumption 6.1. There exists a solution (z, r) ∈ R n(d+1) × R nd of (35) with µ = 0, such that z + J ⊤ r ∈ int(L 2n ).

The strict complementarity hypothesis can be equivalently reformulated relatively to the original problems ( 9) and ( 11), thanks to the orthogonality condition (12). Let (z, r) be a solution of ( 35) with µ = 0. Then (v, u = Jz) is the unique optimal solution of (9). Assumption 6.1 is satisfied if and only if for all i ∈ {1, . . . , n}, one of the following three assertions is satisfied:

• r i0 > max{∥r i ∥, ∥r i ∥} and u i = 0,

• u i0 > ∥ū i ∥ + ∥ũ i ∥ and r i = 0,

• r i0 = max{∥r i ∥, ∥r i ∥} > 0 and u i0 = ∥ū i ∥ + ∥ũ i ∥ > 0.

As in Section 4, under the strict complementarity assumption, the index set {1, . . . , n} can be partitioned into two partitions ( B, N, R) and ( B, Ñ, R), whose definition is a direct extension of the previous one to the current framework. Let Z and R be the primal and dual optimal sets of problems ( 33) and (34). The first partition is defined as follows, the second in a similar way:

B = {i : ∃r ∈ R, (r i,0 , ri ) ∈ int(L)}, N = {i : ∃z ∈ Z, ( ti , ūi ) ∈ int(L)}, R = {i : ∃(z, r) ∈ Z × R, (( ti , ūi ), (r i,0 , ri )) ∈ bd(L 2 )\{0}}.
We can then define the analytic center of the optimal set R as follows. If R is reduced to a singleton, then it is this point, otherwise it is the minimum of the problem

min r∈ri( R) - i∈ B log det(r i,0 , ri ) - i∈ B log det(r i,0 , ri ) - i∈ R,i∈ R log r i,0 . (37)
The analytic center can be characterized like in Proposition 4.4, by which it can be shown that Theorem 4.5 still holds for the rolling friction framework.

Algorithm 1 is modified in order to solve a RF problem. This is described by Algorithm 2.

Algorithm 2 One iteration of Mehrotra primal-dual algorithm for solving a RF problem Parameters: η 1 ∈ (0, 1), η 2 ≥ 1, η 3 ≥ 1, τ 1 ∈ (0, 1), τ 2 ∈ (0, 1 -τ 1 ), tol > 0 1: if the stopping criterion (32) with u = Jz is satisfied then return (v, z, r) as primal-dual solution of (33);

2: Set µ ← r ⊤ Jz/n; 3: Compute (d v a , d r a , d z a
) solution of (38) (resp. ( 39)-( 40)) without the square bracketed term; 4: Find the greatest α a ∈ (0, 1] such that (z,

J ⊤ r) + α a (d z a , J ⊤ d r a ) ∈ L 4n ; 5: Set µ a ← (r + α a d r a ) ⊤ J(z + α a d z a )/n; 6: if µ > η 1 then set e ← max{1, η 2 α 2
a } else set e ← η 3 ; 7: Set σ ← min{1, (µ a /µ) e }; 8: Compute (d v , d r , d z ) solution of (38) (resp. ( 39)-( 40)); 9: Set τ ← τ 1 + α a τ 2 ; 10: Find the greatest α ∈ (0, 1] such that (z, J ⊤ r) + α(d z , J ⊤ d r ) ∈ (1 -τ )(u, r) + L 4n ; 11: Set (v, z, r) ← (v, z, r) + α(d v , d z , d r ) and goto 1.

The linear system solved at each iteration is obtained by linearizing the system (35). It is reformulated under the form of the following augmented system

  M -H ⊤ 0 -H 0 J 0 ZJ ⊤ R     d v d r d z   =   -M v -f + H ⊤ r Hv + w -Jz -z • J ⊤ r -[d z a • (J ⊤ d r a ) -2σµe]   , (38)
where Z = Arw(z), R = Arw(J ⊤ r) and µ = r ⊤ Jz n . As in Algorithm 1, the affine scaling direction (d v a , d r a , d z a ) is the solution of (38) without the square bracketed term in the right-hand side, while the full step (d v , d r , d z ) is the solution of the complete system.

The scaling strategy is similar to that described in Section 5. The NT direction p is defined by the formula (25) where u and r are respectively replaced by z and J ⊤ r. The change of variables is done by setting ẑ := Q p z and y := Q p -1 J ⊤ r.

Recall that ẑ = y, which allows to symmetrize the linear system under the form Because of the ill-conditioning of the matrix Q p 2 , an equivalent form of (39) has been considered:

  M -H ⊤ 0 -H 0 J 0 J ⊤ Q p 2     d v d r d z   =   -M v -f + H ⊤ r Hv + w -Jz -J ⊤ r -[Q p (ẑ -1 • ( dz a • ďy a )) -2σµz -1 ]   . (39)
  M -H ⊤ 0 -H 0 JQ p -1 0 Q p -1 J ⊤ I     d v d r dz   =   -M v -f + H ⊤ r Hv + w -Jz -y -[ẑ -1 • ( dz a • ďy a -2σµe)]   , (40)
where dz = Q p d z . Figure 3 shows the condition number of the three matrices ( 38)-( 40) along the iterations of the numerical resolution of a RF problem. It can be seen that the system ( 40) is better conditioned than (39). We also tried several reductions to a 2 × 2 form as in ( 28) or (30), leading to matrices of the form

M -H ⊤ -H -JQ p -2 J ⊤ or M Ĥ⊤ -Ĥ -I , (41)
where Ĥ = P -1 H and P P ⊤ = JQ p -2 J ⊤ . We also tried several ways to compute the matrix P , by performing a Cholesky factorization or by directly exploiting the structure of a quadratic representation matrix. Despite such a reduction and in contrast to the results obtained with the RF problems, the numerical performance has not been improved. Moreover, the computation of the matrix P increases the overall computational cost, without any real improvement. We also observe the same for the 1 × 1 system like (31). The numerical tests were carried out on 526 RF problems of the FCLIB family [START_REF] Acary | FCLIB: a collection of discrete 3D Frictional Contact problems[END_REF], whose characteristics are in Table 4. The numerical results and performances of Algorithm 2 with the three linear systems described previously, are reported in Table 5 andFigure 4. These results show that, with a tolerance tol = 10 -10 , the choice of system (40) gives the best performance. But the performance gap between the systems with the NT scaling is smaller than those observed for the FC problems. It can also be observed that, as for the RF problems, without scaling the failures are of type 1, while with NT scaling the failures only occur when a NaN is detected. Moreover, in the latter case the maximum value of the residual (32) (right column of Table 5), shows that the stopping point of the algorithm is nearly optimal.

Conclusion

At the beginning of this work, the solution of the relaxed convex form of the friction contact models seemed very simple. However, after our initial experience with existing software and our first implementation of Mehrotra's

Figure 4: Performance profiles for the three linear systems for the RF model primal-dual algorithm, we were not satisfied with its efficiency and robustness. The Nesterov-Todd scaling strategy is a wonderful theoretical tool, but numerically very painful. As the iterates approach the boundary of the second order cones, the conditioning of the linear system explodes, the iterations get stuck on the boundary and divisions by zero occur, producing NaN and thus an emergency stop. We have therefore examined a large number of equivalent formulations of the linear system and found that the one which gives the best results, or shall we say the least bad, is the one in which the quadratic representation matrix which allows scaling is not a direct component of the matrix of the linear system. In addition, special attention must be paid to the way in which the matrix-vector products are performed in order to construct the system to be factorized. Unfortunately, we have not been able to find a formulation for rolling friction problems that is as efficient and robust as for friction cones. The accuracy we have achieved is slightly lower. Nevertheless, the accuracies and computation times we have achieved for both models seem to us to be quite adequate and can be used for real applications.. The next step in this research work is to extend the primal-dual algorithm to solve the original model (1). A natural idea is to solve a sequence of systems parametrized by µ > 0:

M v + f = H ⊤ r, Hv + w + ϕ(s) = u, s i = ∥ū i ∥, i ∈ {1, . . . n}, u • r = 2µe.
The main issues are that this system is not the optimality system of an optimization problem and that there is a non-smooth equation.

One theoretical question remains from this study. It concerns the characterisation of the limit point of the central path, like the formula (17) which defines the analytic center, but without the assumption of strict complementarity. Wriggers, P. (2006). Computational contact mechanics. Vol. 2. Springer Berlin Heidelberg. [doi], [oa].

A Euclidean Jordan algebra

Let us consider the set K = n i=1 K i where K i is an n i -dimensional Lorentz cone defined by

K i = x = (x i0 , x⊤ i ) ⊤ ∈ R × R ni-1 : ∥x i ∥ ≤ x i0 .
Let N = n i=1 n i . For x ∈ R N , we denote x = (x 1 , . . . , x n ), where x i = (x i0 , xi ). For two vectors x and y in R N , the Jordan product is defined by

x • y =    x 1 • y 1 . . ., x n • y n    , where x i • y i = x ⊤ i y i x i0 ȳi + y i0 xi
, for i = 1, . . . , n.

Let x ∈ R N and x 2 = x • x. A fundamental property of the Jordan algebra for interior-point algorithms, is that the Lorentz cone is the cone of squares, that is K = {x 2 : x ∈ R N }, see (Alizadeh and Goldfarb, 2003, pp. 18-19).

For matrices A and B, we define For i ∈ {1, . . . , n}, let det(x i ) = x 2 i0 -∥x i ∥ 2 be the determinant of x i ∈ R ni . If x i is nonsingular, i.e., det(x i ) ̸ = 0, the inverse of x i is the unique vector of R ni such that x i • x -1 i = e i and is given by x -1 i = det(x i ) -1 R ni x i , where R ni is the reflexion matrix defined by

A ⊕ B := A 0 0 B .
R ni = 1 0 ⊤ 0 -I ∈ R ni×ni .
If for all i ∈ {1, . . . , n} x i is nonsingular, we have

x -1 =    x -1 i . . . x -1 n    = Arw(x) -1 e.
For i ∈ {1, . . . , n}, the spectral decomposition of a vector x i ∈ R ni is defined by x i = λ i c i + λ ′ i c ′ i , where

λ i = x i0 + ∥x i ∥, λ ′ i = x i0 -∥x i ∥, c i = 1 2 1 xi ∥xi∥ and c ′ i = 1 2 1 -xi ∥xi∥ .
The scalars λ i and λ ′ i are eigenvalues of Arw(x i ), with the corresponding eigenvectors c i and c ′ i . This pair of eigenvectors is called the Jordan frame of x i . It is said that two vectors commute if they share the same Jordan frame. In that case the corresponding arrow matrices commute. From these definitions, it follows that x i ∈ K i (resp. x i ∈ int(K i )) if and only if λ ′ i ≥ 0 (resp. λ ′ i > 0). If x i is nonsingular, then

x -1 i = λ -1 i c i + λ ′ -1 i c ′ i . If x i ∈ int(K i ), then x 1/2 i = λ 1/2 i c i + λ ′ 1/2 i c ′ i .
More generally, for a continuous function f we can define

f (x i ) = f (λ i )c i + f (λ ′ i )c ′ i .
At last, the quadratic representation of a vector x i ∈ K i is defined as

Q xi = Arw 2 (x i ) -Arw(x 2 i ) = 2x i x ⊤ i -det(x i )R ni .
The scalars λ 2 i and λ ′ 2 i are eigenvalues of Q xi with the eigenvectors c i and c ′ i . In particular, two vectors commute if and only if their quadratic representation matrices commute. For x ∈ K, we define Q x = Q x1 ⊕ . . . ⊕ Q xn . For i ∈ {1, . . . , n} and x ∈ int(K i ), we have ∇ x (log det(x i )) = 2x -1 i and ∇ 2

x (log det(x i )) = -2Q x -1 i . These operators Arw and Q • are useful for the design of interior algorithms. See (Alizadeh and Goldfarb, 2003, §4) for a review of their interesting properties.

B Proof of Proposition 3.2

Before giving a proof of Proposition 3.2, we propose some equivalent formulations of Assumption 3.1. They are proved thanks to the following well-known lemma (Skarpness and Sposito, 1982, Lemma 2), called Tucker's theorem of the alternative in the case where K is the nonnegative orthant, see (Mangasarian, 1969, p. 29). The proof can be made by applying a separation theorem of convex sets, see, e.g., (Rockafellar, 1970, Theorem 11.3).

Lemma B.1. Let K ⊂ R n be a closed pointed convex cone, A ∈ R m×n and B ∈ R p×n . One and only one of the following statements is true.

(i) There exists a non-zero x ∈ K such that Ax = 0 and Bx ≤ 0.

(ii) There exists (y, z) ∈ R m × R p + such that A ⊤ y + B ⊤ z ∈ int(K * ).

The following lemma provides equivalent formulations of Assumption 3.1.

Lemma B.2. Let F be the product of second order cones of the form (5) or (7). Let M ∈ R m×m be symmetric and positive-definite, f ∈ R m , H ∈ R nd×m , w ∈ R nd and W = HM -1 H ⊤ . The following four assertions are equivalent.

(i) There exists v ∈ R m such that Hv + w ∈ int(F * ).

(ii) There exists (v, t) ∈ R m × R + such that Hv + tw ∈ int(F * ).

(iii) There does not exist a nonzero d ∈ F, such that H ⊤ d = 0 and w ⊤ d ≤ 0.

(iv) There does not exist a nonzero d ∈ F, such that W d = 0 and q ⊤ d ≤ 0.

Proof. The implication (i) ⇒ (ii) is obvious. The equivalence (ii) ⇔ (iii) is a direct consequence of the fact that F * is a closed pointed convex cone and of Lemma B.1. The equivalence (iii) ⇔ (iv) follows from the definitions of W and q in (10), and of the positive definiteness of M . It remains to prove that (ii) ⇒ (i). Let (v, t) ∈ R m × R + such that Hv + tw ∈ int(F * ). If t > 0, set v ′ = 1 t v. Since Hv + tw = t(Hv ′ + w) and int(F * ) is a cone, Hv ′ + w ∈ int(F * ). Suppose now that t = 0. There exists ϵ > 0 such that Hv + B ϵ ⊂ F * , where B ϵ is the open ball centered at 0 with radius ϵ. Since F * is a cone, for all t > 0, tHv + B tϵ ⊂ F * . Let us choose t > 1 ϵ ∥w∥. We then have tHv + w ∈ int(F * ), which proves (i). We can prove now that the Slater's assumption is equivalent to the non-emptiness and compactness of the optimal set of the reduced problem (11).

Proof of Proposition 3.2. The assertions (i) and (ii) are direct consequences of the weak duality and of the strong convexity of the objective function of (9). The outcome (iii) can be proved by means of some useful tools from asymptotic analysis in convex optimization [START_REF] Auslender | Asymptotic cones and functions in optimization and variational inequalities[END_REF]. Let us define the function

Figure 1 :

 1 Figure 1: Condition number (κ) of the matrix of the linear systems along the iterations when solving two different problems. The left figure is related to a problem with only one contact point, (nd, m) = (3, 3), while the right figure is with sixty-nine contact points, (nd, m) = (207, 450).

Figure 3 :

 3 Figure 3: Condition number (κ) of the matrix of the linear systems along the iterations of the Algorithm 2 when solving the problem PrimitiveSoup-ndof-6000-nc-37-0 with n = 37 contact points.

For

  x ∈ R N , the arrow-shaped matrix is defined byArw(x) = Arw(x 1 ) ⊕ . . . ⊕ Arw(x n ), where Arw(x i ) = x i0 x⊤ i xi x i0 I for i = 1, . . . , n. For i ∈ {1, . . . , n}, let e i = (1, 0) ∈ R × R ni-1and let e = (e 1 , . . . , e n ) ∈ R N be the neutral element of the Jordan product. For vectors x and y in R N , one hasx • y = Arw(x)y = Arw(x)Arw(y)e

Table 1 :

 1 Sizes of problems of the FCLIB collectionAlgorithm 1 One iteration of Mehrotra primal-dual algorithm for solving a FC problem Parameters:

		# Problems	n	m
	BoxStacks	28	[31, 557]	450
	Capsules	200	[15, 314]	600
	Chute	182	[3, 3224]	[1002, 12672]
	KaplasTower	240	[48, 899]	792
	PrimitveSoup	200	[37, 3123]	6000
	Spheres	200	[344, 4290]	12000
	Spheres1mm	41	[715, 4213]	12000

Table 2 :

 2 Number of successes and failures when solving the FCLIB problems with tol = 10 -10 .

	Linear system # Success # Failure 1 # Failure 2
	(26)	1069	22	0
	(27)	947	0	144
	(28)	852	5	234
	(29)	308	1	782
	(29) + refin	1085	0	6
	(30)	1091	0	0
	(31)	1083	1	7

Table 4 :

 4 Sizes of rolling friction problems of the FCLIB collection

	Family	# problems	n	m
	Chute	155	[4, 1372] [768, 6528]
	PrimitiveSoup	171	[37, 2269]	6000
	SpherePile	200	[2, 542]	[24, 1500]

Table 5 :

 5 Number of successes and failures when solving the RF problems with tol = 10 -10

				linear system # success # failure 1 # failure 2 max res
				(38)		391	135		0		3e-3
				(39)		451	0		75		4e-9
				(40)		508	0		18		5e-9
				Iterations						CPU
		1						1
		0.9						0.8
	ρ s (τ )	0.8					ρ s (τ )	0.4 0.6
		0.7						0.2
		0.6	0	0.5	1	1.5	2	0	0	0.5	1	1.5	2
					τ						τ
						(38)	(39)

https://frictionalcontactlibrary.github.io

g(r) = 1 2 r ⊤ W r + q ⊤ r + δ F (r), where δ F is the indicator function of the set F. The asymptotic function of g is defined by g ∞ (s) = sup t>0 g(ts) -g(0) t = q ⊤ s if s ∈ F and W s = 0, ∞ otherwise, for s ∈ R nd , see (Auslender and Teboulle, 2003, Proposition 2.5.2). The optimal set R coincides with the set of minima of g. It is non-empty and compact if and only if g is coercive, that is g ∞ (s) > 0 for all nonzero s ∈ R nd (Auslender and Teboulle, 2003, Proposition 3.1.3). It follows that R is nonempty and compact if and only if there is no nonzero s ∈ F, such that W s = 0 and q ⊤ s ≤ 0. By Lemma B.2, this is equivalent to the Slater hypothesis.

C Proofs of Propositions 4.1 and 4.2

Proof of Proposition 4.1. Let µ > 0. Let us show that the barrier function defined in ( 16) is coercive. As in the proof of Proposition 3.2, we show that for all nonzero s ∈ R nd , (φ µ ) ∞ (s) > 0. Let us write φ µ = p + ψ µ , where p is the quadratic part. We have (φ µ ) ∞ = p ∞ + (ψ µ ) ∞ (Auslender and Teboulle, 2003, Proposition 2.6.1.). We have

To compute the asymptotic derivative of the function ψ µ , let us look at one element of the sum. Let s ∈ R d and let us define ρ(s) := -log det s. We have

We then have (ψ µ ) ∞ = δ L n . It follows that (φ µ ) ∞ = g ∞ , where g is defined in the proof of Proposition 3.2. Therefore, the same conclusion holds, that is: under Assumption 3.1, (φ µ ) ∞ (s) > 0 for all nonzero s ∈ L n . This implies that the set of minima of φ µ is nonempty and compact. Since φ µ is strictly convex, this set is reduced to a singleton.

To show that the central path is bounded, let us proceed to a reasoning by contradiction. Suppose that there exists a positive sequence {µ k }, converging to zero and such that {(r k , u k )} is unbounded, where we denote

Dividing the first equation by t k and the second one by t 2 k , then taking the limit k → ∞, we obtain W r * = u * and r * • u * = 0. Multiplying the first equation by r * on the left, we obtain r * ⊤ W r * = r * ⊤ u * = 0. Since W = HM -1 H ⊤ and M is positive definite, it follows that H ⊤ r * = 0. Multiplying the first equation of ( 42) by

Since W is positive semi-definite, we then get q ⊤ r k ≤ 2nµ k . Dividing by t k and passing through the limit, we obtain q ⊤ r * ≤ 0. Since for all k ∈ N, r k ∈ int(L n ), we also have r * ∈ L n . Finally, we have shown that there exists a non-zero vector r * ∈ L n , such that H ⊤ r * = 0 and w ⊤ r * = q ⊤ r * ≤ 0. In view of Lemma B.2, this is equivalent to the fact that the Slater hypothesis does not hold, which is in contradiction with the assumption.

The proof of Proposition 4.2 is done in two steps. The convergence is proved first. It follows the one of (Halická et al., 2002, Theorem A.3). It is based on the use of the curve selection lemma from algebraic geometry, see (Milnor, 1968, Lemma 3.1). In the second part of the proof, it is shown that the limit is maximally complementary. It is inspired by [START_REF] Klerk | Initialization in semidefinite programming via a self-dual skewsymmetric embedding[END_REF].

Proof of Proposition 4.2. Let (r, û) be a limit point of the central path. Since û is unique, it remains to show that r(µ) → r as µ → 0. Let us define the subsets of R nd+1 :

The set V is a real algebraic set and U is an open set defined by a finite number of polynomial inequalities. Moreover, 0 ∈ R nd+1 is in the closure of U ∩ V. By the curve selection lemma, there exists ε > 0 and a real analytic curve p : [0, ε) → R nd+1 , with p(0) = 0 and for all t ∈ (0, ε), p(t) ∈ U ∩ V. By setting p(t) = (ρ(t), ν(t)), we have (W (r + ρ(t)) + q) • (r + ρ(t)) = 2ν(t)e, for all t ∈ (0, ε). Since the central path is uniquely defined by ( 15), for all t ∈ (0, ε) we have r(ν(t)) = r + ρ(t).

Since ν(0) = 0, ν(t) > 0 for all t ∈ (0, ε) and ν is analytic on [0, ε), there exists ε ′ > 0 such that for all t ∈ (0, ε ′ ), ν ′ (t) > 0. It follows that ν is invertible on (0, ε ′ ). Moreover, for all t ∈ (0, ν(ε ′ )), ν -1 (t) > 0 and lim t→0 ν -1 (t) = 0. We then have, for µ > 0 small enough,

It follows that lim µ→0 r(µ) = r.

Let us show now that r is maximally complementary. That is, for all i ∈ B, ri ∈ int(L) and for all i ∈ R ∪ T 1 , ri ∈ bd(L)\{0}. Let r ∈ R be a maximally complementary optimal solution. For all µ > 0, by using ( 14) and ( 15) and the positive semi-definiteness of W , we have

By dividing by 2µ and by using the partition of the index set, we obtain

Since each term in these sums is nonnegative, we have that for all i ∈ B∪R∪T 1 and for all µ > 0, r ⊤ i r -1 i (µ) ≤ n. By using the inequality (20) in the proof of Theorem 4.5, we deduce that for for all i ∈ B ∪ R ∪ T 1 and all µ > 0, we have 0 < 1 n r i,0 ≤ r i,0 (µ) + ∥r i (µ)∥ < 2r i,0 (µ).

Taking the limit as µ → 0, we have ri,0 > 0. It remains to show that for all i ∈ B, det(r i ) ̸ = 0. For all i and µ > 0, we have

det ri(µ) . It follows that for all i ∈ B and all µ > 0, 0 ≤ r ⊤ i Rr i (µ) ≤ n det(r i (µ)).

Passing through the limit µ → 0, we obtain 0 ≤ r i Rr i ≤ n det(r i ). If det(r i ) = 0, we then have r i,0 ri,0 = r⊤ i ri ≤ ∥r i ∥∥ ri ∥ < r i,0 ∥ ri ∥, from which we deduce that ri / ∈ L, which is in contradiction with the fact that r ∈ K and allows to conclude that r is maximally complementary.