domain: math info

text:
Jean Charles Gilbert
email: jean-charles.gilbert@inria.fr.

Claude Lemaréchal
email: claude.lemarechal@inria.fr.

The module M1QN3 Version 3.3 (

Goal of the software

The module m1qn3 is a Fortran-77 software that can solve a large scale unconstrained minimization problem of the form:

min {f (x) : x ∈ R n },
where the function f : R n → R is supposed to be smooth (continuously differentiable). M1qn3 has solved problems with n as large as 10 8 . The solver requires evaluations of the function f and its gradient g, computed for a given inner product, denoted by •, • . The gradients are used to form a limited memory quasi-Newton approximation of the Hessian of f . The decrease in f is enforced by the Wolfe line-search.

The software can be used either with direct or reverse communication (these terms are used throughout the presentation of the software; the reader not familiar with this terminology can find clarification in sections 4.1 and 4.2).

In principle, m1qn3 is intended for problems with a large number n of variables, but it can also be suitable for small or medium scale problems. The method used to solve the problem is a limited memory (or variable-storage) quasi-Newton method. This type of methods allows m1qn3 to take advantage of the storage that is declared available in core memory. However, it cannot really be ascertained that, when more storage is used, the algorithm becomes more efficient and the decrease in the objective function f more effective (see section 4.3). The software has two scaling modes: the Scalar Initial Scaling mode (or SIS mode) and the Diagonal Initial Scaling mode (or DIS mode). The latter should normally give better results than the former, but it requires one additional vector of storage.

The module m1qn3 includes of the following subroutines:

-m1qn3: main subroutine that has to be called and that calls m1qn3a after having structured the available memory, -m1qn3a: optimization subroutine, dd, dds: subroutines that compute the descent direction, -mlis3, ecube: subroutines realizing the line-search, mupdts, ystbl: special subroutines for very large problems.

On his part, the user has to supply:

-a program that calls the optimization module m1qn3, -(used only in direct communication) a subroutine, called the simulator , and named simul by m1qn3, that computes the value of f and its gradient g (associated with the inner product •, • ) at a given point, -a subroutine, named prosca by m1qn3, that computes the inner product u, v of two vectors of u, v ∈ R n , -(used only in DIS mode) a subroutine, named ctonb by m1qn3, that gives the coordinates y in an orthonormal basis (for the inner product prosca) of a vector whose coordinates in the canonical basis of R n are x, -(used only in DIS mode) a subroutine, named ctcab by m1qn3, that does the operation converse to ctonb: it gives the coordinates x in the canonical basis of R n of a vector whose coordinates in the orthonormal basis are y.

The standard distribution also provides the subroutines simul rc, euclid, ctonbe, and ctcabe. The subroutine simul rc plays the role of simul in reverse communication (it is an empty subroutine). The subroutines euclid, ctonbe, and ctcabe can play the role of prosca, ctonb, and ctcab inside m1qn3, when the Euclidean (or dot) product is used to compute the gradient of f , i.e., when u, v = u ⊤ v.

2 Description of the method 2

Brief description

In SIS mode, m1qn3 is an implementation of an algorithm proposed by Nocedal [START_REF]Updating quasi-Newton matrices with limited storage[END_REF], while the DIS mode uses a diagonal preconditioner improving the algorithm. Both scaling modes are based on the same principles, which we now describe.

At each iteration k, k ≥ 1, a descent direction d k of f at the current iterate x k is determined. This direction has the form d k = -W k g k , where W k is the current approximation of the inverse Hessian ∇ 2 f (x k ) -1 of f at x k and g k := ∇f (x k ) is the gradient of f at x k , computed for the inner product •, • . Then, a step-size α k is determined along the direction d k by a line-search procedure (mlis3). The next iterate has the form:

x k+1 = x k + α k d k .
When k ≥ m + 1, the approximation W k of the inverse Hessian is obtained by doing m updates of a positive definite diagonal matrix D k , using the inverse BFGS formula. For this, the following m pairs of vectors are used in order:

{(y i , s i ) : k -m ≤ i ≤ k -1}, (2.1)
where s i := x i+1 -x i is the step and y i := g i+1 -g i is the change in the gradient.

The matrix W k is not stored in memory but instead the product W k g k is formed from the pairs (2.1) by an efficient appropriate algorithm. The integer parameter m is determined by m1qn3 in order to use as much as possible the storage declared available in core memory (see the meaning of the argument ndz in section 3).

The user has the possibility to improve the conditioning of the problem by using an adapted inner product •, • , the one realized by the subroutine prosca. Let us denote by • := •, • 1/2 the norm associated with that inner product. Let us insist on the fact that this inner product must be the one used by the simulator to compute the gradient g(x) of f at x. Indeed, for m1qn3, the derivative of f at x along the direction h ∈ R n is given by:

f ′ (x) • h = g(x), h .
The BFGS formula used by m1qn3 is adapted to this inner product (see section 2.2). Very often, the Euclidean (or dot) inner product, defined by

x, y := x ⊤ y = n i=1 x i y i , (2.2)
is used. In that case, the gradient of f is formed of the partial derivatives ∂f /∂x i , i = 1, . . . , n.

The way of computing the starting matrix D k depends on and characterizes the scaling mode of m1qn3. In the Scalar Initial Scaling (SIS) mode, D k has the form δ k-1 I, where I is the identity matrix and

δ k-1 := y k-1 , s k-1 y k-1 2 (2.3)
aims at giving W k a good scaling. In the Diagonal Initial Scaling (DIS) mode, D k is a diagonal matrix obtained from D k-1 by using a diagonal update formula (see section 2.2 for the details). The first update (to obtain D 2 ) uses δ 1 I as starting matrix. The notion of diagonal matrix is defined with respect to the given inner product and the orthonormal basis (e i ) 1≤i≤n .

In DIS mode, an implicit reference is made to an orthonormal basis of R n for the inner product •, • . Let (e i ) 1≤i≤n denote this basis. If x, y = x ⊤ L ⊤ Ly with L nonsingular and if (ê i ) 1≤i≤n is the canonical basis of R n , i.e., êi = (0, . . ., 0, 1, 0, . . ., 0) with 1 in the ith position, you can take e i = L -1 êi , for 1 ≤ i ≤ n. Luckily, this basis needs not be stored, but you have to provide subroutines that give the coordinates in the (e i ) 1≤i≤n basis of a vector whose coordinates in the (ê i ) 1≤i≤n basis are given (subroutine ctonb) and vice-versa (subroutine ctcab). With the inner product and the orthonormal basis given above, ctonb must compute y := Lx and ctcab must compute x := L -1 y. If •, • is the Euclidean inner product, ctonb has just to copy x in y and ctcab will copy y in x.

If N and L denote respectively the number of iterations and the total number of step-size trials done during a particular run, then the subroutine prosca is called 4 + 7N + L times in the two scaling modes and the subroutines ctonb and ctcab are called respectively 3N and N times in DIS mode.

The methods and numerical results are given and discussed by Gilbert and Lemaréchal [START_REF] Gilbert | Some numerical experiments with variable storage quasi-Newton algorithms[END_REF]. See also the numerical experiments by Liu and Nocedal [START_REF] Liu | On the limited memory BFGS method for large scale optimization[END_REF].

Algorithmic details

The BFGS formula used in the software updates a matrix W to Ŵ using two vectors y and s in R n by: Ŵ = I -s ⊗ y y, s W I -y ⊗ s y, s + s ⊗ s y, s .

In this formula, u ⊗ v is the linear operator that to d ∈ R n associates v, d u ∈ R n . This makes the formula adapted to the inner product •, • . We denote it by Ŵ = bfgs(W, y, s).

When k ≥ m + 1, given the diagonal matrix D k (its computation is discussed below), the matrix W k used to form the search direction d k = -W k g k is obtained as follows:

W 0 k := D k , W i+1 k := bfgs(W i k , y k-m+i , s k-m+i ), for 0 ≤ i ≤ m -1, W k := W m k .
Let us summarize this scheme by: W k := bfgs k-1 k-m (D k ). When 1 ≤ k ≤ m, the matrix W k depends on how the optimizer is tuned by imode(2). For a cold start (imode(2) = 0), W 1 = I and, for 2 ≤ k ≤ m, W k is obtained by W k := bfgs k-1 1 (D k ). For a warm restart (imode(2) = 1), W k is obtained by: W k := bfgs k-1 k-m (D k ), in which (y i , s i ) with index i ≤ 0 refers to a pair of a previous run, stored in the working zone dz.

As mentioned in section 2.1, SIS and DIS modes differ by the way the diagonal matrix D k is computed. In SIS mode, D k = δ k-1 I for k ≥ 2, while in DIS mode, D k is a diagonal matrix updated by the following formula:

D (i) k+1 = D k y k , y k y k , s k D (i) k + y k , e i 2 y k , s k - D k y k , y k s k , e i 2 y k , s k D -1 k s k , s k (D (i) k ) 2 -1
.

The subroutine m1qn3a uses the subroutine mlis3 from modulopt for the linesearch. The step-size α k is determined to satisfy the following two Wolfe's conditions:

f (x k +α k d k ) ≤ f (x k ) + ω 1 α k g k , d k , g(x k +α k d k ), d k ≥ ω 2 g k , d k .
It is necessary to have 0 < ω 1 < 1/2 and ω 1 < ω 2 < 1. We took ω 1 = 0.0001 and ω 2 = 0.9. These constants are set as parameter (rm1 for ω 1 and rm2 for ω 2 ) in m1qn3a.

Usage

In the description of the subroutines below, an argument flagged with (I) means that it is an input variable, which has to be initialized before calling m1qn3; an argument flagged with (O) means that it is an output variable, which only has a meaning on return from m1qn3; and an argument flagged with (IO) is an inputoutput argument, which has to be initialized and that has a meaning after the call to m1qn3. Arguments of the type (O) and (IO) are generally modified by m1qn3 and therefore should not be Fortran-77 constants!

The solver

Here is the declaration of the solver m1qn3. subroutine m1qn3 (simul, prosca, ctonb, ctcab, n, x, f, g, dxmin, df1, epsg, normtype, impres, io, imode, omode, niter, nsim, iz, dz, ndz, reverse, indic, izs, rzs, dzs).

simul: In direct communication, simul is the name of the simulator inside m1qn3.

The simulator is the subroutine that computes the value of the function f and the value of its gradient g at the current iterate. When m1qn3 needs these values, it executes the instruction call simul (...)

Hence the actual name of the simulator must be declared external in the program calling m1qn3 and given as actual argument when calling m1qn3. For example, if the actual name of the simulator is mysimul, m1qn3 will be called by external mysimul ... call m1qn3 (mysimul, ...)

See section 3.2, for more information on the simulator, in particular on its arguments.

In reverse communication, a simulator is not needed, since the user can organize the computation of f and g in an arbitrary manner (this is one of the advantage of reverse communication, see section 4.2). The empty subroutine simul rc, provided in the standard distribution, should be used as actual argument in reverse communication (it is never used actually). The calling sequence is then external simul_rc ... call m1qn3 (simul_rc, ...) Section 4.2 tells you more on the use of reverse communication.

prosca: Calling name of the subroutine that computes the inner product u, v of two vectors u and v of R n . This subroutine is supposed to have the following declaration statement:

subroutine prosca (n, u, v, ps, izs, rzs, dzs).

Its actual name has to be declared external in the program calling m1qn3. The arguments n, izs, rzs and dzs have the same meaning as below. The arguments u and v are two double precision arrays of dimension n. The inner product of these vectors is supposed to be returned in the double precision variable ps.

ctonb (used only in DIS mode): Calling name of the subroutine that gives the coordinates (y i ) 1≤i≤n in the basis (e i ) 1≤i≤n of a vector whose coordinates in the basis (ê i ) 1≤i≤n are (x i ) 1≤i≤n (see section 2.1). This subroutine is supposed to have the following declaration statement:

subroutine ctonb (n, x, y, izs, rzs, dzs).

Its actual name has to be declared external in the program calling m1qn3.

The arguments n, izs, rzs and dzs have the same meaning as below. The arguments x and y are two double precision arrays of dimension n.

ctcab (used only in DIS mode): Calling name of the subroutine that does the operation reverse to the one done by ctonb, so that calling ctonb and next ctcab will not change the value of x: it gives the coordinates (x i ) 1≤i≤n in the basis (ê i ) 1≤i≤n of a vector whose coordinates in the basis (e i ) 1≤i≤n are (y i ) 1≤i≤n . It is supposed to have the following declaration statement:

subroutine ctcab (n, y, x, izs, rzs, dzs).

Its actual name has to be declared external in the program calling m1qn3. The arguments are similar to those of ctonb.

n (I): Positive integer variable. It gives the dimension n of the problem.

x (IO): Double precision array of dimension n. On entry, it is supposed to be the value x 1 of the initial point. On return, it is the value of the final point calculated by m1qn3.

f (IO): Double precision variable. On entry, it is supposed to be the value of the function f at the initial point x 1 . This value will be obtained by calling the simulator simul before calling m1qn3. On return, it is the value of f at the final point.

g (IO): Double precision array of dimension n. On entry, it is supposed to be the value of the gradient of f at the initial point x 1 . This value will be obtained by calling the simulator simul before calling m1qn3. On return with omode = 1 (see below), it is the value of the gradient of f at the final point. For other output modes, the value of g is undetermined.

dxmin (I): Positive double precision variable. This argument gives the resolution in x for the l ∞ norm: two points whose distance in the sup-norm is less than dxmin will be considered as indistinguishable by the optimizer (in fact by mlis3).

df1 (I): Positive double precision variable. This argument gives an estimation of the expected decrease in f during the first iteration. See section 4.3 for more details.

epsg (IO): Double precision variable, with value in ]0, 1[ on entry. On entry, epsg gives the precision of the stopping criterion that is based on the norm of the gradient. Indeed, the optimizer considers that the convergence has been obtained at x k and stops with omode = 1, if

ǫ k := |||g k ||| |||g 1 ||| < epsg,
where, g 1 and g k are the gradients of f at the initial and current points respectively and

= ( i v 2 i ) 1/2 ), -'sup' denotes the sup or ℓ ∞ norm (|||v||| = max i |v i |),
-'dfn' denotes the norm • associated with the scalar product defined in the user-supplied subroutine prosca (|||v||| = v, v 1/2 ). It is this norm that is used to print the gradient norms. impres (I): Integer variable that controls the outputs on channel io.

