Jean Charles Gilbert

Numerical methods for large-scale minimization *

Keywords: Automatic differentiation, conjugate gradient, large-scale optimization, limited memory quasi-Newton method, partitioned quasi-Newton method, unconstrained optimization

This paper recalls some basic tools for unconstrained minimization problems. It includes a comparison of performance of automatic differentiation codes and methods when they are used to compute a gradient. Some algorithms for minimizing a function with a large number of variables are reviewed and the impact of automatic differentiation on them is discussed.

Basic notions

Let us equip R n with an inner product •, • and its associated norm | • |. We shall use the following matrix norm on R n : ||A|| := sup{|Au| : |u| = 1}.

Given a smooth function f : Ω → R, defined on an open set Ω of R n , the unconstrained minimization problem consists in finding a point x * ∈ Ω, such that

f (x * ) f (x), ∀ x ∈ N (x * ), (1.1)
where N (x * ) denotes some neighborhood of x * . A point x * satisfying the above inequalities is called a local minimizer of f . If (1.1) holds for N (x * ) = Ω, then x * is called a global minimizer of f on Ω. Finding a global minimizer is a hard task and for this subject, we refer the reader to the collections of papers [START_REF]Towards Global Optimization[END_REF][START_REF]Towards Global Optimization[END_REF][START_REF]Handbooks in Operations Research and Management Science 1: Optimization[END_REF]. The algorithms considered in this paper are only able to find local minimizers.

If x * is a local minimizer, it satisfies the first order optimality conditions:

