Chebyshev methods for optimization of large systems on stiff target functionals with weakly filled structured hessian matrices

Abstract

We analyze matrix gradient methods for optimizing large systems of arbitrary nature according to functions with a special character of Hessian matrices. It is assumed that the Hessian matrix is sparse and structured (nonzero elements occupy fixed positions). The last condition is for large systems optimization consisting of interconnected subsystems of a lower dimension. We suggest that the use of standard Newton-type methods is difficult because of the possible sign uncertainty of the Hessian matrices (non-convexity) and the need to store full-size high-order matrices. In the procedures under consideration, memory savings are achieved by storing only sparse matrices in packaged form.

Keywords

• gradient methods
• relaxation functions
• non-convex problems
• stiff functionals

1 Introduction

To solve the problem of unconditional minimization

a class of matrix gradient methods is considered, k Gk , hk  J   xk  , hk  R1, (1) where Gk  J   xk  , H k  matrix function Gk .

This class of methods includes some particular cases: classical procedures as gradient descent methods, Levenberg

Marquardt methods, Newton’s methods, methods with exponent relaxation (ER-methods).

We present approach based on the concept of the relaxation function [ 1 ] for constructing and analyzing nontraditional gradient methods with the Chebyshev relaxation function. The main assumptions used in the construction of this class of methods are:  high dimension of the argument vector x (n>>1) and the undesirability of storing a full-size matrix (nxn) in the computer's memory;  a high stiffness degree [ 1, 2 ] of the functional J(x) in a wide range of the argument variation;  convexity of the minimized functional J(x) is guaranteed only in the neighborhood of the minimum point.

If the stiffness degree of the optimality criterion J(x) is sufficiently high, then standard computational tools are not very effective due to the reasons stated in [ 1, 2 ]. Newton-type methods, as it is known [ 3-5 ], are not intended for solving nonconvex problems. Moreover, they lose their effectiveness in conditions of high stiffness. The use of algorithms from the class of well-known Markuardt-Levenberg methods with a high stiffness degree of the functional J(x) is associated with the complexity of the algorithmic selection of the corresponding regularization parameter.

Most often in this situation, it is recommended to use different non-matrix forms of the conjugate gradient method (SG). However, further it will be shown that in the class of matrix gradient schemes (1) there are algorithms that are more efficient for the considered problems than the SG methods. 2 Multiparameter Optimization Tasks Large systems we define as systems described by models with a large number of controlled parameters.

An optimized system may be represented as a set of smaller interconnected subsystems (Fig. 1). Figure 1 – The complex of interrelated subsystems

Y is the output characteristics of the subsystems. The requirements for the output parameters of the system (specifications) are given in the form of a system of inequalities:

yj(xj, xq)  tj, j  [1: q – 1]; yq(xq)  tq, where j is nj-dimensional local vector of controlled parameters; xq – nq-dimensional vector of controlled parameters, affecting all q output parameters and interconnecting individual subsystems of the system being optimized. The dimension of the full vector of controlled parameters x  [x1, x2, , xq] is equal to n   n .

i q i1 q

Using the known technique of reducing the solution of inequality systems to optimization tasks, we obtain the target functional of the form

J  x   U j  x j , xq   min, x  Rn .

If we have the system of inequalities yi ( x)  ti , it is equivalent to this system:

zi ( x)  ti  yi  x   0

Corresponding optimization task:

J2 ( x)  max(exp( zi ( x)))  min

i x The last task with a sufficiently large  is equivalent to the following:

m J3 ( x)   exp( zi ( x))  min

i1 x This is the main target functional in solving systems of inequalities. It is non-convex and has the form (4).

Functionals (4) also arise in other formulations of optimization tasks of real systems, therefore the problem (4) has a rather general character.

Further, we will consider methods for problem solving (4) with the following additional assumptions:  solving analyzing task of the optimized system requires significant computational costs. Therefore, in the optimization process, it is required to minimize the number of references to the calculation of values J(x);  fill coefficient  of matrix G(x) = J(x) is small. We can usually assume  ~ 1/q.

