Newton-type Methods for Unconstrained Optimization

Reference: Section 2.4.2 in Biegler (2010)

Newton models for equations and optimization¶

A row-oriented Jacobian for equations¶

Let $D\subseteq\mathbb R^n$ be open and $\mathbf F:D\to\mathbb R^n$ differentiable. Define $(\mathbf J_F(\mathbf x))_{ij}=\partial F_i(\mathbf x)/\partial x_j$ for $i,j\in\{1,\ldots,n\}$. Linearizing at fixed $\mathbf x^k$ gives

with the quadratic remainder when the Jacobian is locally Lipschitz. Newton sets the linear model to zero and solves

The actual residual at the update need not be zero. The row-oriented Jacobian requires no extra transpose. An inverse is useful in analysis, but a numerical implementation solves the linear system.

A quadratic model for a scalar objective¶

For $f\in C^2(D)$ on an open $D$, write $\mathbf g^k=\nabla f(\mathbf x^k)$ and $\mathbf G^k=\nabla^2f(\mathbf x^k)$. With $\mathbf x^k$ fixed,

Differentiating the model with respect to the chosen displacement gives

This computes a model stationary point. It is the model minimum if $\mathbf G^k\succ0$. The model equals a quadratic objective exactly; for a nonlinear objective with locally Lipschitz Hessian its value remainder is $O(\|\mathbf p\|^3)$.

Sources: Biegler (2010), Eqs. (2.37)–(2.39), pp. 34–35; Nocedal and Wright (2006), §11.1, pp. 274–277. The gradient in the numerical examples below is a column vector.

The existing nonlinear-model sketch uses $f(x)=e^x-2x$ at $x^k=1.8$. The Newton displacement is approximately -0.6694, giving $x^{k+1}=1.1306$; the exact minimizer is the natural logarithm of 2. The objective drops from approximately 2.45 to 0.84, but the model stationary point still differs from the true minimizer.

A quadratic Newton step in Python¶

For $f(\mathbf x)=(x_1-1)^2+4(x_2+2)^2$ at $\mathbf x^0=[3,0]^T$, the gradient is $[4,16]^T$ and the Hessian is $\operatorname{diag}(2,8)$. Solve the Newton equation, then add the displacement to the current iterate. This cell runs on its own; no helper functions are needed.

The step is $[-2,-2]^T$ and the updated point is $[1,-2]^T$. A quadratic objective equals its quadratic model everywhere, so one exact Newton step reaches its stationary point when the Hessian is nonsingular.

Checking a gradient norm in Python¶

Compute the gradient at the updated point using the same quadratic objective. A small gradient norm checks approximate stationarity; classification still requires curvature information. This cell also runs independently.

A local model of the scalar website example¶

The scalar function from Descent and Globalization is

with $f'(x)=6-8x-3x^2+2x^3$ and $f''(x)=-8-6x+6x^2$. At $x^k=-3$, the values are $f=15$, $f'=-51$, and $f''=64$. The quadratic model is $m_k(p)=15-51p+32p^2$, whose stationary displacement is $p^k=51/64$. The true objective at $x^{k+1}=-2.203125$ is approximately -8.66, and the true slope is approximately -12.32. The model slope is zero at this point; the objective slope is not.

The globalization notebook compares Newton, steepest descent, and modified curvature on this same function.

Curvature, descent, and step length¶

A strict local minimum need not have a positive-definite Hessian: $f(x)=x^4$ at zero is a counterexample. The sufficient second-order conditions explicitly require a zero gradient and positive-definite Hessian. Continuity then preserves positive definiteness in a small neighborhood.

For $\mathbf g\ne0$, a direction is descent when $\mathbf g^T\mathbf p<0$. If $\mathbf G\succ0$ and $\mathbf G\mathbf p=-\mathbf g$, then

This guarantees decrease for sufficiently small positive step lengths. It does not guarantee decrease for every $\alpha\in(0,1]$. For a $C^2$ objective, the exact Taylor identity is

At fixed $\mathbf x,\mathbf p$, the negative $O(\alpha)$ term dominates the bounded $O(\alpha^2)$ remainder for small $\alpha$. The intermediate point includes $\alpha$.

A singular Hessian does not specify a unique basic Newton step. An indefinite Hessian can produce an ascent direction; an inaccurate positive-definite model can still overshoot. These are distinct reasons to modify curvature or control the accepted step length.

Sources: Biegler (2010), Eq. (2.40), p. 34, and §3.1, pp. 39–40; Nocedal and Wright (2006), pp. 22 and 46–47.