= 0: No print. ≥ 1: Initial and final printouts, error messages.

≥ 3: One line of printout per iteration that gives: the index k of the iteration going from the point x k to the point x k+1 , the number of time the simulator has been called, the value f (x k ) of the objective function and the directional derivative g k , d k . ≥ 4: Printouts from mlis3 during the line-search: see the write-up on mlis3 in modulopt library. ≥ 5: Some more printouts at the end of iteration k: information on the matrix update (the factor δ k for the SIS mode, see (2.3); the fitting factor and the look of the diagonal initial scaling matrix D k for the DIS mode), the value ǫ k+1 of the stopping criterion, and the angle between -g k and the search direction d k .

io (I): Integer variable that will be taken as the channel number for the outputs, i.e., outputs are written by: write (io,...) ... imode (I): Integer array of dimension 3 that specifies the input mode of m1qn3 and tunes its behavior. The following values are meaningful.

imode(1) determines the scaling mode of m1qn3. M1qn3 will run in DIS mode if imode(1) = 0, meaning that the matrix W k is formed by updating m times a dynamically updated diagonal matrix D k .

M1qn3 will run in SIS mode if imode(1) = 1, meaning that the matrix W k is formed by updating m times a variable multiple of the identity matrix.

The DIS mode is generally more efficient than the SIS mode and is recommended.