It is easy to establish that the structure of the matrix G (x) does not depend on the point x: 0

G  x    0 G11     Gq1

G22 Gq2

G33 Gq3  

G1q 

 G2q 

 G3q  .

 

 Gqq  The submatrices Gij have dimensions ninj, and the total number of nonzero elements is q q1 i1 i1 Thus, taking into account the symmetry of the matrix G (x), it is necessary to store in the computer’s memory q1  ni2  2nq  ni . q  ni2  ni  2  nq  ni i1 i1 nonzero elements. The necessary information about the storage schemes of sparse matrices are widely presented in the literature describing large data [ 6 ]. 3 Methods with Chebyshev relaxation functions We take the eigenvalues i(Gk)  [– m, M], M  m  0. According to above assumptions and requirements for relaxation functions formulated in [ 1 ], the most rational method should have a relaxation function R (), the values of which sharply decrease from R = 1 at  = 0, remaining small throughout the entire range [0, M]. And opposite, if <0, the function R () should increase intensively. In addition, the matrix function H corresponding to R () must be constructed without matrix multiplications in order to preserve the sparsity property of the matrix Gk  J(xk). Основной текст должен быть написан в одну колонку. Размер шрифта – 10, межстрочный интервал – одинарный. Рекомендуется использовать шрифт Times New Roman.

We show that as such R () with accuracy up to multiplier, offset Chebyshev polynomials of the 2nd kind Ps() can be used, satisfying the known recurrence relations:

P1()  1, P2()  2(1 – 2); Ps + 1()  2(1 – 2)Ps() – Ps – 1().

(5)

Dependency graph Ps() /s for some s is shown in Fig. 2. Supposing R()  PL() /L with a sufficiently large value of L, we get an arbitrarily fast relaxation of any addend in the presentation of approximating paraboloid [ 1 ] f  xk 1   1 n i2,k i R2  i , (6) where

2 i1 xk   i,k ui n i1 is a decomposition of the current vector into eigenvectors of the matrix Gk .

This statement follows from the known fact of uniform convergence of the sequence {Ps() /s} to zero when s   in the open space (0, 1). Further, we assume that the Gk matrix eigenvalues are normalized to the interval (0, 1). In this case, it suffices to consider the matrix G / || G || instead of the matrix G, and the vector g / || G || instead of the gradient vector g.

Correlating to the adopted R (), H () dependence has the form

Н()  [1 – R()] /  [1 – PL() /L] /.

The methods design (1) directly with function (7) is possible, but it makes it necessary to solve at each step k large linear systems of equations with sparse matrix. Below it is shown that there are more effective implementation techniques.

According to (7) it follows that H () is a polynomial of L–2 degree, while R () has L–1 degree. Therefore, to implement the matrix gradient method with the indicated function H (), generally speaking, there is no need to solve linear systems. The method will be as follows.

xk 1  xk   1E  2Gk  ...  L1GkL2  g k  xk  H Gk  g k . (7) (8) from (5) we can get the recurrence relation

Consequently, we have

The implementation of method (8) can be based on the techniques of calculating the coefficients i for various L degrees. In this case, the number L should be chosen from the condition of the most rapid decrease of J(x). An alternative, more economical approach based on other considerations is discussed below.

For function:

H s     1   2  ...   s 1 s 2 , s  2, 3, ... (s + 1)Hs + 1  2s(1 – 2)Hs – (s –1)Hs – 1 + 4s; H1  0, H2  2, s  [2: L – 1]. (9) or xk 1  s  1  xk  H s1g k  xk  s  2 : L  1 2s s  1

s  1  E  2G  H s g k  k

H s1g k  4s s  1

g k , In the left part of the spectrum (  0) we have

R     1  R  0 ,

s s therefore, the values of the derivatives R′s(0) in the last row of the table characterize the relaxation multipliers for the negative terms in (6). Calculation of derivatives R′s(0) can be performed on the basis of the following recurrence relations:

P1  0, P2  4; Ps1  2Ps  4s  Ps1; RL  0  PL L .

