*********************************************************************** * * * PNUL - A DISCRETE NEWTON ALGORITHM FOR LARGE-SCALE * * EQUALITY CONSTRAINED OPTIMIZATION BASED ON THE * * NULL SPACE APPROACH. * * * *********************************************************************** 1. Introduction: ---------------- The double-precision FORTRAN 77 basic subroutine PNUL is designed to find a close approximation to a local minimum of a twice continuously differentiable function F(X) on the set defined by the nonlinear equality constraints FC_1(X) = 0, FC_2(X) = 0, ..., FC_NC(X)=0. Here X is a vector of NF variables and FC_I(X), 1 <= I <= NC, are twice continuously differentiable functions. We assume that the sparsity pattern of the objective Hessian matrix is known and that NF and NC are large, but constraint functions FC_I(X), 1 <= I <= NC depend on a small number of variables. The latter assumption implies that the mapping CF(X) = [FC_1(X), FC_2(X), ..., FC_NC(X)] has a sparse Jacobian matrix, which will be denoted by CG(X) (it has NC rows and NF columns). The sparsity pattern of the objective Hessian matrix (only the upper part) is stored in the standard compressed row format using arrays IH and JH. For example, if the Hessian matrix has the following pattern H = | * * * 0 * | | * * 0 * 0 | | * 0 * 0 * | | 0 * 0 * 0 | | * 0 * 0 * | (asterisks denote nonzero elements) then arrays IH and JH contain elements IH(1)=1, IH(2)=5, IH(3)=7, IH(4)=9, IH(5)=10, IH(6)=11, JH(1)=1, JH(2)=2, JH(3)=3, JH(4)=5, JH(5)=2, JH(6)=4, JH(7)=3, JH(8)=5, JH(9)=4, JH(10)=5, i.e., IH contains pointers of the diagonal elements in the upper part of the Hessian matrix and JH contains column indices of the nonzero elements. Note that IH has NF+1 elements, the last element is equal MH+1, where MH is the number of nonzero elements stored. The sparsity pattern of the constraint Jacobian matrix is stored in the standard compressed row format using arrays ICG and JCG. For example, if the Jacobian matrix has the following pattern CG = | * * 0 * | | * * * 0 | | * 0 0 * | | 0 * * 0 | | * 0 * 0 | (asterisks denote nonzero elements) then arrays ICG and JCG contain elements ICG(1)=1, ICG(2)=4, ICG(3)=7, ICG(4)=9, ICG(5)=11, ICG(6)=13, JCG(1)=1, JCG(2)=2, JCG(3)=4, JCG(4)=1, JCG(5)=2, JCG(6)=3, JCG(7)=1, JCG(8)=4, JCG(9)=2, JCG(10)=3, JCG(11)=1, JCG(12)=3, i.e., ICG contains pointers of the first elements in rows of the Jacobian matrix and JCG contains column indices of the nonzero elements. Note that ICG has NC+1 elements, the last element is equal MC+1, where MC is the number of nonzero elements. To simplify user's work, an additional easy to use subroutine PNULU is added. It calls the basic general subroutine PNUL. All subroutines contain a description of formal parameters and extensive comments. Furthermore, test program TNULU is included, which contains several test problems (see e.g. [2]). This test program serve as an example for using the subroutine PNULU, verify its correctness and demonstrate its efficiency. In this short guide, we describe all subroutines which can be called from the user's program. A detailed description of the method is given in [1]. In the description of formal parameters, we introduce a type of the argument that specifies whether the argument must have a value defined on entry to the subroutine (I), whether it is a value which will be returned (O), or both (U), or whether it is an auxiliary value (A). Note that the arguments of the type I can be changed on output under some circumstances, especially if improper input values were given. Besides formal parameters, we can use a COMMON /STAT/ block containing statistical information. This block, used in each subroutine has the following form: COMMON /STAT/ NRES,NDEC,NIIT,NIT,NFV,NFG,NFH The arguments have the following meaning: Argument Type Significance ---------------------------------------------------------------------- NRES O Positive INTEGER variable that indicates the number of restarts. NDEC O Positive INTEGER variable that indicates the number of matrix decompositions. NIIT O Positive INTEGER variable that indicates the number of inner iterations (for solving linear systems). NIT O Positive INTEGER variable that indicates the number of iterations. NFV O Positive INTEGER variable that indicates the number of function evaluations. NFG O Positive INTEGER variable that specifies the number of gradient evaluations. NFH O Positive INTEGER variable that specifies the number of Hessian evaluations. 2. Subroutine PNULU: -------------------- The calling sequence is CALL PNULU(NF,NC,X,IH,JH,CF,ICG,JCG,IPAR,RPAR,F,GMAX,CMAX, & IPRNT,ITERM) The arguments have the following meaning. Argument Type Significance ---------------------------------------------------------------------- NF I Positive INTEGER variable that specifies the number of variables of the partially separable function. NC I Positive INTEGER variable that specifies the number of constraints. X(NF) U On input, DOUBLE PRECISION vector with the initial estimate to the solution. On output, the approximation to the minimum. IH(NF+1) I INTEGER array which contains pointers of the diagonal elements in the upper part of the Lagrangian function Hessian matrix. JH(NH) I INTEGER array which contains column indices of the nonzero elements and additional working space for the Choleski factor. CF(NC+1) O DOUBLE PRECISION vector which contains values of constraint functions. ICG(NC+1) I INTEGER array which contains pointers of the first elements in rows of the constraint Jacobian matrix. JCG(MC) I INTEGER array which contains column indices of the nonzero elements. IPAR(7) I INTEGER parameters: IPAR(1)=MIT, IPAR(2)=MFV, IPAR(3)=MFG, IPAR(4)=MOS, IPAR(5)=MD, IPAR(6)=MDE, IPAR(7)=IFIL. Parameters MIT, MFV, MFG, MOS, MD, MDE are described in Section 3 together with other parameters of the subroutine PIND. Parameter IFIL specifies a relative size of the space reserved for fill-in. The choice IFIL=0 causes that the default value IFIL=1 will be taken. RPAR(5) I DOUBLE PRECISION parameters: RPAR(1)=XMAX, RPAR(2)=TOLX, RPAR(3)=TOLC, RPAR(4)=TOLG, RPAR(5)=RPF1. Parameters XMAX, TOLX, TOLC, TOLG, RPF1 are described in Section 3 together with other parameters of the subroutine PIND. F O DOUBLE PRECISION value of the objective function at the solution X. GMAX O DOUBLE PRECISION maximum absolute value of a partial derivative of the Lagrangian function. CMAX O maximum constraint violation. IPRNT I INTEGER variable that specifies PRINT: IPRNT= 0 - print is suppressed, IPRNT= 1 - basic print of final results, IPRNT=-1 - extended print of final results, IPRNT= 2 - basic print of intermediate and final results, IPRNT=-2 - extended print of intermediate and final results. ITERM O INTEGER variable that indicates the cause of termination: ITERM= 1 - if |X - XO| was less than or equal to TOLX in two subsequent iterations, ITERM= 2 - if |F - FO| was less than or equal to TOLF in two subsequent iterations, ITERM= 3 - if F is less than or equal to TOLB, ITERM= 4 - if GMAX is less than or equal to TOLG, ITERM=11 - if NIT exceeded MIT, ITERM=12 - if NFV exceeded MFV, ITERM=13 - if NFG exceeded MFG, ITERM< 0 - if the method failed. If ITERM=-6, then the termination criterion has not been satisfied, but the point obtained is usually acceptable. Values ITERM<=-40 detect a lack of space. In this case, parameter IPAR(7)=IFIL has to be increased (IFIL=2, IFIL=3, etc.). The subroutine PNULU requires the user supplied subroutines OBJ, DOBJ that define the objective function and its gradient and CON, DCON that define constraint functions and their gradients. These subroutines have the form SUBROUTINE OBJ(NF,X,F) SUBROUTINE DOBJ(NF,X,G) SUBROUTINE CON(NF,KC,X,FC) SUBROUTINE DCON(NF,KC,X,GC) The arguments of the user supplied subroutines have the following meaning. Argument Type Significance ---------------------------------------------------------------------- NF I Positive INTEGER variable that specifies the number of variables of the objective function. X(NF) I DOUBLE PRECISION an estimate to the solution. F O DOUBLE PRECISION value of the objective function at the point X. G(NF) O DOUBLE PRECISION gradient of the objective function at the point X. KC I INTEGER index of the partial function. FC O DOUBLE PRECISION value of the KC-th partial function at the point X. GC(NF) O DOUBLE PRECISION gradient of the KC-th partial function at the point X. 3. Subroutine PNUL: ------------------- This general subroutine is called from all subroutines described in Section 2. The calling sequence is CALL PNUL(NF,NC,MAXH,MAXB,X,GF,H,IH,JH,CG,ICG,JCG,KCG,CF,CFO, & GC,G,S,XO,GO,XS,CZ,CZO,CP,B,IB,JB,CGD,DD,BD,COL,PSR,PERC,INVP, & WN11,WN12,WN13,WN14,XMAX,TOLX,TOLC,TOLG,RFF1,F,GMAX,CMAX,MIT, & MFV,MFG,MOS,MD,MDE,IPRNT,ITERM) The arguments NF, NC, X, IH, JH, CF, ICG, JCG, F, GMAX, CMAX, IPRNT, ITERM have the same meaning as in Section 2. Other arguments have the following meaning: Argument Type Significance ---------------------------------------------------------------------- MAXH I INTEGER size of array H. MAXB I INTEGER size of array B. GF(NF) A DOUBLE PRECISION gradient of the objective function. H(MAXH) A DOUBLE PRECISION nonzero elements of the approximation of the Hessian matrix and nonzero elements of the Choleski factor. CG(MC) A DOUBLE PRECISION nonzero elements of the Jacobian matrix. KCG A INTEGER auxiliary array. CFO(NC+1) O DOUBLE PRECISION vector which contains old values of constraint functions. GC(NF) A DOUBLE PRECISION gradient of the constraint function. G(NF) A DOUBLE PRECISION gradient of the merit function. S(NF+NC) A DOUBLE PRECISION direction vector. XO(NF) A DOUBLE PRECISION array which contains increments of variables. GO(NF+1) A DOUBLE PRECISION array which contains increments of gradients. XS(NF) A DOUBLE PRECISION auxiliary array. CZ(NC) A DOUBLE PRECISION vector of Lagrange multipliers. CZO(NC) A DOUBLE PRECISION old vector of Lagrange multipliers. CP(NF) A DOUBLE PRECISION auxiliary array. B(MAXB) A nonzero elements of the preconditioner together with the nonzero elements of its Choleski factor. IB(NC+1) A INTEGER array which contains pointers of the diagonal elements in the upper part of the preconditioner. JB(MAXB) A INTEGER array which contains column indices of the nonzero elements. GCD(NF*ND) A DOUBLE PRECISION array which contains dense columns. DD(ND) A DOUBLE PRECISION auxiliary array. BD(ND*ND) A DOUBLE PRECISION auxiliary array. COL(NF) A INTEGER auxiliary array. PSR(NC+1) A INTEGER pointer vector in the compact form of the Choleski factor. PERC(NC) A INTEGER permutation vector INVP(NC) A INTEGER inverse permutation vector WN11(NC+1) A INTEGER auxiliary array. WN12(NC+1) A INTEGER auxiliary array. WN13(NC+1) A INTEGER auxiliary array. WN14(NC+1) A INTEGER auxiliary array. XMAX I DOUBLE PRECISION maximum stepsize; the choice XMAX=0 causes that the default value 1.0D+16 will be taken. TOLX I DOUBLE PRECISION tolerance for the change of the coordinate vector X; the choice TOLX=0 causes that the default value TOLX=1.0D-16 will be taken. TOLC I DOUBLE PRECISION tolerance for the constraint violation; the choice TOLC=0 causes that the default value TOLC=1.0D-6 will be taken. TOLG I DOUBLE PRECISION tolerance for the Lagrangian function gradient; the choice TOLG=0 causes that the default value TOLG=1.0D-6 will be taken. RPF1 I DOUBLE PRECISION value of the penalty parameter in the merit function; the choice RPF1=0 causes that the default value 1.0D-4 will be taken. MIT I INTEGER variable that specifies the maximum number of iterations; the choice MIT=0 causes that the default value 1000 will be taken. MFV I INTEGER variable that specifies the maximum number of function evaluations; the choice |MFV|=0 causes that the default value 1000 will be taken. MFG I INTEGER variable that specifies the maximum number of gradient evaluations; the choice |MFV|=0 causes that the default value 10000 will be taken. MOS I INTEGER choice of preconditioning strategy: MOS= 1 - preconditioning by the constraint preconditioner with complete Gill-Murray decomposition. MOS=-1 - preconditioning by the constraint preconditioner with incomplete Gill-Murray decomposition. The choice MOS=0 causes that the default value MOS=1 will be taken. MD I INTEGER maximum number of the dense rows; the choice MD=0 causes that the default value 10 will be taken. MDE I INTEGER maximum number of nonzero elements in a sparse rows; the choice MDE=0 causes that the default value 50 will be taken. The choice of parameter XMAX can be sensitive in many cases. First, the objective function can be evaluated only in a relatively small region (if it contains exponentials) so that the maximum stepsize is necessary. Secondly, the problem can be very ill-conditioned far from the solution point so that large steps can be unsuitable. Finally, if the problem has more local solutions, a suitably chosen maximum stepsize can lead to obtaining a better local solution. A suitable choice of parameter RPF1 can improve the efficiency of the method. The default value RPF1=0.0001 is relatively small. Larger value (RPF1=1 say) can sometimes give better result. The subroutine PNUL requires the user supplied subroutines OBJ DOBJ, CON, DCON, which are described in Section 2. 4. Verification of the subroutines: ----------------------------------- Subroutine PNULU can be verified and tested using the program TNULU. This program calls the subroutines TIUS20, TINS20 (initiation), TFFU20 (objective function evaluation), TFGU20 (objective gradient evaluation), TCFU20 (constraint functions evaluation) and TFGU20 (constraint gradients evaluation) containing 18 equality constrained test problems with at most 1000 variables [2]. The results obtained by the program TNULU on a PC computer with Microsoft Power Station Fortran compiler have the following form. NIT= 6 NFV= 7 NFG= 42 C= 0.536016E-12 G= 0.763704D-08 ITERM= 4 NIT= 11 NFV= 12 NFG= 168 C= 0.351663E-13 G= 0.478977D-10 ITERM= 4 NIT= 11 NFV= 12 NFG= 72 C= 0.315126E-11 G= 0.434427D-09 ITERM= 4 NIT= 20 NFV= 21 NFG= 126 C= 0.285079E-10 G= 0.555964D-08 ITERM= 4 NIT= 15 NFV= 16 NFG= 160 C= 0.626330E-08 G= 0.199604D-06 ITERM= 4 NIT= 16 NFV= 19 NFG= 238 C= 0.524025E-13 G= 0.320156D-09 ITERM= 4 NIT= 11 NFV= 13 NFG= 84 C= 0.652065E-10 G= 0.358497D-09 ITERM= 4 NIT= 40 NFV= 54 NFG= 287 C= 0.236033E-12 G= 0.569894D-06 ITERM= 4 NIT= 20 NFV= 37 NFG= 147 C= 0.323282E-06 G= 0.393792D-06 ITERM= 4 NIT= 13 NFV= 20 NFG= 84 C= 0.120160E-06 G= 0.580382D-06 ITERM= 4 NIT= 7 NFV= 8 NFG= 48 C= 0.838722E-09 G= 0.263766D-09 ITERM= 4 NIT= 6 NFV= 7 NFG= 49 C= 0.303063E-07 G= 0.569950D-06 ITERM= 4 NIT= 15 NFV= 25 NFG= 128 C= 0.278610E-07 G= 0.977278D-06 ITERM= 4 NIT= 13 NFV= 14 NFG= 98 C= 0.529941E-10 G= 0.346411D-10 ITERM= 4 NIT= 19 NFV= 22 NFG= 120 C= 0.124682E-07 G= 0.429103D-06 ITERM= 4 NIT= 7 NFV= 8 NFG= 40 C= 0.110654E-10 G= 0.142114D-10 ITERM= 4 NIT= 8 NFV= 9 NFG= 45 C= 0.121537E-12 G= 0.166850D-11 ITERM= 4 NIT= 11 NFV= 17 NFG= 60 C= 0.753158E-07 G= 0.768419D-06 ITERM= 4 NITER = 249 NFVAL = 321 NSUCC = 18 TIME= 0:00:00.79 The rows corresponding to individual test problems contain the number of iterations NIT, the number of function evaluations NFV, the number of gradient evaluations NFG, the constraint violation C, the norm of the Lagrangian function gradient G and the cause of termination ITERM. References: ----------- [1] Luksan L., Vlcek J.: LSA: Algorithms for large-scale unconstrained and box constrained optimization Technical Report V-896. Prague, ICS AS CR, 2004. [2] Luksan L., Vlcek J.: Sparse and partially separable test problems for unconstrained and equality constrained optimization. Research Report V-767, Institute of Computer Science, Academy of Sciences of the Czech Republic, Prague, Czech Republic, 1998.