imode(2) specifies the starting mode of m1qn3.

A cold start is performed if imode(2) = 0: the first descent direction is then -g 1 . If you have not stored the appropriate information in iz and dz in a previous run, set imode(2) to 0.

A warm start is performed if imode(2) = 1: the first descent direction is -W 1 g 1 , where W 1 is the matrix formed from the (y, s) pairs (and the diagonal matrix, in DIS mode) stored in the double precision working zone dz. The integer working zone iz helps m1qn3 decoding the information in dz. This option can only be used if the appropriate data has been stored in iz and dz in a previous run.

See section 4.5 for more details.

imode(3) specifies in direct communication when the simulator has to be called with indic = 1 or, similarly in reverse communication, when m1qn3 returns to the calling subroutine with indic = 1. As explained in section 3.2, when called with indic = 1, the simulator can do anything but changing the value of its arguments n, x, f , and g (a similar rule must be respected by the calling program in reverse communication). The simulator can take this opportunity to do any regular computation, printout, or plot. When imode(3) = 0, the simulator is never called with indic = 1. When imode(3) > 0, the simulator is called with indic = 1 every imode(3) iteration(s), starting at iteration 1. Hence, when imode(3) = 1, the simulator is called every iteration, which is a way for the simulator of being informed when a new iteration starts (recall that there may be an unpredicted number of simulations per iteration).