Algorithm 2.1: basic unit-step Newton¶

The following reference version adapts Biegler’s Algorithm 2.1 so each quantity is evaluated before it is tested.

1. Choose $\mathbf x^0$, gradient tolerance $\epsilon_g$, and step tolerance $\epsilon_p$.

2. Evaluate $\mathbf g^k$ and $\mathbf G^k$. A small gradient reports an approximately stationary candidate.

3. If the Hessian is singular, report failure of the basic solve; otherwise solve $\mathbf G^k\mathbf p^k=-\mathbf g^k$.

4. A small step reports limited numerical movement, together with the current gradient norm.

5. Otherwise take $\mathbf x^{k+1}=\mathbf x^k+\mathbf p^k$ and repeat.

A candidate must still be classified; neither a small step nor an approximate Hessian proves a minimum. The diagnostic implementation below retains a post-step convergence test and also checks NaN derivatives, severe conditioning, and an iteration limit. Its experiments deliberately include failure cases.

Source: Biegler (2010), Algorithm 2.1, p. 35.

Helper Function¶

Algorithm 1 Python Implementation¶

Quantified rates of convergence¶

The prefix Q means quotient: the ratios $e_{k+1}/e_k$ and $e_{k+1}/e_k^2$ compare successive nonzero errors.

Let $\mathbf x^k\to\mathbf x^*$ and $e_k=\|\mathbf x^k-\mathbf x^*\|$. The convergence is at least Q-linear when

It is at least Q-quadratic when

These are eventual upper bounds; a sequence can satisfy a stronger rate too. For $e_k=10^{-3}$, $r=0.5$ gives $e_{k+1}\le5\times10^{-4}$, while $K=1$ gives $e_{k+1}\le10^{-6}$. Strict minimality by itself supplies neither the rate nor the theorem’s nondegeneracy: unit-step Newton on $f(x)=x^4$ has $x^{k+1}=2x^k/3$ for $x^k\ne0$, so its exact error ratio is $2/3$, and $f''(0)=0$.

Source: Nocedal and Wright (2006), Appendix A.2, Eq. (A.34), p. 619. The following theorem and proofs specify the neighborhood and constants required for Newton’s quadratic bounds.

Quadratic Convergence Proof¶

Theorem 2.20 (Biegler, 2010)

Assume $f\in C^2(D)$ on an open neighborhood $D$ of $\mathbf x^*$, with Lipschitz Hessian there, $\nabla f(\mathbf x^*)=\mathbf0$, and $\nabla^2f(\mathbf x^*)\succ0$. Then, by applying Algorithm 2.1 and with a sufficiently close $\mathbf{x}_0$, there exists a constant $\hat{L} > 0$ such that:

1. The convergence rate for $\{\mathbf{x}^k\}$ is quadratic, i.e., $\| \mathbf{x}^{k+1} - \mathbf{x}^* \| \leq \hat{L} \| \mathbf{x}^k - \mathbf{x}^* \|^2$

2. The convergence rate for $\{\nabla f(\mathbf{x}^k)\}$ is also quadratic

Iterate-error proof and the integral-norm tools¶

The goal is to show that unit-step Newton is well defined near $\mathbf x^*$, that its iterates remain in that neighborhood and converge, and that the resulting error bound is quadratic. Use the local convergence theorem: $f\in C^2(D)$ on an open neighborhood of $\mathbf x^*$, $\nabla f(\mathbf x^*)=\mathbf0$, the Hessian is Lipschitz there, and $\nabla^2f(\mathbf x^*)\succ0$.

Use the Euclidean vector norm and its induced matrix norm. The tools are

The last bound is the triangle inequality for a vector integral. Together these give

1. Fix a neighborhood and its constants. Since $f\in C^2$, the symmetric Hessian is continuous. Positive definiteness at $\mathbf x^*$ gives $\lambda_{\min}(\nabla^2f(\mathbf x^*))>0$. By continuity of the eigenvalues, choose $r>0$ with $B=B(\mathbf x^*,r)=\{\mathbf x:\|\mathbf x-\mathbf x^*\|<r\}\subset D$ and $m>0$ such that $\lambda_{\min}(\nabla^2f(\mathbf x))\ge m$ for every $\mathbf x\in B$. Hence every Hessian in $B$ is invertible and

For a symmetric positive definite matrix, its inverse has Euclidean induced norm $1/\lambda_{\min}$. The Hessian Lipschitz assumption supplies a fixed constant $L>0$ with

