Kazakhstan FAIO 2024 · Correction by IOAI Training

F. SVM

This page gives the full approach. Read it once you have searched.

Metric: accuracy (score 100*acc if acc >= 0.97)
Evaluated on: self-generated tests (official tests are hidden)
Our score: 100 % accuracy on 100/100 self-tests
Expected level: accuracy >= 0.97 per test, max 100 points

Expected level reached

We solved the problem ourselves: downloaded the data, ran the experiments and measured the score below. The statement and data belong to their organisers.

The problem in brief

You get between 1000 and 4000 points in the plane, with coordinates up to $10^5$ in absolute value, each labelled $+1$ or $-1$ , and a promise that some straight line separates the two classes. You must print a line $ax + by + c = 0$ . A point is classified $+1$ when $ax+by+c \ge 0$ and $-1$ otherwise. The statement asks for the line of a support vector machine: among all the lines that separate the classes, the one that stays as far as possible from the nearest points on each side. It suggests finding it by gradient descent on a loss it defines.

How you are scored

On each hidden test the judge measures the accuracy of your line on the points, and the score is

\text{score} = \begin{cases} 100 \cdot \text{accuracy} & \text{if accuracy} \ge 0.97,\\ 0 & \text{otherwise,}\end{cases}

averaged over the tests. Since the data is separable, 100 is reachable on every test. The official tests are hidden, so we use the sample and 100 random separable tests of our own, generated in the allowed ranges with class gaps from 10 to 5000.

A warning about the sample. The sample answer 1 1 0 classifies only 1 of the 10 sample points correctly (for example $(2, 6)$ gives $2+6 = 8 \ge 0$ , so $+1$ , but its label is $-1$ ). The opposite line -1 -1 0 gets 9 of 10, and -1 -1 -3 gets all 10: the printed answer has the wrong sign (and the wrong offset). Do not use it to check your sign convention. The real max-margin line for the sample is about $-3.42x - 7.70y - 12.41 = 0$ , which classifies all 10 points correctly.

The key idea

The margin. Write $\vec n = (a, b)$ . The distance from a point to the line is $\lvert a x + b y + c\rvert / \lvert\vec n\rvert$ , so the geometric margin of a separating line is

\gamma = \min_i \frac{l_i\,(a x_i + b y_i + c)}{\lvert \vec n \rvert} .

The SVM line maximises $\gamma$ . Rescaling $(a,b,c)$ does not change the line, so we can fix the numerator: require $l_i(a x_i + b y_i + c) \ge 1$ for every point, so that $\gamma = 1/\lvert\vec n\rvert$ . Maximising the margin then becomes

\min_{a,b,c}\ \lvert\vec n\rvert^2 \quad\text{subject to}\quad l_i(a x_i + b y_i + c) \ge 1 \ \text{ for all } i .

Why the statement's loss cannot work as written. It only charges the misclassified points, $f = \lvert\vec n\rvert\, d = \lvert a x + b y + c\rvert$ , plus $\lambda\lvert\vec n\rvert^2$ . Take $a=b=c=0$ : every point is classified $+1$ , the misclassified ones cost $\lvert 0 \rvert = 0$ , and the regulariser is 0 too. The loss is exactly 0, its global minimum, and the "classifier" says $+1$ to everyone. The loss must demand a margin, not just the right side. The standard fix is the hinge loss, the soft version of the constraint above:

\text{Loss}(a,b,c) = \lambda \lvert\vec n\rvert^2 + \frac{1}{k}\sum_{i=1}^{k} \max\big(0,\ 1 - l_i(a x_i + b y_i + c)\big).

A point now costs nothing only once it lies beyond the margin.

The geometric shortcut. In two dimensions the hard-margin problem has an exact geometric answer: the SVM line is the perpendicular bisector of the closest pair of points between the convex hulls of the two classes, and the margin is half the distance between the hulls. Intuitively, any separating strip must fit between the two hulls, and the widest one is set by the two hulls' nearest points. Formally, this is the dual of the problem above: the optimal $\vec n$ is proportional to $p^{\star} - q^{\star}$ , where $p^{\star}$ and $q^{\star}$ are the closest points of the two hulls.

A first baseline: gradient descent