omode (O): Integer variable that specifies the output mode of m1qn3. The following values are meaningful. = 0: The simulator asks to stop by returning the value indic = 0. = 1: This is the normal way of stopping for m1qn3: the test on the gradient is satisfied (see the meaning of epsg). = 2: One of the input arguments is not well initialized. This can be:

n ≤ 0, niter ≤ 0, nsim ≤ 0, dxmin ≤ 0.0 or epsg ∈ ]0, 1[, ndz < 5n + 1 (in SIS mode) or ndz < 6n + 1 (in DIS mode): not enough storage in memory, -the contents of iz is not correct for a warm restart, -the starting point is almost optimal (the norm of the initial gradient is less than 10 -20 ). = 3: The line-search is blocked on tmax = 10 20 (see section 4.4 and the documentation on mlis3 in modulopt library). = 4: The maximal number of iterations is reached. = 5: The maximal number of simulations is reached. = 6: Stop on dxmin during the line-search (see section 4.4). = 7: Either g, d is nonnegative or y, s is nonpositive (see section 4.4). For additional information and comments, see section 4. niter (IO): Positive integer variable. On entry, it gives the maximal number of iterations accepted. On return, niter is equal to the number of iterations really done.

nsim (IO): Positive integer variable. On entry, it gives the maximal number of simulations accepted. On return, nsim is equal to the number of simulations really done. In SIS mode, m1qn3 needs a working area formed of at least 3 vectors of dimension n (d k , g k and an auxiliary vector) and it needs for each update one scalar and two vectors: y y, s -1/2 and s y, s -1/2 (see section 2.1). Therefore, if m is the desired number of updates for forming the matrix W k , it is necessary to have: ndz ≥ 3n + m(2n + 1).

In fact, the number m of updates is determined by m1qn3, which takes

m = ndz -3n 2n + 1 , (3.1)
where, ⌊•⌋ refers to the floor operator (⌊x⌋ = i, when i ≤ x < i + 1 and i ∈ N).

This number has to be greater than or equal to 1. Therefore, if ndz is less than 5n + 1, m1qn3 stops with omode = 2.

On the other hand, in DIS mode, m1qn3 needs an additional vector of dimension n for storing D k . So, take ndz ≥ 4n + m(2n + 1) and m ≥ 1. Then, m1qn3 will set:

m = ndz -4n 2n + 1 . (3.2)
The above computation of m from ndz is irrelevant when m1qn3 is used as described in section 4.6. reverse (IO): Integer variable that is used to specify whether direct or reverse communication is desired and, in the latter case, to communicate with the call loop (see section 4.2 and figure 4.2 for the program flow chart on reverse communication). On entry in m1qn3, reverse < 0: implies that m1qn3 will stop immediately using the instruction stop, to prevent entering an infinite call loop in reverse communication, due to the fact that the calling program has not left the call loop when m1qn3 returns a negative value, = 0: indicates that m1qn3 has to work in direct communication, = 1: indicates that m1qn3 has to work in reverse communication. On return from m1qn3 in reverse communication, reverse < 0: when m1qn3 has terminated, in which case the call loop must be interrupted, = 1: the call loop must be pursued. See sections 4.1 and 4.2 for more details. On entry in m1qn3, the calling program can send a message to m1qn3, using the following values of indic: < 0: the computation of f and g required by m1qn3 on its last return was not possible at the given x; in that case, m1qn3 divides the step-size by 10, before calling the simulator again; this feature can be used when implicit constraints are present, i.e., strict inequality constraints or inequalities that are known to be inactive at the solution; by no way this feature can handle inequality constraints that are active (satisfied with equality) at the solution; = 0: m1qn3 has to stop, for example because some events that m1qn3 cannot understand (not in the field of optimization) has occurred; > 0: the required computation has been done. On return from m1qn3, the optimization solver sends a message depending on the values of indic:

= 1: means that the calling program can do anything except changing the values of indic, n, x, f, and g; this value of indic is used by m1qn3 every imode(3) iteration(s), when imode(3) > 0, and never, when imode(3) = 0; = 4: means that the calling program has to compute f and g, to put them in f and g, and to call back m1qn3. Note that these values of indic have the same meanings as those used in the communication with the simulator (see section 3.2). izs (I), rzs (I), dzs (I) (used only in direct communication): Addresses of integer, real, and double precision working areas respectively. M1qn3 does not use them but passes their addresses to the subroutines simul, prosca, ctonb, and ctcab so that they can be used in these subroutines.

The simulator

In direct communication, the simulator is the subroutine that computes the value of the function f and the value of its gradient g at a given arbitrary point x. M1qn3 assumes that the subroutine is defined by the following statement:

subroutine simul (indic, n, x, f, g, izs, rzs, dzs).

The name simul of that subroutine can be changed. Its actual name must be declared external in the program calling m1qn3 and must be given as the first argument of m1qn3 (see section 3.1). Here is a description of the arguments of the simulator.

indic (IO) integer value that monitors the communication between m1qn3 and the simulator.

On entry, m1qn3 sends a message to the simulator by using one of the following values of indic: = 1: means that the simulator can do anything except changing the values of indic, n, x, f, and g; this value of indic is used by m1qn3 every imode(3) iteration(s), when imode(3) > 0, and never, when imode(3) = 0; = 4: means that the simulator has to compute f and g, and to put their values in f and g, respectively. On return, the simulator can send a message to m1qn3 by using one of the following values of indic: < 0: the computation of f and g required by m1qn3 on its last return was not possible at the given x; in that case, m1qn3 divides the step-size by 10, before calling the simulator again; this feature can be used when implicit constraints are present, i.e., strict inequality constraints or inequalities that are known to be inactive at the solution; by no way this feature can handle inequality constraints that are active (satisfied with equality) at the solution; = 0: m1qn3 has to stop, for example because some events that m1qn3 cannot understand (not in the field of optimization) has occurred; > 0: the required computation has been done. The instructions preceding the call to m1qn3 will include the following items:

1. the declaration in external of the following names: the name of the simulator (named simul by m1qn3), the name of the subroutine doing the inner product (named prosca by m1qn3) and the name of the subroutines doing the change of basis (named ctonb and ctcab by m1qn3);

2. the calculation of the starting point x 1 ;

3. the call to the simulator to compute f and g at x 1 ;

4. the calculation of df1 and initialization of dxmin, epsg, impres, io, imode, niter, nsim, ndz, and reverse (set to 0); 5. the reading of values for iz and part of dz from the appropriate area or device, if a warm restart is performed (imode(2) = 1); 6. the call to m1qn3. the dashed blue boxes represent the parts of the program that have to be written by the user of the optimization software, while the single red box represents the m1qn3 solver. In this scheme, the user program must contain a call loop, in a part of the program sometimes called the driver. At each turn of this loop, m1qn3 is called and is asked to do a single iteration. Once the optimization routine needs information on the problem, it returns to the calling program, requiring the computation of the function and gradient by setting indic = 4 (and reverse = 1). Then the user program computes f and g, loops, and calls back m1qn3, using also indic to inform the optimizer on the realization of the required computation or to ask it to stop1 . When the run is finished (for some reason specified by omode), m1qn3 returns with reverse = -1, in which case the call loop must be interrupted (if m1qn3 is called back with reverse < 0, the solver abruptly stops using the instruction stop to avoid a possible infinite loop).

Calling sequence in reverse communication

The two main advantages of reverse communication are that there is no need to write a simulator and that it is easier to call m1qn3 from a language different from Fortran-77.

The instructions preceding the driver include the following items:

1. the declaration in external of the following names: the name of the empty simulator simul rc, the name of the subroutine doing the inner product (named prosca by m1qn3) and the name of the subroutines doing the change of basis (named ctonb and ctcab by m1qn3);

2. the calculation of the starting point x 1 ;

3. the computation of f and g at x 1 ;

4. the calculation of df1 and initialization of dxmin, epsg, impres, io, imode, niter, nsim, and ndz;

5. the reading of values for iz and part of dz from the appropriate area or device, if a warm restart is performed (imode(2) = 1); 6. setting reverse = 1 to indicate that reverse commnunication is required.

In Fortran-77, the call loop in the driver can be written as follows:

external simul_rc, ... ... reverse = 1 1 continue call m1qn3 (simul_rc,...,x,f,g,...,reverse,indic,...) if (reverse.lt.0) goto 2 call simul (indic,n,x,f,g) goto 1 2 continue Above, we have assumed that the computation of f and g can be done in a subroutine called simul.

More on some arguments

The value of df1 is used to obtain an estimate of the step-size at the first iteration, which is taken equal to 2(df1)/ g 1 2 . A good value for df1 can be the total decrease in f up to its minimum or some fraction of this value. As a rule, a too large value is not dangerous: the step-size is reduced by mlis3 in a few internal iterations. On the other hand, an excessively small value of df1 could, due to rounding error, force m1qn3 to stop with omode = 6 at the first iteration.

If the storage is sufficiently large, it is a good idea to choose ndz large enough to have a large number m of updates: see Formulae (3.1) and (3.2). It seems reasonable to take m ≥ 2. However, numerical experiments (see [START_REF] Gilbert | Some numerical experiments with variable storage quasi-Newton algorithms[END_REF] and [START_REF] Liu | On the limited memory BFGS method for large scale optimization[END_REF]) have shown that it is difficult to take full advantage of a large number of updates, especially since the cpu time is increasing with m. A good compromise should be obtained by taking m between 5 and 10.

More on some output modes

The output mode omode = 3 is very unlikely. On the one hand, the matrix W k is usually well scaled when the updates start with the identity matrix multiplied by the factor given by (2.3). And, on the other hand, as soon as the number of updates is larger than 2 or 3 and k ≥ m, the BFGS updates also scale W k . Therefore, the unit step-size is usually accepted by mlis3.

The output mode omode = 6 can have various origins.

• It can come from a mistake in the calculation of the gradient. This is likely to be the case when the solver stops with that output mode after very few iterations.

• If the number of step-size trials at the last iteration is small, this can mean that dxmin has been chosen too large. In that case, decreasing dxmin should have the clear effect of increasing the number of iterations.

• It can also come from rounding error in the simulator. The precision on the optimal solution that can be reached by the solver depends indeed on the precision of the simulator (the amount of rounding errors made in the computation).

It is worth noting that an inaccurate gradient may lead to a direction of search d k that is no longer a descent direction of f at x k . Indeed, m1qn3 allows itself to build search directions d k making with -g k an angle close to 90 • . This fact is essential to have good convergence properties. Now, if the gradient is roughly calculated, the angle between the true negative-gradient and d k could be larger than 90 • , in which case d k will not be a descent direction and an output with omode = 6 will occur. Therefore, it is important to compute the gradient carefully; in particular, a finite-difference gradient is likely to yield less precise solution than a calculated gradient.

The output mode omode = 7 is very unlikely. It can only be due to rounding error. nsim,iz,dz,ndz,reverse,indic,izs,rzs,dzs)