s, s values for s  8 (  0) can be calculated using the asymptotic formula when Rs  0,22.

The relation (11) is obtained from the following representation of Chebyshev polynomials s1  xk 1  s  1  xk 

 E  2Gk  s  1  0, 2  2g k , s  2 : L 1.

2s xk + 1[s] is s-е approximation to vector xk + 1  xk + 1[L].

Thus, with a fixed quadratic approximation f(x) of the functional J(x) in the neighborhood of x  xk, we have the opportunity to move from Рs to Ps + 1 due to one multiplication of the matrix E – 2Gk by the vector, fully using the sparsity property of the matrix Gk and without additional gradient calculations. The efficiency of algorithm (9) with large values of the stiffness coefficient [ 1,2 ] is determined by the relaxation factors for small eigenvalus of the matrix Gk. Consider the positive part of the spectrum (  0), which is especially important in the neighbourhood of the optimum, where the matrix G(x) is positively defined. The main advantage of the method with Rs()  Ps() /s is that already at small s there is a noticeable suppression of the addends from (5) in a wide range of  values . Below there are the Rs values for the internal maximum of Rs() and the boundaries of the ranges s    s, where |Rs()|  Rs: s  3 4 5 6 7 8 Rs s s -R′s (0)     0,333 0,147 0,853 5,30 0,272 0,092 0,908 10,0 0,250 0,061 0,939 16,0 0,239 0,044 0,956 23,3 0,233 0,033 0,967 32,0 0,230 0,025 0,975 42,0 (10) (11) (12) кр  x к2р  6, 523;   кр     x кр   0, 22.

Thus, if we suppose кр  4L2кр, we get the following statement: for the smallest (positive) eigenvalue m the inequality is fulfilled

Indeed, for sufficiently small  we have: If we consider x = sqrt(), we get

We have

where it means, if for all   m we’ll have

From (12) follows (11).

s  1,63/s2, s  1 – s; P    

L sin L L sin  2 

2 ,   sin

, ,   [ 0,1 ].

P        

L sin

 

2 ,   4L .

F()  (x)  sin x /x.

F()  /F(кр) when   кр,

min  4L2m  кр  6,523, m  6,523/ (4L2)  1,63/ L2, |RL()|  0,22.

The enlarged scheme of the algorithm based on the relation (10) can be implemented using the following sequence of steps. It is assumed that all task variables are properly normalized. We also suppose that the variables are numbered in some optimal way, ensuring efficient storage of the sparse matrix (E – Gk ) in the computer's memory.

4 RELCH Algorithm

Step 1. Set the starting point x; calculate J : J(x); set L, which determines the number of recalculations using the formula (10) (about the prior choice of L, see below).

Step 2. Calculate g : J(x), G : J(x); lay g : g/ ||G||, G : G/ ||G||;  :  1.

Step 3. According to the formula (10) build L ; put x t : x  L .

Step 4. Calculate Jt : J(xt). If Jt  J, go to step 5, otherwise go to step 6.

Step 5. Put  : /2, x t : x  L and go to step 4.

Step 6. Put x: xt, J : Jt and go to step 2.

The criterion of the end of the process is not specified here. As a rule, the calculations end when the specified number of calculations of the functional has been exhausted or when the algorithm is explicitly stopped. The number of recalculations L by the formula (10) is a parameter set by the user. According to (11), it is initially advisable to assume L 

1, 63 L  1, 3  where  - assessment of the ravine degree of the minimized functional. With this choice of L, the relaxation multipliers in the positive part of the spectrum will be guaranteed to be less than 0.23. When designing algorithmic methods for L tasks, it is necessary to take into account that the sequence {Js}, where Js  J  x k  s  will not decrease monotonically with s  . In step 5 of the algorithm, the regulation of the advance vector norm was applied in order to prevent the local quadratic model of the functional from leaving the space of validity.

5 Convergence Characteristics