Following the statement, we minimise the hinge loss by full-batch (sub)gradient descent:

\nabla_{\vec n} = 2\lambda \vec n - \frac{1}{k}\sum_{i:\ l_i(\vec n\cdot x_i + c) < 1} l_i\, x_i, \qquad \partial_c = -\frac{1}{k}\sum_{i:\ l_i(\vec n\cdot x_i + c) < 1} l_i .

Scale the data first. With coordinates near $10^5$ , the gradient in $\vec n$ is $10^5$ times larger than the gradient in $c$ , and no single learning rate suits both. We centre the points and divide both axes by the same factor (the largest absolute coordinate). A common factor keeps angles and distance ratios, so the max-margin line is unchanged, and we map the result back at the end: $a = w_1/s$ , $b = w_2/s$ , $c = c' - \vec w\cdot\vec\mu/s$ .

With $\lambda = 10^{-4}$ , 3000 steps and a step size of $1/\sqrt t$ , this gets a score of 99.73 on our 100 tests. Every test clears the 97 % bar, but only 39 are at 100 %, and the worst one stops at 98.9 %. The failures are the tests where the gap between the classes is tiny compared with the spread of the points (10 units out of 200 000). To push every point past the margin, $\lvert\vec n\rvert$ must grow to about $1/\gamma$ in scaled units, which takes a very long time with a decaying step and fights against the regulariser.

The full solution