In case of reverse communication, the statement reverse = 1 must be inserted before the second call to m1qn3.

Observe that, for the second run, the required relative gradient norm, specified by epsg, has been adapted to the one obtained at the end of the first run (we assume there that epsg = 0 at the end of the first run), so that the two runs (with and without a warm restart) will stop at the same iteration. Now, if the iterates generated in the two runs are identical, it is normal, however, to get non identical printouts since the relative gradients differ.

Note that, in the last case, we have made no restrictions on the number m of updates. In particular, if m > 5, the number of pairs stored in dz by the first call is less than the number of updates. M1qn3 can manage this situation correctly. On the other hand, it is imperative that the run of m1qn3 filling the vectors iz and dz should be done with the same number of updates and the same scaling mode (SIS or DIS mode) as the current run. Otherwise, m1qn3 terminates with omode = 2 and the error message: inconsistent warm restart 4.6 Usage for very large scale problems Some dispositions have been taken so that, if you have experience with the code and understand what it does, you can adapt it to the optimization of very large scale problems. The idea is to allow you to store yourself the pairs (2.1) on an auxiliary memory (for instance on a disk or any other device) instead of keeping them in core memory. To achieve this goal, you have to write and use new versions of the subroutines mupdts and ystbl (dystbl in double precision), which are provided in the standard distribution. There are two reasons why we do not make this usage of m1qn3 more friendly. First, it is definitely non standard and it should be useful only in case of a huge number of variables. Secondly, the communication with external device is, in any case, installation dependent, which requires writing appropriate subroutines for this.

Here are the declaration statements of the two subroutines: subroutine mupdts (sscale, inmemo, n, m, ndz), subroutine ystbl (store, ybar, sbar, n, j).

In mupdts, sscale, n and ndz are input arguments and inmemo and m are output arguments. In ystbl, store, n and j are input arguments and ybar and sbar are input/output arguments. The rest of this section describes what these subroutines should look like.

In the standard distribution, subroutine mupdts computes the number m of updates from formula (3.1) (if sscale = .true., i.e., SIS mode) or (3.2) (if sscale = .false., i.e., DIS mode) and informs the optimizer that the pairs (2.1) are to be stored in core memory (inmemo = .true.). If you want to store yourself the pairs (2.1), the length ndz of the working area dz will be at least 6n + m, when calling m1qn3. Note that in this case, the number of updates can no longer be deduced from the value of ndz. Therefore, your version of mupdts will set the value of m and will set inmemo = .false.. Subroutine ystbl is device dependent and will be written accordingly. Its role is to store (when store = .true.) or to fetch (when store = .false.) a pair of vectors (y, s) in or from a certain logical position j in an auxiliary memory containing m (= number of updates specified by the subroutine mupdts) such positions. Position j will always be between 1 and m, inclusive. A pair is always stored before being fetched. If a pair has been stored at position j, m1qn3 expects that a fetch with the same index j will correspond to the same pair. There is however no link between the logical position specified by the index j and the index i of the pair (y i , s i ). Note that subroutine ystbl is not used when inmemo = .true.. Just as an example, suppose that you want to take care yourself of the storage of the pairs (2.1) and to use the core memory as auxiliary memory. Suppose also that m = 10 updates are desired. Then subroutine mupdts should look like this: subroutine mupdts (sscale,inmemo,n,m,ndz) logical sscale,inmemo integer n,m,ndz m = 10 inmemo = .false. return end On the other hand, if n = 1000 and if the 10 pairs of vectors (y, s) are stored in the arrays yarray and sarray, then subroutine ystbl should look like this: subroutine ystbl (store,ybar,sbar,n,j) c logical store integer n,j double precision ybar(n),sbar(n) c integer i double precision yarray,sarray common /ysarray/yarray(1000,10),sarray(1000,10) c if (store) then do i = 1,n yarray(i,j) = ybar(i) sarray(i,j) = sbar(i) end do else do i = 1,n ybar(i) = yarray(i,j) sbar(i) = sarray(i,j) end do end if c return end

This subroutine stores (when store = .true.) or fetches (when store = .false.) a pair of vectors (y, s) = (ybar, sbar) of dimension n at or from position j in the "auxiliary" memory represented by the common block /ysarray/. Variables store, n and j are input variables and ybar and sbar are input/output variables, depending on the value of store.

For a warm restart (imode(2) = 1, see section 4.5), you need to give m1qn3 correct values in the 5 locations of iz and, for the DIS mode, in the first n locations in dz, which is the normal place for the diagonal preconditioner (for the SIS mode, the values in dz are irrelevant). Appropriate saves will be done accordingly at the end of the previous run.

5 Frequently asked questions 5.1 Preconditioning the problem

Dealing with fixed variables

Suppose that I have a simulator for an optimization problem defined in the variables x = (x 1 , . . . , x n ) and that I want to solve the same problem, but now with some of the variables fixed. How can I proceed?

We see two possibilities. Below, we have denoted by simul the name of the original simulator, defined in all the variables x 1 , . . . , x n .

1. A solution that always works is to write a cover simulator, say simul cover, which is called by m1qn3 and that calls the original simulator simul.

To be specific, if nfree is the number of free variables and if their indices in {1, . . . , n} are given by free(i), for i = 1, . . . , nfree, the cover simulator can look like below.

