<jupyter_start><jupyter_text>Digit Recognition 1. IntroductionIn this analysis, the handwritten digits are identified using support vector machines and radial basis functions. 1.1 LibrariesThe essential libraries used here are numpy, matplotlib, and scikit-learn.  For convenience, pandas and IPython.display are used for displaying tables, and tqdm is used for progress bars.<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from itertools import product
from sklearn.svm import SVC
from sklearn.model_selection import cross_val_score, cross_val_predict, ShuffleSplit, KFold
from tqdm import tqdm
from IPython.display import display, Math, Latex, HTML
%matplotlib inline
np.set_printoptions(precision=4,threshold=200)
tqdm_bar_fmt='{percentage:3.0f}%|{bar}|'<jupyter_output><empty_output><jupyter_text>1.2 DatasetThe US Postal Service Zip Code dataset is used, which contains handwritten digits zero to nine.  The data has been preprocessed, whereby features of intensity and symmetry are extracted.<jupyter_code>def download_data():
    train_url = "http://www.amlbook.com/data/zip/features.train"
    test_url  = "http://www.amlbook.com/data/zip/features.test"
    column_names = ['digit','intensity','symmetry']
    train = pd.read_table(train_url,names=column_names,header=None,delim_whitespace=True)
    test = pd.read_table(test_url,names=column_names,header=None,delim_whitespace=True)
    train.digit = train.digit.astype(int)
    test.digit = test.digit.astype(int)
    return train,test

def process_data(train,test):
    X_train = train.iloc[:,1:].values
    y_train = train.iloc[:,0].values
    X_test = test.iloc[:,1:].values
    y_test = test.iloc[:,0].values
    return X_train,y_train,X_test,y_test
train,test = download_data()
X_train,y_train,X_test,y_test = process_data(train,test)<jupyter_output><empty_output><jupyter_text>2. Support Vector Machines for Digit Recognition  2.1 Polynomial KernelsWe wish to implement the following polynomial kernel for our support vector machine:$$K\left(\mathbf{x_n,x_m}\right) = \left(1+\mathbf{x_n^Tx_m}\right)^Q$$This is implemented in scikit-learn in the subroutine [sklearn.svm.SVC](http://scikit-learn.org/stable/modules/svm.html), where the kernel function takes the form:$$\left(\gamma \langle x,x' \rangle + r\right)^d$$where $d$ is specified by the keyword `degree`, and $r$ by `coef0`.  2.1.1 One vs Rest ClassificationIn the following subroutine, the data is split into "one-vs-rest", where $y=1$ corresponds to a match to the digit, and $y=0$ corresponds to all the other digits.  The training step is implemented in the call to `clf.fit()`.<jupyter_code>def get_misclassification_ovr(X_train,y_train,X_test,y_test,digit,
                              Q=2,r=1.0,C=0.01,kernel='poly',verbose=False):
    clf = SVC(C=C, kernel=kernel, degree=Q, coef0 = r, gamma = 1.0,
          decision_function_shape='ovr', verbose=False)
    y_in  = (y_train==digit).astype(int)
    y_out = (y_test==digit).astype(int)
    model = clf.fit(X_train,y_in)  # print(model)
    E_in  = np.mean(y_in  != clf.predict(X_train))
    E_out = np.mean(y_out != clf.predict(X_test))
    n_support_vectors = len(clf.support_vectors_)
    if verbose is True:
        print()
        print("Q = {}, C = {}: Support vectors: {}".format(Q, C, n_support_vectors))
        print("{} vs all: E_in  = {}".format(digit,E_in))
        print("{} vs all: E_out = {}".format(digit,E_out))
    return E_in,E_out,n_support_vectors<jupyter_output><empty_output><jupyter_text>The following code trains on the data for the cases: 0 vs all, 1 vs all, ..., 9 vs all.  For each of the digits, 0 to 9, the errors $E_{in}, E_{out}$ and the number of support vectors are recorded and stored in a pandas dataframe.<jupyter_code>results = pd.DataFrame()