If the Hessian is constant, any positive $L$ is a valid upper-bound constant. Fix this ball and these constants before choosing the smaller initial radius in Step 6.

2. Derive the exact gradient identity. Write $\mathbf e_k=\mathbf x^k-\mathbf x^*$, $\mathbf g^k=\nabla f(\mathbf x^k)$, and $\mathbf G^k=\nabla^2f(\mathbf x^k)$. Suppose $\mathbf x^k\in B$ and define

A ball is convex, so this entire segment stays in $B$. For $\boldsymbol\psi(t)=\nabla f(\mathbf z(t))$, the chain rule gives $\boldsymbol\psi'(t)=\nabla^2f(\mathbf z(t))\mathbf e_k$. Apply the one-variable fundamental theorem of calculus (FTC) to each gradient component:

Stationarity makes $\nabla f(\mathbf x^*)=\mathbf0$. More generally, the same chain-rule/FTC argument gives

In one dimension, $u=x+tp$, $du=p\,dt$ turns the integral into $\int_x^{x+p}f''(u)\,du=f'(x+p)-f'(x)$; its oriented limits include negative $p$.

3. Rewrite the next iterate error as a Hessian difference. Unit-step Newton solves $\mathbf G^k\mathbf p^k=-\mathbf g^k$ and updates $\mathbf x^{k+1}=\mathbf x^k+\mathbf p^k$. The inverse exists by Step 1. Thus

The second line factors the inverse using $\mathbf e_k=(\mathbf G^k)^{-1}\mathbf G^k\mathbf e_k$. For the last line, write the constant $\mathbf G^k\mathbf e_k$ as its integral over $[0,1]$, then substitute Step 2. This is an exact identity before any inequality is applied.

4. Use the remaining segment length in the Lipschitz bound. Subtracting the intermediate point gives

Both Hessian arguments are in $B$ by Step 2, which is why the local Lipschitz bound applies.

5. Take norms and integrate. Apply the induced-norm bound and the integral triangle inequality to Step 3, then use Steps 1 and 4:

The first inequality uses induced norms and the integral triangle inequality, not a scalar Cauchy-Schwarz argument. The factor $1/2$ comes from integrating the segment distance. Set $K=CL/2$; using the looser constant $CL$ is also valid.

6. Close the local induction and prove convergence. Keep $B,C,L$ fixed. Choose $0<\delta<r$ such that $q=K\delta<1$, and initialize $\|\mathbf e_0\|<\delta$. If $\|\mathbf e_k\|<\delta$, Step 5 gives

Induction therefore keeps every iterate in $B(\mathbf x^*,\delta)\subset B$, so all the solves and segment estimates used above remain valid. Repeated contraction gives $\|\mathbf e_k\|\le q^k\|\mathbf e_0\|\to0$. The stronger inequality $\|\mathbf e_{k+1}\|\le K\|\mathbf e_k\|^2$ proves at least Q-quadratic convergence. If an iterate equals $\mathbf x^*$, Newton has already reached the solution and the same bound holds with zero error.

The gradient-residual proof uses the same FTC, cancellation, and Lipschitz-remainder mechanics along a different segment.

Sources: Biegler (2010), Theorem 2.14, pp. 27-28, and Theorem 2.20, p. 35, with the corrected segment in the errata. The segment here runs from $\mathbf x^*$ to $\mathbf x^k$, the reverse parametrization of the corrected book segment.

Gradient-residual proof¶

Under the same hypotheses and sufficiently close initialization as the local convergence theorem, the gradient residual obeys

where $\mathbf g^k=\nabla f(\mathbf x^k)$, $\mathbf G^k=\nabla^2f(\mathbf x^k)$, $L$ is the Hessian Lipschitz constant, and $C$ bounds the inverse Hessian on the fixed local ball. This is the second result stated in the handout. The full derivation follows the same mechanics as the iterate-error proof: choose a segment, use the chain rule and FTC, cancel the linear part, and bound the remaining Hessian difference.

1. Check that the step segment stays in the neighborhood. The iterate-error induction keeps both $\mathbf x^k$ and $\mathbf x^{k+1}$ in $B(\mathbf x^*,\delta)\subset B$. Since a ball is convex,

The Hessian Lipschitz bound therefore applies on the entire step, and $\|(\mathbf G^k)^{-1}\|\le C$ at its starting point.

2. Apply FTC and cancel the Newton model. With $\boldsymbol\phi(t)=\nabla f(\mathbf y(t))$, the chain rule gives $\boldsymbol\phi'(t)=\nabla^2f(\mathbf y(t))\mathbf p^k$. Componentwise FTC gives