subroutine simul_cover & (indic, nfree, xfree, f, gfree, izs, rzs, dzs) ... do i = 1,nfree x(free(i)) = xfree(i) end do call simul (indic, n, x, f, g, izs, rzs, dzs) do i = 1,nfree gfree(i) = g(free(i)) end do return end 2. If the inner product used to compute the gradient (specified by the subroutine prosca) is the Euclidean inner product (2.2), then one can proceed by zeroing the components of the gradient corresponding to the fixed variables. Then, m1qn3 will not modified the fixed variables and the free variables will be modified like in solution 1. This technique is valid in both scaling modes of m1qn3 (imode(1) = 0 or 1). Be careful with a warm restart (imode(2) = 1): the (y, s) pairs stored during the previous run must have been obtained with the same fixed variables.

To be more specific, if nfixed is the number of fixed variables and if their indices in {1, . . . , n} are given by fixed(i), for i = 1, . . . , nfixed, the modifications to bring to the original simulator simul could be as below.

subroutine

  iz (IO): Integer array of dimension 5. It is the address of a working array for m1qn3. dz (IO): Double precision array of dimension ndz. It is the address of a working array for m1qn3. ndz (I): Positive integer variable. It gives the dimension of the working area dz (see also section 4.3).

  indic (IO) In reverse communication, m1qn3 returns to the calling program each time it needs information on the problem (function/gradient values). The calling program must then compute what is required by m1qn3 (sometimes using a simulator, like in direct communication) and call back the solver in what can be viewed as a call loop (see section 4.2 and figure 4.2 for the program flow chart on reverse communication). The parameter indic normalizes the communication between m1qn3 and the part of the program computing function and gradient.

xFigure 4 . 1 :

 41 Figure 4.1: Direct communication flow chart

With a reverse communicationFigure 4 . 2 :

 42 Figure 4.2: Reverse communication flow chart

  simul (indic, n, x, f, g, izs, rzs, dzs) .

	imode, 8	inner product, 1-4
	impres, 7	Euclidean, 3, 20
	inmemo, 18	prosca, 6, 13, 14
	io, 8	IO (input-output argument), 5
	iz, 9 izs, 11, 12	norm, 3, 7
	m, 18 n, 7	O (output argument), 5
	ndz, 9, 18 niter, 9 normtype, 7	scaling mode, 2, 3, 8 diagonal initial scaling (DIS), 2, 3 scalar initial scaling (SIS), 2, 3
	nsim, 9	simulator, 5, 11, 11-12
	omode, 9	SIS, see scaling mode
	prosca, 6	stop, 10
	reverse, 10	subroutine
	rzs, 11, 12	ctcab, 2, 4, 6-7
	simul, 5	ctcabe, 2
	sscale, 18	ctonb, 2, 4, 6
	store, 18 .. x, 7, 12 do i = 1,nfixed g(fixed(i)) = 0.d0 BFGS formula, 3, 4 end do return end call loop, 10, 13, 14 cold start, 4, 8, 16-17 communication direct, 5, 10, 12-13 reverse, 6, 10, 13-14	ctonbe, 2 dd, 2 dds, 2 ecube, 2 euclid, 2 m1qn3a, 2 mlis3, 2 mupdts, 2, 17-19 prosca, 2 simul, 2, 11, 20
	DIS, see scaling mode	simul cover, 19
		simul rc, 2, 6
	fixed variable, 19-20	ystbl, 2, 17-19
	frequently asked questions (FAQ), 19-20	warm restart, 4, 8, 9, 13, 14, 16-17, 19,
	I (input argument), 5	20
	Index
	argument	dz, 9
	ctcab, 6	dzs, 11, 12
	ctonb, 6	epsg, 7
	df1, 7	f, 7, 12
	dxmin, 7	g, 7, 12

When the user decides to stop the iterations, he/she could bypass an additional call to m1qn3 with indic = 0. This, however, will spoil a subsequent call to m1qn3 (for a cold start or a warm restart). It is therefore unadvisable.

Cold start and warm restart

By setting the input variable imode appropriately, you can get a cold start (imode(2) = 0) or a warm restart (imode(2) = 1) of m1qn3. A cold start of the optimizer is the usual way of using it. For example, the first time you use m1qn3, you can only call it with imode(2) = 0. Now, when several optimization runs are done sequentially for similar problems, it can be useful to use the information stored in the pairs (y, s) and possibly in the preconditioning diagonal matrix D of a previous run to improve the efficiency of the first iterations of the current run. This is what we mean by doing a warm restart.

All the information needed by m1qn3 to build the preconditioner W k is contained in the working zones iz and dz. As we have just said, normally, the relevant information is placed in these vectors by a previous call to m1qn3. This can be done in the same run of a main program (sequence of calls to m1qn3) or by storing the data on an auxiliary device between two runs. Note that the 5 integers in iz are crucial for a correct warm restart, but that only a part of the elements in dz are meaningful to this respect, namely the first 2nm elements in SIS mode and the first n(2m + 1) elements in DIS mode.

For example, the following piece of program epsg = 1.e-4 imode(1) = 0 imode(2) = 0 niter = 100 nsim = 100 call m1qn3 (mysimul,euclid,ctonbe,ctcabe,n,x,f,g,dxmin, df1,epsg,normtype,impres,io,imode,omode,niter, nsim,iz,dz,ndz,reverse,indic,izs,rzs,dzs)

will generate the same iterates as the next piece of program, where the optimization phase is split in two parts, 5 iterations for the first call and 95 for the second, which starts at the situation reached by the first call (imode(2) = 1): epsg = 1.e-4 imode(1) = 0 imode(2) = 0 niter = 5 nsim = 100 call m1qn3 (mysimul,euclid,ctonbe,ctcabe,n,x,f,g,dxmin, df1,epsg,normtype,impres,io,imode,omode,niter, nsim,iz,dz,ndz,reverse,indic,izs,rzs,dzs) epsg = 1.e-4/epsg imode(2) = 1 niter = 95 nsim = 100 call m1qn3 (mysimul,euclid,ctonbe,ctcabe,n,x,f,g,dxmin, df1,epsg,normtype,impres,io,imode,omode,niter,


Hitchhiker's Guide to a Practical Automated TFHE Parameter Setup for Custom Applications Jakub Klemsa jakub.klemsa@eurecom.fr Motivation TFHE cipher is instantiated with 8 parameters:

• determine security level and plaintext/evaluation error rate ⇒ parameters are application-specific,

• vast impact on performance ⇒ optimization problem, ⇒ need for a tool for TFHE parameter setup.

Given application specifications:

+ fully automated TFHE parameter setup, + optimized for best performance, + unified approach ⇒ allows to compare different parameters:

• same target security & error rates, aim for best parameters, • e.g., different homomorphization in digit-based arithmetics [START_REF] Klemsa | Parallel Operations over TFHE-Encrypted Multi-Digit Integers[END_REF].

Application Specifications

Specify application needs (security & error rates) by three parameters:

A. bit-security level denoted λ; B. requested cleartext space bit-precision denoted π; and C. bound on the number of homomorphic additions before the sample gets bootstrapped, denoted 2 2∆ , referred quadratic weights;

⇒ input parameters for our TFHE parameter setup tool.

A. Bit-Security Level λ Observation 1. At fixed security λ, the logarithm of stddev of LWE noise (den. α), is roughly linear in the LWE dimension n (with factor den. s λ ):

-log 2 (α) ≈ s λ • n; (1)
cf. Figure 1. Due to the collision attack, the relation is limited to n ≥ 2λ, also the behavior changes forlog 2 (α) > τ with τ the torus precision. [START_REF] Albrecht | LWE Estimator[END_REF]. For λ = 128 bits, s λ ≈ 0.0235.

B. Cleartext Space Bit-Precision π

Before selecting an appropriate π, bare in mind: i. complexity of TFHE bootstrapping is roughly exponential in π,

• practical times for up to π ≈ 6 bits (cf. Figure 2), ii. bootstrapping Look-Up Table (LUT) is inherently negacyclic:

LU T (2 π-1 + m) = -LU T (m), m ∈ 0, 2 π-1 . (2)
C. Quadratic weights 2 2∆

Observation 2. LWE noises accumulate with each homomorphic addition: for indep. TLWE samples c i with equal noise variance denoted V 0 :

Var Err

w i • c i = w 2 i 2 2∆
• Var(Err(c i ))

V 0 , w i ∈ Z. (3)
To ensure correct LUT evaluation during bootstrapping:

• bound on 2 2∆ before a sample gets bootstrapped (refresh noise),

• 2 2∆ . . . sum of squared weights; ∆ . . . bits of stddev of addit'l noise.

Parameter Restrictions

Goal: bound the noise of a fresh(ly bootstrapped) sample, s.t.

i. limited number of additions can be performed (cf. 2 2∆ ); and ii. the noise can be refreshed correctly during bootstrapping.

2 π cleartext values ⇒ max error ! < 1 /2 π+1 . . . by 3σ-rule: V max ≤ 2 2∆ V 0 + V round ! ≤ 1 3 2 • 2 2π+2 ,
where

V round = n + 1 48N 2 . ( 4
)
V 0 depends on implementation, for plain TFHE [START_REF] Chillotti | TFHE: fast fully homomorphic encryption over the torus[END_REF]:

V 0 ≤ 2nlN 2 2(γ-1) V BK (N ) (♥) + n(1 + N )2 -2(γl+1) (♦) + Var Err(u, v) = 0 + + tN 2 2(κ-1) V KS (n) (♣) + 2 -2(κt+1) N (♠) . (5)
To derive good TFHE parameters, we need to:

i. satisfy the bound (4) (error budget), using (5); and ii. check their quality in terms of bootstrapping time.

Experimental Results

Z/16Z

Demo: π = 6, 2 2∆ = 2
Repo with Orig. param's [START_REF] Zama | Demo Z/8Z[END_REF] New param's exp. code

N, n ; γ, l 2

Figure 2 :

 2 Figure 2: Bootstrapping times of TFHE parameters with η C < 100% for various scenarios, chosen automatically. The width of the bars represents l ∈ 1, 4 . Hatched bars represent incorrect results, presumably due to log(α) < -64 being out of Concrete's v0.1.11 supported range.

  LWE bit-security level λ as estimated by the LWE Estimator by Albrecht et al.

										576
	20 -log 2 (α)	Est. bit-security of LWE							192 256 320 384 448 512
	10									128
	0	256	512	768	1024	1280	1536	1792	2048	64
				LWE dimension n
	Figure 1:

Table 1 :

 1 Original and newly identified TFHE parameters for the Z/16Z demo[START_REF] Zama | Demo Z/8Z[END_REF] with π = 6 bits. 500 runs with Concrete[START_REF] Zama | [END_REF] on Intel Core i7-7800X. η C and η m stand resp. for the usage of the 3σ max error budget as calculated by Concrete and as measured after decryption. Experimental code at https://gitlab.eurecom.fr/fakub/tfhe-param-testing.

	048, 750 ; 7, 3 2 048, 766 ; 21, 1

Acknowledgements

This work was supported by the MESRI-BMBF project UPCARE (ANR-20-CYAL-0003-01).