*********************************************************************** * * * PGAC - HYBRID GAUSS-NEWTON METHOD WITH SECOND-ORDER CORRECTIONS * * AND ITERATIVE CG-BASED TRUST-REGION SUBALGORITHMS FOR * * LARGE-SCALE PARTIALLY SEPARABLE LEAST SQUARES PROBLEMS. * * * *********************************************************************** 1. Introduction: ---------------- The double-precision FORTRAN 77 basic subroutine PGAC is designed to find a close approximation to a local minimum of a sum of squares F(X) = FA_1(X)**2 + FA_2(X)**2 + ... + FA_NA(X)**2 with simple bounds on variables. Here X is a vector of NF variables and FA_I(X), 1 <= I <= NA, are twice continuously differentiable functions. We assume that NF and NA are large, but partial functions FA_I(X), 1 <= I <= NA depend on a small number of variables. This implies that the mapping AF(X) = [FA_1(X), FA_2(X), ..., FA_NA(X)] has a sparse Jacobian matrix, which will be denoted by AG(X) (it has NA rows and NF columns). Simple bounds are assumed in the form X(I) unbounded if IX(I) = 0, XL(I) <= X(I) if IX(I) = 1, X(I) <= XU(I) if IX(I) = 2, XL(I) <= X(I) <= XU(I) if IX(I) = 3, XL(I) = X(I) = XU(I) if IX(I) = 5, where 1 <= I <= NF. The sparsity pattern of the Jacobian matrix is stored in the coordinate form if ISPAS=1 or in the standard compressed row format if ISPAS=2 using arrays IAG and JAG. For example, if the Jacobian matrix has the following pattern AG = | * * 0 * | | * * * 0 | | * 0 0 * | | 0 * * 0 | | * 0 * 0 | (asterisks denote nonzero elements) then arrays IAG and JAG contain elements IAG(1)=1, IAG(2)=1, IAG(3)=1, IAG(4)=2, IAG(5)=2, IAG(6)=2, IAG(7)=3, IAG(8)=3, IAG(9)=4, IAG(10)=4, IAG(11)=5, IAG(12)=5, JAG(1)=1, JAG(2)=2, JAG(3)=4, JAG(4)=1, JAG(5)=2, JAG(6)=3, JAG(7)=1, JAG(8)=4, JAG(9)=2, JAG(10)=3, JAG(11)=1, JAG(12)=3 if ISPAS=1 or IAG(1)=1, IAG(2)=4, IAG(3)=7, IAG(4)=9, IAG(5)=11, IAG(6)=13, JAG(1)=1, JAG(2)=2, JAG(3)=4, JAG(4)=1, JAG(5)=2, JAG(6)=3, JAG(7)=1, JAG(8)=4, JAG(9)=2, JAG(10)=3, JAG(11)=1, JAG(12)=3 if ISPAS=2. In the first case, nonzero elements can be sorted in an arbitrary order (not only by rows as in the above example). Arrays IAG and JAG have to be declared with lengths NA+MA and MA at least, respectively, where MA is the number of nonzero elements. In the second case, nonzero elements can be sorted only by rows. Components of IAG contain total numbers of nonzero elements in all previous rows increased by 1 and elements of JAG contain corresponding column indices (note that IAG has NA+1 elements and the last element is equal to MA+1). Arrays IAG and JAG have to be declared with length NA+1 and MA at least, respectively. To simplify user's work, two additional easy to use subroutines are added. They call the basic general subroutine PGAC: PGACU - unconstrained large-scale optimization, PGACS - large-scale optimization with simple bounds. All subroutines contain a description of formal parameters and extensive comments. Furthermore, two test programs TGACU and TGACS are included, which contain several test problems (see e.g. [2]). These test programs serve as examples for using the subroutines, verify their correctness and demonstrate their 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). 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,NIN,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. NIN 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 indicates the number of gradient evaluations. NFH O Positive INTEGER variable that indicates the number of Hessian evaluations. 2. Subroutines PGACU, PGACS: ---------------------------- The calling sequences are CALL PGACU(NF,NA,MA,X,AF,IAG,JAG,IPAR,RPAR,F,GMAX,IDER,ISPAS,IPRNT, & ITERM) CALL PGACS(NF,NA,MA,X,IX,XL,XU,AF,IAG,JAG,IPAR,RPAR,F,GMAX,IDER, & ISPAS,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. NA I Positive INTEGER variable that specifies the number of partial functions. MA I Number of nonzero elements in the Jacobian matrix. This parameter is used as input only if ISPAS=1 (it defines dimensions of arrays IAG and JAG in this case). X(NF) U On input, DOUBLE PRECISION vector with the initial estimate to the solution. On output, the approximation to the minimum. IX(NF) I On input (significant only for PGACS) INTEGER vector containing the simple bounds types: IX(I)=0 - the variable X(I) is unbounded, IX(I)=1 - the lower bound X(I) >= XL(I), IX(I)=2 - the upper bound X(I) <= XU(I), IX(I)=3 - the two side bound XL(I) <= X(I) <= XU(I), IX(I)=5 - the variable X(I) is fixed (given by its initial estimate). XL(NF) I DOUBLE PRECISION vector with lower bounds for variables (significant only for PGACS). XU(NF) I DOUBLE PRECISION vector with upper bounds for variables (significant only for PGACS). AF(NA) O DOUBLE PRECISION vector which contains values of partial functions. IAG(NA+1) I INTEGER array which contains pointers of the first elements in rows of the Jacobian matrix. JAG(MA) I INTEGER array which contains column indices of the nonzero elements. IPAR(7) U INTEGER parameters: IPAR(1)=MIT, IPAR(2)=MFV, IPAR(3)=MFG, IPAR(4)=MEC, IPAR(5)=MOS1, IPAR(6)=MOS2, IPAR(7)=IFIL. Parameters MIT, MFV, MFG, MEC, MOS1, MOS2 are described in Section 3 together with other parameters of the subroutine PGAC. 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(9) U DOUBLE PRECISION parameters: RPAR(1)=XMAX, RPAR(2)=TOLX, RPAR(3)=TOLF, RPAR(4)=TOLB, RPAR(5)=TOLG, RPAR(6)=FMIN, RPAR(7)=XDEL, RPAR(8)=ETA, RPAR(9)-NONE. Parameters XMAX, TOLX, TOLF, TOLB, TOLG, FMIN, XDEL, ETA are described in Section 3 together with other parameters of the subroutine PGAC. 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 objective function. IDER I INGEGER variable that specifies the degree of analytically computed derivatives (0 OR 1). ISPAS I INTEGER variable that specifies sparse structure of the Jacobian matrix: ISPAS= 1 - the coordinate form is used, ISPAS= 2 - the standard row compresed format is used. 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= 6 - if termination criterion was not satisfied, but the solution is probably acceptable, ITERM=11 - if NIT exceeded MIT, ITERM=12 - if NFV exceeded MFV, ITERM=13 - if NFG exceeded MFG, ITERM< 0 - if the method failed. 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 subroutines PGACU, PGACS require the user supplied subroutines FUN and DFUN that define partial functions and their gradients and have the form SUBROUTINE FUN(NF,KA,X,FA) SUBROUTINE DFUN(NF,KA,X,GA) If IDER=0, the subroutine DFUN can be empty. 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. KA I INTEGER index of the partial function. X(NF) I DOUBLE PRECISION an estimate to the solution. FA O DOUBLE PRECISION value of the KA-th partial function at the point X. GA(NF) O DOUBLE PRECISION gradient of the KA-th partial function at the point X. Note that only nonzero elements of this gradient have to be assigned. 3. Subroutine PGAC: ------------------- This general subroutine is called from all subroutines described in Section 2. The calling sequence is CALL PGAC(NF,NA,NB,MMAX,X,IX,XL,XU,AF,GA,AG,G,HA,AH,H,IH,JH,IAG, & JAG,S,XO,GO,AGO,XS,GS,IW,XMAX,TOLX,TOLF,TOLB,TOLG,FMIN,XDEL,ETA, & GMAX,F,MIT,MFV,MFG,MEC,MOS1,MOS2,IDER,IPRNT,ITERM) The arguments NF, NA, X, IX, XL, XU, AF, IAG, JAG, GMAX, F, IDER, IPRNT, ITERM have the same meaning as in Section 2. Other arguments have the following meaning: Argument Type Significance ---------------------------------------------------------------------- NB I Nonnegative INTEGER variable that specifies whether the simple bounds are suppressed (NB=0) or accepted (NB>0). MMAX I INTEGER size of array H. GA(NF) A DOUBLE PRECISION gradient of the partial function. AG(MA) A DOUBLE PRECISION nonzero elements of the Jacobian matrix. This array is used only if MEC=3. G(NF) A DOUBLE PRECISION gradient of the objective function. HA(ML) A DOUBLE PRECISION Hessian matrix of the partial function. AH(MH) A DOUBLE PRECISION approximation of the partitioned Hessian matrix. H(MMAX) A DOUBLE PRECISION nonzero elements of the approximation of the Hessian matrix and nonzero elements of the Choleski factor. IH(NF+1) I INTEGER array which contains pointers of the diagonal elements in the upper part of the Hessian matrix. JH(MMAX) I INTEGER array which contains column indices of the nonzero elements and additional working space for the Choleski factor. S(NF) A DOUBLE PRECISION direction vector. XO(NF) A DOUBLE PRECISION array which contains increments of variables. GO(NF) A DOUBLE PRECISION array which contains increments of gradients. AGO(MA) A DOUBLE PRECISION difference between the current and the old Jacobian matrices. This array is used only if MEC=3. XS(NF) A DOUBLE PRECISION auxiliary array. GS(NF) A DOUBLE PRECISION auxiliary array. IW(NF+1) A INTEGER auxiliary array. XMAX U DOUBLE PRECISION maximum stepsize; the choice XMAX=0 causes that the default value 1.0D+16 will be taken. TOLX U 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. TOLF U DOUBLE PRECISION tolerance for the change of function values; the choice TOLF=0 causes that the default value TOLF=1.0D-14 will be taken. TOLB U DOUBLE PRECISION minimum acceptable function value; the choice TOLB=0 causes that the default value TOLB=FMIN+1.0D-16 will be taken. TOLG U DOUBLE PRECISION tolerance for the Lagrangian function gradient; the choice TOLG=0 causes that the default value TOLG=1.0D-6 will be taken. FMIN U DOUBLE PRECISION lower bound for the minimum function value; the choice FMIN<0 causes that the default value value TOLG=0.0D 0 will be taken. XDEL U DOUBLE PRECISION trust region stepsize; the choice XDEL=0 causes that a suitable default value is computed. ETA U DOUBLE PRECISION parameter for switch between the Gauss-Newton method and variable metric correction; the choice ETA=0 causes that the default value ETA=1.5D-4 will be taken. MIT U INTEGER variable that specifies the maximum number of iterations; the choice MIT=0 causes that the default value 5000 will be taken. MFV U INTEGER variable that specifies the maximum number of function evaluations; the choice MFV=0 causes that the default value 5000 will be taken. MFG U INTEGER variable that specifies the maximum number of gradient evaluations; the choice MFV=0 causes that the default value 10000 will be taken. MEC U INTEGER method of a second order correction: MEC=1 - correction by the Marwil sparse variable metric update, MEC=2 - correction by differences of gradients (discrete Newton correction). MEC=3 - correction by the Griewank-Toint partitioned variable metric update (symmetric rank-one). This correction uses three additional matrices (arrays AG, AGO and AH). The choice MEC=0 causes that the default value 2 will be taken. MOS1 U INTEGER method for computing trust-region step: MOS1=1 - Steihaug-Toint conjugate gradient method, MOS1=2 - shifted Steihaug-Toint method with five Lanczos steps. MOS1>2 - shifted Steihaug-Toint method with MOS1 Lanczos steps. The choice MOS1=0 causes that the default value 2 will be taken. MOS2 U INTEGER variable defining a type of preconditioning. MOS2=1 - Preconditioning is not used. MOS2=2 - Preconditioning by the incomplete Gill-Murray decomposition. MOS2=3 - Preconditioning by the incomplete Gill-Murray decomposition with a preliminary solution of the preconditioned system which is used if it satisfies the termination criterion. The choice MOS2=0 causes that the default value 2 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. The subroutine PGAC requires the user supplied subroutines FUN and DFUN which are described in Section 2. 4. Verification of the subroutines: ----------------------------------- Subroutine PGACU can be verified and tested using the program TGACU. This program calls the subroutines TIUB15 (initiation), TAFU15 (function evaluation) and TAGU15 (gradient evaluation) containing 22 unconstrained test problems with at most 1000 variables [2]. The results obtained by the program TGACU on a PC computer with Microsoft Power Station Fortran compiler have the following form. NIT= 1108 NFV= 1110 NFG= 1110 F= 0.00000000 G= 0.000E+00 ITERM= 3 NIT= 624 NFV= 640 NFG= 649 F= 66.4089431 G= 0.283E-07 ITERM= 4 NIT= 11 NFV= 12 NFG= 14 F= 0.202412411E-09 G= 0.210E-06 ITERM= 4 NIT= 11 NFV= 13 NFG= 17 F= 134.749772 G= 0.592E-07 ITERM= 4 NIT= 4 NFV= 5 NFG= 7 F= 0.116836300E-10 G= 0.908E-06 ITERM= 4 NIT= 6 NFV= 7 NFG= 13 F= 0.787821743E-26 G= 0.311E-12 ITERM= 3 NIT= 17 NFV= 40 NFG= 29 F= 60734.8551 G= 0.428E-05 ITERM= 6 NIT= 22 NFV= 25 NFG= 25 F= 0.127726626E-07 G= 0.160E-06 ITERM= 4 NIT= 13 NFV= 15 NFG= 38 F= 2216.45871 G= 0.846E-06 ITERM= 4 NIT= 129 NFV= 147 NFG= 176 F= 191.511336 G= 0.104E-07 ITERM= 4 NIT= 3010 NFV= 3016 NFG= 3012 F= 0.402368464E-24 G= 0.902E-11 ITERM= 3 NIT= 205 NFV= 226 NFG= 236 F= 22287.9069 G= 0.449E-08 ITERM= 4 NIT= 123 NFV= 132 NFG= 152 F= 131234.018 G= 0.743E-09 ITERM= 4 NIT= 7 NFV= 8 NFG= 32 F= 108.517888 G= 0.148E-07 ITERM= 4 NIT= 13 NFV= 20 NFG= 42 F= 18.1763146 G= 0.445E-05 ITERM= 2 NIT= 14 NFV= 15 NFG= 35 F= 2.51109677 G= 0.103E-09 ITERM= 4 NIT= 29 NFV= 34 NFG= 33 F= 0.139780007E-09 G= 0.238E-06 ITERM= 4 NIT= 49 NFV= 53 NFG= 52 F= 0.119511868E-21 G= 0.344E-09 ITERM= 3 NIT= 15 NFV= 16 NFG= 23 F= 0.339085307E-13 G= 0.788E-06 ITERM= 4 NIT= 17 NFV= 18 NFG= 32 F= 0.336618309E-11 G= 0.137E-07 ITERM= 4 NIT= 15 NFV= 18 NFG= 23 F= 647.696136 G= 0.262E-06 ITERM= 4 NIT= 47 NFV= 59 NFG= 98 F= 4486.97024 G= 0.663E-08 ITERM= 4 NITER = 5489 NFVAL = 5629 NITCG = 8282 NSUCC = 22 TIME= 0:00:05.01 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 final value of the objective function F, the norm of gradient G and the cause of termination ITERM. Subroutine PGACS can be verified and tested using the program TGACS. This program calls the subroutines TIUB15 (initiation), TAFU15 (function evaluation), TAGU15 (gradient evaluation) containing 21 box constrained test problems with at most 1000 variables [2]. The results obtained by the program TGACS on a PC computer with Microsoft Power Station Fortran compiler have the following form. NIT= 1028 NFV= 1031 NFG= 1030 F= 0.00000000 G= 0.000E+00 ITERM= 3 NIT= 263 NFV= 262 NFG= 511 F= 1959.28649 G= 0.819E-10 ITERM= 4 NIT= 10 NFV= 12 NFG= 13 F= 0.814133878E-09 G= 0.897E-06 ITERM= 4 NIT= 11 NFV= 14 NFG= 17 F= 134.761343 G= 0.206E-10 ITERM= 4 NIT= 4 NFV= 5 NFG= 7 F= 0.438081882E-11 G= 0.697E-06 ITERM= 4 NIT= 6 NFV= 7 NFG= 13 F= 0.791460667E-17 G= 0.934E-08 ITERM= 3 NIT= 15 NFV= 16 NFG= 42 F= 145814.000 G= 0.000E+00 ITERM= 4 NIT= 17 NFV= 18 NFG= 20 F= 0.201167216E-07 G= 0.162E-06 ITERM= 4 NIT= 54 NFV= 55 NFG= 203 F= 2220.17880 G= 0.614E-10 ITERM= 4 NIT= 95 NFV= 105 NFG= 126 F= 191.511336 G= 0.527E-06 ITERM= 4 NIT= 4577 NFV= 3591 NFG= 3587 F= 0.00000000 G= 0.000E+00 ITERM= 3 NIT= 49 NFV= 50 NFG= 84 F= 67887.2385 G= 0.351E-07 ITERM= 4 NIT= 18 NFV= 19 NFG= 41 F= 147906.000 G= 0.000E+00 ITERM= 4 NIT= 1 NFV= 2 NFG= 6 F= 126.690556 G= 0.000E+00 ITERM= 4 NIT= 33 NFV= 72 NFG= 85 F= 18.1763146 G= 0.182E-04 ITERM= 6 NIT= 8 NFV= 12 NFG= 29 F= 3.59074140 G= 0.542E-06 ITERM= 4 NIT= 25 NFV= 29 NFG= 29 F= 0.716457302E-10 G= 0.922E-07 ITERM= 4 NIT= 0 NFV= 1 NFG= 3 F= 0.00000000 G= 0.000E+00 ITERM= 3 NIT= 28 NFV= 32 NFG= 36 F= 0.209831733E-13 G= 0.620E-06 ITERM= 4 NIT= 937 NFV= 938 NFG= 2806 F= 498.800124 G= 0.572E-12 ITERM= 4 NIT= 21 NFV= 22 NFG= 34 F= 649.598077 G= 0.324E-08 ITERM= 4 NIT= 44 NFV= 51 NFG= 89 F= 4488.96148 G= 0.645E-10 ITERM= 4 NITER = 7244 NFVAL = 6344 NITCG =11206 NSUCC = 22 TIME= 0:00:06.56 References: ----------- [1] Luksan L., Matonoha C., Vlcek J.: LSA: Algorithms for large-scale unconstrained and box constrained optimization. Research Report V-896, Institute of Computer Science, Academy of Sciences of the Czech Republic, Prague, Czech Republic, 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.