Practical aspects of primal-dual nonlinear rescaling method with dynamic scaling parameter update

MATHEMATICAL ANALYSIS AND NUMERICAL MATHEMATICS VOLUME: 12 |NUMBER: 2 |2014 |JUNE Practical Aspects of Primal-Dual Nonlinear Rescaling Method with Dynamic Scaling Parameter Update Richard ANDRASIK Department of Mathematical Analysis and Applications of Mathematics, Faculty of Science, Palacky University Olomouc, st. 17. listopadu 12, 771 46 Olomouc, Czech Republic andrasik.ric[email protected] Abstract. Primal-dual nonlinear rescaling method with dynamic scaling parameter update (PDNRD) is an optimization method from a class of nonlinear rescaling techniques. Previous work does not discuss practical aspects of PDNRD method such as the explanation and the setting of the parameters. To complete this framework, the parameters were described. Moreover, PDNRD method was applied on two quadratic programming problems with quadratic constraints and recommendations about the setting of the parameters were made. Keywords Convex optimization, dynamic scaling parameter update, nonlinear rescaling method, parameters setting. 1. Introduction The optimization theory is developing in parallel with the appearance of real-life problems and with the need to solve them. Thanks to that, some types of problems gained a special status among the others and many tools to solve these problems have been developed. An example is the linear programming (LP). In 1947, George Dantzig introduced the simplex method for solving LP [3]. However, the simplex method has a strong competitor – interior-point methods, especially primal-dual interior-point methods [15]. In the field of the nonlinear programming (NLP), interior-point methods have been also applied. However, the situation in NLP is more complicated in comparison with LP calculations and interior point methods are sometimes experiencing numerical difficulties. This fact motivated Roman Polyak and Igor Griva to design an alternative method based on the nonlinear rescaling (NR) theory. Nowadays, NR methods can solve large-scale NLP problems with thousands of variables and constraints. They were successfully used to the radiotherapy treatment planning and are applied at some hospitals in USA and Europe [1]. The basic idea of NR methods is a nonlinear transformation of constraint functions to improve the properties of Lagrangian. Originally, the modified barrier methods [8] were introduced along with few modified barrier functions. Afterwards, the log-sigmoid function was also considered usable for NR [9], [10]. Consequently, the pieces of knowledge were refined, put together and a generalization of these techniques led to the concept of NR methods and NR functions. Similar to progress with interior-point methods, the primaldual nonlinear rescaling (PDNR) method was developed [11]. PDNR method is locally convergent with the Qlinear convergence rate. To improve these properties, PDNR method can be combined with another optimization method (e.g. the primal-dual path-following method) to obtain the global convergence [11]. Another way how to improve the convergence of PDNR method is a dynamic scaling parameter update [5] together with some globalization strategy (e.g. a step length computation). Recently, generalizations and other improvements [12], [13] were developed to improve the asymptotic convergence rate and to reduce the computational effort. The main purposes of this paper are to describe in detail parameters of PDNRD method and to give recommendations about their setting. The paper is organized as follows. First, the convex optimization problem is stated and basic assumptions c 2014 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 144 MATHEMATICAL ANALYSIS AND NUMERICAL MATHEMATICS VOLUME: 12 |NUMBER: 2 |2014 |JUNE are discussed. Then, NR functions are defined and the key idea of NR method is explained. Afterwards, basic primal-dual variant of NR method is presented. Next, PDNRD method is explained and its parameters are described. Finally, numerical experiments with different parameter settings were made and the results are presented in Section 7. 2. Statement of the Problem We consider the convex optimization problem (minimize f(x),x∈ <n, subject to ci(x)≥0,i= 1, . . . , r. (1) Function fis convex and functions ciare concave, ∀i= 1,. . .,r. Let S⊆ <nbe the admissible set of problem Eq. (1). For simplicity we define mapping c:<n→ <ras c(x)=(c1(x), c2(x), . . . , cr(x))T,∀x∈ <n. We suppose that: •Functions f,ci,∀i= 1, . . . , r, are at least twice continuously differentiable on the set <n. •The optimal set X∗=Argmin {f(x); x∈S}is bounded and not empty. •The Slater condition holds. For problem Eq. (1) we define the Lagrangian L(x;λ) = f(x)− r X i=1 λici(x).(2) Due to assumption (C), Karush-Kuhn-Tucker’s (KKT) conditions can be used to test the optimality. If ˆx∈X∗then there is a vector ˆ λ∈ <rsuch that ∇xL(ˆx;ˆ λ)=0, ˆ λ≥0, ˆ λici(ˆx)=0,∀i= 1, . . . , r. (3) Conversely, if a pair (ˆx, ˆ λ)∈S× <rsatisfies Eq. (3) then ˆx∈X∗. 3. Nonlinear Rescaling Functions First, we define functions that will be used to transform constraints of problem Eq. (1). Definition 1. Twice continuously differentiable function ψ: (t0; +∞)→ <, where −∞ < t0<0, satisfying conditions (i)ψ(0) = 0, ψ0(0) = 1, (ii)ψ0(t)>0,∀t∈(t0; +∞), (iii)ψ00(t)<0,∀t∈(t0; +∞), (iv)∃a > 0 : ψ(t)≤ −at2,∀t∈(t0; 0), (v)∃b > 0 : ψ0(t)≤bt−1,∀t > 0, (vi)∃c > 0 : ψ00(t)≥ −ct−2,∀t > 0 is called NR function. Remark 1. It follows from (ii) and (iii) that NR function ψis increasing and concave on the whole domain. Because of (ii) and (v) it is true that lim t→+∞ψ0(t) = 0. Similarly from (iii) and (vi) it holds that lim t→+∞ψ00(t) = 0. For example, the exponential transformation, the modified logarithmic function and the hyperbolic barrier function defined by the following formulas ψ1(t)=1−e−t, ψ2(t) = ln (t+ 1), ψ3(t) = t t+1 . are nonlinear rescaling functions [11]. Fig. 1: Graphs of functions ψ2(blue) and ψq2(red), τ=−1 2. Functions ψi,i= 2,3, can be modified so that ψi∈ C2(<). The function ψ1is already twice continuously differentiable and the following modification is not necessary, of course it can be done. For a given parameter τ∈(−1; 0) the quadratic extrapolation is defined by the relations ψqi(t) = ψi(t),for t≥τ, qi(t) = Ait2+Bit+Ci,for t < τ, whereas coefficients of the function qican be determined from the formulas ψi(τ) = qi(τ), ψ0 i(τ) = q0 i(τ), ψ00 i(τ) = q00 i(τ). c 2014 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 145 MATHEMATICAL ANALYSIS AND NUMERICAL MATHEMATICS VOLUME: 12 |NUMBER: 2 |2014 |JUNE From this equations we obtain Ai=1 2ψ00 i(τ), Bi=ψ0 i(τ)−τψ00 i(τ), Ci=ψi(τ)−τψ0 i(τ) + 1 2τ2ψ00 i(τ). 4. Nonlinear Rescaling Approach NR methods are based on the idea to transform problem Eq. (1) using a nonlinear rescaling function ψdefined on the whole real axis to the equivalent problem (minimize f(x), x ∈ <n, subject to k−1ψ(kci(x)) ≥0,i= 1, . . . , r. (4) From the definition of NR function, it is obvious that problems Eq. (1) and Eq. (4) have the same admissible sets, because ψ(kci(x)) ≥0⇔ci(x)≥0,∀i= 1, . . . , r. Hence, the optimal sets are also the same. Positive real number kis the scaling parameter, which can be fixed or dynamically enlarged during the iteration process, see Algorithm 3. The Lagrangian for equivalent problem Eq. (4) is given by the following formula: L(x;λ, k) = f(x)−k−1 r X i=1 λiψ(kci(x)).(5) For any r∈Nwe denote <r +={v∈ <r;vi≥0}, <r ++ ={v∈ <r;vi>0}. Suppose for a while that we know the solution of the dual problem λ∗∈ <r +. Then it is sufficient to minimize the function L(x;λ∗, k)in the primal variable x. The constrained optimization problem would be transformed to an unconstrained optimization problem. Since the Lagrange multipliers λ∗are not known, we estimate them and update them in every step of the method - just like the primal problem. In consequence, the constrained optimization problem is converted to a sequence of unconstrained optimization problems. Newton’s method or its variant is applied in each step to minimize the Lagrangian Lin the primal variable. Algorithm 1. (The basic concept of NR methods) Let k > 0be a scaling parameter. Initial approximations x0∈ <nand λ0∈ <r ++ are given. We suppose that an approximation (xs, λs)∈ <n× <r ++,s∈N0, is known already. We find the next primal-dual pair (xs+1, λs+1)using the following formulas xs+1 :∇xL(xs+1;λs, k) = 0, λs+1 i=ψ0kci(xs+1)λs i, i = 1, . . . , r. (6) The NR method converges for any fixed but large enough k > 0under the standard second order optimality conditions [5]. Remark 2. If λs∈ <r ++, then also λs+1 ∈ <r ++. This property follows from condition (ii)of NR function ψ. In other words, NR methods are interior-point methods in the dual variable. The Algorithm 1 is well defined due to the following theorem. Theorem 1. Suppose that X∗is bounded. For any given (λ, k)∈ <r ++ × <++ there exists one and only one ˆx∈ <nsuch that L(ˆx;λ, k) = min x∈<nL(x;λ, k). Proof: [10], page 206. The main purpose of NR is to improve properties of the Lagrangian. The classical Lagrangian L– as a connection between the constrained and the unconstrained optimization – does not always work, because the existence of the unconstrained Lagrange minimizer is unknown in general. On the other hand, the unconstrained minimizer of the Lagrangian Lalways exists (according to Theorem 1). Moreover, NR dramatically sharpens the reaction of the Lagrangian to the constraint violation, which has an impact on the computations. In comparison to interior-point methods, there is no "infinite wall". NR methods are exterior point methods. Also, the unbounded increase of the scaling parameter is not needed to guarantee the convergence of the method. We illustrate the influence of NR on the Lagrangian on the following simple example. Example 1. The Lagrangian of the following problem (minimize x, x ∈ <1, subject to x+ 1 ≥0,(7) and the Lagrangian of the equivalent problem are considered. (minimize x, x ∈ <1, subject to k−1ψq2(k(x+ 1)) ≥0,(8) c 2014 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 146 MATHEMATICAL ANALYSIS AND NUMERICAL MATHEMATICS VOLUME: 12 |NUMBER: 2 |2014 |JUNE We take a look on the graphs of both Lagrangians, which are defined by the formulas L(x;λ) = x−λ(x+ 1), L(x;λ, k) = x−k−1λψq2(k(x+ 1)) .(9) NR fundamentally changes the shape of the Lagrangian. The sharp reaction of the Lagrangian Lto the constraint violation is obvious in Fig. 2. Fig. 2: The Lagrangian for the problem Eq. (7) - on the top - and for the equivalent problem Eq. (8) with k= 1 - on the bottom. We project the graphs of the both Lagrangians to the phase plane in Fig. 3. Suppose that we want to find Lagrange minimizer for λ= 2. After substitution, formulas Eq. (9) have the form: L(x; 2) = −x−2, L(x; 1, k) = x−k−12ψq2(k(x+ 1)) . The function L(x; 2) is a decreasing linear function. We cannot find the Lagrange minimizer of it, because it is unbounded below. If we rearrange the first equation in formulas Eq. (9) in the following way: L(x;λ) = (1 −λ)x−λ, it is apparent that for any λ6= 1 we cannot find the Lagrange minimizer. On the other hand, we know from Theorem 1 that the Lagrangian Lhas a minimizer for any λ > 0. This is the main contribution of NR. Fig. 3: The projection of the Lagrangian for the problem Eq. (7) - on the top - and for the equivalent problem Eq. (8) with k= 1 - on the bottom - to the phase plane. 5. Primal-dual Nonlinear Rescaling Method The update of the Lagrange multipliers is clear from equations Eq. (6). When we want to calculate xs+1, we must solve a nonlinear system of equations using Newton’s method. Although this approach works well far from the solution, primal-dual variant of the NR method performs better in the neighborhood of the solution [5]. In Section 6, we will combine these two techniques to obtain a stable and globally convergent method. Let approximations x∈ <n,λ∈ <r ++ be known. We suppose that ˆx=x+ ∆x, ˆ λ=¯ λ+ ∆λ, (10) c 2014 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 147 MATHEMATICAL ANALYSIS AND NUMERICAL MATHEMATICS VOLUME: 12 |NUMBER: 2 |2014 |JUNE where (∆x, ∆λ)is a primal-dual Newton step and ¯ λis a predictor of the Lagrange multipliers, which is given by the following formula ¯ λ= Ψ0(kc(x)) λ, where Ψ0(kc(x)) = diag (kci(x))r i=1 . The reason, why define ¯ λin this way is hidden in the equality ∇xL(x;λ, k) = ∇xL(x;¯ λ). One step of Newton’s method consists of solving the system of (n+r)linear equations N(·)∆x ∆λ=−∇xL(x;¯ λ) 0,(11) where a pair (∆x, ∆λ)is unknown and N(·) = ∇2 xxL(x;¯ λ)−∇c(x)T −kΨ00 (kc(x)) Λ∇c(x)Ir. Symbols Ψ00 (kc(x)),Λdenotes diagonal matrices defined as Ψ00 (kc(x)) = diag (ψ00(kci(x)))r i=1 , Λ = diag (λi)r i=1 . If N(·)is sparse, we can use numerical algebra techniques for sparse matrices [15] to solve the system Eq. (11). In the opposite case, we express ∆λfrom the second equation of the system Eq. (11) and we substitute it to the first equation. Now, we must solve only a system with nequations, instead of a system of (n+r)equations. The system has the following form M(·)∆x=−∇xL(x;¯ λ),(12) where M(·) = ∇2 xxL(x;¯ λ)−k∇c(x)TΨ00 (kc(x)) Λ∇c(x). Remark 3. It is obvious that M(·)is a symmetric matrix. It can be shown [11], page 120 that if the second order optimality conditions hold true at the point (x∗;λ∗), then M(·)is a positive definite matrix for all (x;λ)sufficiently close to the optimal primal-dual pair, when the scaling parameter is sufficiently large. Algorithm 2. (One step of PDNR method) Let k∈ <++ be fixed. The approximations x∈ <n, λ∈ <r ++ are given. •We compute a dual predictor ¯ λ= Ψ0(kc(x)) λ. •We calculate the Newton step (∆x, ∆λ)from system Eq. (11), or •we find out the primal Newton step ∆xfrom formula Eq. (12) and then we compute the dual Newton step using the following relation ∆λ=kΨ00 (kc(x)) Λ∇c(x)∆x. •We calculate the new primal-dual approximation x:= x+ ∆x, λ := ¯ λ+ ∆λ. The Algorithm 2 describes only one step of PDNR method. Although, the input vector of multipliers belongs to <r ++, it is not guaranteed that the updated vector of multipliers has also all its components positive. In contrast, this statement is true for Algorithm 1, according to Remark 2. Therefore, PDNR method cannot stand alone in a general case and the improvement is needed (see Algorithm 3). 6. Dynamic Scaling Parameter Update To obtain a higher convergence rate of the method, we dynamically change the scaling parameter. Moreover, we use Newton’s method with a step length (e.g. the backtracking line search algorithm) to solve formulas Eq. (6). We introduce a function which measures the distance between the approximation (x, λ)and the solution (x∗, λ∗). Definition 2. The function ν:<n×<r→ <+, defined as follows ν(x, λ) = max { k∇xL(x;λ)k, −min1≤i≤rci(x), −min1≤i≤rλi, Pr i=1 |λici(x)| }, (13) is called the merit function. Formula (13) is motivated by KKT conditions. From the first order optimality conditions it follows that ν(ˆx, ˆ λ) = 0 ⇔(ˆx, ˆ λ)∈X∗. For known primal-dual pair (x, λ)we set the parameter kaccording to the relation k=ν(x, λ)−1/2.(14) It is obvious that if a primal-dual sequence {(xs, λs)}+∞ 0tends to (x∗, λ∗), then ν(xs, λs)→0+ and also ks→+∞when s→+∞. c 2014 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 148 MATHEMATICAL ANALYSIS AND NUMERICAL MATHEMATICS VOLUME: 12 |NUMBER: 2 |2014 |JUNE To prevent from a singularity in system Eq. (11), we solve the following regularized system Nk(·)∆x ∆λ=−∇xL(x;¯ λ) 0,(15) where the primal-dual pair (x, λ)is an approximation of solution, kis the scaling parameter, ¯ λ= Ψ0(kc(x)) λ is the dual predictor of the Lagrange multipliers and Nk(·) = ∇2 xxL(x;¯ λ) + 1 k2In−∇c(x)T −kΨ00 (kc(x)) Λ∇c(x)Ir. It is important to remark that the regularization has no effect to the 1.5-superlinear convergence rate of PDNRD method [5]. In a comparison to interior-point methods, the inverse matrix to Nk(·)exists at the optimal point for every k∈ <++. Assume that the approximation (x, λ)of the point (x∗, λ∗)is known. First, we use PDNRD method. If there is a superlinear decrease of the merit function, we have found the next approximation. In the opposite case, we use the primal Newton step ∆x(calculated during the use of PDNRD method) to minimize the function L(x;λ, k), where λand kare fixed. We apply the backtracking line search method to guarantee the global convergence in the minimization process. In this way, we obtained the globally convergent PDNRD method with the 1.5-superlinear convergence rate. Algorithm 3. (The globally convergent PDNRD method) An initial approximation x0∈ <nis given. An accuracy parameter ε > 0and an initial scaling parameter k∈ <++ are given. Parameters q∈(0; 1),η∈(0; 0.5), ω > 1,σ > 0and θ > 0are also given. Set x:= x0, λ:= (1,1,...,1) ∈ <r,λg:= λand H:= ν(x, λ). •If H≤ε, then stop, output (x, λ). •Find ¯ λand (∆x, ∆λ)from primal-dual system Eq. (15) with known (x, λ)and set ˆx:= x+ ∆x, ˆ λ:= ¯ λ+ ∆λ, ˆ H:= ν(ˆx, ˆ λ). •If ˆ H≤min H3/2−θ,1−θ, then set x:= ˆx, λ := ˆ λ, H := ˆ H, k := max 1 √H, k and go to step 1. •Find α∈(0; 1iso that it holds L(x+α∆x;λg, k)−L(x;λg, k) ≤ηα∆xT∇xL(x;λg, k), using the backtracking line search algorithm. •Set x:= x+α∆x, ˆ λ:= Ψ0(kc(x))λg. •If k∇xL(x;λg, k)k ≤ σ kkˆ λ−λgk, then go to step 8. •Find (∆x, ∆λ)from primal-dual system Eq. (15) with known (x, λg)and go to step 4. •If ν(x, ˆ λ)≤qH, then set λ:= ˆ λ, λg:= ˆ λ, H := ν(x, ˆ λ), k := max 1 √H, k and go to step 1. •Set k:= ωk and go to step 7. 7. Numerical Experiments From the finite element approximation of contact problems of the linear elasticity with the friction in three space dimensions arise a minimization problem      minimize 1 2xTAx−xTb,x∈ <n subject to g2 i−x2 i+m−x2 i+2m≥0,i∈ I, xi−li≥0,i∈ I, (16) where n= 3mis the number of variables, A∈ <n×n is a symmetric and positive definite matrix, b∈ <n, g∈ <m +,l∈ <mand I={1,2, . . . , m}. This is a convex programming problem so we can use NR approach to solve it. We use the function ψq2with τ=1 2to rescale the conditions and subsequently we obtain the equivalent problem. The reason for this choice is that for a sufficiently large class of functions c(x)is the function −ln c(x)self-concordant. If we apply Newton’s method to a self-concordant function, we can say something more about the convergence of Newton’s method [2]). PDNRD method was tested on two model problems - the chord problem, which in contrary to Eq. (16) contains also unconstrained variables, and the steel brick problem, whose finite element approximation leads exactly to Eq. (16). All computations were performed in MATLAB on the PC Intel Core i7 (2.4GHz) with 8GB RAM. In tables below we reported the number of iterations (it) , the number of solutions of the primal-dual system (nS) and the solution time in seconds (time). If some data is missing, it means that the solution time was too long in the comparison to the other cases in the table. c 2014 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 149 MATHEMATICAL ANALYSIS AND NUMERICAL MATHEMATICS VOLUME: 12 |NUMBER: 2 |2014 |JUNE 7.1. Chord Problem We consider a problem min u∈K J(u),(17) where J(u) = 1 2Z1 0ku0(t)k2dt−Z1 0 u(t)Tf(t)dt, K={u∈H1 0(0; 1)2:u2(t)≥0,∀t∈(0; 0.5), ku(t)k ≤ 1.4,∀t∈(0.5; 1) }, f(t) = 36π2sin 6πt, −4π2sin 2πtT. Fig. 4: The chord deformation. Minimization problem Eq. (17) describes a loaded chord fixed at the endpoints. The chord is partially above a plain and partially inside a cylindrical tube Fig. 4. The function u(t)is the chord deflection. The chord problem was presented as a model problem in [6]. The objective function 1 2xTAx−xTb,x∈ <n matches the convex quadratic functional J(u), linear constraints in problem Eq. (16) matches the constraint u2(t)≥0,∀t∈(0; 0.5) from the definition of the set Kand quadratic constraints matches ku(t)k ≤ 1,4, ∀t∈(0.5; 1). First, we decided if the primal-dual system has a sparse matrix. We say that the matrix is sparse if it has at most 10 % nonzero elements. In Fig. 5, the structure of the matrix Nk(·)∈ <192×192 for n= 128 is depicted. This matrix has only 636 nonzero elements (1.7 %). Hence, Nk(·)is the sparse matrix. In the same way we could argue for other choices of n. We solved the chord problem for different settings of parameters of PDNRD method. The main result from Tab. 1 is a non-increasing number of iterations and also a non-increasing number of solutions of the primal-dual system while increasing the number of variables. The best choice of the scaling parameter was kinit = 10 in Fig. 5: The structure of the matrix Nk(·)∈ <192×192 (the chord problem). Black colour indicates 636 nonzero elements. Tab. 1: The chord problem. PDNRD method with parameters ω= 10,σ=1 2kinit,θ= 0.4,q= 0.5,η= 0.01, ε= 10−6. kinit n r 10 102103 64 32 10/31/0.173 4/23/0.118 7/7/0.021 128 64 7/50/0.238 6/33/0.126 3/79/0.351 256 128 8/60/0.515 8/37/0.454 3/143/1.099 512 256 11/65/1.989 5/49/2.037 7/150/5.036 1024 512 12/70/10.330 6/72/13.088 8/100/22.229 2048 1024 12/58/37.640 6/70/53.256 7/66/50.650 4096 2048 13/74/221.834 – – it/nS/time it/nS/time it/nS/time Tab. 2: The chord problem (n= 1024). PDNRD with parameters ω= 10,σ=1 2kinit,kinit = 10,ε= 10−6. η= 0.01 q θ = 0.1θ= 0.4 0.1 6/49/8.015 6/49/8.054 0.5 10/65/11.554 12/70/10.518 0.9 30/99/16.347 31/77/13.548 η= 0.3 q θ = 0.1θ= 0.4 0.1 6/49/8.133 6/49/8.036 0.5 10/68/11.503 12/71/10.671 0.9 30/92/15.094 31/78/12.885 it/nS/time it/nS/time this case, because of the best performance for larger problems. From Tab. 2 it is obvious that the choice q= 0.1was the best one. Parameters ηand θhad no significant impact on the computation in this case. 7.2. Steel Brick Problem Let us consider a steel brick lying on the rigid obstacle. The brick occupies the domain S= (0; 3) ×(0; 1) × c 2014 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 150 MATHEMATICAL ANALYSIS AND NUMERICAL MATHEMATICS VOLUME: 12 |NUMBER: 2 |2014 |JUNE (0; 1). The boundary ∂Sis divided into three parts Γu={0}×(0; 1) ×(0; 1), Γc= (0; 3) ×(0; 1) ×{0}, Γp=∂S\¯ Γu∪¯ Γc, on which different boundary conditions are given (see Fig. 6). The problem is described in the detail in [6]. Fig. 6: The steel brick. The progress was similar to the previous one considering the chord problem. First, the matrix Nk(·) was tested whether it is sparse or not. In Fig. 7, the structure of the matrix Nk(·)∈ <150×150 for n= 90 is depicted. This matrix contains 8340 nonzero elements (37 %). Hence, the matrix Nk(·)is not sparse and we solved the primal-dual system in a different way (see Section 5). Fig. 7: The structure of matrix Nk(·)∈ <150×150 (steel brick problem). Black colour indicates 8340 nonzero elements. The number of variables was denoted as n= 3 m and the number of conditions as r= 2 m. The steel brick problem was solved for different choices of n. According to the data in Tab. 3, the best setting for the initial values of the scaling parameter was kinit = 10 in this problem. Number nSrose with increasing kinit. The results from Tab. 4 supports that the choice η= 0.01 is better than η= 0.3. The parameter θhas no significant impact on the computation. Tab. 3: The steel brick problem. PDNRD method with parameters ω= 10,σ=1 2kinit,θ= 0.4,q= 0.5,η= 0.01, ε= 10−6. kinit n r 10 102103 54 36 6/39/0.153 5/54/0.268 4/119/0.570 90 60 6/38/0.218 6/45/0.239 4/102/0.743 180 120 8/40/0.257 4/90/0.700 4/136/1.653 324 216 7/37/0.663 10/35/0.619 4/130/3.146 648 432 12/67/7.258 10/77/10.003 6/123/15.696 it/nS/time it/nS/time it/nS/time Tab. 4: The steel brick problem (n= 648). PDNRD method with parameters ω= 10,σ=1 2kinit,kinit = 10,ε= 10−6. η= 0.01 q θ = 0.1θ= 0.4 0.1 5/46/5.562 5/46/5.767 0.5 12/45/4.460 12/67/7.135 0.9 12/67/7.545 12/67/7.105 η= 0.3 q θ = 0.1θ= 0.4 0.1 5/44/5.593 5/44/5.482 0.5 11/41/4.398 12/62/8.820 0.9 12/62/8.761 12/62/8.718 it/nS/time it/nS/time The characteristic property of PDNRD method is the "hot" start of the method (as it was observed and proved in [11]). "Hot" start means that from some point only one Newton step is needed to sufficiently shrink the distance between a current approximation and the solution. We can expect that a computational effort to solve the problem with the accuracy 10−12 will not be dramatically higher than with an accuracy 10−6. Our results support this statement (compare Tab. 3 and Tab. 5). Tab. 5: The steel brick problem. PDNRD method with parameters ω= 10,σ=1 2kinit,kinit = 10,θ= 0.4,q= 0.5, η= 0.01,ε= 10−12. n r kinit = 10 54 36 8/45/0.165 90 60 10/51/0.252 180 120 12/44/0.269 324 216 10/47/0.806 648 432 14/71/7.371 it/nS/time 7.3. Parameters When using a numerical method the setting of the parameters is very important. Experience with numerical experiments helps us to find the optimal setting. Based c 2014 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 151 MATHEMATICAL ANALYSIS AND NUMERICAL MATHEMATICS VOLUME: 12 |NUMBER: 2 |2014 |JUNE on the tests of PDNRD method on examples from Sections 7.1 and 7.2, we made considerations about a suitable setting of the parameters. Factor ω The factor ωaffects the rate of the increase of the scaling parameter. For ω∈ h5; 20iwe obtain almost the same results. So we set ω= 10. However, even for the other choices ω > 1there are not any significant changes. At most, it may happen that it takes a few extra steps of the method. Parameter σ The choice of the parameter σis less important than the choice of the ratio between σand the initial choice of the scaling parameter kinit. It is the fraction σ k which decides about the number of inner iterations (and thus about the number of Newton steps in the damped phase of Newton’s method). It is clear that for σ >> k there are too little inner steps. On the other hand, for σ << k there is too many of them. According to the results of the numerical experiments, the choice σ=1 2kinit is suitable. Parameter θ The parameter θaffects whether the inner solver runs or not. If we set θ= 0, the inner solver runs whenever the 1.5-superlinear decrease of the merit function was not achieved. Choosing θ= 0.5, we are saying that we are satisfied with only linear decrease of the merit function. Therefore, it is wise to set θ∈ h0; 0.5i. It appears (see Tab. 2 and Tab. 4) that the method does not depend on this parameter. In the condition, which is related to the parameter θ, the term min H3/2−θ,1−θis calculated. We can ask why do not simply use the term H3/2−θinstead of the previous one. As we know, the classical Newton method is effective only in the neighborhood of the solution, so if we are "far" from the solution the damped Newton method is better suited to use. The expression "far" means that the merit function is greater than 1−θ. Factor q The number of inner steps is influenced by the fraction σ ktogether with the factor q. This factor also affects how often is the scaling parameter increased. Due to the way in which is this parameter used in PDNRD method, it is needed to set q∈(0; 1). Moreover, we must choose the factor qso that the method does not use too many inner steps. Based on the data in Tab. 2 and Tab. 4, any factor q∈(0; 1) was suitable. However, the best results were obtained by setting q= 0.1. Parameter η The backtracking line search parameter ηis usually chosen in the range from 0.01 to 0.3(see [2] page 466). According to data in Tab. 2 and Tab. 4 the choice of the parameter η∈ h0.01; 0.3iis arbitrary. Initial value of the scaling parameter Based on the results shown in Tab. 1 and Tab. 3, increasing the initial value of the scaling led to instabilities in computations. The best results in the both model examples were obtained with kinit = 10. Therefore, we recommend setting kinit around ten or less. Further reason to set kinit rather lower is that the scaling parameter can be increased during the computation, but it cannot be decreased. Initial approximation The computation also depends on the initial approximation. The number of steps can be decreased, if "lucky" initial approximation is chosen. In problems like a beam deflection or contact problems we choose an initial state of the system as the initial approximation, because the shape changes of a body are usually very small. Thus we set x0= (0,0,...,0). 8. Conclusion The practical aspects of PDNRD method were described. Especially, the meaning of the parameters of PDNRD method were explained in detail. PDNRD method was tested on two quadratic programming problems with quadratic constraints (the chord problem and the steel brick problem). Based on the tests, the recommendations about setting the parameters were made. The solution of the primal-dual system is the most expensive operation, thus the number of solutions of the primal-dual system determines the total complexity of computations. It was found out that the increasing number of variables in both presented problems has not a consequence in the increasing number of solutions of the primal-dual system (see Tab. 1 and Tab. 3). This fact supports the applicability of PDNRD method on problems of an arbitrary size. Acknowledgment This research was supported by the grant PrF 2013, 013 Mathematical models of the Internal Grant Agency of Palacky University Olomouc. c 2014 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 152