Pedro Gajardo
email: pedro.gajardo@usm.cl

Alberto Seeger
email: alberto.seeger@univ-avignon.fr

Equilibrium problems involving the Lorentz cone

Keywords: Lorentz cone, Euclidean Jordan algebra, Complementarity problem, System of nonlinear equations

We study a general equilibrium model formulated as a smooth system of equations coupled with complementarity conditions relative to the n-dimensional Lorentz cone. For the purpose of analysis, as well as for the design of algorithms, we exploit the fact that the Lorentz cone is representable as a cone of squares in a suitable Euclidean Jordan algebra.

Introduction

A large variety of equilibrium models arising in applications can be formulated as a nonlinear system of equations coupled with complementarity conditions relative to a closed convex cone. The cone under consideration is usually the nonnegative orthant of R n , but sometimes one needs to consider the Lorentz cone also known as the second-order cone or the ice-cream cone. Geometrically speaking, L n is a revolution cone whose central axis is the half-line generated by the n-dimensional vector e n = (0 n-1 , 1). Here and in the sequel the symbol 0 d refers to the d-dimensional zero vector. This paper addresses the issue of analyzing and solving a nonlinear system f (x, y, λ) = 0 n (1)

g(x) = 0 m , (2)
where x, y ∈ R n are unknown vectors linked by means of complementary conditions of the type

0 n x ⊥ y 0 n . ( 3
)
Here ⊥ indicates orthogonality with respect to the usual inner product of R n and is the order relation induced by the Lorentz cone.

The entries of x and y are called primal and dual variables, respectively. The Eq. ( 2) reflects a normalization condition and/or a structural constraint imposed on the primal variables. For achieving more generality in the modeling, the Eq. ( 1) includes an unknown vector of exogenous variables λ 1 , . . . , λ p . Exogenous variables are not restricted by complementarity, so they are different in nature from the pair of primal and dual variables.

Although [START_REF] Andreani | Box-constrained minimization reformulations of complementarity problems in second-order cones[END_REF] cannot be written in a componentwise manner as in the classical theory of complementarity problems, such inconvenience is largely compensated by the rich algebraic structure of the Lorentz cone. We shall elaborate on this issue in Sect. [START_REF] Andreani | Box-constrained minimization reformulations of complementarity problems in second-order cones[END_REF].

The purpose of this work is to study and design algorithms for solving the Lorentz Equilibrium Model LEM( f, g) ⎧ ⎨ ⎩ 0 n x ⊥ y 0 n f (x, y, λ) = 0 n g(x) = 0 m , [START_REF] Aubin | Set-Valued Analysis[END_REF] as well as the particular version LEM( f ) 0 n x ⊥ y 0 n f (x, y, λ) = 0 n [START_REF] Chen | Analysis of nonsmooth vector-valued functions associated with secondorder cones[END_REF] in which the normalization or constraint function g is missing. The name attributed to this kind of equilibrium model emphasizes the role played by the Lorentz cone. Some comments on hypotheses are appropriate:

-For avoiding a situation of overdetermination in (4) one assumes that p ≥ m ≥ 1. The choice p = 0, which by convention means that f is free of exogenous variables, is tolerated in the particular model [START_REF] Chen | Analysis of nonsmooth vector-valued functions associated with secondorder cones[END_REF]. -For simplicity in the presentation one supposes that f : R 2n+ p → R n and g : R n → R m are smooth functions (i.e., continuously differentiable). Much of our analysis extends however to a nonsmooth setting.

The next definition distinguishes between different sorts of solutions to a LEM. The interior and boundary of the Lorentz cone are given respectively by

int(L n ) = x ∈ R n : x n > x 2 1 + • • • + x 2 n-1 1/2 , bd(L n ) = x ∈ R n : x n = x 2 1 + • • • + x 2 n-1 1/2 .
Definition 1.1 A solution (x * , y * , λ * ) to a LEM is called trivial if x * = 0 n , otherwise it is nontrivial. Furthermore, (i) A nontrivial solution is of interior type (respectively, boundary type) if x * belongs to the interior (respectively, boundary) of L n . (ii) An interior type solution is central if x * lies on the central axis of L n , otherwise it is eccentric.

Remark 1.2 As a consequence of Moreau's decomposition theorem [START_REF] Moreau | Décomposition orthogonale d'un espace hilbertien selon deux cônes mutuellement polaires[END_REF], the component y * of an interior type solution is necessarily equal to 0 n . By contrast, for a boundary type solution one may have y * = 0 n or not.

The remaining part of this introductory section reviews some interesting examples of equilibrium problems that fit into the framework of a LEM.

LCP relative to the Lorentz cone

Let M n be the space of square matrices of order n. The Linear Complementarity Problem relative to the Lorentz cone reads as follows:

(LCP) ⎧ ⎪ ⎨ ⎪ ⎩
given a matrix M ∈ M n and a vector q ∈ R n , find a vector x ∈ R n such that

0 n x ⊥ (M x + q) 0 n . (6)
In this example the Eq. ( 1) is free of exogenous variables and the Eq. ( 2) does not show up, because neither normalization nor structural constraints are imposed on x. So, one needs to solve the particular LEM

0 n x ⊥ y 0 n M x + q -y = 0 n .
LCP's relative to the Lorentz cone are considered in [START_REF] Kong | The solution set structure of monotone linear complementarity problems over second-order cone[END_REF][START_REF] Malik | On Q and R 0 properties of a quadratic representation in linear complementarity problems over the second-order cone[END_REF]. Existence results for (6) can be derived by adapting [START_REF] Gowda | Some P-properties for linear transformations on Euclidean Jordan algebras[END_REF]Theorem 12] or various abstract theorems stated in [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF]Chapter 2].

As a mathematical curiosity we mention a variant of [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF] in which the matrix M depends on a certain number of exogenous variables and the solution x is sought for instance on a unit sphere:

⎧ ⎨ ⎩ given a smooth function M : R p → M n and a vector q ∈ R n , find λ ∈ R p and a unit vector

x ∈ R n such that 0 n x ⊥ (M(λ)x + q) 0 n .
The later equilibrium problem corresponds to the particular LEM ⎧ ⎨ ⎩