The Newton equation is $\mathbf g^k+\mathbf G^k\mathbf p^k=\mathbf0$. Subtract its zero linear-model residual and integrate the constant matrix term over $[0,1]$:

This identity isolates the exact gradient-model remainder; there is no unestimated Taylor remainder.

3. Bound the Hessian difference by the step distance. The intermediate point is $t\|\mathbf p^k\|$ from $\mathbf x^k$, so

Applying the integral triangle inequality and the induced-norm bound to Step 2 yields

Here the integral of $t$ supplies $1/2$, whereas the iterate-error proof integrates $1-t$.

4. Express the step bound in terms of the current gradient. The Newton solve and the uniform inverse bound give

Substitution into Step 3 proves

The iterate-error result already establishes $\mathbf x^k\to\mathbf x^*$; continuity of the gradient and stationarity imply $\mathbf g^k\to\mathbf0$. Together with this squared-residual bound, that gives at least Q-quadratic convergence of the gradient norm. Small steps alone do not imply the rate: the uniform inverse bound and Hessian Lipschitz assumption do the work. The local positive-definiteness hypothesis is sufficient here; for the analogous root-finding theorem a nonsingular limiting Jacobian suffices and does not by itself classify a stationary point as a minimum.

Source: Biegler (2010), Theorem 2.20 and Eq. (2.43), p. 35. The handout states this result and links here for the full second derivation.

Test Problem¶

Example 2.19 Consider the following two-variable unconstrained optimization problem:

with

The book reports the rounded candidate $\mathbf{x}^* = [0.7395, 0.3144]$ with $f(\mathbf{x}^*) = -5.0893$.

The nearby exact stationary point has zero gradient; rounding its coordinates does not preserve exact stationarity. The book reports the following approximate Hessian there:

which has eigenvalues $\lambda = 43.417$ and $\lambda = 426.362$.

Source: Biegler (2010), Example 2.19, p. 30. The canonical optimality notebook compares rounded and refined points and reports actual derivative residuals; the Newton experiments below use the same objective for a distinct iteration task.

Checking Hessian eigenvalues in Python¶

Use the reported symmetric Hessian from Example 2.19. scipy.linalg.eigvalsh returns its real eigenvalues. Both are positive, so this matrix is positive definite. The local-minimum conclusion also requires a zero gradient at the point under consideration. The displayed matrix and coordinates are rounded data, not an exact stationarity certificate.

Two starts: a minimum and a saddle¶

Biegler’s Example 2.21 applies unit-step Newton to the Example 2.19 objective defined above. Starting from $[0.8,0.3]^T$ gives the minimum near $[0.7395,0.31435]^T$; starting from $[1.0,0.5]^T$ gives a saddle near $[1.0768,0.9504]^T$. Both endpoint gradient norms are below 10^-12. The saddle Hessian has one positive and one negative eigenvalue, approximately 11.752 and -3.034.

The prepared figure uses the analytic gradient and a complex-step Hessian. The complete plotting implementation contains the objective, both derivatives, the Newton loop, and the checks against Biegler’s Tables 2.2 and 2.3. The finite histories illustrate fast convergence; they do not prove an asymptotic rate. Newton applied to $\nabla f=\mathbf0$ can converge quadratically to a saddle with a nonsingular Hessian when the Hessian is locally Lipschitz and the initial point is sufficiently close (Nocedal and Wright, Theorem 11.2, p. 276).

The following experiments use the same objective with finite-difference derivatives and additional starting points. They retain their deliberate failure cases.

Starting Point Near Optimal Solution¶

Activity: A Different Starting Point¶

Try with $\mathbf{x}_0 = [0,0]^T$

Activity: Let’s break it¶

Activity: Try $\mathbf{x}_0 = [-0.2, -0.2]^T$. Why does the gradient or Hessian return NaN? Hint: redefine calc_grad and create my_f_verbose using verbose=True.

Activity: Use $50\mathbf{I}$ in place of Hessian¶

Use the displayed constant model $50\mathbf I$ in place of the Hessian. This computes a gradient step with scale $1/50$; it does not compute the objective’s actual curvature.

Test with starting point $\mathbf{x}=[0.7,0.3]^T$

Test with starting point $\mathbf{x}=[0.0,0.0]^T$

Adjust Hessian with Levenberg-Marquardt Correction¶

Enhancement Ideas¶

Here are some ideas to improve this notebook:

• Plot the test function

• On the plot of the test function, visualize the solution path