6 · Why Line Search Matters#

A Newton direction \(-H^{-1}\nabla E\) points the right way, but the step length still matters. Too large a step overshoots, inverts elements, and the energy explodes instead of decreasing.

Backtracking line search starts from a full step and shrinks it until the energy actually decreases (the Armijo condition). To make the failure easy to watch, we move the handle slowly right and back and take a few Newton iterations per frame — once with a fixed, too-large step, once with line search.

%matplotlib inline
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import simkit
from simkit import backtracking_line_search
import utils

A thick, stiff, nearly-incompressible beam#

Same objective as tutorial 5 but thicker (so bending stores more energy) and nearly incompressible (\(\nu = 0.49\) — resisting volume change makes the energy much stiffer), so an over-eager step bites fast and the value of line search is obvious. The handle (the full right edge) follows a there-and-back ramp.

BeamObjective composes the now-familiar terms — the Neo-Hookean material, a fixed pin spring on the left edge, and a movable handle spring on the right edge — and sums them one per line in energy / gradient / hessian.

class BeamObjective:
    """Minimize  elastic(x) + pin spring + handle spring  for a thick beam."""

    def __init__(self, X, T, K=1e6):
        mu, lam = simkit.ympr_to_lame(50.0, 0.49)            # high Poisson ratio
        self.p   = utils.precompute(X, T, mu=mu, lam=lam)    # J, vol, masses, ...
        self.psi = utils.make_material("Neo-Hookean")        # elastic term
        self.pin    = utils.PenaltySpring(self.p.n, self.p.dim, K)
        self.handle = utils.PenaltySpring(self.p.n, self.p.dim, K)

    def set_pin(self, idx, targets):    self.pin.set(idx, targets);    return self
    def set_handle(self, idx, targets): self.handle.set(idx, targets); return self

    def energy(self, x):
        E_elastic = self.psi.energy(x, self.p)
        E_pin     = self.pin.energy(x)
        E_handle  = self.handle.energy(x)
        return E_elastic + E_pin + E_handle

    def gradient(self, x):
        g_elastic = self.psi.gradient(x, self.p)
        g_pin     = self.pin.gradient(x)
        g_handle  = self.handle.gradient(x)
        return g_elastic + g_pin + g_handle

    def hessian(self, x):
        H_elastic = self.psi.hessian(x, self.p)
        H_pin     = self.pin.hessian(x)
        H_handle  = self.handle.hessian(x)
        return H_elastic + H_pin + H_handle

X, T = utils.triangulated_grid(nx=16, ny=7, width=2.0, height=0.6)   # thick beam
n, dim = X.shape

pin_idx   = np.where(X[:, 0] <= X[:, 0].min() + 1e-6)[0]
right_idx = np.where(X[:, 0] >= X[:, 0].max() - 1e-6)[0]

ramp = np.concatenate([np.linspace(0, 1, 22), np.linspace(1, 0, 22)])    # right then back
offsets = [np.array([1.0 * r, 0.0]) for r in ramp]

The quasi-static driver#

At each handle position we take a few Newton iterations. The only difference between the two runs is the step: a fixed STEP (deliberately too big) versus the backtracking step that guarantees the energy decreases.

We solve the linear system \(H\,\Delta x = -g\) for the Newton direction, then either take a fixed multiple of it or let backtracking_line_search pick a safe step length \(\alpha\).

DIVERGE = 1e3

def drive(use_line_search, STEP=2.0, iters_per_frame=2):
    beam = BeamObjective(X, T).set_pin(pin_idx, X[pin_idx])
    x = X.flatten().reshape(-1, 1).copy()
    states, blew_at = [], None
    for k, off in enumerate(offsets):
        beam.set_handle(right_idx, X[right_idx] + off)
        for _ in range(iters_per_frame):
            g = beam.gradient(x)
            H = beam.hessian(x)
            dx = sp.sparse.linalg.spsolve(H.tocsc(), -g).reshape(-1, 1)
            if use_line_search:
                alpha = backtracking_line_search(beam.energy, x, g, dx)[0]
            else:
                alpha = STEP
            x = x + alpha * dx
        states.append(x.reshape(n, dim).copy())
        if not np.isfinite(np.abs(x).max()) or np.abs(x).max() > DIVERGE:
            blew_at = k
            break
    return states, blew_at

fixed_states, blew = drive(use_line_search=False, STEP=2.0)
ls_states,    _    = drive(use_line_search=True)
print("fixed step (no line search): blew up at frame", blew, "of", len(offsets))
print("line search: completed all", len(ls_states), "frames")

fixed step (no line search): blew up at frame 9 of 44
line search: completed all 44 frames

Without line search: the solver explodes (fixed step = 2x the Newton step)#

lims = ((-1.3, 2.2), (-1.6, 1.6))
fig, anim = utils.animate_mesh(fixed_states, T, lims=lims,
    title="Fixed step, NO line search  ->  explodes", pin_pts=X[pin_idx],
    handle_traj=[s[right_idx] for s in fixed_states], fps=12)
utils.show_video(fig, anim, "media/06_no_line_search.mp4", fps=12)

With line search: smooth and stable#

fig, anim = utils.animate_mesh(ls_states, T, lims=lims,
    title="Backtracking line search  ->  stable", pin_pts=X[pin_idx],
    handle_traj=[s[right_idx] for s in ls_states], fps=20)
utils.show_video(fig, anim, "media/06_line_search.mp4", fps=20)

Takeaways#

A good Newton direction isn’t enough — an unchecked step length can diverge.
Backtracking line search shrinks the step until the energy decreases.
Cheap insurance: a few extra energy evaluations buy a solver that won’t blow up.