Checks a user-supplied Hessian of an analytic function.
GRAD —
User-supplied subroutine to compute the gradient at the point X. The usage is
CALL GRAD
(N, X,
G), where
N – Length of X and G. (Input)
X – The point at which the gradient is evaluated. X should not be changed by GRAD. (Input)
G – The gradient evaluated at the point X. (Output)
GRAD must be declared EXTERNAL in the calling program.
HESS —
User-supplied subroutine to compute the Hessian at the point X. The usage is
CALL HESS
(N, X,
H, LDH), where
N – Length of X. (Input)
X – The point at which
the Hessian is evaluated. (Input)
X should not be changed by
HESS.
H – The Hessian evaluated at the point X. (Output)
LDH – Leading dimension of H exactly as specified in in the dimension statement of the calling program. (Input)
HESS must be declared EXTERNAL in the calling program.
X — Vector of length N containing the point at which the Hessian is to be checked. (Input)
INFO — Integer matrix of dimension N by N. (Output)
INFO(I, J) = 0 means the Hessian is a poor estimate for function I at the point X(J).
INFO(I, J) = 1 means the Hessian is a good estimate for function I at the point X(J).
INFO(I, J) = 2 means the Hessian disagrees with the numerical Hessian for function I at the point X(J), but it might be impossible to calculate the numerical Hessian.
INFO(I, J) = 3 means the Hessian for function I at the point X(J) and the numerical Hessian are both zero, and, therefore, the gradient should be rechecked at a different point.
N — Dimension of
the problem. (Input)
Default: N = size
(X,1).
LDINFO — Leading
dimension of INFO exactly as
specified in the dimension statement of the calling program.
(Input)
Default: LDINFO = size
(INFO,1).
Generic: CALL CHHES (GRAD, HESS, X, INFO [,…])
Specific: The specific interface names are S_CHHES and D_CHHES.
Single: CALL CHHES (GRAD, HESS, N, X, INFO, LDINFO)
Double: The double precision name is DCHHES.
The routine CHHES uses the following finite-difference formula to estimate the Hessian of a function of n variables at x:
where hj = ɛ1/2 max{|xj|, 1/sj} sign(xj), ɛ is the machine epsilon, ej is the j-th unit vector, sj is the scaling factor of the j-th variable, and gi(x) is the gradient of the function with respect to the i-th variable.
Next, CHHES checks the user-supplied Hessian H(x) by comparing it with the finite difference approximation B(x). If
|Bij(x) - Hij(x)| < τ |Hij(x)|
where τ = ɛ1/4, then Hij(x) is declared correct; otherwise, CHHES computes the bounds of calculation error and approximation error. When both bounds are too small to account for the difference, Hij(x) is reported as incorrect. In the case of a large error bound, CHHES uses a nearly optimal stepsize to recompute Bij(x) and reports that Bij(x) is correct if
|Bij(x) - Hij(x)| < 2τ |Hij(x)|
Otherwise, Hij(x) is considered incorrect unless the error bound for the optimal step is greater than τ |Hij(x)|. In this case, the numeric approximation may be impossible to compute correctly. For more details, see Schnabel (1985).
Workspace may be explicitly provided, if desired, by use of C2HES/DC2HES. The reference is
CALL C2HES (GRAD, HESS, N, X, INFO, LDINFO, G, HX, HS, XSCALE, EPSFCN, INFT, NEWX)
The additional arguments are as follows:
G — Vector of length N containing the value of the gradient GRD at X.
HX — Real matrix of dimension N by N containing the Hessian evaluated at X.
HS — Real work vector of length N.
XSCALE — Vector of length N used to store the diagonal scaling matrix for the variables.
EPSFCN — Estimate of the relative noise in the function.
INFT — Vector of length N. For I = 1 through N, INFT contains information about the Jacobian.
NEWX — Real work array of length N.
The user-supplied Hessian of
at (-1.2, 1.0) is checked, and the error is found.
USE CHHES_INT
IMPLICIT
NONE
INTEGER LDINFO, N
PARAMETER (N=2, LDINFO=N)
!
INTEGER INFO(LDINFO,N)
REAL X(N)
EXTERNAL GRD, HES
!
! Input values for X
! X = (-1.2, 1.0)
!
DATA X/-1.2, 1.0/
!
CALL CHHES (GRD, HES, X, INFO)
!
END
!
SUBROUTINE GRD (N, X, UG)
INTEGER N
REAL X(N), UG(N)
!
UG(1) = -400.0*X(1)*(X(2)-X(1)*X(1)) + 2.0*X(1) - 2.0
UG(2) = 200.0*X(2) - 200.0*X(1)*X(1)
RETURN
END
!
SUBROUTINE HES (N, X, HX, LDHS)
INTEGER N, LDHS
REAL X(N), HX(LDHS,N)
!
HX(1,1) = -400.0*X(2) + 1200.0*X(1)*X(1) + 2.0
HX(1,2) = -400.0*X(1)
HX(2,1) = -400.0*X(1)
! A sign change is made to HX(2,2)
!
HX(2,2) = -200.0
RETURN
END
*** FATAL ERROR 1 from CHHES. The
Hessian evaluation with respect
to
*** X(2) and X(2) is
a poor estimate.
Visual Numerics, Inc. PHONE: 713.784.3131 FAX:713.781.9260 |