Indefinite systems - Executable Engineering: Interactive Numerical Methods in Python

The quadratic $f(x, y) = x^2 - 6xy + y^2$ is an example of a problem with an indefinite Hessian matrix. Such functions are neither convex nor concave, and the stationary point where $\nabla f = \mathbf{0}$ is neither a minimum nor a maximum, but a saddle point.

Saddle point problems form an important subclass of numerical optimization—especially when solving constrained problems via Lagrange multipliers.

Newton’s method and gradient-based solvers readily identify points where $\nabla f = \mathbf{0}$ , but direct methods often struggle unless the objective function is reformulated (e.g., minimizing $\|\nabla f\|^2$ ).

import numpy as np
import plotly.graph_objects as go

Eigenvalues and matrix definiteness¶

The fact that $f(x, y) = x^2 - 6xy + y^2$ is indefinite may come as a surprise considering the slices along the axes, $f(x, 0) = x^2$ and $f(0, y) = y^2$ , are strictly positive parabolas curving upward.

The key is coordinate rotation. If we rotate coordinates by $45^\circ$ ( $x' = \frac{x - y}{\sqrt{2}}$ , $y' = \frac{x + y}{\sqrt{2}}$ ), the cross term vanishes:

f(x', y') = 8 (x')^2 - 4 (y')^2

(1)

Here the saddle nature is obvious: positive curvature (+8) along the anti-diagonal ( $y = -x$ ), and negative curvature (-4) along the diagonal ( $y = x$ ).

These principal axes are the eigenvectors of the Hessian $\mathbf{H}$ , and the curvatures along them are the eigenvalues $\lambda_i$ :

Definiteness	Eigenvalues $\lambda_i$	Geometric Shape
Positive Definite	All $\lambda_i > 0$	Local minimum (bowl curves upward)
Negative Definite	All $\lambda_i < 0$	Local maximum (dome curves downward)
Indefinite	Mixed signs (some $\lambda_i > 0$ , some $\lambda_j < 0$ )	Saddle point (curves up along some directions, down along others)
Positive Semi-Definite	All $\lambda_i \ge 0$	Valley / trough (flat along zero-eigenvalue directions)
Negative Semi-Definite	All $\lambda_i \le 0$	Ridge (flat along zero-eigenvalue directions)

def hessian_f(x):
    return np.array([[2.0, -6.0], [-6.0, 2.0]])

hessian = hessian_f([0.0, 0.0])
eigenvalues, eigenvectors = np.linalg.eigh(hessian)

print("Eigenvalues:", eigenvalues)
print("Eigenvectors:\n", eigenvectors)

x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)
Z = X**2 - 6 * X * Y + Y**2

fig = go.Figure(data=[go.Surface(x=X, y=Y, z=Z, colorscale='Viridis', opacity=0.85)])
fig.update_layout(
    title='Saddle Surface: f(x, y) = x^2 - 6xy + y^2',
    scene=dict(xaxis_title='x', yaxis_title='y', zaxis_title='f(x, y)'),
    width=700,
    height=500
)
fig.show()

Limitations of pure Newton optimization¶

Newton’s optimization routine is highly sensitive to initial guesses. When $\mathbf{H}$ is indefinite:

The standard step $\mathbf{H} \Delta \mathbf{x} = -\nabla f(\mathbf{x})$ can point directly toward a saddle point or even ascend toward a local maximum rather than a minimum.
Inversion of the Hessian is computationally demanding ( $O(N^3)$ ) and can be ill-conditioned or singular near inflection regions.

Rather than taking an unconstrained Newton step, we need methods that constrain step size to a region where our local quadratic model can be trusted—motivating Trust Region Methods.