i=0
for digit in tqdm(range(10),bar_format=tqdm_bar_fmt):
    ei, eo, n = get_misclassification_ovr(X_train,y_train,X_test,y_test,digit)
    df = pd.DataFrame({'digit': digit, 'E_in': ei, 'E_out': eo, 'n': n}, index=[i])
    results = results.append(df)
    i += 1
display(HTML(results[['digit','E_in','E_out','n']].iloc[::2].to_html(index=False)))
display(HTML(results[['digit','E_in','E_out','n']].iloc[1::2].to_html(index=False)))
from tabulate import tabulate
print(tabulate(results, headers='keys', tablefmt='simple'))<jupyter_output>E_in      E_out    digit     n
--  ---------  ---------  -------  ----
 0  0.105884   0.111609         0  2179
 1  0.0144013  0.0219233        1   386
 2  0.100261   0.0986547        2  1970
 3  0.0902483  0.0827105        3  1950
 4  0.0894253  0.0996512        4  1856
 5  0.0762584  0.079721         5  1585
 6  0.0910712  0.0847035        6  1893
 7  0.0884652  0.0732436        7  1704
 8  0.0743382  0.0827105        8  1776
 9  0.0883281  0.0881913        9  1978<jupyter_text>2.1.2 One vs One ClassificationOne vs one classification makes better use of the data, but is more computatationally expensive.  The following subroutine splits the data so that $y=0$ for the first digit, and $y=1$ for the second digit.  The rows of data corresponding to all other digits are removed.<jupyter_code>def get_misclassification_ovo(X_train,y_train,X_test,y_test,digit1,digit2,
                              Q=2,r=1.0,C=0.01,kernel='poly'):
    clf = SVC(C=C, kernel=kernel, degree=Q, coef0 = r, gamma = 1.0,
          decision_function_shape='ovo', verbose=False)
    select_in = np.logical_or(y_train==digit1,y_train==digit2)
    y_in  = (y_train[select_in]==digit1).astype(int)
    X_in  =  X_train[select_in]
    select_out = np.logical_or(y_test==digit1,y_test==digit2)
    y_out = (y_test[select_out]==digit1).astype(int)
    X_out =  X_test[select_out]
    model = clf.fit(X_in,y_in)

    E_in  = np.mean(y_in  != clf.predict(X_in))
    E_out = np.mean(y_out != clf.predict(X_out))
    n_support_vectors = len(clf.support_vectors_)
    return E_in,E_out,n_support_vectors<jupyter_output><empty_output><jupyter_text>In the following code, a 1-vs-5 classifier is tested for $Q=2,5$ and $C=0.001,0.01,0.1,1$.<jupyter_code>C_arr = [0.0001, 0.001, 0.01, 1]
Q_arr = [2, 5]
CQ_arr = list(product(C_arr,Q_arr))
results = pd.DataFrame()
i=0
for C, Q in tqdm(CQ_arr,bar_format=tqdm_bar_fmt):
    ei, eo, n = get_misclassification_ovo(X_train,y_train,X_test,y_test,
                                          digit1=1,digit2=5,Q=Q,r=1.0,C=C)
    df = pd.DataFrame({'C': C, 'Q': Q, 'E_in': ei, 'E_out': eo, 'n': n}, index=[i])
    results = results.append(df)
    i += 1
