The method to calculate gradient in this case is Calculus (analytically, NOT numerically!). So we differentiate loss function with respect to W(yi) like this:

and with respect to W(j) when j!=yi is:

The 1 is just indicator function so we can ignore the middle form when condition is true. And when you write in code, the example you provided is the answer.

Since you are using cs231n example, you should definitely check note and videos if needed.

Hope this helps!

To answer your question, unless you have a very good idea of why you want to define a custom kernel, I'd stick with the built-ins. They are very fast, flexible, and powerful, and are well-suited to most applications.

That being said, let's go into a bit more detail:

A Kernel Function is a special kind of measure of similarity between two points. Basically a larger value of the similarity means the points are more similar. The scikit-learn SVM is designed to be able to work with any kernel function. Several kernels built-in (e.g. linear, radial basis function, polynomial, sigmoid) but you can also define your own.

Your custom kernel function should look something like this:

def my_kernel(x, y):
    """Compute My Kernel

    Parameters
    ----------
    x : array, shape=(N, D)
    y : array, shape=(M, D)
        input vectors for kernel similarity

    Returns
    -------
    K : array, shape=(N, M)
        matrix of similarities between x and y
    """
    # ... compute something here ...
    return similarity_matrix

The most basic kernel, a linear kernel, would look like this:

def linear_kernel(x, y):
    return np.dot(x, y.T)

Equivalently, you can write

def linear_kernel_2(x, y):
    M = np.array([[1, 0],
                  [0, 1]])
    return np.dot(x, np.dot(M, y.T))

The matrix M here defines the so-called inner product space in which the kernel acts. This matrix can be modified to define a new inner product space; the custom function from the example you linked to just modifies M to effectively double the importance of the first dimension in determining the similarity.

More complicated non-linear modifications are possible as well, but you have to be careful: kernel functions must meet certain requirements (they must satisfy the properties of an inner-product space) or the SVM algorithm will not work correctly.

I don't have enough reputation to comment, so I am answering here. Whenever you compute loss vector for x[i], ith training example and get some nonzero loss, that means you should move your weight vector for the incorrect class (j != y[i]) away by x[i], and at the same time, move the weights or hyperplane for the correct class (j==y[i]) near x[i]. By parallelogram law, w + x lies in between w and x. So this way w[y[i]] tries to come nearer to x[i] each time it finds loss>0.

Thus, dW[:,y[i]] += -X[i] and dW[:,j] += X[i] is done in the loop, but while update, we will do in direction of decreasing gradient, so we are essentially adding X[i] to correct class weights and going away by X[i] from weights that miss classify.

Here's a vectorized approach -

delta = 1
N = X.shape[0]
M = W.shape[1]
scoresv = X.dot(W)
marginv = scoresv - scoresv[np.arange(N), y][:,None] + delta

mask0 = np.zeros((N,M),dtype=bool)
mask0[np.arange(N),y] = 1
mask = (marginv<0) | mask0
marginv[mask] = 0

loss_out = marginv.sum()/num_train # mean
loss_out += 0.5 * reg * np.sum(W * W) # l2 regularization

Additionally, we could optimize np.sum(W * W) with np.tensordot, like so -

float(np.tensordot(W,W,axes=((0,1),(0,1))))

Runtime test

Proposed approach as function -

def svm_loss_vectorized_v2(W, X, y, reg):
    delta = 1
    N = X.shape[0]
    M = W.shape[1]
    scoresv = X.dot(W)
    marginv = scoresv - scoresv[np.arange(N), y][:,None] + delta

    mask0 = np.zeros((N,M),dtype=bool)
    mask0[np.arange(N),y] = 1
    mask = (marginv<=0) | mask0
    marginv[mask] = 0

    loss_out = marginv.sum()/num_train # mean
    loss_out += 0.5 * reg * float(np.tensordot(W,W,axes=((0,1),(0,1))))
    return loss_out

Timings -

In [86]: W= np.random.randn(3073,10)
    ...: X= np.random.randn(500,3073)
    ...: y= np.random.randint(0,10,(500))
    ...: reg = 4.56
    ...: 

In [87]: svm_loss_naive(W, X, y, reg)
Out[87]: 70380.938069371899

In [88]: svm_loss_vectorized_v2(W, X, y, reg)
Out[88]: 70380.938069371914

In [89]: %timeit svm_loss_naive(W, X, y, reg)
100 loops, best of 3: 10.2 ms per loop

In [90]: %timeit svm_loss_vectorized_v2(W, X, y, reg)
100 loops, best of 3: 2.94 ms per loop

Let's use the example of the SVM loss function for a single datapoint:

$L_i = \sum_{j\neq y_i} \left[ \max(0, w_j^Tx_i - w_{y_i}^Tx_i + \Delta) \right]$

Where $\Delta$ is the desired margin.

We can differentiate the function with respect to the weights. For example, taking the gradient with respect to $w_{yi}$ we obtain:

$\nabla_{w_{y_i}} L_i = - \left( \sum_{j\neq y_i} \mathbb{1}(w_j^Tx_i - w_{y_i}^Tx_i + \Delta > 0) \right) x_i$

Where 1 is the indicator function that is one if the condition inside is true or zero otherwise. While the expression may look scary when it is written out, when you're implementing this in code you'd simply count the number of classes that didn't meet the desired margin (and hence contributed to the loss function) and then the data vector $x_i$ scaled by this number is the gradient. Notice that this is the gradient only with respect to the row of $W$ that corresponds to the correct class. For the other rows where $j≠{{y}_{i}}$ the gradient is:

$\nabla_{w_j} L_i = \mathbb{1}(w_j^Tx_i - w_{y_i}^Tx_i + \Delta > 0) x_i$

Once you derive the expression for the gradient it is straight-forward to implement the expressions and use them to perform the gradient update.

Taken from Stanford CS231N optimization notes posted on github.