We give an estimate of the efficiency of the method (10) in comparison with the methods of conjugate gradients (SG methods). For tasks of large dimension (when the number of iterations is less than the dimension), one can guarantee the convergence of SG-methods only with the speed of a geometric progression even for strongly convex quadratic functionals.

We’ll consider the case

f(x)  1/2Gx, x, G  0 and estimate the rate of convergence of the SG method to the extreme point x = 0.

The iteration xk, obtained by the SG method can be represented as

xk  (E + c1G + c2G2 +  + ckGk)x0  Pk(G)x0, where Pk(G) - matrix polynomial of degree k. Moreover, it follows from the properties of the SG method that the coefficients с1, , ck of the polynomial Рk(G) at each iteration take such values to minimize the value f(xk), which differs only by a multiplier from the error function. In other words k-е approximation minimizes f(xk) among vectors x0 + V, where vector V is the element of a subspace stretched by vectors Gx0, G2x0, , Gkx0.

We suggest

i 1 where {ui} - orthonormal basis of eigenvectors of the G matrix, therefore we get

n n x k  Pk G   i,0u i   i,0 Pk    u i , Pk 0   1,

i1 Hence, we have f  x k   1 2 Gx k , x k  1 2  i,0 Pk2  i  i .

2 (13) x 0   i,0 u i , n i1 n i1 x 0 2

n   2i,0 ,

i1

As a polynomial Pk() we choose the closest to optimal polynomial, the least deviating from zero on the interval [m, M], containing all the eigenvalues of the positive-definite G matrix and normalized so that Pk(0)  1.

By linear change of variables

the task is reduced to constructing a polynomial least deviating from zero on the interval t  [ – 1, 1 ] and taking at the point t0  (M + m) / (M – m), (corresponding   0) value 1. The solution of the last problem is given by the polynomial. n cos  k arccos t0  

Tk t  Tk  t0 

, m1at x1 Tk  t  

1

Tk  t0  m1at x1 Tk  t  . m1at x1 Tk  t   1,

1

Tk  t0  Lk  max Pk     mtax Tk  t  

,   m, M  , t   1,1 .  M  m  M  m 2 k  M  m   M  m   1      MM  mm   MM  mm 2  1     k    M  M 

k m  

   m  

M  M  m k 

  .

m     2i,0 Pk2  i   max Pk2  i  x 0 2 i .

(14) (15) (16) where Tk(t)  cos(karccos t) is the Chebyshev polynomial. At the same time

It is obvious that is why

Since the conception is fair

With sufficiently large k (k  k0) we have

 Lk  2  

M  M 

k m   m   2 

k  1     1   2 1 

k 2 

 ,    

M m . or

From (13) and (14) we get where q  1  2

  . Thus, the convergence of the SG method with the speed of a geometric progression is proved. The exact value of Lk, valid for any k, will be equal to

||xk||  Lk||x0|| ||xk||  2qk||x0||, k  k0, > 0, not exceeding the value of 0.23.

We will consider the problem of minimizing a quadratic functional f(x)  1/2Gx, x with a positively defined matrix G. Let us estimate the number of calculations f(x) required to achieve the control vector x with the norm ||x||  0,23 by the SG method and the RELCH algorithm from the starting point x0 c ||x0||  1. When reaching the x point the whole situation repeats, therefore, the comparative efficiency estimates obtained below are of a rather general character.

We will assume that two-sided finite difference relations are used to calculate the derivatives, which in the following analysis provides additional advantages to the SG method.

To achieve the vector x, the RELCH algorithm needs to calculate at the point x0 a weakly filled Hesse matrix and a gradient vector f(x0). With fill coefficient , it would require about 2n2 calculations of f. Further, we iterate L  1, 3  using formula (10), which do not require additional calculations of the objective functional f .

To obtain the vector x the SG method will need N iterations, where N is calculated this way:

||xN||  2qN  0,23, it means N  – 2,2/ ln q. To perform each iteration, it is necessary to update the gradient vector, that in the general case of the application of two-sided finite-difference approximations of derivatives is associated with 2n calculations of f(x). The total number of calculations f is – 4,4n / ln q. The relative benefit in the number of f calculations by the RELCH method compared to the SG method is given by the function ()  – 2,2/ (nln q). It is obvious that when    we have q()  1 and ()  . Typical values  for   0,01 and n  1000 are given below:   100 1000 1500 104 105   1,0 3,4 4,0 11,0 35,0

