1"""Exact weighted linear-probe coreset by preserving normal equations."""
 2import numpy as np
 3
 4
 5def numerical_rank(x, tol=1e-6):
 6    s = np.linalg.svd(x, compute_uv=False)
 7    return int(np.sum(s > (tol * s[0] if len(s) and s[0] else 0)))
 8
 9
10def fit_min_norm(x, y, weights=None):
11    if weights is None:
12        weights = np.ones(len(x))
13    a = x.T @ (weights[:, None] * x)
14    b = y.T @ (weights[:, None] * x)
15    return b @ np.linalg.pinv(a)
16
17
18def residual_coreset(x, y, tol=1e-10, max_iter=None):
19    """Return <= (m+1)r examples preserving the full optimum's normal equations.
20
21    For W*, each atom is residual_i outer x_i.  A nonnegative nullspace move
22    among k+1 atoms preserves sum w_i residual_i outer x_i, where k=m*d.
23    Local elimination avoids large full nullspace SVDs.
24    """
25    n, d = x.shape
26    m = y.shape[1]
27    weights = np.ones(n)
28    W = fit_min_norm(x, y)
29    r = numerical_rank(x)
30    R = ((y - x @ W.T)[:, :, None] * x[:, None, :]).reshape(n, m*d)
31    active = list(range(n))
32    k = R.shape[1]
33    target = (m+1) * max(r, 1)
34    if max_iter is None: max_iter = n + 5
35    it = 0
36    while len(active) > target and it < max_iter:
37        # Any k+1 atoms are linearly dependent. If this particular block has
38        # rank deficiency, its null vector still gives a valid elimination.
39        block = np.asarray(active[:min(k+1, len(active))], dtype=int)
40        _, s, vh = np.linalg.svd(R[block].T, full_matrices=True)
41        c = vh[-1]
42        if np.linalg.norm(c) < tol:
43            break
44        if np.all(c <= tol): c = -c
45        pos = c > tol
46        if not np.any(pos): break
47        wb = weights[block]
48        t = np.min(wb[pos] / c[pos])
49        weights[block] -= t*c
50        weights[np.abs(weights) < 1e-12] = 0.0
51        active = [i for i in active if weights[i] > 1e-12]
52        it += 1
53    active = np.asarray([i for i in active if weights[i] > 1e-12], dtype=int)
54    return active, weights[active], W, r, it
55
56
57def weighted_loss(x, y, idx, weights, W):
58    e = x[idx] @ W.T - y[idx]
59    return float(np.sum(weights[:, None] * e*e))