***********************************************************************
* *
* PHYB - HYBRID METHODS FOR UNCONCTRAINED AND LINEARLY *
* CONSTRAINED NONLINEAR LEAST SQUARES. *
* *
***********************************************************************
1. Introduction:
----------------
The double-precision FORTRAN 77 basic subroutine PHYB is designed
to find a close approximation to a local minimum of the sum of squares
F(X) = FA_1(X)**2 + FA_2(X)**2 + ... + FA_NA(X)**2
with simple bounds on variables and general linear constraints. Here
X is a vector of N variables. Here X is a vector of NF variables and
FA_I(X), 1 <= I <= NA, are twice continuously differentiable functions.
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 <= N. General linear constraints are assumed in the form
C(I) unbounded if IC(I) = 0,
CL(I) <= C(I) if IC(I) = 1,
C(I) <= CU(I) if IC(I) = 2,
CL(I) <= C(I) <= CU(I) if IC(I) = 3,
CL(I) = C(I) = CU(I) if IC(I) = 5,
where C(I) = A_I*X, 1 <= I <= NC, are linear functions. To simplify
user's work, three additional easy to use subroutines are added. They
call the basic general subroutine PHYB:
PHYBU - unconstrained optimization,
PHYBS - optimization with simple bounds,
PHYBL - optimization with simple bounds and general linear
constraints.
All subroutines contain a description of formal parameters and
extensive comments. Furthermore, two test programs THYBU and THYBL 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 methods 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,NREM,NADD,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.
NRED O Positive INTEGER variable that indicates the number of
reductions.
NREM O Positive INTEGER variable that indicates the number of
constraint deletions during the QP solutions.
NADD O Positive INTEGER variable that indicates the number of
constraint additions during the QP solutions.
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. Subroutines PHYBU, PHYBS, PHYBL:
-----------------------------------
The calling sequences are
CALL PHYBU(NF,NA,X,AF,IPAR,RPAR,F,GMAX,IPRNT,ITERM)
CALL PHYBS(NF,NA,NB,X,IX,XL,XU,AF,IPAR,RPAR,F,GMAX,IPRNT,ITERM)
CALL PHYBL(NF,NA,NB,NC,X,IX,XL,XU,CF,IC,CL,CU,CG,AF,IPAR,RPAR,
& F,GMAX,IPRNT,ITERM)
The arguments have the following meaning.
Argument Type Significance
----------------------------------------------------------------------
NF I Positive INTEGER variable that specifies the number of
variables of the objective function.
NA I Positive INTEGER variable that specifies the number of
partial functions in the sum of squares.
NB I Nonnegative INTEGER variable that specifies whether the
simple bounds are suppressed (NB=0) or accepted (NB>0).
NC I Nonnegative INTEGER variable that specifies the number
of linear constraints; if NC=0 the linear constraints
are suppressed.
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 if NB>0) 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 if NB>0).
XU(NF) I DOUBLE PRECISION vector with upper bounds for variables
(significant only if NB>0).
CF(NC) A DOUBLE PRECISION vector which contains values of
constraint functions (only if NC>0).
IC(NC) I On input (significant only if NC>0) INTEGER vector which
contains constraint types:
IC(K)=0 - the constraint CF(K) is not used,
IC(K)=1 - the lower constraint CF(K) >= CL(K),
IC(K)=2 - the upper constraint CF(K) <= CU(K),
IC(K)=3 - the two side constraint
CL(K) <= CF(K) <= CU(K),
IC(K)=5 - the equality constraint CF(K) = CL(K).
CL(NC) I DOUBLE PRECISION vector with lower bounds for constraint
functions (significant only if NC>0).
CU(NC) I DOUBLE PRECISION vector with upper bounds for constraint
functions (significant only if NC>0).
CG(NF*NC) I DOUBLE PRECISION matrix whose columns are normals of the
linear constraints (significant only if NC>0).
AF(NA) O DOUBLE PRECISION vector which contains values of partial
functions.
IPAR(6) A INTEGER parameters:
IPAR(1)=MIT, IPAR(2)=MFV, IPAR(3)-NONE,
IPAR(4)-NONE, IPAR(5)=MET, IPAR(6)=MET1.
Parameters MIT, MFV, MET, MET1 are described in Section 3
together with other parameters of the subroutine PHYB.
RPAR(7) A DOUBLE PRECISION parameters:
RPAR(1)=XMAX, RPAR(2)=TOLX, RPAR(3)=TOLF,
RPAR(4)=TOLB, RPAR(5)=TOLG, RPAR(6)-NONE,
RPAR(7)=XDEL.
Parameters XMAX, TOLX, TOLF, TOLB, TOLG, XDEL are
described in Section 3 together with other parameters of
the subroutine PHYB.
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.
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< 0 - if the method failed.
The subroutines PHYBU, PHYBS, PHYBL require the user supplied
subroutines FUN and DFUN that defines partial functions and their
gradients and have the form
SUBROUTINE FUN(NF,KA,X,FA)
SUBROUTINE DFUN(NF,KA,X,GA)
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 PHYB:
-------------------
This general subroutine is called from all the subroutines described
in Section 2. The calling sequence is
CALL PHYB(NF,NA,NB,NC,X,IX,XL,XU,CF,IC,CL,CU,CG,ICA,CFD,CR,CZ,AF,
& GA,G,GN,H,HH,B,S,SN,XO,GO,XMAX,TOLX,TOLF,TOLB,TOLG,XDEL,GMAX,F,
& MIT,MFV,MET,MET1,IPRNT,ITERM).
The arguments NF, NA, NB, NC, X, IX, XL, XU, CF, IC, CL, CU, CG, AF,
GMAX, F, IPRNT, ITERM, have the same meaning as in Section 2. Other
arguments have the following meaning:
Argument Type Significance
----------------------------------------------------------------------
ICA(NC) A INTEGER vector containing indices of active constraints.
CFD(NC) A DOUBLE PRECISION vector of constraint function
increments.
CR(NCR) A DOUBLE PRECISION matrix containing triangular
decomposition of the orthogonal projection kernel
(NCR is equal to NF*(NF+1)/2).
CZ(NF*NF) A DOUBLE PRECISION matrix containing orthogonal basis
of linear manifold defined by active constraints.
GA(NF) A DOUBLE PRECISION gradient of the partial function.
G(NF) A DOUBLE PRECISION gradient of the objective function.
GN(NF) A DOUBLE PRECISION reduced gradient of the objective
function.
H(NH) A DOUBLE PRECISION Gauss-Newton approximation of the
Hessian matrix.
HH(NH) A DOUBLE PRECISION variable metric approximation of the
second order term.
B(NH) A DOUBLE PRECISION auxiliary matrix which is decomposed.
S(NF) A DOUBLE PRECISION direction vector.
SN(NF) A DOUBLE PRECISION reduced direction vector.
XO(NF) A DOUBLE PRECISION vector which contains increments of
variables.
GO(NF) A DOUBLE PRECISION vector which contains increments of
gradients.
XMAX I DOUBLE PRECISION maximum stepsize; the choice XMAX=0
causes that the default value 1.0D+3 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.
TOLF I DOUBLE PRECISION tolerance for the change of function
values; the choice TOLF=0 causes that the default
value TOLF=1.0D-16 will be taken.
TOLB I DOUBLE PRECISION minimum acceptable function value;
the choice TOLB=0 causes that the default value
TOLB=FMIN+1.0D-16 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.
XDEL I DOUBLE PRECISION trust region radius.
MIT I INTEGER variable that specifies the maximum number of
iterations; the choice MIT=0 causes that the default
value 200 will be taken.
MFV I INTEGER variable that specifies the maximum number of
function evaluations; the choice |MFV|=0 causes that
the default value 500 will be taken.
MET I INTEGER variable that specifies the method used:
MET=1 - the Gauss-Newton method is used.
MET=2 - the structured rank one method is used.
MET=3 - the Huschens modification of the structured
rank one method is used.
MET=4 - a combination of the Gauss-Newton and the
BFGS nethods is used.
The choice MET=0 causes that the default value MET=4
will be taken.
MET1 I INTEGER variable that specifies scaling strategy:
MET1=1 - scaling is not used.
MET1=2 - scaling in the first iteration is used.
MET1=3 - controled scaling is used.
MET1=4 - permanent scaling is used.
The choice MET1=0 causes that the default value MET1=3
will be taken.
The subroutine PHYB requires the user supplied subroutines FUN and DFUN
which are described in Section 2.
4. Verification of the subroutines:
-----------------------------------
Subroutine PHYBU can be verified and tested using the program THYBU.
This program calls the subroutines TIUD15 (initiation), TAFU15 (partial
function evaluation), TAGU15 (partial gradient evaluation) containing 22
academic test problems with at most 100 variables [2]. The results obtained
by the program THYBU on a PC computer with Microsoft Power Station Fortran
compiler have the following form.
NIT= 125 NFV= 129 NFG= 129 F= 0.228324632E-21 G= 0.235E-09 ITERM= 3
NIT= 69 NFV= 81 NFG= 81 F= 0.154234698E-13 G= 0.272E-05 ITERM= 4
NIT= 9 NFV= 10 NFG= 10 F= 0.779710653E-08 G= 0.393E-05 ITERM= 4
NIT= 20 NFV= 24 NFG= 24 F= 12.6030647 G= 0.528E-05 ITERM= 4
NIT= 4 NFV= 5 NFG= 5 F= 0.619992637E-16 G= 0.679E-08 ITERM= 3
NIT= 5 NFV= 6 NFG= 6 F= 0.182897736E-17 G= 0.453E-08 ITERM= 3
NIT= 14 NFV= 21 NFG= 21 F= 5982.28867 G= 0.832E-05 ITERM= 4
NIT= 12 NFV= 15 NFG= 15 F= 0.108254394E-05 G= 0.823E-05 ITERM= 4
NIT= 22 NFV= 26 NFG= 26 F= 217.459747 G= 0.568E-05 ITERM= 4
NIT= 18 NFV= 26 NFG= 26 F= 19.0247088 G= 0.461E-05 ITERM= 4
NIT= 356 NFV= 363 NFG= 363 F= 0.761416163E-24 G= 0.124E-10 ITERM= 3
NIT= 26 NFV= 35 NFG= 35 F= 2146.09845 G= 0.414E-04 ITERM= 2
NIT= 28 NFV= 35 NFG= 35 F= 12594.4298 G= 0.375E-05 ITERM= 4
NIT= 15 NFV= 16 NFG= 16 F= 10.6560652 G= 0.289E-05 ITERM= 4
NIT= 6 NFV= 7 NFG= 7 F= 1.78484853 G= 0.236E-07 ITERM= 4
NIT= 16 NFV= 21 NFG= 21 F= 0.246580645 G= 0.636E-05 ITERM= 4
NIT= 7 NFV= 8 NFG= 8 F= 0.159973696E-09 G= 0.747E-05 ITERM= 4
NIT= 17 NFV= 18 NFG= 18 F= 0.367899931E-15 G= 0.290E-06 ITERM= 3
NIT= 12 NFV= 13 NFG= 13 F= 0.252033641E-12 G= 0.215E-05 ITERM= 4
NIT= 14 NFV= 15 NFG= 15 F= 0.144405321E-12 G= 0.406E-07 ITERM= 4
NIT= 20 NFV= 25 NFG= 25 F= 61.9617704 G= 0.522E-05 ITERM= 4
NIT= 61 NFV= 72 NFG= 72 F= 436.970239 G= 0.736E-05 ITERM= 4
NITER = 876 NFVAL = 971 NSUCC = 22
TIME= 0:00:03.14
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 PHYBL can be verified and tested using the program
THYBL. This program calls the subroutines TILD22 (initiation), TAFU22
(partial function evaluation), TAGU22 (partial gradient evaluation)
containing 7 academic test problems with at most 16 variables. The results
obtained by the program THYBL on a PC computer with Microsoft Power
Station Fortran compiler have the following form.
NIT= 6 NFV= 7 NFG= 7 F= 0.109474579 G= 0.199E-05 ITERM= 4
NIT= 6 NFV= 7 NFG= 7 F= 0.472606902 G= 0.778E-07 ITERM= 4
NIT= 17 NFV= 18 NFG= 18 F= 1.61649353 G= 0.222E-06 ITERM= 4
NIT= 8 NFV= 9 NFG= 9 F= 1.90396173 G= 0.255E-05 ITERM= 4
NIT= 17 NFV= 28 NFG= 28 F= 0.678270395E-11 G= 0.431E-05 ITERM= 4
NIT= 14 NFV= 16 NFG= 16 F= 0.234493933 G= 0.114E-05 ITERM= 4
NIT= 6 NFV= 7 NFG= 7 F= 0.113241803E-20 G= 0.435E-08 ITERM= 3
NITER = 74 NFVAL = 92 NSUCC = 7
TIME= 0:00:00.00
References:
-----------
[1] Luksan L., Spedicato E.: Variable metric methods for unconstrained
optimization and nonlinear least squares. Journal of Computational
and Applied Mathematics 124 (2000) 61-93.
[2] Luksan L., Vlcek J.: Test problems for unconstrained optimization.
Research Report V-897, Institute of Computer Science, Academy of
Sciences of the Czech Republic, Prague, Czech Republic, 2003.