display(HTML(results[['C','Q','E_in','E_out','n']].to_html(index=False)))<jupyter_output>100%|██████████|<jupyter_text>2.1.3 Polynomial Kernel with Cross-ValidationFor k-fold cross-validation, the subroutine [`sklearn.model_selection.KFold`](http://scikit-learn.org/stable/modules/generated/sklearn.model_selection.KFold.htmlsklearn.model_selection.KFold)`()` is employed to split the data into folds.  Random shuffling is enabled with the `shuffle=True` option.<jupyter_code>def select_digits_ovo(X,y,digit1,digit2):
    subset = np.logical_or(y==digit1,y==digit2)
    X_in  =  X[subset].copy()
    y_in  = (y[subset]==digit1).astype(int).copy()
    return X_in,y_in
def get_misclassification_ovo_cv(X_in,y_in,Q=2,r=1.0,C=0.01,kernel='poly'):
    kf = KFold(n_splits=10, shuffle=True)
    kf.get_n_splits(X_in)
    E_cv = []
    for train_index, test_index in kf.split(X_in):
        X_trn, X_tst = X_in[train_index], X_in[test_index]
        y_trn, y_tst = y_in[train_index], y_in[test_index]
        clf = SVC(C=C, kernel=kernel, degree=Q, coef0 = r, gamma = 1.0,
              decision_function_shape='ovo', verbose=False)
        model = clf.fit(X_trn,y_trn)
        E_cv.append(np.mean(y_tst != clf.predict(X_tst)))
    return E_cv<jupyter_output><empty_output><jupyter_text>In this example, our parameters are $Q=2, C \in \{0.0001,0.001,0.01,0.1,1\}$ and we are considering the 1-vs-5 classifier.<jupyter_code>C_arr = [0.0001, 0.001, 0.01, 0.1, 1]; d1 = 1; d2 = 5
Q=2
X_in,y_in = select_digits_ovo(X_train,y_train,d1,d2)<jupyter_output><empty_output><jupyter_text>Due to the effect of random shuffling, 100 runs are carried out, and the results tabulated.<jupyter_code>E_cv_arr = []
count_arr = []
for n in tqdm(range(100),bar_format=tqdm_bar_fmt):
    counts = np.zeros(len(C_arr),int)
    kf = KFold(n_splits=10, shuffle=True)
    kf.get_n_splits(X_in)
    for train_index, test_index in kf.split(X_in):
        X_trn, X_tst = X_in[train_index], X_in[test_index]
        y_trn, y_tst = y_in[train_index], y_in[test_index]
        E_cv = []
        for C in C_arr:
            clf = SVC(C=C, kernel='poly', degree=Q, coef0=1.0, gamma=1.0,
                  decision_function_shape='ovo', verbose=False)
            model = clf.fit(X_trn,y_trn)
            E_cv.append(np.mean(y_tst != clf.predict(X_tst)))
        counts[np.argmin(E_cv)] += 1
        E_cv_arr.append(np.array(E_cv))
    count_arr.append(counts)<jupyter_output>100%|██████████|<jupyter_text>The number of times each particular value of $C$ is picked (for having the smallest $E_{cv}$) is calculated, as well as the mean cross validation error.  The data is tabulated as follows:<jupyter_code>df = pd.DataFrame({'C': C_arr})
df['count'] = np.sum(np.array(count_arr),axis=0)
df['E_cv']  = np.mean(np.array(E_cv_arr),axis=0)
display(HTML(df.to_html(index=False)))
fig = plt.figure(figsize=(10,5))
ax1 = fig.add_subplot(121)
_ = df.plot(y=['count'],ax=ax1,marker='x',color='red',markersize=10,grid=True)
ax2 = fig.add_subplot(122)
_ = df.plot(y=['E_cv'],ax=ax2,marker='o',color='green',markersize=7,grid=True)<jupyter_output><empty_output><jupyter_text>2.2 Scaling Law for Support Vector MachinesThis is based on [scikit-learn example code](http://scikit-learn.org/stable/tutorial/basic/tutorial.html), but modified to demonstrate how the training time scales with the size of the data.  The dataset in this example is the original MNIST data.<jupyter_code>from sklearn.datasets import fetch_mldata
from sklearn import svm
import timeit
digits = fetch_mldata('MNIST original')  # stores data in ~/scikit_learn_data by default
n_max = len(digits.target)
selection = np.random.permutation(n_max)
n_arr = [500, 1000, 1500, 2000, 5000, 10000]
t_arr = []
for n in n_arr:
    sel = selection[:n]
    X = digits.data[sel]
    y = digits.target[sel]
    clf = svm.SVC(gamma=0.001, C=100.)
    t0 = timeit.default_timer()
    clf.fit(X,y)
    t_arr.append(timeit.default_timer() - t0)
    print("n = {}, time = {}".format(n, t_arr[-1]))
plt.plot(np.array(n_arr),np.array(t_arr),'bx-')
plt.xlabel('no. of samples')
plt.ylabel('time (s)')
plt.grid()<jupyter_output><empty_output><jupyter_text>3. Radial Basis Functions 3.1 BackgroundThe hypothesis is given by:$$h\left(\mathbf{x}\right) = \sum\limits_{k=1}^K w_k \exp\left(-\gamma \Vert \mathbf{x} - \mu_k \Vert^2\right)$$This is implemented in the subroutine [`sklearn.svm.SVC`](http://scikit-learn.org/stable/modules/svm.html)`(..., kernel='rbf', ...)` as shown in the code below.<jupyter_code>def get_misclassification_rbf(X_train,y_train,X_test,y_test,C,digit1=1,digit2=5):
    clf = SVC(C=C, kernel='rbf', gamma = 1.0,
          decision_function_shape='ovo', verbose=False)
    select_in = np.logical_or(y_train==digit1,y_train==digit2)
    y_in  = (y_train[select_in]==digit1).astype(float)
    X_in  =  X_train[select_in]
    select_out = np.logical_or(y_test==digit1,y_test==digit2)
    y_out = (y_test[select_out]==digit1).astype(float)
    X_out =  X_test[select_out]
    model = clf.fit(X_in,y_in)
    E_in  = np.mean(y_in  != clf.predict(X_in))
    E_out = np.mean(y_out != clf.predict(X_out))
    n_support_vectors = len(clf.support_vectors_)
    return E_in,E_out,n_support_vectors<jupyter_output><empty_output><jupyter_text>For $C \in \{0.01, 1, 100, 10^4, 10^6\}$, the in-sample and out-of-sample errors, as well as the number of support vectors are tabulated as follows:<jupyter_code>C_arr = [0.01, 1., 100., 1e4, 1e6]
results = pd.DataFrame()
i=0
for C in tqdm(C_arr,bar_format=tqdm_bar_fmt):
    ei, eo, n = get_misclassification_rbf(X_train,y_train,X_test,y_test,
                                          C,digit1=1,digit2=5)
    df = pd.DataFrame({'C': C, 'E_in': ei, 'E_out': eo, 'n': n}, index=[i])
    results = results.append(df)
    i += 1
display(HTML(results.to_html(index=False)))<jupyter_output>100%|██████████|<jupyter_text>The table is also plotted graphically below:<jupyter_code>def plot_Ein_Eout(df):
    fig = plt.figure(figsize=(12,5))
    ax1 = fig.add_subplot(121)
    df.plot(ax=ax1, x='C', y='E_in', kind='line', marker='o', markersize=7, logx=True)
    df.plot(ax=ax1, x='C', y='E_out', kind='line', marker='x', markersize=7, logx=True)
    ax1.legend(loc='best',frameon=False)
    ax1.grid(True)
    ax1.set_title('In-Sample vs Out-of-Sample Errors')

    ax2 = fig.add_subplot(122)
    df.plot(ax=ax2, x='C', y='n', kind='line', marker='+', markersize=10, logx=True)
    ax2.legend(loc='best',frameon=False)
    ax2.grid(True)
    ax2.set_title('Number of Support Vectors')

plot_Ein_Eout(results)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14742071073748567, 'fraction_numerical': 0.040275124188001526, 'mean_word_length': 4.345588235294118, 'pattern_counts': {'":': 0, '<': 53, '<?xml version=': 0, '>': 53, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '45953', 'license': 'Apache-2.0', 'n_tokens_mistral': 5454, 'n_tokens_neox': 5066, 'n_words': 1023, 'path': 'perceptron/digit-recognition.ipynb', 'repo_name': 'nathanielng/machine-learning', 'size': '13085'}