0 n x ⊥ y 0 n M(λ)x + q -y = 0 n x 2 -1 = 0. ( 7
)
1.2 AVE relative to the Lorentz cone Given a vector b ∈ R n and matrices A, B ∈ M n , one must solve the Absolute Value Equation

Az + B|z| = b. ( 8
)
The absolute value of z ∈ R n is taken with respect to the Lorentz cone, i.e., |z| = z + + z - with z + denoting the projection of z onto L n and z -= z +z. AVE's relative to the nonnegative orthant are studied in [START_REF] Mangasarian | Absolute value programming[END_REF][START_REF] Mangasarian | A generalized Newton method for absolute value equations[END_REF][START_REF] Mangasarian | Absolute value equations[END_REF][START_REF] Prokopyev | On equivalent reformulations for absolute value equations[END_REF][START_REF] Rohn | A theorem of the alternatives for the equation Ax + B|x| = b[END_REF] and in many other places. AVE's relative to the Lorentz cone are considered in [START_REF] Hu | A generalized Newton method for absolute value equations associated with second order cones[END_REF]. As pointed out in [9, Theorem 1.1], the Eq. ( 8) can be reformulated as

0 n x ⊥ y 0 n (B + A)x + (B -A)y -b = 0 n ,
i.e., as a particular instance of the model (5).

Second-order cone programming

Consider a conic optimization problem of the form inf{c(x) :

x 0 n , Ax = b}, ( 9
)
where c : R n → R is a twice continuously differentiable function, b ∈ R m , and A is a real matrix of size m × n. This type of optimization problem has been widely studied in the last decade, specially when c is a linear function (cf. [START_REF] Alizadeh | Second-order cone programming[END_REF]). If one introduces a Lagrange multiplier vector λ ∈ R m for the equality constraint Ax = b and a Karush-Kuhn-Tucker multiplier vector y ∈ R n for the nonnegativity constraint x 0 n , then the optimality conditions for ( 9) can be expressed in the form

x 0 n , Ax = b primal feasibility y 0 n dual feasibility

x, y = 0 complementarity slackness ∇c(x) -A T λy = 0 n stationarity, where the superscript "T" indicates transposition. One gets in this way a LEM ⎧ ⎨ ⎩

0 n x ⊥ y 0 n ∇c(x) -A T λ -y = 0 n Ax -b = 0 m (10)
in which p is equal to m. If the triplet (x * , y * , λ * ) solves [START_REF] Kong | The solution set structure of monotone linear complementarity problems over second-order cone[END_REF], then under some conditions, like for instance the convexity of the cost function c, the primal component x * is a solution to the optimization problem (9).

Lorentz eigenvalue problem

The Lorentz eigenvalue problem is formulated in [START_REF] Seeger | On eigenvalues induced by a cone constraint[END_REF] in the following terms:

⎧ ⎨ ⎩ given a matrix A ∈ M n , find λ ∈ R and a nonzero vector x ∈ R n such that 0 n x ⊥ (Ax -λx) 0 n . ( 11
)
By Corollary 2.1 in [START_REF] Seeger | Eigenvalue analysis of equilibrium processes defined by linear complementarity conditions[END_REF] one knows that (11) admits always a solution. If (x * , λ * ) solves [START_REF] Lancaster | Lambda-Matrices and Vibrating Systems[END_REF], then x * is called a Lorentz eigenvector of A and λ * is an associated Lorentz eigenvalue. By a positive homogeneity argument, there is no loss of generality in assuming that x * satisfies the equation e n , x = 1. Such normalization condition guarantees that x * is a nonzero vector. So, the Lorentz eigenvalue problem [START_REF] Lancaster | Lambda-Matrices and Vibrating Systems[END_REF] is nothing but the particular LEM ⎧ ⎨ ⎩

0 n x ⊥ y 0 n Ax -λx -y = 0 n e n , x -1 = 0. ( 12
)
Theoretical aspects of the Lorentz eigenvalue problem has been discussed in [START_REF] Seeger | On eigenvalues induced by a cone constraint[END_REF]. Numerical methods for solving [START_REF] Lasserre | Moments and sums of squares for polynomial optimization and related problems[END_REF] has been proposed in [START_REF] Adly | Nonsmooth algorithms for cone-constrained eigenvalue problems[END_REF] and [20, Section 2.2].

Lorentz quadratic eigenvalue problem

The general theory of cone-constrained quadratic eigenvalue problems is developed in [START_REF] Seeger | Quadratic eigenvalue problems under conic constraints[END_REF].

The quadratic version of (11) involves a matrix

Q(λ) := λ 2 A + λB + C ( 13
)
that depends quadratically on the variable λ ∈ R. One refers to the function Q : R → M n as the quadratic pencil associated to the triplet

(A, B, C) ∈ T n := M n × M n × M n .
The vector space of quadratic pencils is identified with T n . The formulation of the Lorentz quadratic eigenvalue problem reads as follows:

⎧ ⎨ ⎩ given a quadratic pencil Q ∈ T n , find λ ∈ R and a nonzero vector x ∈ R n such that 0 n x ⊥ Q(λ)x 0 n . (14)
If (x * , λ * ) is a solution to the above problem, then x * is called a Lorentz eigenvector of Q and λ * is an associated Lorentz eigenvalue. Of course, one is concerned here with the particular LEM ⎧ ⎨ ⎩

0 n x ⊥ y 0 n Q(λ)x -y = 0 n e n , x -1 = 0. (15)
The problems [START_REF] Lancaster | Lambda-Matrices and Vibrating Systems[END_REF] and ( 14) look similar, but they differ in a number of aspects. For instance, the problem (11) admits always a solution, whereas [START_REF] Mangasarian | A generalized Newton method for absolute value equations[END_REF] may not be solvable.

An existence result for LEMs of explicit type

A LEM is of explicit type if f has the special structure f (x, y, λ) = h(x, λ)y, where h : R n+ p → R n is a smooth function. In such case the dual variables can be expressed explicitly in terms of the primal and exogenous variables. The next proposition is an existence result for the explicit LEM ⎧ ⎨ ⎩

0 n x ⊥ y 0 n h(x, λ) -y = 0 n g(x) = 0 m . ( 16
)
Its proof is based on a celebrated theorem by Ky Fan on the existence of solutions to variational inequalities.

Proposition 2.1

The combination of the following two hypotheses ensures that [START_REF] Markus | Introduction to the Spectral Theory of Polynomial Operator Pencils[END_REF] is solvable:

(i) := {x ∈ L n : g(x) = 0 m } is a convex compact set such that L n = {tw : t ≥ 0, w ∈ }.
(ii) For all x ∈ , the equation x, h(x, λ) = 0 has a solution λ x that depends continuously on x.

Proof Recall that h is assumed to be smooth. In particular, h is continuous. Consider the function : × → R given by

(x, w) = w, h(x, λ x ) .
Since is a convex compact set and ⎧ ⎨ ⎩ for all w ∈ , ( • , w) is continuous, for all x ∈ , (x, • ) is linear, for all x ∈ , (x, x) = 0, Ky Fan's theorem (cf. [4, Theorem 3.1.1]) ensures the existence of a point x * ∈ such that (x * , w) ≥ 0 for all w ∈ . If one sets λ * = λ x * and y * = h(x * , λ * ), then it is not difficult to check that (x * , y * , λ * ) is a solution to [START_REF] Markus | Introduction to the Spectral Theory of Polynomial Operator Pencils[END_REF]. The details are omitted.

The hypothesis (i) is strong in general. For instance, it does not hold for the explicit LEM [START_REF] Gatermann | Symmetry groups, semidefinite programs, and sums of squares[END_REF] because the set

{x ∈ L n : x 2 -1 = 0}
is not convex. However, the hypothesis (i) is in force for the explicit LEMs [START_REF] Lasserre | Moments and sums of squares for polynomial optimization and related problems[END_REF] and [START_REF] Mangasarian | Absolute value equations[END_REF]. The next corollary is obtained as a consequence of Proposition 2.1 or by specializing [START_REF] Seeger | Quadratic eigenvalue problems under conic constraints[END_REF]Theorem 3.3] to the case of the Lorentz cone.

Corollary 2.2

Let Q ∈ T n be a quadratic pencil as in [START_REF] Mangasarian | Absolute value programming[END_REF]. Suppose that x, Ax = 0 for all x 0 n , x = 0 n and that Q is Lorentz hyperbolic in the sense that x, Bx 2 ≥ 4 x, Ax x, C x for all x 0 n . Then (14) admits a solution.

The companion system associated to a LEM

The main source of difficulties in a LEM is the presence of the nonnegativity constraints x 0 n and y 0 n . Fortunately, it is possible to get rid of these constraints by introducing the change of variables

x = u • u = 2r ξ, ξ 2 + r 2 , ( 17
)
y = v • v = 2sη, η 2 + s 2 , ( 18
)
where • is a vector product on R n defined by

u • v := (r η + sξ, ξ, η + rs) .
In the above line and in the sequel one uses the decomposition

u = (ξ, r ), v = (η, s) with ξ, η ∈ R n-1 and r, s ∈ R. It is well known (cf. [8, Example 2.0]) that the triplet (R n , •, • , •)
satisfies the axioms of an Euclidean Jordan algebra and that

L n = {w [2] :

w ∈ R n } (19)
with w [2] := w • w. Hence, LEM( f, g) admits the equivalent formulation ⎧ ⎪ ⎨ ⎪ ⎩ u [2] , v [2] = 0 f (u [2] , v [2] , λ) = 0 n g(u [2] )

= 0 m . ( 20
)
A result from the theory of Euclidean Jordan algebras (cf. [START_REF] Gowda | Some P-properties for linear transformations on Euclidean Jordan algebras[END_REF]Proposition 6]) asserts that the scalar equation u [2] , v [2] = 0 can be written in the vector form u [2] • v [2] = 0 n . Hence, [START_REF] Pinto Da Costa | Numerical resolution of cone-constrained eigenvalue problems[END_REF] can be reformulated as

⎧ ⎪ ⎨ ⎪ ⎩
u [2] • v [2] = 0 n f (u [2] , v [2] , λ) = 0 n g(u [2] )

= 0 m . ( 21
)
By the way, if p = m (as occurs quite often in practice), then ( 21) is a square system, i.e., the number of equations is equal to the number of unknown variables. Note that ( 21) fits into the abstract framework of a system of nonlinear equations

(z) = 0, ( 22
)
where 0 is a zero vector of appropriate size and : Z → W is a smooth function between Euclidean spaces of possible different dimensions. When it comes to solve [START_REF] Qi | A nonsmooth version of Newton's method[END_REF] numerically, there are two cases for consideration:

-Square case. If dim(W) = dim(Z), then ( 22) is a square system, which can be solved by initializing Newton's method

z τ +1 = z τ -[ (z τ )] -1 (z τ )
at a suitable starting point z 0 . Here (z τ ) : 22) is underdetermined, i.e., the number of unknown variables is greater than the number of equations. In such a situation one uses the Normal Flow Algorithm (NFA)

Z → W is the differential of at z τ . -Underdetermined case. If dim(W) < dim(Z), then (
z τ +1 = z τ -[ (z τ )] † (z τ ),
where L † : W → Z denotes the Moore-Penrose inverse of a linear map L : Z → W. A detailed convergence analysis of the NFA can be found in [START_REF] Walker | Newton-like methods for underdetermined systems[END_REF][START_REF] Walker | Least-change secant update methods for underdetermined systems[END_REF].

The general theory of nonlinear equations tells us that an element z * of the solution set

-1 (0) = {z ∈ Z : (z) = 0}
is difficult to detect numerically if the differential map (z * ) is nonsurjective. In particular, for the system (21) it will be difficult to detect a solution z * = (u * , v * , λ * ) with v * = 0 n . Indeed, the corresponding differential map (u * , 0 n , λ * ) is clearly nonsurjective.

The phenomenon of lack of surjectivity in ( 21) can be remediated by shifting the attention to a slightly modified system, namely [2] , v [2] , λ) = 0 n g(u [2] ) = 0 m .

⎧ ⎨ ⎩ u • v = 0 n f (u
(

) 23
One refers to (23) as the companion system associated to LEM( f, g). As we shall see in Lemma 3.3, the equality u 2 , v 2 = 0 implies that u •v = 0 n . Since the reverse implication is false, a solution to [START_REF] Rohn | A theorem of the alternatives for the equation Ax + B|x| = b[END_REF] does not solve necessarily [START_REF] Pinto Da Costa | Numerical resolution of cone-constrained eigenvalue problems[END_REF]. 23), but x = u [2] and y = v [2] are not orthogonal, then one refers to the triplet (x, y, λ) as a fake solution to LEM( f, g). Fake solutions to LEM( f ) are defined in a similar way.

