fminNCG

fminNCG(func, X0, fprime, fhess, [xtol=1e-5], [maxIter], [c1=1e-4], [c2=0.9])

func is the function to minimize. The return value of the function must be numeric type.

X0 is a numeric scalar or vector indicating the initial guess.

fprime is the gradient of func.

fhess is the function to compute the Hessian matrix of func.

xtol (optional) is a positive number. Convergence is assumed when the average relative error in the minimizer falls below this amount. The default value is 1e-5.

maxIter (optional) is a non-negative integer indicating the maximum number of iterations. The default value is 15000.

c1 (optional) is a number in (0,1) indicating the parameter for Armijo condition rule. The default value is 1e-4.

c2 (optional) is a number in (0,1) indicating the parameter for curvature condition rule. The default value is 0.9. Note that c2 must be greater than c1.

Perform unconstrained minimization of a function using the Newton-CG method.

Return value: A dictionary with the following members:

• xopt: A floating-point vector indicating the parameters of the minimum.

• fopt: A floating-point scalar indicating the value of func at the minimum, i.e., fopt=func(xopt).

• iterations: The number of iterations.

• fcalls: The number of function calls made.

• hcalls: The number of Hessian calls made.

• warnFlag: An integer, which can be

  • 0: Minimization performed.

  • 1: Maximum number of iterations exceeded.

  • 2: Line search failure (precision loss).

  • 3: NULL result encountered.

To minimize function rosen:

def rosen(x) { 
	N = size(x);
	return sum(100.0*power(x[1:N]-power(x[0:(N-1)], 2.0), 2.0)+power(1-x[0:(N-1)], 2.0));
}

def rosen_der(x) {
	N = size(x);
	xm = x[1:(N-1)]
	xm_m1 = x[0:(N-2)]
	xm_p1 = x[2:N]
	der = array(double, N)
	der[1:(N-1)] = (200 * (xm - xm_m1*xm_m1) - 400 * (xm_p1 - xm*xm) * xm - 2 * (1 - xm))
	der[0] = -400 * x[0] * (x[1] - x[0]*x[0]) - 2 * (1 - x[0])
	der[N-1] = 200 * (x[N-1] - x[N-2]*x[N-2])
	return der
}

def diag1(x, k) {
	N = size(x)
	m = matrix(type(x), N+1,N+1)
	if (k == 1) {
		for(i in 0:N){
			m[i, i+1] = x[i]
		}
	} else {
		for(i in 0:N){
			m[i+1, i] = x[i]
		}
	}

	return m
}

def rosen_hess(x) {
	N = size(x);
	x1= x[0:(N-1)] * 400
	H = diag1(-x1, 1) - diag1(x1, -1)
	diagonal = array(type(x), N)
	diagonal[0] = 1200 * x[0]*x[0] - 400 * x[1] + 2
	diagonal[N-1] = 200
	diagonal[1:(N-1)] = 202 + 1200 * x[1:(N-1)]*x[1:(N-1)] - 400 * x[2:N]
	H = H + diag(diagonal)
	return H
}

X0 = [4, -2.5]
fminNCG(rosen, X0, rosen_der, rosen_hess)

Output:

xopt->[0.999999966120496,0.999999932105584]
fopt->1.149654357653714E-15
iterations->34
fcalls->45
gcalls->45
hcalls->34
warnFlag->0