f ′ (x * ) • h = 0, ∀ h ∈ R n ,
where f ′ (x * ) is the Fréchet derivative of f at x * . These conditions can be rewritten as follows: g(x * ) = 0, (

where g(x) := ∇f (x) denotes the gradient of f at x. It is important to keep in mind that the gradient depends on the inner product chosen on R n . It is indeed defined by means of the Riesz theorem:

∃ g(x) ∈ R n , unique : g(x), h = f ′ (x) • h, ∀ h ∈ R n .
If the Euclidean inner product is used, u, v = n i=1 u i v i , then g i (x) is the i-th partial derivative of f at x. A point x * satisfying (1.2) is called a stationary point of f . A minimizer is a stationary point and the converse is true for convex function.

Here, we have to be still more modest as to what the algorithms described in this paper can do. In general, they can only find stationary points. This is due to the fact that they will use only the value of f and the value of its gradient. If for some reasons, a local minimizer has to be found (and not only a stationary point) and if there are many stationary points that are not minimizers, something more has to be done. A good idea is to use second order information to locate directions of negative curvature. Along this line, see the papers by Gill and Murray [START_REF] Gill | Newton-type methods for unconstrained and linearly constrained optimization[END_REF], McCormick [START_REF] Mccormick | A modification of Armijo's step-size rule for negative curvature[END_REF], Moré and Sorensen [START_REF] Moré | On the use of directions of negative curvature in a modified Newton method[END_REF] and Goldfarb [START_REF] Goldfarb | Curvilinear path steplength algorithms for minimization which use directions of negative curvature[END_REF].

An important class of methods generating sequences of approximations of a solution is formed by the descent-direction methods. They proceed as follows. It is assumed that some approximation x k is known at the beginning of iteration k, k 1 (x 1 has to be given at the beginning of the run). Then, an iteration consists of calculating a new approximation x k+1 that is better than the previous one, in the sense that f (x k+1 ) < f (x k ).

(1.3)

To obtain this descent property, the first stage consists in choosing a descent direction at x k , i.e., a direction d k verifying

f ′ (x k ) • d k < 0. (1.4)
The second stage consists in selecting a positive step-size ρ k along d k and taking

x k+1 = x k + ρ k d k
such that x k+1 satisfies, at least, the descent property (1.3). This stage is called the line-search. Let us summarize this method as follows.

The descent-direction algorithm: 1. choose x 1 ∈ Ω; 2.

k := 1; 3.

until convergence do { 3.1.

choose a descent direction d k ; 3.2.

select a positive step-size ρ k along d k such that (1.3) holds; 3.3.

x k+1 := x k + ρ k d k ; 3. 4. k := k + 1

}

The way of calculating the descent direction characterizes a descent-direction method. For example, when g k := g(x k ) = 0, an obvious descent direction is:

d k = -g k .
(1.5)

The method that uses this descent direction at each iteration is called the gradient method or the steepest descent method. More generally, if H k is a positive definite matrix of order n,

d k = -H k g k (1.6)
will be a descent direction. Many practical algorithms for minimizing a function have this form of descent direction: the Newton method (near a solution), the bfgs method (or more generally, algorithms in the family of quasi-Newton methods) and even conjugate gradient methods, as we shall see in Section 5.

The descent condition (1.4) implies that for sufficiently small positive step-size ρ, we have f (x k + ρd k ) < f (x k ). Therefore, to satisfy (1.3), a step-size selection procedure could simply choose a small value for ρ. This is not satisfactory, however, because examples can be given for which taking arbitrarily small positive step-sizes forces a descent algorithm to generate a sequence {x k } k 1 that converges to a nonstationary point (false convergence). This phenomenon is dangerous and hopefully, it can be avoided in several ways. We mention two of them.

A first idea is to minimize the function ξ k (ρ) := f (x k +ρd k ) for positive ρ. We shall say that a step-size selection procedure realizes an optimal line-search if it selects an optimal step-size, that is to say a positive step-size ρ k such that

x k+1 = x k + ρ k d k satisfies f (x k+1 ) f (x k + ρk d k ), (1.7)
where ρk is the smallest positive stationary point of the function ξ k . A step-size that is a stationary point of ξ k will be said to be exact and is obtained by a exact line-search.

Actually, the optimal and exact line-search strategies have two drawbacks. First, the stepsize ρk may not exist at some iteration although problem (1.1) has a solution. Secondly, finding a stationary point of ξ k is a nonlinear minimization problem in itself and this cannot be solved exactly, except in particular cases (for instance, when f is quadratic). Therefore, although this procedure is theoretically suitable (see Theorem 1.1 below), it is useless in practice.

A better choice is to satisfy some inequalities that will guarantee convergence if the descent directions are reasonable. A first inequality assures sufficient decrease of the objective function: for

ω 1 ∈ ] 0, 1[ choose ρ k such that f (x k+1 ) f (x k ) + ρ k ω 1 g k , d k . (1.8)
This inequality is always satisfied for small positive ρ k . In order to avoid too small stepsizes one have to satisfy another inequality like the following:

g(x k+1 ), d k ω 2 g k , d k (1.9)
with ω 2 ∈ ] ω 1 , 1[. Conditions (1.8) and (1.9) are called the Wolfe conditions and a point x k+1 satisfying them is called a Wolfe point. It can be proved (see Wolfe [START_REF] Wolfe | Convergence conditions for ascent methods[END_REF]) that if f is bounded below and if 0 < ω 1 < ω 2 < 1, it is always possible to find a Wolfe point. Typical values for the parameters are ω 1 = 10 -4 . . . 10 -2 and ω 2 = 0.5 . . . 0.99. In any case, it is healthy to take ω 1 < 0.5, so that the step-size minimizing the quadratic approximation of ξ k will be accepted by the inequality (1.8). Calculating a Wolfe point requires a specific algorithm. The one proposed by Lemaréchal [START_REF]A view of line-searches[END_REF] has been proved to find a Wolfe point in a finite number of trials. For some algorithms, we shall need a step-size satisfying the strong Wolfe conditions, that is to say, satisfying (1.8) and an inequality stronger than (1.9), namely

| g(x k+1 ), d k | ω 2 | g k , d k |. (1.10)
An algorithm for finding a strong Wolfe point can be found in [START_REF] Al-Baali | An efficient line search for nonlinear least squares[END_REF] and [START_REF] Moré | Line search algorithms with guaranteed sufficient decrease[END_REF].

The preceding two line-search techniques allow us to obtain the following result that can be derived from the papers by Zoutendijk [START_REF] Zoutendijk | Nonlinear programming, computational methods[END_REF] and Wolfe [START_REF] Wolfe | Convergence conditions for ascent methods[END_REF][START_REF] Wolfe | Convergence conditions for ascent methods II: some corrections[END_REF]. Before stating the theorem, let us specify some notation. We shall denote by θ k the angle between d k and -g k , that is

cos θ k := - g k , d k |g k | |d k | .
The notation "lim inf ǫ k " means "lim k→∞ inf 1 i k ǫ i ".

Theorem 1.1

Let f be a differentiable function defined on the open set Ω, with a Lipschitz continuous derivative. Suppose that a descent-direction method generates a sequence of points {x k } k 1 in Ω and a sequence of directions {d k } k 1 verifying (1.4). Suppose also that at each iteration the line-search procedure finds either an optimal point satisfying (1.7) or a Wolfe point satisfying (1.8) and (1.9). Suppose finally that {f (x k )} k 1 is bounded below. Then, the following conclusions hold:

(i) k 1 |g k | 2 cos 2 θ k is convergent; (ii) |g k | cos θ k → 0; (iii) if k 1 cos 2 θ k is divergent, then lim inf |g k | = 0.
The theorem says nothing on the convergence of the sequence of points {x k } k 1 . Only the sequence of gradients is concerned, although it is not claimed that g k converges to zero. Actually, Theorem 1.1 expresses the contribution of the line-search procedure to force convergence. Now, to make its conclusions useful, some more information concerning the way the descent directions d k are chosen is needed. For instance, if the angle θ k remains bounded away from ±π/2, then cos θ k remains bounded away from 0 and the conclusion (ii) of Theorem 1.1 implies that the sequence of gradients converges to zero and therefore that any limit point of {x k } k 1 is stationary. This is the case for the gradient method, because θ k = 0, ∀ k 1. This is also the case for the method (1.6), if the condition numbers κ(H k ) of the positive definite matrices H k remain bounded; that is, if there exists a positive constant κ such that

κ(H k ) := ||H k || ||H -1 k || κ, ∀ k 1.
Indeed, we have cos θ k κ(H k ) -1 .

These remarks indicate that to be sure to have a convergent algorithm, it is enough to use method (1.6) with constant matrices H k (the identity matrix, for instance) or with matrices having bounded condition number. Our task is not finished, however. Indeed, such directions can be very poor descent directions in the sense that the decrease of the objective function may be small along them. Therefore, in the subsequent sections, we shall play the dangerous game of trying to define "good" descent directions by taking information from the function with the risk of violating the conditions of convergence stated in Theorem 1.1. For some of the methods below, it is still not clear whether they can fail to converge when they aim at minimizing general functions.

Automatic differentiation implementations

As we have seen, the gradient ∇f (x) of the function f to minimize gives precious information on the behavior of f around x. Efficient algorithms always use this information.

Sometimes, the cost of the computation of the gradient has been invoked to motivate users' preference for methods only using the values of the function or to use divided differences to evaluate the gradient. The latter choice is very expensive, since it usually requires n more evaluations of the function. It has also low accuracy. These two reasons should be sufficiently deterrent factors for not using it.

Fortunately, for about a decade, the technique known as the reverse mode of automatic differentiation has provided a way of computing the gradient exactly (up to the precision of the arithmetic of the computer) at a relative cost independent of the number of variables, meaning that

κ ≡ T (f, ∇f ) T (f ) C. (2.1)
Here, T (f ) and T (f, ∇f ) measure the time to compute f and the pair (f, ∇f ), respectively, and C is a constant independent of n. This constant is small, between 3 and 5, depending on the operations taken in consideration [START_REF] Baur | The complexity of partial derivatives[END_REF][START_REF]Mathematical Programming[END_REF]23]. The method requires that the function to differentiate be described by a computer program. It is also a computer program that will represent the gradient function. This technique is opposed to the direct mode of automatic differentiation, which is, in spirit and cost, very similar to the divided difference gradient evaluation. For review papers on the subject, see [23,[START_REF] Griewank | The chain rule revisited in scientific computing[END_REF][START_REF] Gilbert | La différentiation automatique de fonctions représentées par des programmes[END_REF].

The reverse mode goes back at least to Linnainmaa [START_REF] Linnainmaa | Taylor expansion of the accumulated rounding error[END_REF], who introduced it for error propagation analysis, and to Speelpenning [START_REF] Speelpenning | Compiling fast partial derivatives of functions given by algorithms[END_REF], whose interest was precisely the automatic differentiation of functions represented by program. For short (but with a little inaccuracy), the method can be seen as an adjoint method to compute the gradient of f , the lines of the program evaluating f being considered as constraints. More specifically, the method basically deals with programs computing f by a sequence of instructions:

x µ k := ϕ k (x P k ), for k = 1, . . . , K, f := x N .
(2.2)

We have denoted by x 1 , . . . , x n the independent variables, i.e., the variables defining the function f (x) = f (x 1 , . . . , x n ). The other intermediate N -n variables appearing in the program are denoted by x n+1 , . . . , x N . Instruction (2.2) modifies the variable x µ k , for some 1 µ k N , by means of an intermediate function ϕ k , which depends on x P k ≡ {x j : j ∈ P k }, P k being some subset of {1, . . . , N }. The gradient can be evaluated by introducing variables p k , 1 k N (p k is called the dual -or adjoint -variable to x k ), initialized to

p k := 0 for k = 1, . . . , N -1 1 for k = N, (2.3)
and evaluated backwards by the following adjoint program:

p j := p j + ∂ϕ k ∂x j p µ k , ∀ j ∈ P k \{µ k } p µ k := ∂ϕ k ∂x µ k p µ k        , for k = K, . . . , 1. (2.4)
Then, the gradient is obtained by

∂f ∂x i = p i , for i = 1, . . . , n. (2.5)
The partial derivatives of ϕ k used in (2.4) have to be evaluated during the "forward sweep" (2.2), and stored. This delay between the evaluation of ∂ϕ k /∂x j and its use in (2.4) gives rise to different implementations of the reverse mode.

In the first implementation, the calculations made in (2.2) are stored in a data structure representing the computational graph of the program. The size of this graph is, as a result, proportional to the number of elementary operations executed to obtain f (x). In other words, the size is proportional to the computation time T (f ). This is not a very desirable feature and constitutes the main drawback of the approach, limiting its use to small or medium size problems. On this point, some promising remedies have been proposed to control the memory requirement. The price to pay is an increase in computation time: see Griewank [START_REF] Griewank | Achieving logarithmic growth of temporal and spatial complexity in reverse automatic differentiation[END_REF], for instance. On the other hand, this approach has a number of nice features: the implementation can be easily done by operator overloading [START_REF] Griewank | ADOL-C, a package for the automatic differentiation of algorithms written in C/C++[END_REF], the technique can face programs with intricated structures (recursive subroutines, goto's, . . . ) and higher order derivatives can also be obtained, using a hardly more complicated algorithm than (2.3)-(2.5). Representatives of this approach are jakef of Speelpenning [START_REF] Speelpenning | Compiling fast partial derivatives of functions given by algorithms[END_REF] and Hillstrom [START_REF] Hillstrom | User guide for JAKEF[END_REF], padre2 of Iri and Kubota [START_REF] Iri | PADRE 2, version 1 -User's manual[END_REF] and adol-c of Griewank et al. [START_REF] Griewank | ADOL-C, a package for the automatic differentiation of algorithms written in C/C++[END_REF].

The other approach is the adjoint code implementation, as used in the meteorology world since its introduction by Ph. Courtier and O. Talagrand -see [START_REF] Talagrand | The use of adjoint equations in numerical modelling of the atmospheric circulation[END_REF]. The idea is to write a program, following the rules given by formula (2.4). Here, the structure of the calculation is stored in a program, not in a computational graph. The main advantage of this approach is its sobriety in memory requirement. Indeed, the only pieces to memorize are the values of the variables entering non linearly in the intermediate functions ϕ k , since those entering linearly no longer appear in (2.4). From these values, the adjoint code will evaluate the partials ∂ϕ k /∂x j during the backward sweep. Interestingly enough, this implies that, for the differentiation of linear functions, nothing has to be stored during the direct calculation (2.2). This is in complete opposition with the computational graph approach. The approach also offers the possibility to recalculate some variables in the adjoint code, instead of memorizing them. Unfortunately, there is presently no automatic generation of adjoint codes, so that they have to be written by hand.

We now give the results of some tests made in [START_REF] Gilbert | La différentiation automatique de fonctions représentées par des programmes[END_REF], comparing the efficiency of the automatic differentiators jakef, padre2, adol-c with the divided difference approach (dd), as well as the hand-written adjoint code approach (meteo). All tests have been made on a SPARCstation 1. In Tables 2.1 and 2.2, label "Kb" refers to the additional memory requirement in Kbytes, label "sec" gives the CPU time (in seconds) for the computation of f and/or ∇f and label "κ" is the quotient in (2.1). The additional memory spaces used by jakef, padre2 and adol-c are approximately identical: it is mainly used to store the computational graph. This is much more than what is necessary in the original program or for the adjoint code method (meteo). For the latter, we chose to store only one additional vector and to recompute in the backward phase (2.4) all the other variables appearing nonlinearly in the code.

The computing times were obtained by averaging on 10 trials (except for dd!). The CPU time used by the divided difference method is very important: it exceeds n times the one of f , since for some components, the increment dx i is readapted and the difference recomputed. The times realized by jakef and padre2 are similar and a bit less than the one of adol-c (for the latter, "f with graph" is the time to build the computational graph, "f without graph" and "∇f without graph" are the times to compute f and ∇f using the graph previously built). Note that with adol-c, the code is translated from fortran to c++ (using f2c from [START_REF] Feldman | A Fortran-to-C converter[END_REF]) and is executed in part in double precision, while only single precision arithmetic is used with the other codes. This may explain the difference in computing time between adol-c and the pair jakef-padre2. The time realized by meteo is really amazing and encouraging. It is less than half the one used by the original code and less than 10 times the one used by jakef or padre2, despite our choice to recalculate most of the intermediate variables. Note finally that only the method meteo computes the pair (f , ∇f ) at a relative cost not exceeding 5, as expected by the theory.

Another example will highlight the fact that the increase in memory space required by the computational graph implementation, due to the increase in the number of elementary operations for the computation of f , results in a deterioration of the performance of the reverse mode. We consider the following subroutine loops, which computes f (x) = (ntimes) 3 x 2 in a laborious way: three nested do-loops ranging from 1 to ntimes are used to obtain (ntimes) 3 and x 2 is calculated in the innermost do-loop.

subroutine loops (n,ntimes,x,f) integer n,ntimes real x(n),f integer i,j,k,l real r f=0. do Results are given in Table 2.2 for n = 10 and ntimes ranging from 5 to 25. We observe that for jakef and padre2, the memory requirement (in Kb) increases as (ntimes) 3 . We also see that, for jakef and padre2, the quotient κ deteriorates as ntimes increases. This has to be attributed to the computational graph implementation of the reverse mode.

On the other hand, it can be shown that, for meteo (at least for our way of writing the adjoint code), κ is less than or equal to 3 when ntimes and n are large. This is what we observe in Table 2.2. These numerical results are very encouraging for the adjoint code approach of the reverse mode of automatic differentiation.

Suppose now that we are given a vector-valued function F : R n → R n . As we have seen, the reverse mode of automatic differentiation will compute cheaply the gradient

f
∇F i (x) of a component of F , i.e.
, a line of the Jacobian matrix ∇F (x). On the other hand, the direct mode can computes cheaply the directional derivative F ′ (x) • e j , i.e., a column of ∇F (x). Therefore, it makes sense to ask whether there exists a method, DIAG say, that would compute the n elements of the diagonal of the Jacobian at approximately the same price as a line or a column. The question is important since the diagonal can be a good preconditioner for many numerical algorithms. Apparently, there is no definite answer to this question but it can be shown ( [START_REF] Gilbert | La différentiation automatique de fonctions représentées par des programmes[END_REF]) that if such an algorithm existed, it would be possible to compute the product of 2 matrices in O(n 2 ) operations. As the best presently known method computes this product in O(n 2.50 ) operations, our result puts some doubt on the existence of DIAG. It also suggests that it is more valuable to know the diagonal of the Jacobian rather than one of its lines or columns.

Conjugate gradient methods without optimal linesearches

The conjugate gradient (cg) method was introduced to solve linear systems of equations [START_REF] Hestenes | Methods of conjugate gradients for solving linear systems[END_REF] and later to minimize strongly convex quadratic functions [START_REF] Fletcher | Function minimization by conjugate gradients[END_REF]. The basic idea is to generate a set of conjugate directions, i.e., directions that are orthogonal with respect to the inner product induced by the Hessian matrix A of f . As a result, f is minimized on the larger and larger affine subspace spanned by the successive directions, the Krylov subspace. When the first direction is -g 1 , the conjugacy of the directions can be obtained by the following simple formula:

d k = -g k + β HS k d k-1 ,
where β HS k is the Hestenes-Stiefel beta (see [START_REF] Hestenes | Methods of conjugate gradients for solving linear systems[END_REF]):

β HS k := y k-1 , g k y k-1 , d k-1 .
We note y k := g k+1 -g k and s k := x k+1 -x k . To have conjugacy, it is necessary to do optimal line-searches, which does not pose any difficulty in case of quadratic functions.

The linear CG algorithm: 1. choose a point x 1 ∈ Ω; 2.

d 1 := -g 1 ; k := 1; 3. repeat { 3.1. ρ k := -g k , d k / Ad k , d k ;
/ * the optimal step-size * / 3.2.

x k+1 :=

x k + ρ k d k ; 3.3.
if convergence then stop; 3.4.

d k+1 := -g k+1 + β HS k+1 d k ; 3.5. k := k + 1 3. }
Usually, real-life large-scale models are not quadratic. The question is therefore to know whether, for minimizing a general nonlinear function, it still makes sense to use a cg-like method defined for some scalar β k by

d 1 = -g 1 d k = -g k + β k d k-1 , for k 2.
(3.1)

Note that for non quadratic functions, it is still possible to define a local notion of conjugacy by using the average Hessian of f :

Bk = 1 0 ∇ 2 f (x k + ts k ) dt.
Indeed, because y k = Bk s k , equation (3.1) with β k = β HS k implies that Bk d k , d k+1 = 0 for k 1. Therefore, d k and d k+1 are conjugate with respect to Bk . However, it does not seem that this concept is very helpful in the non quadratic case and therefore, on this basis only, it is hard to claim that the HS beta is better than any other beta.

Another point is to know whether the directions defined by (3.1) are downhill: do we have g k , d k negative? As

g k , d k = -|g k | 2 + β k d k-1 , g k ,
an obvious way of getting descent directions is to do exact line-searches: d k-1 , g k = 0. We have said that such strategy is not realistic. Hopefully, it can be avoided. The following theorem (see [START_REF] Gilbert | Global convergence properties of conjugate gradient methods for optimization[END_REF]Theorem 3.2]) is an extension of a result of Al-Baali [START_REF] Al-Baali | Descent property and global convergence of the Fletcher-Reeves methods with inexact line search[END_REF]. It shows that provided |β k | is not too large, a line-search procedure realizing the strong Wolfe conditions automatically ensures descent of the direction d k . Furthermore, global convergence can be proved. Note that this result strongly relies on the fact that the first search direction is -g 1 . In this theorem, the Fletcher-Reeves beta [START_REF] Fletcher | Function minimization by conjugate gradients[END_REF] is used to control the magnitude of β k . It is defined by (see [START_REF] Fletcher | Function minimization by conjugate gradients[END_REF]):

β FR k := |g k | 2 |g k-1 | 2 .
For quadratic functions and optimal line-searches, β HS k = β FR k .

Theorem 3.1 Let f be a twice continuously differentiable function defined on the open set Ω = R n and let x 1 be a point in Ω. Suppose that the level set { x ∈ Ω : f (x) f (x 1 )} is bounded. If the nonlinear cg method (3.1) with |β k | β FR k starts from x 1 and uses the strong Wolfe line-search with 0 < ω 1 < ω 2 < 1/2, then the directions d k are downhill and lim inf |g k | = 0. This result could let think that the cg method (3.1) with β k = β FR k and the strong Wolfe line-search is a good algorithm. Unfortunately, it is not the case. The point is that this theorem only ensures the global convergence of the method, i.e., the convergence from any starting point, but not its fast convergence. In fact, numerical experiments with the Fletcher-Reeves [START_REF] Fletcher | Function minimization by conjugate gradients[END_REF] formula have shown that the method can converge very slowly. A clarification of this phenomenon has been given by Powell [START_REF] Powell | Restart procedures of the conjugate gradient method[END_REF][START_REF] Powell | Convergence properties of algorithms for nonlinear optimization[END_REF] who gave a reason why the Polak-Ribière beta, defined by (see [START_REF] Polak | Note sur la convergence de méthodes de directions conjuguées[END_REF]):

β PR k := y k-1 , g k |g k-1 | 2 ,
should by used instead of β FR k (for exact line-searches, β HS k = β PR k ). Numerical experiments have confirmed the superiority of the PR beta over the FR beta.

The misfortune is that the PR method is not as safe as the FR method! If it can be proved that the method is convergent for convex functions (see [START_REF] Polak | Note sur la convergence de méthodes de directions conjuguées[END_REF]), Powell [55] has given an example of a nonconvex function where the PR method diverges: the generated sequence cycles and each of its limit points is nonstationary! Two globally convergent remedies to this situation have been proposed in [START_REF] Gilbert | Global convergence properties of conjugate gradient methods for optimization[END_REF]. The first one is inspired by Theorem 3.1. It consists of using the PR beta only when it lies in the interval [-β FR k , +β FR k ]. Specifically,

Method PR|FR : β k =    -β FR k if β PR k < -β FR k β PR k if |β PR k | β FR k β FR k if β PR k > β FR k .
This strategy avoids one of the main disadvantages of the FR method, as discussed in [START_REF] Gilbert | Global convergence properties of conjugate gradient methods for optimization[END_REF].

Another globally convergent remedy consists of restarting the PR method in the opposite direction of the gradient, each time β PR k is negative. Synthetically, this corresponds to the following choice of β k :

Method PR + : β k = max{0, β PR k }.
For this method, a specific line-search assuring descent of the directions has to be implemented.

FR

PR|FR PR PR + > 4.07 1.55 1.02 1.00 [START_REF] Gilbert | Global convergence properties of conjugate gradient methods for optimization[END_REF] on a collection of 26 largescale problems with a number of variables between 100 and 10000. The table gives the total number of function/gradient calls used for all the test-problems divided by the 18792 calls required by method PR + . The sign > means that the method was stopped on some test-problems because the number of function/gradient calls exceeds 9999. The actual number should therefore be greater than the one given in the table. Methods PR and PR + appear to be the best one, while the FR method is on average more than four times slower than method PR + , but this behavior of the FR method is well known. [START_REF] Baur | The complexity of partial derivatives[END_REF] The Newton and truncated Newton algorithms

The Newton method

The Newton method for solving nonlinear equations (here ∇f (x) = 0) is well known. It consists in linearizing the equations and solving them. Here, this gives the following linear system to solve in d:

∇ 2 f (x k ) d + g k = 0,
where ∇ 2 f (x k ) denotes the Hessian matrix of f at x k . This matrix also depends on the inner product. If ∇ 2 f (x k ) is nonsingular, the solution is

d N k = -∇ 2 f (x k ) -1 g k . (4.1)
The direction has the form (1.6) and is a descent direction when ∇ 2 f (x k ) is positive definite, which is the case when x k is close to "normal" solution.

An attractive feature of Newton's method is its rate of convergence. If the function f is sufficiently smooth and if x 1 is close enough to x * , then the method generates a quadratically convergent sequence {x k } k 1 , that is

||x k+1 -x * || C ||x k -x * || 2 , ∀ k 1,
for some positive constant C. This property implies that if x k is close enough to x * , its number of correct significant digits grows very quickly. Quadratic convergence is typical of algorithms using second order information, i.e., second derivatives.

If convergence occurs with the Newton method when the initial point x 1 is close to a solution, divergence is often observed if the guess is not good enough. This drawback can be tempered by line-search techniques as described in Section 1 or by a trust region approach (see the excellent survey by Moré [START_REF] Moré | Recent developments in algorithms and software for trust region methods[END_REF]). Other difficulties may occur. For example, when the Hessian is not positive definite, the Newton direction may no longer be downhill. One remedy to this situation is to modify the Hessian in equation (4.1) to obtain a positive definite matrix: one of the most effective way of doing this is due to Gill and Murray [START_REF] Gill | Newton-type methods for unconstrained and linearly constrained optimization[END_REF]. Another remedy is to solve the linear system (4.1) within a trust region. A last possibility, which can be extended to large-scale problems, is the truncated Newton method.

The truncated Newton method

Newton's method is difficult to use for large-scale problems because the storage of the full Hessian requires O(n 2 ) memory locations and the solution of Newton's equation (4.1) demands O(n 3 ) operations. The truncated Newton method attempts to overcome these two difficuties.

The basic idea is to solve only partially the linear system (4.1), using an iterative procedure (formed of inner iterations, as opposed to the outer iterations of Newton's method). For instance, the cg method can be used as long as the possible non definite positiveness of the Hessian is not detected [START_REF] Dembo | Truncated-Newton algorithms for large-scale unconstrained optimization[END_REF]. Then, it can be shown that the resulting direction is a descent direction of f (and a negative curvature direction if the process has been interrupted). Another possibility would be to use a Lanczos tridiagonalization process [START_REF] O'leary | A discrete Newton algorithm for minimizing a function of many variables[END_REF].

These two methods are suitable for large-scale problems, since the only information used from the Hessian is its action on a vector v. There are different possibilities to evaluate the product Hessian * vector (see, for instance, the discussion in [START_REF] Overton | A powerful truncated Newton method for potential energy minimization[END_REF]), but when the Hessian is dense, it is generally approximated by divided differences:

∇ 2 f (x k ) v ≃ ∇f (x k + ǫv) -∇f (x k ) ǫ ,
where ǫ is a small positive number. This formula costs one gradient evaluation. Now, automatic differentiation offers the possibility to obtain this product exactly and cheaply. Indeed, ∇ 2 f (x) v is the directional derivative of x → ∇f (x) in the direction v. If the gradient has been represented by an adjoint code (reverse mode), the product ∇ 2 f (x) v can be evaluated by the direct mode at a cost similar to one gradient evaluation.

5 The quasi-Newton, LBFGS and partitioned quasi-Newton methods

Quasi-Newton methods

The Newton method requires the calculation of the Hessian matrix of f , which can be in certain cases, either too expensive or impossible. This is precisely the domain where the quasi-Newton (qN) methods can be useful. See the review paper by Dennis and Moré [START_REF] Dennis | Quasi-Newton methods, motivation and theory[END_REF].

The descent direction of a qN method has the form (1.6), where H k is a matrix that is updated at each iteration by some formula. Hence, the method generates two sequences: a sequence of points {x k } k 1 and a sequence of matrices {H k } k 1 . At the beginning of the run, the algorithm is given a point x 1 and a matrix H 1 . Often, H 1 is the identity matrix (we shall give better choices later) or a diagonal matrix, like the diagonal of the Hessian ∇ 2 f (x 1 ), when it is available.

The general scheme of a qN method is the following.

The QN algorithm: 1.

choose x 1 ∈ Ω and a positive definite matrix H 1 ; 2.

k := 1; 3.

until convergence do { 3.1.

d k := -H k g k ; 3.2.
select a step-size ρ k by the Wolfe procedure (1.8) and (1.9); 3.3.

x k+1 := x k + ρ k d k ; 3.4. calculate H k+1 ; 3.5. k := k + 1 3. }
The only thing to examine in this algorithm is the way the matrix H k is updated. If we set B k = H -1 k (we shall see how to maintain H k positive definite, hence nonsingular), we see that B k plays the part of ∇ 2 f (x k ) in Newton's method. On the other hand, with s k := x k+1 -x k and y k := g k+1 -g k , we have by Taylor's theorem:

y k = 1 0 ∇ 2 f (x k + ts k ) dt s k .
These two remarks show that it makes sense to impose the matrix B k+1 , which has to be calculated at the end of iteration k, to satisfy:

y k = B k+1 s k .
This equation is called the secant or quasi-Newton equation. If n 2, the equation does not determine B k+1 completely and some other conditions can be imposed. Generally, they make B k+1 keeping information from B k . There is today a wide consensus that one of the best update formulae for minimization is the bfgs update. To write this formula, it is convenient to use the following tensor product of two vectors u and v of R n : u ⊗ v is the rank one matrix defined by (u ⊗ v)h = v, h u, ∀ h ∈ R n . If the Euclidean inner product is used, u ⊗ v is the matrix uv T . See the appendix in [START_REF] Gilbert | Some numerical experiments with variable-storage quasi-Newton algorithms[END_REF] for a derivation of the bfgs formula in this setting, which differs from the usual one by the use of a general inner product. The bfgs formula can be written like this:

B k+1 = B k + y k ⊗ y k y k , s k - (B k s k ) ⊗ (B k s k ) B k s k , s k =: U (B k , y k , s k ). (5.1)
Note that the inverse matrices H k can also be updated by using the inverse bfgs formula:

H k+1 = I - s k ⊗ y k y k , s k H k I - y k ⊗ s k y k , s k + s k ⊗ s k y k , s k =: Ū (H k , y k , s k ). (5.2) Precisely, if H k = B -1 k , B k+1 = U (B k , y k , s k ) and H k+1 = Ū (H k , y k , s k ), then H k+1 = B -1
k+1 . Formulae (5.1) and (5.2) have the nice property of transmitting the positive definiteness of B k to B k+1 (respectively of H k to H k+1 ), if and only if y k , s k is positive. Therefore to have H k positive definite, it will suffice to take H 1 positive definite, to update the matrices by the inverse bfgs formula and to ensure the positivity of y k , s k at each iteration. Note that for strictly convex functions, y k , s k is automatically positive and that for nonconvex functions, the second inequality (1.9) of the Wolfe line-search implies the positivity of y k , s k . This important remark makes the Wolfe step-size selection procedure the most suitable for the qN framework in optimization. It justifies also Statement 3.2 of the qN algorithm.

For problems having small or medium dimension, the bfgs method is very efficient. It is globally convergent for convex problems (see [START_REF] Powell | Some global convergence properties of a variable metric algorithm for minimization without exact line searches[END_REF]) and for general functions, when convergence occurs, it is fast, namely superlinearly convergent. This means that

||x k+1 -x * ||/||x k -x * || → 0.
Although, this is less good than the quadratic rate of convergence of the Newton method, it is generally satisfactory in practice.

The limited memory BFGS method

There are many proposals to adapt qN methods to large-scale problems and many of them fit into the following framework.

Although n may be large, the initial matrix H 1 takes generally little space in memory: it is most commonly a positive multiple of the identity matrix. On the other hand, H k is formed from H 1 and k -1 couples {(y i , s i )} 1 i<k by using formula (5.2). Therefore, to take less storage, one can think of storing these elements instead of H k and computing H k g k by an appropriate algorithm. Of course, when the number of iterations increases, these pieces of information become more and more cumbersome in memory and we must get rid of some of the couples (y i , s i ). A method will be called an m-storage qN method if only m of these couples are used to form H k from a starting matrix H 0 k . In this type of method, the inverse update formula (5.2) is preferable to the direct formula (5.1), because the inversion of B k may be problematic.

The most famous limited memory methods can be seen from the preceding point of view: see [START_REF] Gilbert | Some numerical experiments with variable-storage quasi-Newton algorithms[END_REF] for the details. For example, the algorithm conmin of Shanno and Phua [START_REF] Shanno | Algorithm 500, minimization of unconstrained multivariate functions[END_REF] is a 2-storage qN method (see also [START_REF] Shanno | Conjugate gradient methods with inexact searches[END_REF]). The algorithm of Buckley and LeNir [5] is an mstorage qN method where m varies with the iteration index. Finally, the method proposed by Nocedal [START_REF]Updating quasi-Newton matrices with limited storage[END_REF] is also of this type. We describe it in more detail because we feel it is important.

Nocedal proposed to use always the last m couples (y, s), i.e., {(y i , s i )} k-m i k-1 , to form the matrix H k : going from iteration k to iteration k + 1, the couple (y k-m , s k-m ) is discarded and cyclically replaced by the new couple (y k , s k ). Therefore, H k is always formed from the most recent information and, as a result, it should be more appropriate.

Description of Nocedal's proposal to calculate H

k : 1.
choose a positive definite matrix H 0 k ; 2.

for i := 0 to m -1 do

H i+1 k := Ū (H i k , y k-m+i , s k-m+i ); 3. H k := H m k ;
Of course, this scheme is formal: in practice, the matrices H i k and H k are not stored (they are too large). Instead, just H 0 k and the couples {(y i , s i ) : k -m i k -1} are kept in memory and an elegant algorithm, described in [START_REF] Shanno | Conjugate gradient methods with inexact searches[END_REF], is provided to compute H k g k .

Note that H k g k can be considered as the gradient in g k of the function g k → ϕ(g k ) := 1 2 H k g k , g k . Interestingly enough, this observation has been used in [START_REF] Gilbert | Automatic differentiation and the step computation in the limited memory BFGS method[END_REF] to show that the algorithm computing H k g k in [START_REF]Updating quasi-Newton matrices with limited storage[END_REF] is in fact an adjoint code (in the sense of automatic differentiation, see Section 2) of a program computing ϕ(g k ). This shows the usefulness of automatic differentiation, in particular its reverse mode by adjoint coding, not only to compute the gradient of a numerical function represented by a program, but also to find new algorithms.

The above description of Nocedal's proposal outlines the need of a choice for the starting matrices H 0 k . Contrary to the standard qN methods where this choice has to be made only once at the beginning of a run, here it has to be made at each iteration. Its role is therefore crucial. One of the aims of the work [START_REF] Gilbert | Some numerical experiments with variable-storage quasi-Newton algorithms[END_REF] is to examine several possible starting matrices. From this study, two choices emerge.

A first choice is simply to take a judicious multiple of the identity matrix:

H 0 k = δ k-1 I , δ k-1 := y k-1 , s k-1 |y k-1 | 2 .
(5.3)

The factor δ k-1 is known as the Oren-Spedicato factor (see [START_REF] Oren | Optimal conditioning of self-scaling variable metric algorithms[END_REF]) and aims at giving H 0 k a good scaling. Some justifications for using this factor are given in [START_REF] Gilbert | Some numerical experiments with variable-storage quasi-Newton algorithms[END_REF].

The other possibility is to take a diagonal matrix and to update it at each iteration. In this case, the algorithm generates and stores two sequences, {x k } k 1 and {H 0 k } k 1 , and a circular order array for the couples (y i , s i ). Then H k is formed by updating m times H 0 k , using the m couples {(y i , s i )} k-m i k-1 . At the beginning, H 0 2 = δ 1 I. Then, for k 2, the formula for updating the i-th diagonal element of H 0 k is the following:

(H 0 k+1 ) (i) = H 0 k y k , y k y k , s k (H 0 k ) (i) + y k , e i 2 y k , s k - H 0 k y k , y k ( s k , e i /(H 0 k ) (i) ) 2 y k , s k (H 0 k ) -1 s k , s k -1 , (5.4)
where {e i } 1 i n is an orthonormal basis for the given inner product on R n . It was obtained by scaling a diagonalized version of the bfgs formula and proved to be very effective.

Table 5.1 gathers the numerical results from [START_REF] Gilbert | Some numerical experiments with variable-storage quasi-Newton algorithms[END_REF], which were obtained with 8 real-life test-problems having a dimension between 34 and 1559 (the average value of n is 692). A comparison is done between the standard bfgs method implemented by Buckley and LeNir [5], the cg method of Hestenes and Stiefel [START_REF] Hestenes | Methods of conjugate gradients for solving linear systems[END_REF] implemented in the code va14 from the Harwell library and three algorithms from the modulopt library (INRIA): m1gc3 is the code of Buckley and LeNir [5], m1qn2 and m1qn3 use Nocedal's proposal with the scaling (5.3) and (5.4) respectively. Table t:lbfgs gives the number of function/gradient calls divided by the one required by the standard bfgs algorithm. The number m of updates ranges from 2 to 50. The sign > used for va14 means that the true number should be larger than the one given because the algorithm failed on one of the testproblems. Observe that the performance of the algorithms increases with m, except for va14, which is not a variable storage method. Observe also that m1qn2 and m1qn3 are better than the bfgs method as soon as m is greater than 10. . . 20 and 2. . . 5, respectively. This improvement seems mainly due to the scaling made by the initial matrix in these algorithms, which acts as a dynamic preconditioner. An algorithm similar to m1qn2 has been tested by Liu and Nocedal [35] who report similar results in their comparison with the method of Buckley and LeNir.

The convergence of m1qn2 on strongly convex functions has been proved in [35]. However, because of the counter-examples of Powell [55], there is some doubt about the possibility of extending this result for nonconvex functions. To make this remark clear, observe that, with exact line-searches, the limited memory methods can be seen as a generalization of the cg method of Hestenes-Stiefel or Polak-Ribière. Indeed, let m = 1 and H 0 k = I in Nocedal's proposal and suppose that exact line-searches are done, i.e., s k-1 , g k = 0. Then, y k-1 , d k-1 = -g k-1 , d k-1 = |g k-1 | 2 and using formula (5.2), we obtain:

d k = -H k g k = -Ū (I, y k-1 , s k-1 ) g k = -g k + y k-1 , g k y k-1 , d k-1 d k-1 = -g k + β PR k d k-1 .
As the PR method with exact line-searches does not converge on nonconvex functions, the limited memory method (at least for m = 1 and without scaling) does not converge either.

The partitioned quasi-Newton methods

The algorithms we have seen so far are suitable for any objective function f . The partitioned quasi-Newton (pqN) method proposed by Griewank and Toint [START_REF] Griewank | On the unconstrained optimization of partially separable functions[END_REF] is dedicated to the minimization of functions having a particular structure. A function f is said partially separable if it can be written as a sum of p elements functions f i :

f (x) = p i=1 f i (x) , ∀ x ∈ Ω
and if each of the element functions f i is invariant on some subspace E i ⊂ R n :

f i (x + e) = f i (x) , ∀ e ∈ E i .
A typical case is when the subspaces E i are spanned by some basic vectors of R n , which means that the f i 's only depend on few of the n variables. Note that the E i 's need not be disjoint, hence the name of partial separability.

In the pqN method, the Hessian of f , ∇ 2 f (x) = p i=1 ∇ 2 f i (x), is approximated at iteration k by a matrix of the form

B k = p i=1 B i k .
The matrix B i k approximates the Hessian of the i-th element function and is updated to satisfy an element-wise secant condition:

B i k+1 s k = y i k := ∇f i (x k+1 ) -∇f i (x k ).
This approach has some weak points:

• If the element functions are nonconvex, their Hessian is not positive definite and y i k , s k is not necessarily positive; therefore the overall approximation is not necessarily positive definite either and may even be singular although the true Hessian is positive definite.

• The preceding point implies that different update formulae (bfgs and sr1) are used, depending on the element function, which is a lack of homogeneity and an additional complexity of the method. • It requires more effort from the user who has to supply separately the gradient of the element functions.

But it has also important nice features:

• if the dimension of the subspaces E i is large, the element Hessian are identified very quickly because the B i 's have low rank. This implies that the overall approximation is rapidly correct; • the overall approximation is updated by a consistent high rank update;

• the notion of partial separability is basis independent, in contrast with the notion of sparsity.

Liu and Nocedal [35] have made some numerical comparison between the pqN approach and m1qn2. They report excellent results for the pqN method when the dimension n i of the subspaces E i is small. If the n i 's are not very small, m1qn2 becomes competitive. The comparison was made on the computing time and not on the number of function/gradient calls, which generally favors the pqN code. This seems due to the iteration time of the pqN method that can be 2. . . 5, sometimes 10. . . 15, times larger than the one of m1qn2.

To conclude this section, let us mention that the pqN method has been implemented by Toint [START_REF] Ph | VE08AD, a routine for partially separable optimization with bounded variables[END_REF] in the Harwell routine ve08.

m = 2 m

 2 = 5 m = 10 m = 20 m = 50 va14 > 1.

Table 2 .

 2 The first test-problem we consider is u1mt1, a program determining atmospheric pressure and velocity using a meteorological model and measurements collected at various points. The least-square function evaluated in this code depends on n = 1875 variables and requires 34.3 Kbytes of memory. Results are presented in Table2.1.

	Type of calculation	Kb	sec	κ
	f alone	-	0.04	-
	original f and ∇f	7.5	0.23	5.8
	f and ∇f by dd	-83.86	-
	f and ∇f by jakef	1142.4	1.19 29.8
	f and ∇f by padre2 1348.3	1.14 28.5
	by adol-c:	1362.4	-	-
	f with graph	-	4.75	-
	f without graph	-	2.24	-
	∇f with graph	-	5.50	-
	∇f without graph	-	3.14 79.5
	∇f by meteo	14.8	0.09	3.3

1: Comparison of performance on u1mt1

Table 2 .

 2 2: Comparison of performance on loops

			jakef			padre2		meteo
	ntimes sec	Kb	sec	κ	Kb	sec	κ	sec	κ
	5 0.01	78.	0.08	8.0	49. 0.04	4.0 0.01 2.0
	10 0.03	625.	0.64 21.3	390. 0.38 12.7 0.04 2.3
	15 0.10 2109.	2.19 21.9 1318. 1.25 12.5 0.13 2.3
	20 0.21 5000.	5.24 25.0 3125. 3.20 15.2 0.32 2.5
	25 0.41 9765. 10.41 25.4 6103. 7.21 17.6 0.62 2.5

Table 3 .

 3

1: Comparison of the methods FR, PR|FR, PR and PR + . Table 3.1 allows us to compare the performance of several cg methods: FR, PR, PR|FR and PR + . The results come from tests made in

Table 5 .

 5 92 > 1.92 > 1.92 > 1.92 > 1.92

	m1gc3	2.07	1.60	1.46	1.37	1.23
	m1qn2	1.49	1.14	1.07	0.97	0.92
	m1qn3	1.00	0.85	0.85	0.76	0.74

1: Comparison between the limited memory, cg and bfgs methods.

Further references

We conclude by mentioning some further references that can be useful to deepen the topics addressed in this paper. The books by Gill, Murray and Wright [21], McCormick [38], Minoux [39], also available in English, Dennis and Schnabel [8] and Fletcher [14] are classical textbooks, see also Lemaréchal [33]. The following references are collections of papers giving a state of the art: [9], [36], [54], [3], [31], [44], [50], and [12].


Luciano Da Fontoura
email: luciano@ifsc.usp.br

Luciano Da

F Costa

Similarity Means: Seeking for Robustness to Outliers

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The center of mass (or centroid), as well as the directly related concept of mean (e.g. [START_REF] Da | By means of the means: Arithmetic, geometric, harmonic[END_REF]), play a key role in many scientific and technologic areas as a consequence of their special important properties. More specifically, from the Newtonian mechanics perspective, the center of mass of an object corresponds to the point such that respective application of forces will only induce only linear, but no angular acceleration (e.g. [START_REF] Serway | Physics for scientists and engineers[END_REF]).

When taken from the statistical point of view, the center of mass becomes the arithmetic mean (as well as its weighted version), which is analogous to the concept of statistical expectance of a given random variable (e.g. [START_REF] Da | Statistical modeling[END_REF]). It should be observed that though the arithmetic mean is a mean, there are other types of means including geometric and harmonic (e.g. [START_REF] Da | By means of the means: Arithmetic, geometric, harmonic[END_REF]).

The arithmetic mean of a discrete set of values x i , i = 1, 2, . . . , N , sampled from a respective random variable X can be simply understood as corresponding to the sum of these values divided by the cardinality of the data set. An analogous operation applies for continuous data, with the sum being replaced by the integration. In the case of a random variable described by a respective probability density function, the respective arithmetic mean can be obtained in terms of the respective expectance, corresponding to the sum of the product between the random variable values and respective relative frequencies.

Other types of means, including the geometric, harmonic, and contraharmonic, also rely on adding terms de-rived from the samples values x i . The weighted arithmetic mean represents an immediate extension of the arithmetic mean which assigns weights to of each of the possible random variable values.

Though there are several types of means -including geometric, harmonic, contraharmonic, and root mean square, the arithmetic and weighted arithmetic mean have the special, and often desirable, property of being invariant to translations of the original random variable scores, in the sense of adding a constant value k to them. However, all the above mentioned means, including the arithmetic, are particularly susceptible to the presence of outliers and skewness (e.g. [START_REF] Da | By means of the means: Arithmetic, geometric, harmonic[END_REF]).

Though many of the types of means rely on adding terms, which are then normalized, the specific case of arithmetic value of a single variable can also be understood in terms of distances between the values of the sample and the respective mean. Given that distances are intrinsically related, through the concept of dissimilarity, to the concept of similarity (e.g. [START_REF] Da | On similarity[END_REF])), it becomes an interesting prospect to consider types of centers and means of a data set that to be obtained not from the distance between the values of the data elements and the mean, but in terms of respective similarities.

The present work aims at developing these ideas, more specifically the consideration of similarity indices (e.g. [START_REF] Cha | Comprehensive survey on distance/similarity measures between probability density functions[END_REF][START_REF] Vijaymeena | A survey on similarity measures in text mining[END_REF][START_REF] Brusco | A comparison of 71 binary similarity coefficients: The effect of base rates[END_REF][START_REF] Hamers | Similarity measures in scientometric research: The jaccard index versus salton's cosine formula[END_REF][START_REF] Cha | Comprehensive survey on distance/similarity measures between probability density functions[END_REF][START_REF]Jaccard index[END_REF][START_REF] Mirkin | Mathematical Classification and Clustering[END_REF][START_REF] Akbas | L1 norm based multiplicationfree cosine similiarity measures for big data analysis[END_REF][START_REF] Akbas | Multiplication-free neural networks[END_REF][START_REF] Da | On similarity[END_REF][START_REF] Da | Further generalizations of the Jaccard index[END_REF][START_REF] Da | Multisets[END_REF][START_REF] Da | Coincidence complex networks[END_REF][START_REF] Wolda | Similarity indices, sample size and diversity[END_REF][START_REF] Hamers | Similarity measures in scientometric research: The Jaccard index versus Salton's cosine formula[END_REF]) as subsidy for defining respective means.

For simplicity's sake, we shall be restricted to discrete 1 random variables taking exclusively positive values, but the extension to non-positive values is straightforward to be achieved by considering the real-valued version of the Jaccard and coincidence similarity indices (e.g. [START_REF]Jaccard index[END_REF][START_REF] Mirkin | Mathematical Classification and Clustering[END_REF][START_REF] Akbas | L1 norm based multiplicationfree cosine similiarity measures for big data analysis[END_REF][START_REF] Akbas | Multiplication-free neural networks[END_REF][START_REF] Da | On similarity[END_REF][START_REF] Da | Further generalizations of the Jaccard index[END_REF][START_REF] Da | Multisets[END_REF][START_REF] Da | Coincidence complex networks[END_REF]). Though there are many types of similarity indices (e.g. [START_REF] Cha | Comprehensive survey on distance/similarity measures between probability density functions[END_REF][START_REF] Vijaymeena | A survey on similarity measures in text mining[END_REF][START_REF] Brusco | A comparison of 71 binary similarity coefficients: The effect of base rates[END_REF][START_REF] Hamers | Similarity measures in scientometric research: The jaccard index versus salton's cosine formula[END_REF][START_REF] Cha | Comprehensive survey on distance/similarity measures between probability density functions[END_REF][START_REF]Jaccard index[END_REF][START_REF] Mirkin | Mathematical Classification and Clustering[END_REF][START_REF] Akbas | L1 norm based multiplicationfree cosine similiarity measures for big data analysis[END_REF][START_REF] Akbas | Multiplication-free neural networks[END_REF][START_REF] Da | On similarity[END_REF][START_REF] Da | Further generalizations of the Jaccard index[END_REF][START_REF] Da | Multisets[END_REF][START_REF] Da | Coincidence complex networks[END_REF][START_REF] Wolda | Similarity indices, sample size and diversity[END_REF][START_REF] Hamers | Similarity measures in scientometric research: The Jaccard index versus Salton's cosine formula[END_REF]), in the present work we shall focus on the Jaccard similarity index (e.g. [START_REF] Da | Further generalizations of the Jaccard index[END_REF][START_REF] Jaccard | Distribution de la flore alpine dans le bassin des dranses et dans quelques régions voisines[END_REF][START_REF]Jaccard index[END_REF]).

As it will be illustrated in the present work, the extension of the concept of mean to consider similarities allows the sensitivity to outliers to be substantially reduced.

We start by presenting the Jaccard similarity index for positive values, and then present the respectively obtained similarity mean, followed by an illustration of its substantial robustness to data translations and presence of outliers.

Similarity Indices

Given two N -dimensional non-zero vectors x and y containing only positive elements, the respective Jaccard similarity (e.g. [START_REF]Jaccard index[END_REF][START_REF] Da | Further generalizations of the Jaccard index[END_REF]) can be quantified as:

J(x, y) = N i=1 min(x i , y i ) N i=1 max(x i , y i ) (1)
Interestingly, we necessarily have that 0 < J(x, y) ≤ 1. We also have that:

J(x, y) = J(y, x) (2)
3 Toward Similarity Means

The rationale henceforth adopted for the consideration of a similarity index as a subsidy for defining a respective means corresponds to:

The similarity mean of a set of samples of a random variable X corresponds to the scalar value x that is most similar to the values in that set.

In the case of a set of N observations x i of a random variable X, the similarity between that set and a value x can be quantified by using the following expression:

J(x, y) = N i=1 min(x, y i ) N i=1 max(x, y i ) (3)
which corresponds to the Jaccard index above considering the et of available scores as one of its vectors and another vector with the same dimension containing all entries identical to x.

Thus, the similarity mean in the above situation can be obtained by searching for the value x that maximizes Equation 3. This maximization can be done in several ways, including considering a sequence of equally spaced values along the range of the original scores, or an optimization method as gradient descent.

The approach described above for defining the similarity mean given a set of samples can be adapted to considering a respective discrete probability density function, involving respective relative frequencies f i in the following manner:

S = N i=1 min(x, x i ) f i N i=1 max(x, x i ) f i (4)
where the relative frequencies f i act as respective weights.

Figure 1 illustrates the arithmetic, as well as the Jaccard means, respectively to a Gaussian-like probability density function. In this particular case, we have that the arithmetic and Jaccard means coincided, taking place right at the middle of the distribution. However, observe that the profile of Jaccard similarities obtained in terms of x (shown in green) resulted moderately skewed.

The effect of scaling a density probability function on the respective arithmetic and Jaccard means are illustrated in Figure 2, which indicates that directly analogous results are obtained.

Figure 3 presents the same distribution as above, but after being translated by 2 to the right-hand side. Again, the arithmetic and Jaccard means coincided, while the similarity profile became less skewed. In order to illustrate the enhanced robustness of the Jaccard mean, we now add an outlier to the previous distribution, as illustrated in Figure 4. Interestingly, the Jaccard mean resulted substantially less affected by the presence of the outlier, while the re-spective similarity profiled underwent a relatively moderate modification when compared with the original distribution in Figure 1.

Figure 5 illustrates the effect of an even stronger outlier, in which case the arithmetic mean resulted substantially. more displaced, while the Jaccard undergoing a more moderate shift.

Figure 5: The effect of adding an even stronger outlier to the distribution in Fig. 1. The Jaccard mean (in red) suffered a substantially smaller effect than the arithmetic mean (in black). The similarity profile (in green) underwent a relatively moderate modification.

Concluding Remarks

The concept of mean, and particularly its arithmetic and weighted arithmetic versions, are frequently adopted in a wide range of areas in science and technology as a means to convey centrality and summarize properties of data sets.

Though the arithmetic and weighted arithmetic means distinguish themselves from other types of means (e.g. geometric, harmonic, and root mean square) by being translational invariant, these two types of means are also particularly susceptible to the presence of outliers or skewness.

In the present work, we developed the concept of similarity mea, in which some type of similarity such as the Jaccard index, are employed to define a respective mean. This is done by searching for the value x that is maximally similar to the set S of values in the supplied set of observations (or scores) of a respective random variable X.

In the present work, the similarity between each candidate value x has its overall similarity with the values in the set S quantified in terms of the Jaccard similarity between the vector containing the N original values in S and an N -dimensional vector with all entries identical to the candidate value x.

This approach can be readily extended to data sets described by probability density functions instead of a set of sampled observations. This has been done by using the relative frequencies of each entry as a weights while calculating the similarity between the two vectors as described above.

The reported approach has been illustrated, at least for the considered data sets, to retain its translational invariance while providing, when compared to the arithmetic mean, substantially enhanced robustness to the presence of outliers.

Future works could involve using other similarity approaches such as the Sørensen-Dice and coincidence indices. In addition, it would be interesting to study the properties of similarity means respectively to more types of data sets and in presence of noise and interference between distributions. Another particularly promising perspective concerns the extension of the suggested concepts and methods to higher dimensional spaces. Yet another interesting prospect regards the identification, in complex networks in which the nodes have been characterized in terms of respective topological measurements, of the node that is most similar to the other nodes, therefore corresponding to its similarity mean node.

The described concepts and developments are believed to have particular potential for applications in data characterization, pattern recognition and deep learning, among many other scientific and technological areas.

Figure 1 :

 1 Figure1: A gaussian-like discrete probability density (in blue) and the respectively estimated arithmetic mean (in black) and Jaccard mean (in red). In this particular case, these two means can be verified to coincide one another. Also shown (in green) is the profile of Jaccard similarity values obtained for a fine-grained sequence of values x ranging along the domain of the density function (i.e. from 1 to 7).

Figure 2 :

 2 Figure 2: The arithmetic and Jaccard means obtained for a scaled (by a factor of 4) version of the distribution in Fig. 1. The obtained results are directly analogous to those shown in Fig. 1.

Figure 3 :

 3 Figure 3: The arithmetic and Jaccard means obtained for a righthand translated (shifted) version of the distribution in Fig. 1. Though these two mean again coincided one another, the similarity profile became less skewed as a consequence of the implemented translation.

Figure 4 :

 4 Figure 4: The effect of adding an outlier to the distribution in Fig. 1. The Jaccard mean (in red) suffered a substantially smaller effect than the arithmetic mean (in black). The similarity profile (in green) underwent a relatively moderate modification.

Acknowledgments

Luciano da F. Costa thanks CNPq (grant no. 307085/2018-0) and FAPESP (grant 15/22308-2).

Observations

As all other preprints by the author, the present work contains preliminary work subject to further revision and validation. Respective modification, commercial use, or distribution of any of its parts are not possible, as this work has author copyright. Many of the preprints by the author are also available in Hal and/or arXiv. This work can also be cited by using the DOI number or article identification link. It is also interesting to check for new versions of this work, as they may become available. Thanks for reading.


Espen Gaarder Haug
email: espenhaug@mac.com

Reissner-Nordström Kerr Kerr-Newman

Keywords: general relativity theory, cosmological constant extremal solution, Reissner-Nordström, Kerr, Kerr-Newman

Einstein's field equation is a highly general equation consisting of sixteen equations. However, the equation itself provides limited information about the universe unless it is solved with di↵erent boundary conditions. Multiple solutions have been utilized to predict cosmic scales, and among them, the Friedmann-Lemaître-Robertson-Walker solution that is the back-bone of the development into todays standard model of modern cosmology: the ⇤-CDM model. However this is naturally not the only solution to Einsteins field equation.

We will explore the extremal solutions of the Reissner-Nordström, Kerr, and Kerr-Newman metrics. Interestingly, in their extremal cases, these solutions yield identical predictions for horizons and escape velocity. These solutions can be used to formulate a new cosmological model that resembles the Friedmann equation. However, a significant distinction lies in the fact that the extremal universe solution does not require the ad hoc insertion of the cosmological constant; instead, it naturally emerges from the derivation itself. To the best of our knowledge, all other solutions relying on the cosmological constant do so by initially ad hoc inserting it into Einstein's field equation. This clarification unveils the true nature of the cosmological constant, suggesting that it serves as a correction factor for strong gravitational fields, accurately predicting real-world cosmological phenomena only within the extremal solutions of the discussed metrics, all derived strictly from Einstein's field equation.

Moreover, in the extremal model, the cosmological redshift seems to be explained simply as a special type of cosmic gravitational redshift arising from the extremal solution. This eliminates the need for the hypothesis of an expanding space. Although this view may challenge the prevailing consensus in cosmology today, we must question the validity of the Big Bang hypothesis based on pure science and derivations from Einstein's field equation, even if it should make us unpopular.

The Extremal Universe Exact Solution from Einstein's Field Equation gives the Cosmological Constant directly and explains the Cosmological Redshift without the need for the Big Bang Hypothesis

1 Extremal solutions to Einsteins field equation

The Reissner-Nordström [START_REF] Reissner | Über die eigengravitation des elektrischen feldes nach der einsteinschen theorie[END_REF][START_REF] Nordström | On the energy of the gravitation field in Einstein's theory[END_REF] metric for a spherical charged gravitational object is an exact solution to Einsteins [START_REF] Einstein | Die grundlage der allgemeinen relativitätstheorie[END_REF] field equation and is given by:

ds 2 = ✓ 1 2GM c 2 r + r 2 Q r 2 ◆ c 2 dt 2 ✓ 1 2GM c 2 r + r 2 Q r 2 ◆ 1 dr 2 r 2 (d✓ 2 + sin 2 ✓d 2 ) (1)
Here, in SI units, we have

r 2 Q = k e qq G c 4 ,
where k e is the Coulomb constant and q is the charge. The special case when r Q = GM c 2 is well known as the extremal solution of the Reissner-Nordström metric, seen for example Zee [START_REF] Zee | Einstein gravity in a nutshell[END_REF] and Aretakis [START_REF] Aretakis | Dynamics of Extremal Black Holes[END_REF]. Furthermore, the Kerr [START_REF] Kerr | Gravitational field of a spinning mass as an example of algebraically special metrics[END_REF] metric is given by:

ds 2 = ⌃ (cdt a sin 2 ✓d ) 2 ⌃ dr 2 ⌃d✓ 2 sin 2 ✓ ⌃ (adt (r 2 + a 2 )d ) 2 (2)
where = r 2 r s r + a 2 and ⌃ = r 2 + a 2 cos 2 ✓, and r s = 2GM c 2 . The Kerr metric also has an extremal solution when a = GM c 2 . The Kerr-Newman [START_REF] Newman | Note on the kerr spinning-particle metric[END_REF][START_REF] Newman | Metric of a rotating, charged mass[END_REF] metric extends the Kerr metric to include charge, and it is given by:

ds 2 = ⌃ (cdt a sin 2 ✓d ) 2 ⌃ dr 2 ⌃d✓ 2 sin 2 ✓ ⌃ (adt (r 2 + a 2 )d ) 2 (3)
Here, = r 2 r s r + a 2 + r 2 Q and ⌃ = r 2 + a 2 cos 2 ✓, where r s = 2GM c 2 . In the special case of a = 0, it simplifies to the Reissner-Nordström metric and in the special case of r Q = 0, it simplifies to the Kerr metric.

The extremal solutions of the Reissner-Nordström metric (r Q = GM/c 2 ), as well as the extremal solution of the Kerr metric (a = GM/c 2 ) and the extremal solution of the Kerr-Newman metric [START_REF] Zee | Einstein gravity in a nutshell[END_REF][START_REF] Guidry | Modern General Relativity[END_REF]), all have one and the same horizon given by:

(a 2 + r 2 Q = G 2 M 2 /c 4 ) (see
r h = GM c 2 (4)
This is half the Schwarzschild radius. In the Schwarzschild [START_REF] Schwarzschild | über das gravitationsfeld einer kugel aus inkompressibler flussigkeit nach der einsteinschen theorie[END_REF] metric, the escape velocity is

v e = q 2GM
r , see [START_REF] Augousti | An observation on the congruence of the escape velocity in classical mechanics and general relativity in a Schwarzschild metric[END_REF]. However, in the extremal solutions of the Reissner-Nordström, Kerr, and Kerr-Newman metrics, the escape velocity is given by:

v e = r 2GM r G 2 M 2 c 2 r 2 (5)
We will soon see that this higher-order term G 2 M 2 c 2 r 2 , which di↵erentiates it from, for example, the Schwarzschild metric, could play a critical role in understanding the cosmos.

Table 1 summarizes the known key results from the extremal solutions. However, relatively few people have shown interest in these extremal solutions, resulting in a limited number of predictions being discussed based on them. From our perspective, these predictions are unrelated to cosmological models, which we will explore in the next section.

G 2 M 2 c 2 r 2 v e = q 2GM r G 2 M 2 c 2 r 2
Table 1: The table sumarize the extremal solutions of the metrics we will be looking at.

Cosmological Model

The horizon and escape velocity play a central role in predicting black holes. Interestingly, the Hubble sphere also exhibits several mathematical aspects of a black hole, including a horizon known as the Hubble radius (r H = c H 0 ) as we now will discuss. The Friedmann critical mass-equivalent for the universe is given by M c = c 3 2GH 0 . When considering a Hubble sphere with the Friedmann critical mass, the Schwarzschild radius is given by

r s = 2GM c c 2 = 2G c 3 2GH 0 c 2 = c H 0 (6)
This implies that the Schwarzschild radius is exactly identical to the Hubble radius if the Hubble sphere were filled with the Friedmann critical mass-energy. The mathematical similarities between Hubble spheres and black holes have led several researchers in prominent journals such as Nature to suggest and speculate that the observable universe could be inside a black hole (see Pathria [START_REF] Pathria | The universe as a black hole[END_REF] and Stuckey [START_REF] Stuckey | The observable universe inside a black hole[END_REF]).

It is important to note that our intention is not to claim that we live inside a black hole, but rather to highlight the mathematical properties shared by the Hubble sphere and black holes. However, in the ⇤-CMD model, the universe has expanded well beyond the Hubble radius due to the assumption of space expansion, including an accelerating expansion attributed to dark energy. The hypothesis of dark energy appeared necessary to reconcile the model with high-redshift supernova observations.

The escape velocity and the horizon derived from the metric of interest have practical applications, even in cosmology. For instance, Schutz [START_REF] Schutz | Gravity from the Ground Up[END_REF] derived the critical density of the universe based on simply the escape velocity formula from the Schwarzschild metric, resulting in an equation that depends on the Hubble constant. However, in this discussion, we will not explore cosmology through the escape velocity derived from the Schwarzschild metric. Instead, we will focus on the predictions of the escape velocity using the extremal solutions of the Reissner-Nordström, Kerr, and Kerr-Newman metrics.

Subsequently, we will derive an equation analogous to the Friedmann [START_REF] Friedmann | Über die krüng des raumes[END_REF] equation from the escape velocity in the extremal solutions. Our interest lies initially in the special case where the escape velocity is equal to the speed of light (c). Let us rephrase the escape velocity formula as follows:

r 2GM r G 2 M 2 c 2 r = c 2GM r G 2 M 2 c 2 r 2 = c 2 8⇡G⇢r 2 3 G 2 M 2 c 2 r 2 = c 2 8⇡G⇢ 3 G 2 M 2 c 2 r 4 = c 2 r 2 8⇡G⇢ 3 G 2 M 2 c 2 r 4 3 = c 2 r 2 (7)
Here, ⇢ = M V = M 4 3 ⇡r 3 represents the volumetric mass density of a sphere with radius r. Next, we replace r with the Hubble radius r H = c H 0 , and furthermore, it must be equal to r h = GM c 2 since the escape velocity is c at the horizon r h = GM c 2 . Therefore, we assume that

r H = r h = GM u c 2 = c H 0 (8)
Here, M u is the mass equivalent of all mass and energy in the universe we are considering. Solving for

M u gives M u = c 3 GH 0 ⇡ 1.77 ⇥ 10 53 kg (when assuming H 0 = 70 km • s 1 • Mpc 1
), which means the mass is exactly twice that of the critical mass (mass equivalent) of the Friedmann universe, which is

M c = c 3 2GH 0 . Next, let's replace r = r H = c
H 0 in equation 7, and we obtain:

8⇡G⇢ 3 G 2 M 2 u c 2 r 4 H 3 = c 2 r 2 H 8⇡G⇢ 3 G 2 M 2 u c 4 r 4 H c 2 3 = c 2 r 2 H 8⇡G⇢ 3c 2 r 2 H 3 = c 2 c 2 H 2 0 8⇡G⇢ 3 r 2 H 3 = H 2 0 8⇡G⇢ 3 H 0 c 2 c 2 3 = H 2 0 (9)
Next, it's important to note that the ad hoc inserted cosmological constant in general relativity theory is identical to

⇤ = 3 H 0 c 2 ⌦ ⇤ . When ⌦ ⇤ = 1, we have ⇤ = 3 H 0 c
2 , which means we can rewrite the equation as follows:

H 2 0 = 8⇡G⇢ ⇤c 2 3 ( 10
)
This equation is identical to the Friedmann equation for homogeneous, isotropic universe, except our cosmological constant has been derived from the extremal solutions of Einstein's field equations rather than being ad hoc inserted in the field equation (or later). To our knowledge, the extremal solutions have not been previously utilized to construct a cosmological model. In 1917, Einstein [START_REF] Einstein | Cosmological considerations in the general theory of relativity[END_REF] ad hoc inserted a cosmological constant into his field equation, which he referred to as an extended field equation. It was based on sound reasoning and was actually aimed at achieving a steady-state universe model. However, after Hubble's discovery of cosmological redshift in 1929, Einstein removed the cosmological constant, allegedly referring to it as his biggest blunder (although this statement is uncertain, as it comes from a single source, Gamow [START_REF] Gamow | My World Line : An Informal Autobiography[END_REF]). Later, in 1998, a astrophysicists team led by Saul Perlmutter [START_REF] Perlmutter | Measurements of ! and from 42 high-redshift supernovae[END_REF] and another led by Brian Schmidt and Adam Riess [START_REF] Reiss | Observational evidence from supernovae for an accelerating universe and a cosmological constant[END_REF] that observed high-redshift supernovae that did not conform to the model, the cosmological constant was again reintroduced and praised along with the hypothesis of dark energy.

However, for the first time, we have a cosmological model that is identical to the Friedmann model, except the cosmological constant is derived and likely carries a considerably di↵erent interpretation than it is traditionally given. If we use energy density rather than mass density we get

H 2 0 = 8⇡G⇢ E ⇤c 4 3 ( 11
)
where:

⇢ E = E V = E 4 3 ⇡r 3 H = Mc 2 4 3 ⇡r 3 H
, represents the volumetric energy density of a sphere with radius equal to the Hubble radius.

Cosmological Redshift from Einstein's Extremal Universe

The cosmological redshift in the extremal solution is as follow:

z = 1 2 2 = q 1 2GMu r 1 c 2 + G 2 M 2 u c 4 r 2 q 1 2GMu r 2 c 2 + G 2 M 2 u c 4 r 2 2 1 ( 12
)
Here, r 1 and r 2 represent the distance from the emitter (for example a supernova or galaxy) to the Hubble sphere horizon and r 2 the distance from the Hubble sphere horizon to the observer (in our case basically Earth based observatories). The first term of the Taylor series expansion is given by:

z ⇡ (r 1 r 2 )GM u c 2 r 1 r 2 (13) If r 1 = r H = GMu c 2
, which corresponds to the observer's distance to the Hubble sphere horizon, then we have:

z ⇡ (r H r 2 )GM u c 2 r H r 2 Furthermore, since M u = c 3
GH 0 , we can substitute GMu c 2 with r H , resulting in:

z ⇡ (r H r 2 ) r 2 (14)
Moreover, when the object emitting the photons (galaxies, quasars 1 , supernovas) is not too far away from us, we can approximate also r 2 in the denominator as r H . Substituting this approximation into the denominator, we obtain:

z ⇡ (r H r 2 ) c H 0 z ⇡ (r H r 2 )H 0 c (15)
We define the distance d as the di↵erence between r H and r 2 . This distance represents the distance from us to the object that emits the photons, such as stars, galaxies, supernovas, and quasars. Consequently, the expression becomes:

z ⇡ dH 0 c (16)
This corresponds to the well-known prediction of cosmological redshift approximation also used in the standard model (see for example [START_REF] Schutz | Gravity from the Ground Up[END_REF][START_REF] Wald | General Relativity[END_REF][START_REF] Schneider | Extragalactic Astronomy and Cosmology[END_REF]). However, in the extremal solutions of Einstein's field equations that we have derived, this redshift does not seem to be related to the expansion of space. Instead, it appears to be a pure gravitational redshift caused by the mass (energy) within the entire Hubble sphere. Therefore, we must, based on pure scientific principles, assert that the universe is likely not expanding and even question the entire Big Bang hypothesis.

To our knowledge, the extremal solution is the only solution in Einstein's field equations that directly yields the cosmological constant without any ad hoc insertion, as we demonstrated in the previous section, so our model seems simpler and likely more robust than models that need ad hoc inserted constants. Our new cosmological model is an exact solution to Einsteins original field equation with no ad hoc inserted constant. This alone should at least hopefully be enough to make other researchers curious and investigate the model in more depth than we can do in a single paper. Only rigorous investigation by a series of researchers over time can determine whether this extreme universe is a better model than the current standard model. However, this possibility should not be prematurely excluded, despite going against the consensus. After all, both solutions come from Einstein's field equation.

Implications

We will now summarize some of the most important implications in terms of cosmological predictions of the extremal universe:

• The cosmological constant does not need to be ad hoc inserted as done today; it automatically arises from the extremal solutions of Einsteins field equation and is given by

⇤ = 3 H 0 c 2 .
• The cosmological redshift prediction in the non-extremal solutions will incorrectly predict cosmological redshift for objects very far away, i.e., those significantly close to the Hubble sphere horizon (the Hubble horizon). If the extremal universe is the correct model of the real universe, then other solutions of the Einstein field equation will give incorrect predictions of the cosmological constant for objects significantly close to the Hubble horizon. To address 1 quasars are assumed to be early forming glaxies.

this, hypotheses of dark energy have been added, but this is likely unnecessary and requires careful investigation over some time before any firm conclusions can be made, see also [START_REF] Haug | Lorentz relativistic mass makes dark energy superfluous?[END_REF]. We also note that the extremal solutions of Einsteins field equation that we have discussed gives exactly the same escape velocity and horizon as a relativistic modified Newtonian theory [START_REF] Haug | A new full relativistic escape velocity and a new Hubble related equation for the universe[END_REF], it should also be investigated if this is a coincidence or not, we think not.

• The amount of energy (mass) in the Hubble sphere is twice the Friedmann critical mass. Therefore, it is

M u = c 3 GH 0 instead of M c = c 3 2GH 0 .

Deeper philosophical aspects of the Extremal solutions

The extremal solutions have received relatively little attention, especially regarding their predictions in cosmology. Although all three metrics studied yield the same horizon and escape velocity (in the extremal solutions), they di↵er in their interpretation. The extremal solution of the Reissner-Nordström metric lacks rotation but possesses charge, while the extremal solution of the Kerr metric has rotation but no charge. The Kerr-Newman metric's extremal solution possesses both charge and rotation. Nonetheless, all of these solutions lead to the same cosmological equation, as shown in Equation 11. Still from a deeper philosophical aspect they have di↵erent interpretations.

Our findings should not be considered in isolation be considered in conjunction with recent research by Haug and Spavieri [START_REF] Haug | Micro black hole candidates and the Planck scale: Schwarzschild micro black holes can only match a few properties of the Planck scale, while a Reissner-Mordström and Kerr micro black hole matches all the properties of the planck scale[END_REF], who have recently demonstrated that only the extremal solutions of the Reissner-Nordström, Kerr, and Kerr-Newman metrics provide a perfect match between the Planck mass micro black hole and all the properties of the Planck scale. In contrast, the Schwarzschild metric can only match a few aspects of the Planck scale.

It is also considered a mystery why the electromagnetic force is enormous compared to the gravitational force. If we compare the Coulomb force between a proton and an electron to the theoretical gravitational force between a proton and an electron, we obtain:

|F c | |F G | = k e |e||e| r 2 G Mprme r 2 ⇡ 2.26 ⇥ 10 39 (17)
Where F c represents the Coulomb force [START_REF] Coulomb | Premier mémoire sur l'électricité et le magnétisme[END_REF], F G denotes the Newtonian gravitational force [START_REF] Newton | Philosophiae Naturalis Principia Mathematica[END_REF], and e represents the elementary charge. Additionally, M pr and m e respectively refer to the proton and electron masses. The significant disparity in strength between the electrostatic and gravitational forces is well-documented in the literature. However, despite this knowledge, the gravitational force between a proton and an electron has never been measured. Thus, there is still clearly room for us to gain a deeper understanding of gravity at the atomic and subatomic scales.

On the other hand for two Planck [START_REF] Planck | Der Königlich Preussischen Akademie Der Wissenschaften[END_REF][START_REF] Planck | Vorlesungen über die Theorie der Wärmestrahlung[END_REF] masses (m p = q ~c G ) the electrostatic force is identical to the gravitational force as we have

|F c | |F G | = k e |qp||qp| r 2 G mpmp r 2 = 1 (18)
Here, q p represents the Planck charge: q p = e p ↵ . The fact that these forces are equal at the Planck scale indicates the potential unification of electromagnetic and gravitational forces at the Planck scale as expected by multiple researchers. However, this is based on Newton's theory, and we need to move beyond it. In the extremal solution of the Reissner-Nordström metric, we have:

r 2 Q = nk e |q p ||q p | G c 4 = nGm p m p G c 4 (19)
Here, n represents the number of Planck masses in the large gravitational mass M , so we have n = M mp . The extremal solution of Reissner-Nordström is consistent with the electrostatic force being identical to the gravitational force at the Planck scale. The term G c 4 is identical to part of Einsteins gravitational constant and is needed to convert the units to the right form needed for predicting gravity phenomena.

If gravity is ultimately caused at the Planck scale as first suggested by Eddington [START_REF] Eddington | Report On The Relativity Theory Of Gravitation[END_REF] in 1918 and assumed by most researchers working on quantum gravity theory today (see for example [START_REF] Adler | Six easy roads to the Planck scale[END_REF][START_REF] Hossenfelder | Can we measure structures to a precision better than the Planck length? Classical and Quantum Gravity[END_REF][START_REF] Hossenfelder | Minimal length scale scenarios for quantum gravity[END_REF][START_REF] Haug | Di↵erent mass definitions and their pluses and minuses related to gravity[END_REF]), then the extremal solution could be the only truly valid exact solution for real phenomena. This suggests that the extremal solution of the Reissner-Nordström metric could be the most realistic model for the universe. This possibility could explain why no ad hoc inserted constants are needed in this specific solution to fit cosmological observations. Naturally, this hypothesis needs to be carefully investigated and, at this stage, can be seen as plausible.

Conclusion

We have demonstrated that the extremal solutions of the Reissner-Nordström, Kerr, and Kerr-Newman metrics all yield the same cosmological model, which is analogous to the Friedmann equation. However, there is a significant distinction in that the cosmological constant is now derived rather than being ad hoc inserted into the field equation. Additionally, the cosmological redshift observed in the extremal solution appears to be a specific instance of gravitational redshift caused by the mass and energy within the Hubble sphere.

This discovery implies that the universe may ultimately exist in a steady state, suggesting that the Big Bang never actually occurred. Naturally, this notion requires careful study and investigation. Nonetheless, considering that this new cosmological model is an exact solution to Einstein's field equation and the only solution known to us where the cosmological constants emerge automatically without the need for ad hoc insertion, we believe it warrants careful consideration by the research community over an extended period before being rejected or accepted.

Acknowledgment

Thanks to (AB-anonymous) and William Bliss for useful comments.

Data Availability Statement:

No data has been used for this study.

Conflicts of Interest:

The author declares no conflict of interest.

Appendix

c 4 r 4 and solve with respect to M this gives

For the Hubble sphere the mass is

Alternatively set the Cosmolocal constant ⇤ = 0 and the mass of the universe is

GH 0 is the first taylor series approximation