Thus, to obtain comparable results when   104 , the RELCH algorithm will require approximately 11 times less f calculations than with the SG method. However, it should be taken into account that as  increases, the number of L recalculations by the formula (10) grows. This can lead to an increase in the influence of computational errors in the calculation of s with large s numbers.

Example. We will consider the model task of the quadratic functional minimization f(x) с n = 200, h = 1500, g = 0,025. For definiteness, we assume that the time of a single calculation of f(x) is equivalent to performing 102n multiplication operations with floating point. The execution time of a single multiplication operation for some device for definiteness conditionally we will assume equal to ty  310– 5sec. The calculation of the f(x) value takes at the same time tf  0,6 sec of processor time. To calculate f и f using the general finite-difference formulas, we need, respectively, t  2ntf  4 min, t  2n2t  20 min. The number of recalculations by the formula (10) is equal to L  1, 3   50 . At each recalculation, a slightly filled matrix E – 2Gk is multiplied by a vector 50 , that requires n2ty  310– 2 sec of computer time. The vector construction time 50 without taking into account the calculation of f, f will be about 50310– 2  1,5 seconds and may not be taken into account.

The result is that to build the control vector x when ||x||  0,23 by using the RELCH method, it will take about t  20 min of machine time. The SG method will cost, respectively, (1500)20  1,3 hours. When the RELCH algorithm is re-applied to the constructed vector x, we will get the vector x с ||x||  0,23||x|| etc. Therefore, if we denote the corresponding sequence of vectors by{xm}, the norm of the vector x will decrease according to the law of geometric progression ||xm||  dm||x0||, where d <0.23 regardless of the  and n values.

An important additional advantage of the RELCH algorithm in comparison with the SG method is its rather high efficiency in the non-convex case, since the relaxation function of the method in the left half-plane is entirely located in the allowed area and the relaxation multipliers for   0 grow rapidly by absolute value when going from s to s1 . Growth characteristics were given before.

6 Conclusion

The described class of matrix gradient methods has shown in practice a sufficiently high performance in conditions of high stiffness and non-convexity of objective functionals. When optimizing large systems, it is possible to use efficient “packed” storage forms for matrices of second derivatives, that significantly reduces the requirements for the necessary computer memory. However, the method retains its main characteristics for small-sized systems, competing with the main optimization procedures of nonlinear programming.

Acknowledgements The work was financially supported by the Ministry of Education and Science of the Russian Federation in the framework of the Federal Targeted Program for Research and Development in Priority Areas of Advancement of the Russian Scientific and Technological Complex for 2014-2020 (No. 14.584.21.0022, ID RFMEFI58417X0022).

References

[1] I. Chernorutskiy , P. Drobintsev , V. Kotlyarov and N. Voinov . A New Approach to Generation and Analysis of Gradient Methods Based on Relaxation Function . Proceedings - 2017 UKSim-AMSS 19th International Conference on Modelling and Simulation , UKSim 2017 : 83 - 88 , May 2018 .

[2] I. G. Chernorutsky . Algorithmic problems of stiff optimization . St . Petersburg Polytechnic University Journal of Engineering Science and Technology, № 6 : 141 - 152 , 2012 .

[3] Yu. E. Nesterov and B. T. Polyak . Cubic regularization of Newton method and its global performance . Mathematical Programming , 108 ( 1 ): 177 - 205 , August 2006 .

[4] M. Avriel . Nonlinear programming: analysis and methods - Mineola , NY: Dover Publishing, 2003 .

[5] P. Deuflhard . Newton methods for nonlinear problems. Affine invariance and adaptive algorithms - Volume 35 / Springer Series in Computational Mathematics, Springer, 2004 .

[6] S. Pissanetski . Technology of sparse matrices. - Mir , 1988 .