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