Formula: x_{t+1} = x_t - H(x_t)^-1 g(x_t).
Usage
Newraf(f, grad_f, hess_f, x0, n_iter = 50, tol = 1e-12)
Arguments
- f
Objective, f(x) returning a scalar.
- grad_f
Gradient, grad_f(x) returning length p.
- hess_f
Hessian, hess_f(x) returning p by p.
- x0
Starting point, length p.
- n_iter
Maximum iterations.
- tol
Stop when the Euclidean gradient norm falls below this.
Value
List with x, estimate, fval,
grad_norm, iterations, converged, p.
Details
The step is obtained by SOLVING H d = -g, never by forming
H^-1: inverting costs three times the work and loses the
conditioning the solve keeps. A general solve is used rather than a
Cholesky factorisation, because H is only guaranteed positive
definite at a minimum, not on the way to one – on a log-likelihood
being maximised the caller passes the Hessian of the log-likelihood
and the method descends on its negative.
Determinism: a fixed iteration cap with a gradient-norm stopping
rule; no line search, no random restarts.
References
Newton, I. (1669/1711). De analysi per aequationes
numero terminorum infinitas; Raphson, J. (1690). Analysis
aequationum universalis. The modern statement is Nocedal, J. &
Wright, S. J. (2006), Numerical Optimization, 2nd ed., Springer,
algorithm 3.2.
Examples
f <- function(x) sum((x - c(1, 2))^2); g <- function(x) 2 * (x - c(1, 2)); H <- function(x) diag(2, 2)
Newraf(f, g, H, x0 = c(0, 0))$x # one step on a quadratic
#> [1] 1 2
