This type of SVM is often implemented with the SMO algorithm. You may want check for the original published version (Platt, John. Fast Training of Support Vector Machines using Sequential Minimal Optimization, in Advances in Kernel Methods Support Vector Learning, B. Scholkopf, C. Burges, A. Smola, eds., MIT Press (1998)), but it is quite complicated as for me.

A bit simplified version is presented in Stanford Lecture Notes, but derivation of all the formulas should be found somewhere else (e.g. this random notes I found on the Internet).

As an alternative I can propose you my own variation of the SMO algorithm. It is highly simplified, implementation contains a bit more than 30 lines of code

class SVM:
  def __init__(self, kernel='linear', C=10000.0, max_iter=100000, degree=3, gamma=1):
    self.kernel = {'poly':lambda x,y: np.dot(x, y.T)**degree,
                   'rbf':lambda x,y:np.exp(-gamma*np.sum((y-x[:,np.newaxis])**2,axis=-1)),
                   'linear':lambda x,y: np.dot(x, y.T)}[kernel]
    self.C = C
    self.max_iter = max_iter

  def restrict_to_square(self, t, v0, u):
    t = (np.clip(v0 + t*u, 0, self.C) - v0)[1]/u[1]
    return (np.clip(v0 + t*u, 0, self.C) - v0)[0]/u[0]

  def fit(self, X, y):
    self.X = X.copy()
    self.y = y * 2 - 1
    self.lambdas = np.zeros_like(self.y, dtype=float)
    self.K = self.kernel(self.X, self.X) * self.y[:,np.newaxis] * self.y
    
    for _ in range(self.max_iter):
      for idxM in range(len(self.lambdas)):
        idxL = np.random.randint(0, len(self.lambdas))
        Q = self.K[[[idxM, idxM], [idxL, idxL]], [[idxM, idxL], [idxM, idxL]]]
        v0 = self.lambdas[[idxM, idxL]]
        k0 = 1 - np.sum(self.lambdas * self.K[[idxM, idxL]], axis=1)
        u = np.array([-self.y[idxL], self.y[idxM]])
        t_max = np.dot(k0, u) / (np.dot(np.dot(Q, u), u) + 1E-15)
        self.lambdas[[idxM, idxL]] = v0 + u * self.restrict_to_square(t_max, v0, u)
    
    idx, = np.nonzero(self.lambdas > 1E-15)
    self.b = np.sum((1.0-np.sum(self.K[idx]*self.lambdas, axis=1))*self.y[idx])/len(idx)
  
  def decision_function(self, X):
    return np.sum(self.kernel(X, self.X) * self.y * self.lambdas, axis=1) + self.b

In simple cases it works not much worth than sklearn.svm.SVC, comparison shown below

I have posted this code with some more code producing images for comparison on GitHub. For more elaborate explanation with formulas you may want to refer to my preprint on ResearchGate.

UPDATE: now live version is available, see Github Pages

Say that mat1 is $n \times d$ and mat2 is $m \times d$.

Recall that the Gaussian RBF kernel is defined as $k(x, y) = \exp\left( - \frac{1}{2 \sigma^2} \lVert x - y \rVert^2 \right)$. But we can write $\lVert x - y \rVert^2$ as $(x - y)^T (x - y) = x^T x + y^T y - 2 x^T y$. The code uses this decomposition.

First, the trnorms1 vector stores $x^T x$ for each input $x$ in mat1, and trnorms2 stores $y^T y$ for each $y$ in mat2.

Then, the k1 matrix is obtained by multiplying the $n \times 1$ matrix of $x^T x$ entries by a $1 \times m$ matrix of ones, getting an $n \times m$ matrix with $x^T x$ entries repeated across the rows, so that k1[i, j] is $x_i^T x_i$.

The next line does basically the same thing for the $y$ norms repeated across columns, getting an $n \times m$ matrix with k2[i, j] of $y_j^T y_j$.

k is then their sum, so that k[i, j] is $x_i^T x_i + y_j^T y_j$. The next line then subtracts twice the product of the data matrices, so that k[i, j] becomes $x_i^T x_i + y_j^T y_j - 2 x_i^T y_j = \lVert x_i - y_j \rVert^2$.

Then, the code multiplies by $\frac{-1}{2 \sigma^2}$ and finally takes the elementwise $\exp$, getting out the Gaussian kernel.

If you dig into the scikit-learn implementation, it's exactly the same, except:

• It's parameterized instead with $\gamma = \frac{1}{2 \sigma^2}$.
• It's written in much better Python, not wasting memory all over the place and doing computations in a needlessly slow way.
• It's broken up into helper functions.

But, algorithmically, it's doing the same basic operations.