Specialized Optimization Methods
Topics in Digital Heritage:
Numerical Methods for Digital Reconstruction
Fall 2026
Framing as a Linear Least Squares Problem
Squared Tangents Make a Parabola
- A sum of squared lines is always a quadratic in \(\delta\mathbf{x}\), so its minimizer has a closed form (will see next slide).
- We repeat this linearization at the new point \(\mathbf{x}_1\), and so on, until convergence.
Demo: Gauss-Newton Converges Fast
Back to the Nonlinear Regression Example
Show code
import numpy as np
import matplotlib.pyplot as plt
rng = np.random.default_rng(0)
a_true, c_true = 0.6, 2.0
x_data = np.linspace(0, 4, 25)
y_data = c_true * np.exp(a_true * x_data) + rng.normal(scale=0.4, size=x_data.size)
def residual(theta):
a, c = theta
return y_data - c * np.exp(a * x_data)
def jacobian(theta):
a, c = theta
e = np.exp(a * x_data)
J = np.zeros((x_data.size, 2))
J[:, 0] = -c * x_data * e
J[:, 1] = -e
return J
def loss(theta):
r = residual(theta)
return 0.5 * r @ r
def gradient_descent(theta0, lr, max_iter):
theta = np.array(theta0, dtype=float)
hist = [theta.copy()]
for _ in range(max_iter):
theta = theta - lr * (jacobian(theta).T @ residual(theta))
hist.append(theta.copy())
return np.array(hist)
def gauss_newton(theta0, max_iter):
theta = np.array(theta0, dtype=float)
hist = [theta.copy()]
for _ in range(max_iter):
J, r = jacobian(theta), residual(theta)
delta, *_ = np.linalg.lstsq(J, -r, rcond=None)
theta = theta + delta
hist.append(theta.copy())
return np.array(hist)
theta0_good = [0.3, 1.0]
h_gd = gradient_descent(theta0_good, lr=5e-5, max_iter=20000)
h_gn = gauss_newton(theta0_good, max_iter=8)
fig, axes = plt.subplots(1, 2, figsize=(8.5, 4))
xs = np.linspace(0, 4, 200)
axes[0].plot(x_data, y_data, "o", ms=4, color="gray", label="data")
axes[0].plot(xs, h_gn[-1][1] * np.exp(h_gn[-1][0] * xs), color="steelblue", lw=2, label="Gauss-Newton fit")
axes[0].set_title(f"Fitted $y=ce^{{ax}}$\n($a$={h_gn[-1][0]:.3f}, $c$={h_gn[-1][1]:.3f}, true: 0.6, 2.0)", fontsize=9)
axes[0].legend(fontsize=8)
x_gd = np.arange(1, len(h_gd) + 1)
x_gn = np.arange(1, len(h_gn) + 1)
axes[1].loglog(x_gd, [loss(t) for t in h_gd], color="firebrick", label=f"Gradient descent ({len(h_gd)-1} steps)")
axes[1].loglog(x_gn, [loss(t) for t in h_gn], "-o", ms=4, color="steelblue", label=f"Gauss-Newton ({len(h_gn)-1} steps)")
axes[1].set_xlabel("iteration $t$ (log scale)"); axes[1].set_ylabel("$E(\\theta_t)$")
axes[1].legend(fontsize=8)
axes[1].set_title("Same accuracy, 2500x fewer steps", fontsize=9)
plt.tight_layout()
plt.show()
When Gauss-Newton Misbehaves
Levenberg-Marquardt
Gauss-Newton is not guaranteed to converge:
- far from the solution, the linear approximation of \(F\) can be bad enough that the step makes things worse.
Levenberg-Marquardt
Damp with \(\lambda>0\): \(\mathbf{x} = \mathbf{x}_0 - \big(DF(\mathbf{x}_0)^\top DF(\mathbf{x}_0) + \lambda \mathbf{I}\big)^{-1} DF(\mathbf{x}_0)^\top F(\mathbf{x}_0)\)
- Small \(\lambda\) \(\to\) behaves like Gauss-Newton;
- Large \(\lambda\) \(\to\) behaves like a small gradient-descent step.
- This is exactly the Tikhonov regularization trick applied to the optimization step instead of the least-squares problem.
Demo: Levenberg-Marquardt Rescues a Bad Start
Show code
def levenberg_marquardt(theta0, lam0, max_iter):
theta = np.array(theta0, dtype=float)
lam = lam0
cur = loss(theta)
hist = [theta.copy()]
for _ in range(max_iter):
J, r = jacobian(theta), residual(theta)
delta = np.linalg.solve(J.T @ J + lam * np.eye(2), -J.T @ r)
trial = theta + delta
trial_loss = loss(trial)
if trial_loss < cur:
theta, cur, lam = trial, trial_loss, lam * 0.5
else:
lam *= 2.0
hist.append(theta.copy())
return np.array(hist)
theta0_bad = [0.05, 0.3]
h_gn_bad = gauss_newton(theta0_bad, max_iter=6)
h_lm_bad = levenberg_marquardt(theta0_bad, lam0=1.0, max_iter=30)
losses_gn_bad = np.clip([loss(t) for t in h_gn_bad], None, 1e6)
losses_lm_bad = [loss(t) for t in h_lm_bad]
fig, ax = plt.subplots(figsize=(8.5, 4))
ax.semilogy(range(len(losses_gn_bad)), losses_gn_bad, "-o", color="firebrick", label="Gauss-Newton (undamped)")
ax.semilogy(range(len(losses_lm_bad)), losses_lm_bad, "-o", ms=3, color="steelblue", label="Levenberg-Marquardt")
ax.axhline(loss([a_true, c_true]), color="forestgreen", ls="--", lw=1, label="noise floor")
ax.set_xlabel("iteration $t$"); ax.set_ylabel("$E(\\theta_t)$ (capped at $10^6$ for display)")
ax.set_title(f"Same problem, worse start $(a_0,c_0)=({theta0_bad[0]},{theta0_bad[1]})$", fontsize=10)
ax.legend(fontsize=9)
plt.tight_layout()
plt.show()
Undamped Gauss-Newton’s first step alone overshoots to a loss above \(10^{40}\) (clipped above for display) and then gets stuck; Levenberg-Marquardt spends a few steps growing \(\lambda\) to stay safe, then shrinks it again and converges cleanly.
Example: k-means Clustering as Alternation
k-means objective
Given data \(\vec{x}_1,\dots,\vec{x}_m\), find \(k\) centers \(\vec{y}_1,\dots,\vec{y}_k\) minimizing \[ E(\vec{y}_1,\dots,\vec{y}_k) \equiv \sum_{i=1}^m \min_{c\in\{1,\dots,k\}} \|\vec{x}_i-\vec{y}_c\|_2^2. \]
Introducing the assignment \(c_i\) explicitly splits the variables into two blocks:
\[\text{labels}\quad \{c_i\}\qquad\qquad\text{centers}\quad\{\vec{y}_j\}.\]
- Fix centers \(\{\vec{y}_j\}\), optimize labels \(\{c_i\}\): assign each point to its nearest center (closed form).
- Fix labels \(\{c_i\}\), optimize centers \(\{\vec{y}_j\}\): \(\vec{y}_j \leftarrow\) mean of the points assigned to cluster \(j\) (closed form).
Example: k-means Clustering as Alternation
Example: Shape Deformation (ARAP)
As-rigid-as-possible (ARAP) deformation
For a mesh with vertices \(v\) and edges \((v,w)\), find rotations \(R_v\) and new positions \(\vec{y}_v\) minimizing \[ \sum_{v}\sum_{(v,w)} \|R_v(\vec{x}_v-\vec{x}_w) - (\vec{y}_v-\vec{y}_w)\|_2^2 \qquad \text{subject to} \quad R_v\in SO(3). \]
Fix \(R_v\), optimize \(\vec{y}_v\):
A sparse SPD linear system, exactly what conjugate gradient was built for.
Fix \(\vec{y}_v\), optimize \(R_v\):
Decouples per vertex into a \(2\times2\) closed-form SVD (a rotation-fitting/Procrustes problem).
A hard, nonconvex joint problem, split by alternation into two problems we already know how to solve.
Escaping Bad Local Minima
Recall the image-alignment demo:
- Newton-like local search methods can get stuck in bad local minima due to high-frequency texture in the original images.
- Blurring the photos first, then sharpening, can avoid bad local alignments caused by fine texture.
Graduated optimization generalizes this:
- Solve a sequence \(f_1,f_2,\dots,f_k=f\) of progressively harder problems,
- using each minimizer as the next initial guess.