Definition 3.1 If (u, v, λ) solves (
Example 3.2 By way of illustration we display a fake solution for the LCP stated in [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF] with

M = ⎡ ⎢ ⎣ 1 4 √ 5 2 5 -1 3 6 3 ⎤ ⎥ ⎦ , q = ⎡ ⎢ ⎣ -5 √ 5 - √ 5 ⎤ ⎥ ⎦ .
Here n = 3 and p = 0. The pair (u, v) = ((-1, 2, 0), (4, 2, 0)) solves the companion system u • v = 0 3 Mu [2] + q -v [2] = 0 3 , but u [2] = 0, 0, √ 5 and v [2] = 0, 0, 2 √ 5 are not orthogonal.

The above example has been constructed artificially. A vast majority of solutions to the companion system [START_REF] Rohn | A theorem of the alternatives for the equation Ax + B|x| = b[END_REF] satisfy the orthogonality condition u [2] , v [2] = 0. Said in other words, getting a fake solution to a LEM is an exception and not the general rule. This theme is developed next.

On fake solutions

We start by stating a technical lemma. Lemma 3.3 For all u, v ∈ R n , the orthogonality condition u [2] , v [2] = 0 implies u•v = 0 n . Furthermore, the following statements are equivalent: [2] and v [2] are not orthogonal. (b) ξ, η ∈ R n-1 are nonzero orthogonal vectors and r = s = 0.

(a) u • v = 0 n , but u
Proof The equality u [2] , v [2] = 0 says that

4rs ξ, η + ξ 2 + r 2 η 2 + s 2 = 0.
This can be rearranged as

ξ 2 η 2 -ξ, η 2 + r η + sξ 2 + ( ξ, η + rs) 2 = 0.
Since ξ, η ≤ ξ η , the above equality breaks down into three pieces:

ξ 2 η 2 -ξ, η 2 = 0 ( 2 4
)

r η + sξ = 0 (25) ξ, η + rs = 0. ( 26
)
Note that ( 25)-( 26) says that u • v = 0 n . It is possible to have ( 25)-( 26) without having [START_REF] Seeger | Eigenvalue analysis of equilibrium processes defined by linear complementarity conditions[END_REF]. This happens exactly when ξ, η are nonzero orthogonal vectors and r = s = 0.

Fake solutions are easily recognizable. In fact, one has:

Proposition 3.4 If (x, y, λ
) is a fake solution to a LEM, then x and y are nonzero vectors on the central axis of L n .

Proof Let (x, y, λ) be a fake solution to a LEM. By Lemma 3.3 one has (x, y) = u [2] , v [2] with u = (ξ, 0), v = (η, 0), and ξ, η ∈ R n-1 \{0 n-1 } mutually orthogonal. Hence, x = 0 n-1 , ξ 2 and y = 0 n-1 , η 2 are nonzero vectors on the central axis of L n .

By way of illustration, consider the Lorentz eigenvalue problem as formulated in [START_REF] Lasserre | Moments and sums of squares for polynomial optimization and related problems[END_REF]. The system [START_REF] Prokopyev | On equivalent reformulations for absolute value equations[END_REF] becomes

⎧ ⎪ ⎨ ⎪ ⎩
u [2] • v [2] = 0 n (A -λI n )u [2] -v [2] = 0 n e n , u [2] -

1 = 0, ( 27
)
where I n is the identity matrix of order n. This is a system of 2n + 1 equations in the same number of unknown variables. The companion counterpart [2] -v [2] = 0 n e n , u [2] -1 = 0

⎧ ⎨ ⎩ u • v = 0 n (A -λI n )u
is a system of the same size. Unless the matrix A ∈ M n is constructed in an artificial way, the systems ( 27) and ( 28) have the same solution set. A probabilistic formulation of this statement reads as follows:

Proposition 3.5 Let A ∈ M n be a random matrix with absolutely continuous probability distribution. Then, almost surely, the Lorentz eigenvalue problem [START_REF] Lasserre | Moments and sums of squares for polynomial optimization and related problems[END_REF] is free of fake solutions.

Proof If (x, y, λ) is a fake solution to [START_REF] Lasserre | Moments and sums of squares for polynomial optimization and related problems[END_REF], then x = e n and y = α e n with α > 0. The second equation in [START_REF] Lasserre | Moments and sums of squares for polynomial optimization and related problems[END_REF] becomes Ae n = (α + λ)e n . Hence, the last column of A is a multiple of e n , i.e., a 1,n = 0, a 2,n = 0, . . . , a n-1,n = 0.

Such event occurs with probability zero because A ∈ M n is a random matrix with absolutely continuous probability distribution.

Numerical tests with the companion system

The programming language used for all numerical tests was Scilab. The experiment reported in Table 1 compares Newton's method applied to (27) and Newton's method applied to the companion counterpart [START_REF] Walker | Newton-like methods for underdetermined systems[END_REF]. The idea is measuring the rate of success as function of n and number of initial points. A few practical considerations are in order:

-An iterative algorithm like Newton's method (or the NFA) keeps running until one of the following three termination criteria occurs: τ = 1000 (maximum number of iterations), κ( (z τ )) ≥ 10 5 (ill-conditioning), (z τ ) ≤ 10 -8 (a solution is found). Here κ(M) = σ max (M)/σ min (M) refers to the condition number of a matrix M. -An iterative algorithm is declared successful if a solution is found while working with a certain tolerated Number of Initial Points (for instance, NIP = 1, NIP = 10, or NIP = 10 2 ). -In general, there is no optimal rule for selecting a suitable initial point. For the systems ( 27) and ( 28), the initial point z 0 = (u 0 , v 0 , λ 0 ) is constructed as follows:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ u 0 = uniformly distributed on the unit sphere of R n x 0 = (u 0 ) [2] λ 0 = x 0 , Ax 0 / x 0 2 v 0 = n-dimensional Gaussian vector z 0 = (u 0 , v 0 , λ 0 ).
A Gaussian vector (or matrix) is a vector (or matrix) whose the entries are independent random variables with standard normal distributions.

As one can see from Table 1, Newton's method performs poorly on [START_REF] Tisseur | The quadratic eigenvalue problem[END_REF]. To understand why the percentages of success are so low, recall that Newton's method is unsuccessful if one gets a current point z τ at which (z τ ) is ill-conditioned. This is exactly what is happening while handling the ill-conditioned system [START_REF] Tisseur | The quadratic eigenvalue problem[END_REF]. Newton's method performs much better on the companion system [START_REF] Walker | Newton-like methods for underdetermined systems[END_REF].

Remark 3.6

The performance of Newton's method on [START_REF] Tisseur | The quadratic eigenvalue problem[END_REF] does not improve if one considers higher ill-conditioning conditions. We did a few numerical experiments with κ( (z τ )) ≥ 10 8 and even with κ( (z τ )) ≥ 10 16 , but the percentages of success remain essentially unchanged.

The numerical experiment reported in Table 3 is different in spirit. The purpose of this experiment is to see whether Newton's method on (28) tends to favor a certain type of solution, be it central, eccentric, or boundary type. One considers the test matrix

A = ⎡ ⎢ ⎢ ⎢ ⎣ 1 0 0 0 -5 6 0 0 -2 0 5 0 -6 4 1 4 ⎤ ⎥ ⎥ ⎥ ⎦ , ( 29
)
which has exactly six Lorentz eigenvalues: one coming from a central solution, one coming from an eccentric solution, and four coming from boundary type solutions. The details are The percentages of failure (ill-conditioning or divergence) and the percentages of converge to the different solutions are estimated with a random sample of 10 4 initial points displayed in Table 2. Beware that the classical eigenvalue λ = 1 is not a Lorentz eigenvalue because the associated eigenvector x = (2, 2, 1, 1) is not in the Lorentz cone. Conversely, the Lorentz eigenvalues λ = 2, λ = 3, and λ = 7 are not eigenvalues in the classical sense because the corresponding dual vectors y are nonzero. Table 3 shows that the Lorentz eigenvalue λ = 2 is the most likely to be detected, then comes the Lorentz eigenvalue λ = 3, and so on. In this numerical experiment the particular solution z * = (x * , y * , λ * ) = ((0, 0, 1, 1), (0, 0, 0, 0), 5) was never detected. This fact can be explained by recalling that Newton's method is not suitable for finding z * ∈ -1 (0) such that (z * ) is nonsurjective.

The use of natural coordinates

We come back again to the squaring technique based on the change of variables ( 17)-( 18) and characterize the solutions to a LEM in terms of the "natural" coordinates ξ, r, η, s, λ. For the sake of clarity in the presentation, we distinguish between trivial and nontrivial solutions.

The trivial case

Searching for a trivial solution to a LEM makes sense only if g(0 n ) = 0 m or if the function g is missing altogether. Since the component x of a trivial solution is required to be zero, one must substitute u = 0 n into [START_REF] Pinto Da Costa | Numerical resolution of cone-constrained eigenvalue problems[END_REF]. Such a substitution leads to

f (η, s, λ) = 0 n , ( 30
)
where f : R n+ p → R n is the smooth function given by

f (η, s, λ) = f 0 n , 2sη, η 2 + s 2 , λ .
Note that (30) is a square system if p = 0 and it is an underdetermined system if p ≥ 1.

Example 4.1 By way of illustration, consider the feasibility problem:

given smooth functions M : R → M n and q : R → R n , find λ ∈ R and y 0 n such that M(λ) y = q(λ).

(31)

In other words, one needs to find a value of the parameter λ so that the parametric linear system M(λ) y = q(λ) has a solution y in the Lorentz cone. One may see (31) as the problem of finding a trivial solution to the particular LEM

0 n x ⊥ y 0 n M(λ) y -q(λ) = 0 n .
There is no normalization condition on the primal variables. Finding a solution to the feasibility problem (31) is then a matter of solving (30) with

f (η, s, λ) = M(λ) 2sη η 2 + s 2 -q(λ).
What we are doing, in essence, is to express y in terms of the coordinates η and s. In this example one has p = 1, and therefore one must solve an underdetermined system of nonlinear equations.

The nontrivial case

This case is the most interesting one, but it must be handled with care. For the sake of simplicity, from now on we assume that

g(0 n ) = 0 m . ( 32
)
Such hypothesis ensures that every solution to LEM( f, g) is nontrivial. Note that (32) is in force for the particular LEM's ( 7), [START_REF] Lasserre | Moments and sums of squares for polynomial optimization and related problems[END_REF], and [START_REF] Mangasarian | Absolute value equations[END_REF].

Lemma 4.2 Consider the change of variables (17)-(18)

. Under the assumption (32), the triplet (x, y, λ) is a solution to LEM( f, g) if and only if one of the following situations occur:

(a) r = 0 and (ξ, r, η, s, λ) solves the companion system

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ r η + sξ = 0 n-1 ξ, η + rs = 0 F(ξ, r, η, s, λ) = 0 n G(ξ, r ) = 0 m , ( 33
)
where

F(ξ, r, η, s, λ) = f 2r ξ, ξ 2 + r 2 , 2sη, η 2 + s 2 , λ G(ξ, r ) = g 2r ξ, ξ 2 + r 2 .
(b) r = 0, η = 0 n-1 , s = 0, and (ξ, λ) satisfies

f 0 n-1 , ξ 2 , 0 n , λ = 0 n g 0 n-1 , ξ 2 = 0 m . ( 34
)
Proof As we saw already in the proof of Lemma 3.3, the orthogonality condition in ( 20) is equivalent to ( 24)-( 26). Hence, (x, y, λ) solves ( 4) if and only if (ξ, r, η, s, λ) satisfies ( 24)-( 26) and

F(ξ, r, η, s, λ) = 0 n G(ξ, r ) = 0 m .
Depending on the value of r , there are two cases for consideration. If r = 0, then the Eq. ( 24) becomes redundant because it can be deduced from [START_REF] Seeger | Quadratic eigenvalue problems under conic constraints[END_REF]. One ends up with the system (33).

If r = 0, then the assumption (32) implies that ξ = 0 n-1 . Hence, ( 24)-( 26) is equivalent to saying that η = 0 n-1 and s = 0. One ends up with the system (34).

The new system (33) is the same, of course, as the old companion system [START_REF] Rohn | A theorem of the alternatives for the equation Ax + B|x| = b[END_REF], except that now everything is expressed in terms of the natural coordinates.

Boundary type and interior type solutions

The next theorem is the main result of Sect. 4.

Theorem 4.3 Let the assumption (32) be in force. The triplet (x, y, λ) solves LEM( f, g) if and only if one of the following situations occur:

(α) x = 0 n-1 , ξ 2 , y = 0 n , and (ξ, λ) solves (34). (β) x = 2r 2 (ω, 1), y = 2s 2 (-ω, 1), and (ω, r, s, λ) solves

⎧ ⎨ ⎩ F bd (ω, r, s, λ) = 0 n G bd (ω, r ) = 0 m ω 2 -1 = 0 ( 35
)
where and(ω, r, λ) solves

F bd (ω, r, s, λ) = f 2r 2 (ω, 1) , 2s 2 (-ω, 1) , λ G bd (ω, r ) = g 2r 2 (ω, 1) . (γ ) x = r 2 2ω, 1 + ω 2 , y = 0 n ,
F (ω, r, λ) = 0 n G(ω, r ) = 0 m , ( 36
)
where

F (ω, r, λ) = f r 2 2ω, 1 + ω 2 , 0 n , λ G(ω, r ) = g r 2 2ω, 1 + ω 2 .
Proof Consider again the change of variables ( 17)- [START_REF] Moreau | Décomposition orthogonale d'un espace hilbertien selon deux cônes mutuellement polaires[END_REF]. Let us examine more closely the part (a) of Lemma 4.2. The condition r = 0 and the first equation in (33) yield η = -sω with ω = ξ/r . Hence, Lemma 4.2(a) says that

x = r 2 2ω, 1 + ω 2 , ( 37
)
y = s 2 -2ω, 1 + ω 2 , ( 38
)
where (ω, r, s, λ) solves the following reduced form of the companion system:

s 1 -ω 2 = 0 (39) f r 2 2ω, 1 + ω 2 , s 2 -2ω, 1 + ω 2 , λ = 0 n (40) g r 2 2ω, 1 + ω 2 = 0 m . ( 41
)
Note that we got rid of η and changed ξ by ω. The Eq. (39) says that either ω = 1 or s = 0. By plugging ω = 1 into (37)-( 38) and ( 40)-(41) one gets (β). By plugging s = 0 into (38) and ( 40)-(41) one gets (γ ).

The case (α) is disjoint from (β) and also disjoint from (γ ). However, the cases (β) and (γ ) are not disjoint. Indeed, by setting s = 0 in (β) one recovers the same solutions as those obtained by setting ω = 1 in (γ ). We end this section by giving a full classification of the solutions to LEM( f, g). The main merit of the next corollary is to provide a clear-cut characterization of the boundary type solutions. Proof Part (a). If (α) occurs, then (x, y, λ) is a central solution because x = 0 n-1 , ξ 2 lies on the central axis of the Lorentz cone. Note that ξ = 0 n-1 thanks to (32). Conversely, if (x, y, λ) is a central solution, then y is necessarily equal to 0 n . Furthermore, x = (0 n-1 , t) with t being a positive scalar such that

f ((0 n-1 , t) , 0 n , λ) = 0 n g (0 n-1 , t) = 0 m .
It suffices then to write t = ξ 2 for a suitable ξ ∈ R n-1 . Part (b). If (β) occurs, then (x, y, λ) is a boundary type solution. Indeed,

x n -x 2 1 + • • • + x 2 n-1 1/2 = 2r 2 (1 -ω ) = 0
thanks to the last equation in (35). Conversely, suppose that (x, y, λ) is a boundary type solution. Let x and y be as in ( 17)- [START_REF] Moreau | Décomposition orthogonale d'un espace hilbertien selon deux cônes mutuellement polaires[END_REF]. The case r = 0 must be excluded, so one can rewrite x and y as in (37)-(38). Since x ∈ bd(L n ) and r = 0, one deduces that 2ω, 1

+ ω 2 ∈ bd(L n ).
But this amounts to saying that 2ω = 1 + ω 2 , which in turn is equivalent to ω = 1. One arrives in this way to the condition (β). Part (c). If (γ ) occurs with ω / ∈ {0, 1}, then (x, y, λ) is an eccentric solution. Indeed

2ω, 1 + ω 2 / ∈ bd(L n ),
and x = r 2 2ω, 1 + ω 2 belongs to the interior of L n . Furthermore, x is not on the central axis of L n because ω = 0 n-1 . Conversely, let (x, y, λ) be an eccentric solution. The cases (α) and (β) must be excluded, and therefore (γ ) occurs. Since x is in the interior of the Lorentz cone, but not on its central axis, one necessarily has ω / ∈ {0, 1}.

Analysis of the Lorentz quadratic eigenvalue problem

Testing Lorentz hyperbolicity is not a trivial matter. Recall, from Corollary 2.2, that a quadratic pencil Q is Lorentz hyperbolic if the discriminant function [2] ) is an homogeneous multivariate polynomial of degree 8. The representation [START_REF] Parrilo | Computing sum of squares decompositions with rational coefficients[END_REF] of the cone L n leads straightforwardly to the following characterization of Lorentz hyperbolicity.