Convex hulls of the $+1$ points and of the $-1$ points (Andrew's monotone chain, $O(k\log k)$ ).
Closest pair $p \in H_+$ , $q \in H_-$ : the minimum over all vertex-to-edge distances, in both directions (vectorised with numpy, at most a few million distance computations).
The line: normal $\vec n = p - q$ (pointing towards the $+1$ class) through the midpoint, $c = -\vec n\cdot\frac{p+q}{2}$ .

p, q = closest_points(convex_hull(X[l == 1]), convex_hull(X[l == -1]))
a, b = p - q
c = -(p - q) @ (p + q) / 2

The two classes, the exact max-margin line (solid) with its margin, and the gradient-descent line (dotted)

Both lines separate this example. The exact one has a margin of 8164 (scikit-learn's SVC solver finds 8163), while gradient descent stopped at 6183: close, but not centred in the gap.

Results

100 random separable tests in the allowed ranges, compared with scikit-learn's SVC(kernel="linear", C=1e8) as the reference max-margin solver.

method	min accuracy	tests at 100 %	score	margin / optimal (median)
statement's loss, at its minimum $a=b=c=0$	about 0.5	0	0
gradient descent on the hinge loss (scaled data)	0.989	39	99.73	negative (some points stay inside the margin)
exact: closest points of the convex hulls	1.000	100	100.00	1.000

What to remember

A loss must reward what you actually want: if correct points cost nothing and the regulariser shrinks $\vec n$ , the zero vector wins. The hinge loss $\max(0, 1 - l\,f(x))$ builds the margin into the objective.
Max margin = minimum $\lvert\vec n\rvert$ under $l_i(\vec n\cdot x_i + c) \ge 1$ : fixing the scale of the line turns a max-min problem into a clean quadratic program.
Rescale before gradient descent, with the same factor on every axis when geometry matters, and map the parameters back.
Geometry gives the exact answer in 2-D: the SVM separates the convex hulls through the bisector of their closest points. In higher dimensions the same duality is what SVM solvers exploit.
Do not trust a sample output blindly: check it against the statement's own conventions.

Full code (solution.py)

"""FAIO 2024, problem F: SVM -- a linear max-margin classifier for separable 2-D points, by gradient descent.

stdin : k, then k lines "x y l" with l in {-1, +1} (blank lines allowed)
stdout: a b c  such that sign(a x + b y + c) (>= 0 -> +1) classifies the points

Score per test: 100 * accuracy if accuracy >= 0.97, else 0.

Two fixes to the statement's loss, explained in the correction:
  * the statement penalises only misclassified points, so a = b = c = 0 already has loss 0:
    we use the hinge loss max(0, 1 - l (w.x + c)) which demands a margin;
  * coordinates go up to 1e5: we centre and rescale them (same factor on both axes, which keeps
    angles and therefore the max-margin line) before gradient descent, then map back.

Run:  python solution.py < input.txt
      python solution.py --check    (random separable tests, compared with the exact SVM)
      python solution.py --figure
"""
import sys
import numpy as np


def train_svm(X, l, lam=1e-4, iters=3000, lr0=1.0):
    """Full-batch subgradient descent on  lam |w|^2 + mean(max(0, 1 - l (X w + c)))."""
    mu = X.mean(0)
    s = np.abs(X - mu).max() or 1.0                      # one scale for both axes
    Z = (X - mu) / s
    w = np.zeros(2); c = 0.0
    best = (-1.0, -np.inf, w.copy(), c)
    for t in range(1, iters + 1):
        m = l * (Z @ w + c)                              # functional margins
        act = m < 1                                      # points inside the margin or wrong
        gw = 2 * lam * w - (l[act, None] * Z[act]).sum(0) / len(l)
        gc = -l[act].sum() / len(l)
        lr = lr0 / np.sqrt(t)
        w -= lr * gw; c -= lr * gc
        if t % 10 == 0 or t == iters:
            pred = np.where(Z @ w + c >= 0, 1, -1)
            acc = (pred == l).mean()
            nw = np.linalg.norm(w)
            geo = (l * (Z @ w + c)).min() / nw if nw > 0 else -np.inf
            if (acc, geo) > best[:2]:
                best = (acc, geo, w.copy(), c)
    _, _, w, c = best
    # back to original coordinates: w.(x - mu)/s + c = (w/s).x + (c - w.mu/s)
    a, b = w / s
    return a, b, c - (w @ mu) / s


# ---------------------------------------------------------------------------------------------
# Exact hard-margin SVM by geometry: the max-margin line is the perpendicular bisector of the
# closest pair of points between the two classes' convex hulls.
# ---------------------------------------------------------------------------------------------
def cross(u, v):
    return u[0] * v[1] - u[1] * v[0]


def convex_hull(P):
    """Andrew's monotone chain; returns the hull vertices (counter-clockwise)."""
    P = np.unique(P, axis=0)
    if len(P) <= 2:
        return P
    def half(pts):
        h = []
        for p in pts:
            while len(h) >= 2 and cross(h[-1] - h[-2], p - h[-2]) <= 0:
                h.pop()
            h.append(p)
        return h[:-1]
    return np.array(half(P) + half(P[::-1]))


def closest_points(A, B):
    """Closest pair between polygons A and B (vertex-to-edge distances, both directions)."""
    best = (np.inf, None, None)
    for P, Q, swap in ((A, B, False), (B, A, True)):
        Q0, Q1 = Q, np.roll(Q, -1, axis=0)                     # edges of Q
        E = Q1 - Q0
        L2 = np.maximum((E ** 2).sum(1), 1e-300)
        for i in range(0, len(P), 256):
            p = P[i:i + 256, None, :]
            t = np.clip(((p - Q0) * E).sum(2) / L2, 0, 1)
            proj = Q0 + t[..., None] * E
            d = ((p - proj) ** 2).sum(2)
            j = np.unravel_index(d.argmin(), d.shape)
            if d[j] < best[0]:
                a, b = P[i + j[0]], proj[j]
                best = (d[j], b, a) if swap else (d[j], a, b)
    return best[1], best[2]


def exact_svm(X, l):
    p, q = closest_points(convex_hull(X[l == 1]), convex_hull(X[l == -1]))
    w = p - q                                   # normal points towards the +1 class
    c = -w @ (p + q) / 2                        # the line passes through the midpoint
    return w[0], w[1], c


def parse(text):
    t = text.split()
    k = int(t[0])
    v = np.array(t[1:1 + 3 * k], dtype=float).reshape(k, 3)
    return v[:, :2], v[:, 2].astype(int)


def accuracy(X, l, a, b, c):
    return float((np.where(X @ np.array([a, b]) + c >= 0, 1, -1) == l).mean())


def margin(X, l, a, b, c):
    return float((l * (X @ np.array([a, b]) + c)).min() / np.hypot(a, b))


def solve():
    X, l = parse(sys.stdin.read())
    a, b, c = exact_svm(X, l)
    print(f"{a:.12g} {b:.12g} {c:.12g}")


def random_case(rng, k=None, gap=None):
    k = k or rng.integers(1001, 4000)
    theta = rng.uniform(0, 2 * np.pi); n = np.array([np.cos(theta), np.sin(theta)])
    off = rng.uniform(-3e4, 3e4)
    X = rng.uniform(-99999, 99999, size=(4 * k, 2))
    d = X @ n + off
    g = gap if gap is not None else rng.choice([10.0, 100.0, 1000.0, 5000.0])
    X = X[np.abs(d) > g][:k]
    l = np.where(X @ n + off >= 0, 1, -1)
    return X, l


def check():
    from sklearn.svm import SVC
    X, l = parse("10 2 6 -1 4 -8 1 9 -1 -1 5 2 -1 -9 7 -1 10 4 -1 2 -8 1 -6 -4 1 -10 -4 1 -2 8 -1")
    print("sample: statement's '1 1 0' accuracy", accuracy(X, l, 1, 1, 0),
          "| exact", accuracy(X, l, *exact_svm(X, l)), [round(float(v), 4) for v in exact_svm(X, l)],
          "| GD", accuracy(X, l, *train_svm(X, l)))
    rng = np.random.default_rng(0)
    res = {"gradient descent (hinge)": ([], []), "exact (convex hulls)": ([], [])}
    for _ in range(100):
        X, l = random_case(rng)
        sv = SVC(kernel="linear", C=1e8).fit(X / 1e5, l)             # reference hard-margin SVM
        A, B = sv.coef_[0] / 1e5; C = sv.intercept_[0]
        opt = margin(X, l, A, B, C)
        for name, fn in (("gradient descent (hinge)", train_svm), ("exact (convex hulls)", exact_svm)):
            a, b, c = fn(X, l)
            res[name][0].append(accuracy(X, l, a, b, c)); res[name][1].append(margin(X, l, a, b, c) / opt)
    for name, (acc, rat) in res.items():
        acc = np.array(acc)
        print(f"{name:26s} accuracy min {acc.min():.4f} | 100 % on {(acc == 1).sum()}/100 | "
              f"score {np.mean(np.where(acc >= 0.97, 100 * acc, 0)):.2f} | margin/optimal median {np.median(rat):.4f}"
              f" min {np.min(rat):.4f}")


def figure():
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt
    from sklearn.svm import SVC
    rng = np.random.default_rng(3)
    X, l = random_case(rng, k=1500, gap=8000.0)
    a, b, c = exact_svm(X, l)
    g = train_svm(X, l)
    sv = SVC(kernel="linear", C=1e8).fit(X / 1e5, l)
    A, B = sv.coef_[0] / 1e5; C = sv.intercept_[0]
    fig, ax = plt.subplots(figsize=(7.5, 7))
    ax.scatter(*X[l == 1].T, s=6, color="#2171b5", label="l = +1")
    ax.scatter(*X[l == -1].T, s=6, color="#d95f02", label="l = -1")
    xs = np.array([-1e5, 1e5])
    for (p, q, r), sty, nm in [(g, ":", "gradient descent on the hinge loss"),
                               ((a, b, c), "-", "exact max-margin line (convex hulls)")]:
        ax.plot(xs, -(p * xs + r) / q, "k" + sty, lw=2 if sty == "-" else 2.5, label=nm)
    mg = margin(X, l, a, b, c); nn = np.hypot(a, b)
    for sgn in (-1, 1):
        ax.plot(xs, -(a * xs + c + sgn * mg * nn) / b, color="gray", ls="--", lw=1)
    ax.set_xlim(-1e5, 1e5); ax.set_ylim(-1e5, 1e5); ax.set_aspect("equal")
    ax.set_title(f"Margin: exact {mg:.0f}, gradient descent {margin(X, l, *g):.0f}, sklearn SVC {margin(X, l, A, B, C):.0f}", fontsize=10)
    ax.legend(loc="lower left", fontsize=9)
    fig.tight_layout(); fig.savefig("figures/margin.png", dpi=150, facecolor="white")


if __name__ == "__main__":
    if "--check" in sys.argv:
        check()
    elif "--figure" in sys.argv:
        figure()
    else:
        solve()