Practice set: constrained optimisation

import numpy as np

BOWL = [[2.0, 0.0], [0.0, 4.0]]  # ½xᵀAx = x² + 2y²
TILTED = [[2.0, 1.0], [1.0, 2.0]]  # with b = (3, 0): x² + xy + y² − 3x


def _raises_value_error(*args):
    try:
        minimise_quadratic(*args)
    except ValueError:
        return True
    return False


def _solve(*args):
    out = minimise_quadratic(*args)
    assert out is not None, "minimise_quadratic returned None: return the pair (x, lam)"
    x, lam = out
    return np.asarray(x, dtype=float), float(lam)


def test_binding():
    """x² + 2y² with x + y ≥ 3: the constraint binds at (2, 1), with λ = 4"""
    x, lam = _solve(BOWL, [0, 0], [1, 1], 3)
    expect_close(x, [2.0, 1.0], what="x for x + y ≥ 3")
    expect_close(lam, 4.0, what="λ for x + y ≥ 3")


def test_slack():
    """x² + 2y² with x + y ≥ −1: the origin is allowed, so x = (0, 0) and λ = 0"""
    x, lam = _solve(BOWL, [0, 0], [1, 1], -1)
    expect_close(x, [0.0, 0.0], what="x for x + y ≥ −1", abs_tol=1e-12)
    expect_close(lam, 0.0, what="λ for a slack constraint", abs_tol=1e-12)


def test_on_the_edge():
    """x² + 2y² with x + y ≥ 0: the lowest point lies on the edge itself, and λ = 0"""
    x, lam = _solve(BOWL, [0, 0], [1, 1], 0)
    expect_close(x, [0.0, 0.0], what="x for x + y ≥ 0", abs_tol=1e-12)
    expect_close(lam, 0.0, what="λ for x + y ≥ 0", abs_tol=1e-12)


def test_three_variables():
    """x² + 2y² + 3z² with x + y + z ≥ 1 gives (6/11, 3/11, 2/11) and λ = 12/11"""
    x, lam = _solve(np.diag([2.0, 4.0, 6.0]), [0, 0, 0], [1, 1, 1], 1)
    expect_close(x, [6 / 11, 3 / 11, 2 / 11], what="x in three variables")
    expect_close(lam, 12 / 11, what="λ in three variables")


def test_upper_bound():
    """x² + xy + y² − 3x with x ≤ 1, written as −x ≥ −1, gives (1, −0.5) and λ = 1.5"""
    x, lam = _solve(TILTED, [3, 0], [-1, 0], -1)
    expect_close(x, [1.0, -0.5], what="x for x ≤ 1")
    expect_close(lam, 1.5, what="λ for −x ≥ −1")


def test_tilted_slack():
    """x² + xy + y² − 3x with x ≤ 3: its own minimum (2, −1) is allowed"""
    x, lam = _solve(TILTED, [3, 0], [-1, 0], -3)
    expect_close(x, [2.0, -1.0], what="x for x ≤ 3")
    expect_close(lam, 0.0, what="λ for a slack constraint", abs_tol=1e-12)


def test_kkt_conditions():
    """For a random positive definite A, the answer passes every KKT condition"""
    rng = np.random.default_rng(45)
    M = rng.normal(size=(4, 4))
    A = M.T @ M + np.eye(4)
    b, a = rng.normal(size=4), rng.normal(size=4)
    x0 = np.linalg.solve(A, b)
    c = float(a @ x0) + 2.0  # the unconstrained minimum is not allowed
    x, lam = _solve(A, b, a, c)
    expect_close(a @ x, c, what="a · x (the constraint binds)")
    assert lam >= 0, f"λ must be at least 0 for a ≥ constraint, got {lam}"
    expect_close(A @ x - b, lam * a, what="∇f = Ax − b against λa")


def test_rate():
    """λ is the rate of the best value: for x + y ≥ c at c = 3, a central difference gives 4"""
    def best(c):
        x, _ = _solve(BOWL, [0, 0], [1, 1], c)
        return 0.5 * x @ np.array(BOWL) @ x
    h = 1e-3
    expect_close((best(3 + h) - best(3 - h)) / (2 * h), 4.0, what="(f*(3 + h) − f*(3 − h)) / 2h", rel=1e-6)


def test_not_symmetric():
    """A matrix that is not symmetric raises ValueError"""
    assert _raises_value_error([[2, 1], [0, 2]], [0, 0], [1, 1], 1), "a matrix that is not symmetric should raise ValueError"


def test_not_positive_definite():
    """A matrix with a negative eigenvalue raises ValueError"""
    assert _raises_value_error([[2, 0], [0, -1]], [0, 0], [1, 1], 1), "diag(2, −1) is not positive definite: raise ValueError"


def test_semidefinite():
    """A semidefinite matrix, with a zero eigenvalue, raises ValueError"""
    assert _raises_value_error([[1, 1], [1, 1]], [0, 0], [1, 0], 1), "[[1, 1], [1, 1]] has eigenvalue 0: raise ValueError"


def test_zero_direction():
    """a = (0, 0) raises ValueError"""
    assert _raises_value_error(BOWL, [0, 0], [0, 0], 1), "the zero vector gives no constraint: raise ValueError"


def test_sizes():
    """Vectors of the wrong length raise ValueError"""
    assert _raises_value_error(BOWL, [0, 0, 0], [1, 1], 1), "b has 3 entries for a 2 × 2 matrix: raise ValueError"