x ∈ R n → δ Q (x) = x, Bx 2 -4 x, Ax x, C x is nonnegative on L n . Note that δ Q is an homogeneous multivariate polynomial of degree 4. Hence, u ∈ R n → Q (u) = δ Q (u
Proposition 5.1 A quadratic pencil Q is Lorentz hyperbolic if and only if the homogeneous polynomial Q is nonnegative.

The polynomial Q is certainly nonnegative if it can be written as a sum of squares (SOS) of other polynomials in u. Hence, the existence of a SOS representation for Q is a sufficient condition for Q to be Lorentz hyperbolic. According to the specialized literature [START_REF] Gatermann | Symmetry groups, semidefinite programs, and sums of squares[END_REF][START_REF] Lasserre | Moments and sums of squares for polynomial optimization and related problems[END_REF][START_REF] Parrilo | Computing sum of squares decompositions with rational coefficients[END_REF], checking whether a given homogeneous multivariate polynomial admits a SOS representation reduces to solving a semidefinite program, a convex optimization problem that one may solve efficiently to arbitrary precision, in polynomial time in the input size.

Boundary type solutions to the Lorentz quadratic eigenvalue problem

For the particular LEM given in [START_REF] Mangasarian | Absolute value equations[END_REF] one has

F bd (ω, r, s, λ) = 2r 2 Q(λ) ω 1 -2s 2 -ω 1 G bd (ω, r ) = 2r 2 -1.
The boundary type solutions to [START_REF] Mangasarian | Absolute value equations[END_REF] are therefore found by solving

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 2r 2 Q(λ) ω 1 -2s 2 -ω 1 = 0 n 2r 2 -1 = 0 ω 2 -1 = 0.
By getting rid of r and by making the change of variables μ = √ 2s, one obtains

⎧ ⎨ ⎩ Q(λ) ω 1 -μ 2 -ω 1 = 0 n ω 2 -1 = 0. ( 42
)
Remark 5.2 If one introduces the diagonal matrix

J n := Diag(1, . . . , 1, -1) = I n-1 0 n-1 0 T n-1 -1 ,
⎧ ⎨ ⎩ Q μ (λ) ω 1 = 0 n ω 2 -1 = 0 (43) for the pencil Q μ (λ) = λ 2 A + λB + C + μ 2 J n .
There is a very rich literature devoted to the numerical analysis of classical quadratic eigenvalue problems, see for instance the books [START_REF] Lancaster | Lambda-Matrices and Vibrating Systems[END_REF][START_REF] Markus | Introduction to the Spectral Theory of Polynomial Operator Pencils[END_REF] or the survey paper [START_REF] Tisseur | The quadratic eigenvalue problem[END_REF]. A bothersome aspect of ( 43) is that the quadratic pencil Q μ depends on the unknown parameter μ.

The system (42) is smooth and concerns n + 1 nonlinear equations in the same number of unknown variables. Table 4 illustrates the performance of Newton's method on (42). Some technical comments are in order. If (ω, μ, λ) solves (42), then clearly

μ = ± 1 2 Q (λ) ω 1 , -ω 1 1/2 .
Hence, a natural way to construct an initial point z 0 = (ω 0 , μ 0 , λ 0 ) is as follows:

-one generates a random vector ω 0 with uniform distribution on the unit sphere of R n-1 , a discrete random variable δ with uniform distribution on {-1, 1}, and a random variable λ 0 with uniform distribution on some interval [a n , b n ]. Then one sets

μ 0 = δ 1 2 Q λ 0 ω 0 1 , -ω 0 1 1/2
.

For simplicity we take a n = -10 and b n = 10 for all n, but one may consider more sophisticated options.

Percentages of success are estimated by averaging the results obtained with a sample of 10 3 random quadratic pencils Q ∈ T n . In order to enhance the chances of (42) to have a solution, each random quadratic pencil is constructed as follows: one generates Gaussian symmetric matrices Ã, B, C and takes Q as the pencil associated to

(A, B, C) = PSD ( Ã), B, -PSD (-C) ,
where PSD (•) stands for the projection on the Loëwner cone of positive semidefinite symmetric matrices. Such Q is obviously Lorentz hyperbolic.

Interior type solutions to the Lorentz quadratic eigenvalue problem

Consider now the problem of finding interior type solutions to ( 15), be them central or eccentric. Both cases can be treated simultaneously. The system (36) takes here the form

⎧ ⎨ ⎩ r 2 Q(λ) 2ω 1 + ω 2 = 0 n r 2 1 + ω 2 -1 = 0.
Since r cannot be equal to 0, the above system can we rewritten as

⎧ ⎨ ⎩ Q(λ) 2r 2 ω 1 = 0 n r 2 1 + ω 2 -1 = 0.
Returning to the original variable ξ = r ω, one ends up with the equivalent system

⎧ ⎨ ⎩ Q(λ) 2r ξ 1 = 0 n ξ 2 + r 2 -1 = 0. ( 44
)
The system (44) produces all the central solutions (obtained with r ξ = 0) and all the eccentric solutions (obtained with r ξ = 0, ξ 2 = r 2 ).

Remark 5.3 Beware that (44) may produce also a few boundary type solutions. To see this, just consider the configuration r ξ = 0, ξ 2 = r 2 . Of course, when it comes to search for the whole collection of boundary type solutions it is better to shift the attention to the system (42).

6 Alternative techniques and extensions

Techniques of nonsmooth analysis

Another way of getting rid of the nonnegativity constraints x 0 n and y 0 n is to use a complementarity function associated to the Lorentz cone. This terminology refers to any function

C : R n × R n → R n such that C(x, y) = 0 n ⇔ 0 n x ⊥ y 0 n . ( 45
)
With the help of such a function C one can reformulate LEM( f, g) as a system

⎧ ⎨ ⎩ C(x, y) = 0 n f (x, y, λ) = 0 n g(x) = 0 m (46)
of 2n + m equations and 2n + p unknown variables. The most popular choice of C is the Fischer-Burmeister complementarity function

C FB (x, y) = x + y -x [2] + x [2] [1/2] ,
where the square root operation [ • ] [1/2] is understood in the Jordan algebra sense. Another choice of interest is the natural complementarity function

C nat (x, y) = x -(x -y) + .
A common feature of C FB and C nat is that both functions are locally Lipschitz and semismooth; see [START_REF] Chen | Analysis of nonsmooth vector-valued functions associated with secondorder cones[END_REF] for definitions and proofs. So, with any of the above choices of C, the system (46) has the general structure [START_REF] Qi | A nonsmooth version of Newton's method[END_REF] with : Z → W standing for a locally Lipschitz semismooth function between Euclidean spaces of possible different dimensions. The Newton method and the NFA take now a more general form, to wit:

-Square case. If dim(W) = dim(Z), then ( 22) is a square system. One solves it with the Semismooth Newton Method (SNM)

z τ +1 = z τ -M -1 τ (z τ ),
where M τ : Z → W is a linear map taken from the Clarke generalized differential ∂ (z τ ) of at z τ (cf. [START_REF] Qi | A nonsmooth version of Newton's method[END_REF]). -Underdetermined case. If dim(W) < dim(Z), then ( 22) is underdetermined. In such a situation one uses the Semismooth Normal Flow Algorithm (SNFA)

z τ +1 = z τ -M † τ (z τ ),
where M τ is chosen arbitrarily from ∂ (z τ ).

By way of example, consider again the Lorentz eigenvalue problem as formulated in (12). The system [START_REF] Qi | A nonsmooth version of Newton's method[END_REF] becomes

⎧ ⎨ ⎩ C(x, y) = 0 n (A -λI n )x -y = 0 n e n , x -1 = 0.
This is system of 2n + 1 equations in the same number of unknown variables, so one can solve it by using the SNM. This idea has been suggested in [1, Section 4.1].

Equilibrium models constrained by elliptic cones

The Lorentz equilibrium model (4) can be seen as a particular instance of an equilibrium model

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ K x ⊥ y ∈ K + f (x, y, λ) = 0 n g(x) = 0 m (47)
in which the complementarity conditions are expressed in terms of a closed convex cone K . The notation K + refers to the dual cone of K . The squaring technique developed in this work can be extended to the case in which

K = E(Q) := {(ω, t) ∈ R n : ω, Qω ≤ t} (48)
is an elliptic cone associated to a positive definite matrix Q of order n -1. Elliptic cones have been studied under different angles in a number of references. Complementarity conditions relative to an elliptic cone are considered in [START_REF] Andreani | Box-constrained minimization reformulations of complementarity problems in second-order cones[END_REF]. An elliptic cone can be seen as the image of the Lorentz cone under an invertible linear transformation. Indeed, one has [2] : u ∈ R n }, where S Q : R n → R n is the invertible linear transformation whose matrix representation is

E(Q) = {S Q x : x 0 n } = {S Q u
S Q = ⎡ ⎢ ⎣ Q -1/2 0 n-1 0 T n-1 1 ⎤ ⎥ ⎦ .
One can easily check that the dual of ( 48) is given by [2] : v ∈ R n }.

(E(Q)) + = E(Q -1 ) = {S Q -1 y : y 0 n } = {S Q -1 v
Hence, the equilibrium model (47) can be converted into a system of nonlinear equations, namely ⎧ ⎪ ⎨ ⎪ ⎩ S Q u [2] , S Q -1 v [2] = 0 f (S Q u [2] , S Q -1 v [2] , λ) = 0 n g(S Q u [2] ) = 0 m .

Since S Q is symmetric and S Q -1 = (S Q ) -1 , it is clear that S Q u [2] , S Q -1 v [2] = u [2] , v [2] .

In other words, the seemingly more general model ( 47) is nothing but LEM( f Q , g Q ) with

f Q (x, y, λ) = f (S Q x, S Q -1 y, λ) g Q (x, y, λ) = g(S Q x).
In short, the whole theory of Lorentz equilibrium models extends to the elliptic case.

Final comments

The leading motivation of this work has been to exploit the possibility of getting rid of the nonnegativity constraint x 0 n by expressing the vector x ∈ R n as the "square" of another vector u ∈ R n . Of course, the squaring technique is applied also to the nonnegativity constraint y 0 n imposed on the dual variables. As an alternative to the squaring technique one may consider the use of a complementarity function as in (45). Some preliminary numerical experiments with the Lorentz eigenvalue problem suggest that the combination of the complementarity function technique and the SNM performs slightly better than the combination of the squaring technique and the classical Newton method. Having said this, we would like to point out that the squaring technique has at least two advantages:

• First of all, the squaring technique is well suited to discriminate between different types of solutions to a LEM (boundary type, interior type, et cetera). • Secondly, the squaring technique can be easily extended to the case in which the nonnegative constraints are expressed in terms of an elliptic cone. One could even consider nonnegative constraints expressed in terms of a Cartesian product of several elliptic cones. LEMs with such sort of nonnegative constraints arise for instance in the analysis of elastic systems with unilateral contact and friction.

Corollary 4 . 4

 44 Let the assumption (32) be in force. Then (a) (x, y, λ) is a central solution to LEM( f, g) if and only if (α) occurs. (b) (x, y, λ) is a boundary type solution to LEM( f, g) if and only if (β) occurs. (c) (x, y, λ) is an eccentric solution to LEM( f, g) if and only if (γ ) occurs with ω / ∈ {0, 1}.

Table 1

 1 Solving a Lorentz eigenvalue problem by using Newton's method Percentages of success are estimated with a sample of 10 3 random matrices A with Gaussian distribution

	n	Newton's method applied to (27)	Newton's method applied to (28)
		NIP = 1	NIP = 10	NIP = 10 2	NIP = 1	NIP = 10	NIP = 10 2
		(%)	(%)	(%)	(%)	(%)	(%)
	5	7 . 4	5 8 .0	9 9 .9	4 5 .7	100	100
	10	0.1	5.1	3 3 .5	2 0 .0	8 8 .1	100
	15	0	0.1	1 .0	8 .9	6 3 .2	100
	20	0	0	0	6.0	4 0 .6	9 9 .6

Table 2

 2

	Lorentz eigenvalues of the matrix (29)	Type	λ Primal vector			Dual vector
			x 1	x 2	x 3	x 4	y 1	y 2	y 3	y 4
		Boundary 2

Table 3

 3 Finding the Lorentz eigenvalues of the matrix (29) by applying Newton's method to the companion system[START_REF] Walker | Newton-like methods for underdetermined systems[END_REF]

	Ill-conditioning	Divergence	Central	Eccentric	Boundary
					λ = 2	λ = 3	λ = 5	λ = 7
	55.0 %	0.0 %	3.6 %	2.7 %	21.6 %	15.4 %	0 %	1.7 %

Table 4

 4

	Finding boundary type solutions to the Lorentz quadratic	n	Newton's method applied to (42)
	eigenvalue problem		NIP = 1 (%)	NIP = 10 (%)	NIP = 10 2 (%)
		5	6 3 .0	9 7 .4	9 8 .1
		10	32.4	8 3 .4	9 6 .7
	Percentages of success are estimated with a sample of 10 3 random quadratic pencils Q ∈ T n	15 20	25.3 19.5	7 2 .5 6 5 .7	9 6 .1 9 5 .8
	then one can reformulate (42) as a classical quadratic eigenvalue problem

1/2 , Pedro Gajardo was partially supported by Chilean Fondecyt Grant No. 1120239 and by "Programa de Financiamiento Basal" from the Center of Mathematical Modeling, Universidad de Chile.