Why do minor variations in data scaling produce different SVM class boundaries?

Question

Anti-Closing-Preamble: I know why standardisation ("scaling") is useful and I've read this, this, this, this, this, this, this, this, this, this, this, and other questions. My question is different.

I have a very simple data set---two linearly separable classes, two features, 21 data points---with already reasonable values. I nevertheless standardise it (subtract the mean and divide by the standard deviation), train a linear SVM, and obtain a reasonable class boundary:

However, if I scale the data set (simply multiplying both features by a scalar), the slope of the boundary changes:

Why is it so? There should be nothing magical about the standard deviation of the data being exactly 1, instead of 0.5 or 2. For such small scaling factors numerical instability also shouldn't be the issue. The data are scaled equally along both dimensions, so it's not that one feature suddenly gets the upper hand. The slack variables are unit-less, measured in terms of the margin, and should therefore be invariant to scaling. So what else could be going on here?

Below is the code in case you want to reproduce the results.

import numpy as np
import matplotlib.pyplot as plt
from sklearn import svm


# just for plotting:
def plotSVM(X, y, clf):
  plt.title(f'Scaling = {scaling}')
  w = clf.coef_[0]
  b = clf.intercept_[0]
  k = -w[0] / w[1]
  l = -b    / w[1]
  minX = min(X[:, 0])
  maxX = max(X[:, 0])
  dX   = maxX-minX
  x1 = np.array([minX-.1*dX, maxX+.1*dX])
  x2 = k*x1 + l
  plt.plot(x1, x2, '-', color='gold', linewidth=2)
  dY = max(X[:, 1]) - min(X[:, 1])
  sv  = clf.support_vectors_
  plt.scatter(sv[:, 0], sv[:, 1], s=64, facecolor='none', edgecolor='gray')
  ix = y == +1
  plt.scatter(X[ix, 0], X[ix, 1], s=64, color='seagreen',  marker='+')
  ix = y == -1
  plt.scatter(X[ix, 0], X[ix, 1], s=64, color='firebrick', marker='_')
  plt.xlim(x1)
  plt.ylim([min(X[:, 1])-.1*dY, max(X[:, 1])+.1*dY])
  plt.show()
  

# the original data
X = np.log([
  [29.7 ,  1.68], [29.6 ,  2.13], [32.8 ,  1.46], [24.  ,  2.75], [27.  ,  2.3 ],
  [21.  ,  2.1 ], [29.4 ,  2.34], [12.  ,  1.27], [48.  ,  0.88], [50.  ,  0.97],
  [27.  ,  1.8 ], [ 8.5 ,  0.6 ], [ 5.7 ,  0.5 ], [ 4.3 ,  0.43], [ 3.8 ,  0.52],
  [ 3.2 ,  0.3 ], [ 4.6 ,  0.34], [14.6 ,  0.7 ], [ 9.2 ,  0.54], [ 4.2 ,  0.29],
  [ 4.1 ,  0.32]
])
y = np.array([ 1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,
              -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 ])
### standardise X to mean=0 and sd = 1:
X_std = (X - np.mean(X, axis=0)) / np.std(X, axis=0)


### scale it by a factor:
for scaling in [.5, 1, 2]:
  X_sc = scaling * X_std
  clf = svm.SVC(kernel='linear', C=.1, random_state=1)
  clf.fit(X_sc, y)
  plotSVM(X_sc, y, clf)

Answer 1

The reason lies in the relationship between $\mathbf{w}$, the vector determining the boundary hyperplane, and the slack variables, $\mathbf{\xi}$. The slack variables are unit-less, but the boundary vector $\mathbf{w}$ isn't (but this fact gets somewhat obscured through standardisation).

As we scale up the inputs $\mathbf{x}_i$, $\mathbf{w}$ must scale down in order to keep the functional margin constant ($\pm 1$). So the balance between the slack variables $\xi_i$ and the boundary vector $\mathbf{w}$ shifts as we scale the data.

In more detail: The SVM classifier solves the optimisation problem:

\begin{align} \min_{\mathbf{w}, b, \mathbf{\xi}} \quad & \frac{1}{2}\| \mathbf{w}\|^2 + C\sum_i \mathbf{\xi}_i \\ \text{s.t.} \quad & y_i(\mathbf{w} \cdot \mathbf{x}_i + b) \geq 1-\xi_i \\ \quad & \xi_i \geq 0 \end{align}

So, if $\mathbf{x}_i$'s are small, $\mathbf{w}$ will need to be large in order to satisfy the first condition and it will dominate the objective function. On the other hand, if $\mathbf{x}_i$'s are large, $\mathbf{w}$ will be small and, for a sufficiently large $C$, the objective function will be dominated by the $\mathbf{\xi}_i$'s.

In order to keep the class boundary invariant to scaling, we need to divide $C$ by the square of the scaling factor:

for scaling in [.5, 1, 2]:
  X_sc = scaling * X_std
  clf = svm.SVC(kernel='linear', C=.1/scaling**2, random_state=1)
  clf.fit(X_sc, y)
  plotSVM(X_sc, y, clf)

produces: