{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\nimport plotly.express as px\n\nimport matplotlib.pyplot as plt\n\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.model_selection import GridSearchCV\n\nfrom sklearn.metrics import pairwise_distances\nfrom sklearn.neighbors import KNeighborsClassifier\nfrom sklearn.cluster import KMeans\n\nfrom sklearn.decomposition import PCA\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-06-21T15:03:49.405988Z","iopub.execute_input":"2021-06-21T15:03:49.406691Z","iopub.status.idle":"2021-06-21T15:03:49.422505Z","shell.execute_reply.started":"2021-06-21T15:03:49.406633Z","shell.execute_reply":"2021-06-21T15:03:49.420879Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Hand-written digit recognition with prototype methods\n\nIn this notebook we approach the problem of digit recognition from the point of view of prototype methods, specifically nearest-neighbour methods.\n\nNotice that the abundance of observation in the traning set means we might find it difficult to use K-NN as it is, since it will take a long time for the algorithm to compute the distance between all observations.\n\nFor this reasons, we have:\n\n1. a baseline method where we randomly sample 10k observations and run K-NN with that\n2. an hybrid method where we first cluster the observation and then run K-NN\n3. an hybrid method where we perform PCA and then run K-NN in the reduced vector space\n4. same as method 2, but we use the tangent distance between observations\n\nWe report the accuracy on test set for each of the approaches.\n\n### Function definition\n\nSome useful functions we'll use a couple of times,","metadata":{}},{"cell_type":"code","source":"def show_random_image(df,size = 28):\n    \"\"\" Plot gray image from 1D array of numbers \"\"\"\n\n    example_row = df.sample(1).to_numpy().flatten()\n    example_label = example_row[0]\n    example_digit = example_row[1:]\n\n\n    fig = px.imshow(\n        example_digit.reshape((size,size)),\n        color_continuous_scale='gray',\n        title= 'digit {}'.format(example_label)\n    )\n\n    fig.show()\n    \ndef rescale_image_values(image,new_image):\n    \"\"\" Rescale the image to keep it within [0,255]. \"\"\"\n\n    a0, b0 = image.min(), image.max()\n    a, b = new_image.min(), new_image.max()\n\n    new_image = a0 + (new_image - a) * (b0 - a0) / (b - a)\n\n    return new_image.astype(np.int)\n\ndef image_transform(\n    image : np.ndarray,\n    theta : float,\n    mx : float,\n    my : float,\n    sheer : float,\n    gamma : float) -> np.ndarray:\n    '''\n    Perform a rotation, re-shape, sheer of the image.\n    \n    Input parameters\n        - image : the image aim to transform (1D array)\n        - theta : angle of rotation for the image (0,2*pi)\n        - mx : magnification coefficient for the x-axis\n        - my : magnification coefficient for the y-axis\n        - sheer : sheer parameter\n        - gamma : locality parameter for the radial basis kernel\n    '''\n\n    # Compute the size of the image assuming it is squared\n    n_pixels = image.size\n    size = int(np.sqrt(n_pixels))\n\n    # The centre of the image is the origin of the rotation\n    x0 = (size-1)/2\n    y0 = (size-1)/2\n\n    # We create a grid of pixels that will be transformed\n    pixel_range = np.arange(size,dtype=np.float)\n    x,y = np.meshgrid(pixel_range,pixel_range)\n\n    x = x.flatten()\n    y = y.flatten()\n\n    # Rotate, re-shape, sheer the image\n    x_transformed = x0 + mx * np.cos(theta) * (x - x0) + (my * np.sin(theta) + sheer * np.cos(theta)) * (y - y0)\n    y_transformed = y0 - mx * np.sin(theta) * (x - x0) + (my * np.cos(theta) - sheer * np.sin(theta)) * (y - y0)\n\n    # Compute the distance between any point in the grid and the transformed points\n    x_distance = np.broadcast_to(x,(n_pixels,n_pixels)).T - x_transformed\n    y_distance = np.broadcast_to(y,(n_pixels,n_pixels)).T - y_transformed\n\n    # Use a radial basis kernel to smooth the image\n    kernel = np.exp(-gamma*(x_distance**2 + y_distance**2))\n\n    new_image = np.sum(kernel * image,axis=1).astype(np.int)\n\n    # Rescale the values to get them between 0 and 255\n    new_image = rescale_image_values(image,new_image)\n\n    return new_image\n\ndef first_derivative_rotation(image : np.ndarray, gamma : float) -> np.ndarray:\n    '''\n    Compute the tangent to the curve obtained by rotating the image.\n    \n    Input parameters\n        - image : the array of images we aim to transform (2D array)\n        - gamma : locality parameter for the radial basis kernel\n    '''\n\n    # Compute the size of the image assuming it is squared\n    _,n_pixels = image.shape\n    size = int(np.sqrt(n_pixels))\n\n    # Shape we broadcast the arrays (to compute any possible distance at once)\n    shape_bc = (n_pixels,n_pixels)\n\n    # The centre of the image is the origin of the rotation\n    x0 = (size-1)/2\n    y0 = (size-1)/2\n\n    # We create a grid of pixels that will be transformed\n    pixel_range = np.arange(size,dtype=np.float)\n    x,y = np.meshgrid(pixel_range,pixel_range)\n\n    x = x.flatten()\n    y = y.flatten()\n\n    # Compute the distance between any two points in the grid\n    x_distance = np.broadcast_to(x,shape_bc).T - x\n    y_distance = np.broadcast_to(y,shape_bc).T - y\n\n    # Compute the inner product between \n    inner_product_diff = y_distance * (np.broadcast_to(x,shape_bc) - x0) - x_distance * (np.broadcast_to(y,shape_bc) - y0)\n\n    # Use a radial basis kernel to smooth the image\n    kernel = np.exp(-gamma*(x_distance**2 + y_distance**2))\n\n    # first derivative of the image for a rotation around the centre\n    deriv_image = -2 * gamma * kernel * inner_product_diff @ image.T\n\n    return deriv_image.T\n\ndef optimal_line_distance(line1,line2):\n    ''' The optimal distance between two non-parallel lines. '''\n    \n    if np.all(line1==line2):\n        return 0\n    \n    # Line is defined by a tuple; first element is a point on the line,\n    # the second is the direction of the line\n    x1,n1 = np.split(line1,2)\n    x2,n2 = np.split(line2,2)\n    \n    dist = x1 - x2\n\n    overlap = n1 @ n2\n    d1 = dist @ n1\n    d2 = dist @ n2\n    \n    # We return the squared l2-distance (to avoid having to take a sqrt)\n    return dist @ dist - (d1**2 + d2**2 - 2 * overlap * d1 * d2)/(1 - overlap**2)\n\ndef format_data_tangent_method(X,gamma=3):\n    \"\"\" compute the first dervative of the rotated-image function and add to image\"\"\"\n\n    # compute the first dervative of the rotated-image function\n    X_deriv = first_derivative_rotation(X,gamma)\n\n    N,pixels = X_deriv.shape\n\n    # Compute its inverse norm\n    norm_X_deriv = np.linalg.norm(X_deriv,axis=1)\n    norm_X_deriv = np.full((pixels,N),norm_X_deriv).T\n\n    X_deriv = X_deriv / norm_X_deriv\n    \n    return np.hstack((X,X_deriv))","metadata":{"execution":{"iopub.status.busy":"2021-06-21T15:29:39.14434Z","iopub.execute_input":"2021-06-21T15:29:39.144846Z","iopub.status.idle":"2021-06-21T15:29:39.172051Z","shell.execute_reply.started":"2021-06-21T15:29:39.144798Z","shell.execute_reply":"2021-06-21T15:29:39.170589Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Loading the datasets\n\nLet's load the datasets for training and testing.\n\n### Training dataset","metadata":{}},{"cell_type":"code","source":"df_train = pd.read_csv('/kaggle/input/digit-recognizer/train.csv')\ndf_train.head()","metadata":{"execution":{"iopub.status.busy":"2021-06-21T14:19:17.514848Z","iopub.execute_input":"2021-06-21T14:19:17.515296Z","iopub.status.idle":"2021-06-21T14:19:21.544288Z","shell.execute_reply.started":"2021-06-21T14:19:17.515243Z","shell.execute_reply":"2021-06-21T14:19:21.543188Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's check if the dataset is balanced,","metadata":{}},{"cell_type":"code","source":"df_train['label'].value_counts()","metadata":{"execution":{"iopub.status.busy":"2021-06-21T14:19:21.546386Z","iopub.execute_input":"2021-06-21T14:19:21.546689Z","iopub.status.idle":"2021-06-21T14:19:21.559093Z","shell.execute_reply.started":"2021-06-21T14:19:21.546659Z","shell.execute_reply":"2021-06-21T14:19:21.557841Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It seems to be!\n\nLet's now visualize some random digit in the dataset,","metadata":{}},{"cell_type":"code","source":"show_random_image(df_train)","metadata":{"execution":{"iopub.status.busy":"2021-06-21T14:19:21.561154Z","iopub.execute_input":"2021-06-21T14:19:21.561506Z","iopub.status.idle":"2021-06-21T14:19:22.588834Z","shell.execute_reply.started":"2021-06-21T14:19:21.561471Z","shell.execute_reply":"2021-06-21T14:19:22.587754Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Test dataset","metadata":{}},{"cell_type":"code","source":"df_test = pd.read_csv('/kaggle/input/digit-recognizer/test.csv')\ndf_test.head()","metadata":{"execution":{"iopub.status.busy":"2021-06-21T14:19:22.590612Z","iopub.execute_input":"2021-06-21T14:19:22.590924Z","iopub.status.idle":"2021-06-21T14:19:25.02368Z","shell.execute_reply.started":"2021-06-21T14:19:22.590894Z","shell.execute_reply":"2021-06-21T14:19:25.022661Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_test = df_test.to_numpy()","metadata":{"execution":{"iopub.status.busy":"2021-06-21T14:19:25.025269Z","iopub.execute_input":"2021-06-21T14:19:25.02566Z","iopub.status.idle":"2021-06-21T14:19:25.029885Z","shell.execute_reply.started":"2021-06-21T14:19:25.025626Z","shell.execute_reply":"2021-06-21T14:19:25.029224Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Submission file,","metadata":{}},{"cell_type":"code","source":"submission = pd.read_csv('/kaggle/input/digit-recognizer/sample_submission.csv')\nsubmission.set_index('ImageId',inplace=True)","metadata":{"execution":{"iopub.status.busy":"2021-06-21T14:19:25.030938Z","iopub.execute_input":"2021-06-21T14:19:25.031361Z","iopub.status.idle":"2021-06-21T14:19:25.062917Z","shell.execute_reply.started":"2021-06-21T14:19:25.031303Z","shell.execute_reply":"2021-06-21T14:19:25.061614Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The random seed to use for the models that needs it,","metadata":{}},{"cell_type":"code","source":"rand_seed = 1864","metadata":{"execution":{"iopub.status.busy":"2021-06-21T14:19:25.065683Z","iopub.execute_input":"2021-06-21T14:19:25.066156Z","iopub.status.idle":"2021-06-21T14:19:25.071817Z","shell.execute_reply.started":"2021-06-21T14:19:25.066109Z","shell.execute_reply":"2021-06-21T14:19:25.070134Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Baseline\n\nThe baseline model we use is a K-NN where we sample randomly for the observation to use to train the model (the reason for this is that K-NN with ~40k observations is too much for us (takes too much time).\n\n### Pre-processing\n\nLet's first sample the dataset and create the training set,","metadata":{}},{"cell_type":"code","source":"N_cv = 8000\n\ndf_cv = df_train.sample(n=N_cv, random_state=rand_seed)\ndf_cv['label'].value_counts()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y_train = df_cv.pop('label').to_numpy()\nX_train = df_cv.to_numpy()\n\nstandardizer = StandardScaler()\nX_train = standardizer.fit_transform(X_train)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"And the test set,","metadata":{}},{"cell_type":"code","source":"X_baseline_test = standardizer.transform(X_test)","metadata":{"execution":{"iopub.status.busy":"2021-06-17T11:59:40.379932Z","iopub.execute_input":"2021-06-17T11:59:40.380335Z","iopub.status.idle":"2021-06-17T11:59:40.498422Z","shell.execute_reply.started":"2021-06-17T11:59:40.380292Z","shell.execute_reply":"2021-06-17T11:59:40.497081Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### K-Nearest-Neighbours\n\nWe can now run our baseline model, where we are making a gridsearch for optimizing the hyper-parameters of the problem (and using 5-fold CV to estimate the generalization error).","metadata":{}},{"cell_type":"code","source":"estimator = KNeighborsClassifier()\n\nest_params = {}\nest_params['n_neighbors'] = [1,2,3,4,5,6]\nest_params['weights'] = ['uniform','distance']\nest_params['p'] = [1,2,3]\n\ngrid = GridSearchCV(\n    estimator,\n    est_params,\n    scoring='accuracy',\n    n_jobs=-1,\n    refit=True,\n    cv=5,\n    verbose=2,\n    return_train_score=False\n)\n\ngrid.fit(X_train,y_train)\n\ndf_results = pd.DataFrame(grid.cv_results_)\ndf_results.sort_values('rank_test_score')\n\ndf_results.head()","metadata":{"execution":{"iopub.status.busy":"2021-06-17T09:58:36.757017Z","iopub.execute_input":"2021-06-17T09:58:36.757536Z","iopub.status.idle":"2021-06-17T11:46:28.297707Z","shell.execute_reply.started":"2021-06-17T09:58:36.757495Z","shell.execute_reply":"2021-06-17T11:46:28.296643Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = px.line(\n    df_results,\n    x='param_n_neighbors',\n    y='mean_test_score',\n    error_y='std_test_score',\n    color='param_p',\n    facet_col='param_weights'\n)\n\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2021-06-17T11:46:28.299563Z","iopub.execute_input":"2021-06-17T11:46:28.299906Z","iopub.status.idle":"2021-06-17T11:46:28.498587Z","shell.execute_reply.started":"2021-06-17T11:46:28.299857Z","shell.execute_reply":"2021-06-17T11:46:28.497681Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's now use the trained model to make predictions,","metadata":{}},{"cell_type":"code","source":"opt_model = grid.best_estimator_\ny_test = opt_model.predict(X_baseline_test)","metadata":{"execution":{"iopub.status.busy":"2021-06-17T12:01:14.445175Z","iopub.execute_input":"2021-06-17T12:01:14.445563Z","iopub.status.idle":"2021-06-17T12:04:12.011472Z","shell.execute_reply.started":"2021-06-17T12:01:14.445532Z","shell.execute_reply":"2021-06-17T12:04:12.010397Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's prepare the file for submission,","metadata":{}},{"cell_type":"code","source":"submission['Label'] = y_test\nsubmission.to_csv('baseline_results.csv')","metadata":{"execution":{"iopub.status.busy":"2021-06-17T12:11:22.413191Z","iopub.execute_input":"2021-06-17T12:11:22.413801Z","iopub.status.idle":"2021-06-17T12:11:22.479688Z","shell.execute_reply.started":"2021-06-17T12:11:22.413744Z","shell.execute_reply":"2021-06-17T12:11:22.478915Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Got an accuracy of 0.92985 from the baseline.\n\n## 1. Clustering and NN together\n\nHow about creating a large number of clusters and then running K-NN, instead of selecting the training set at random? This way we might be able to mantain an efficient scoring, but we would not simply forget about the whole information about the digits!\n\n### K-Means\n\nLet's first cluster the observation together. To be sure that we are not mixing digit together, we perform k-Means on each digit separately.","metadata":{}},{"cell_type":"code","source":"clusters = 1000\n\nmodel = KMeans(\n    n_clusters=clusters,\n    n_init=10,\n    max_iter=300,\n    random_state=rand_seed\n)\n\ndigits = df_train['label'].unique()\n\ndf_centroids = pd.DataFrame(columns=df_train.columns)\n_, *pixels = df_train.columns\n\nfor d in digits:\n    \n    # get images of digit g\n    df_digit = df_train.query('label == {}'.format(d))\n    X_digit = df_digit.drop(columns='label').to_numpy()\n    \n    # find centorids using K-means\n    model.fit(X_digit)\n\n    # Fill the dataframe with the centriods\n    column_form = lambda img : {'label':d,**dict(zip(pixels,img))}\n    centroids_columns = np.apply_along_axis(column_form,axis=1,arr=model.cluster_centers_)\n    \n    df_centroids = df_centroids.append(centroids_columns.tolist())","metadata":{"execution":{"iopub.status.busy":"2021-06-21T14:20:56.966267Z","iopub.execute_input":"2021-06-21T14:20:56.966764Z","iopub.status.idle":"2021-06-21T14:44:14.823604Z","shell.execute_reply.started":"2021-06-21T14:20:56.966727Z","shell.execute_reply":"2021-06-21T14:44:14.822041Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's have a look at he centriods (we do not expect them to look much different from standard hand-written digits),","metadata":{}},{"cell_type":"code","source":"show_random_image(df_centroids)","metadata":{"execution":{"iopub.status.busy":"2021-06-21T14:44:14.82607Z","iopub.execute_input":"2021-06-21T14:44:14.826522Z","iopub.status.idle":"2021-06-21T14:44:14.888989Z","shell.execute_reply.started":"2021-06-21T14:44:14.826483Z","shell.execute_reply":"2021-06-21T14:44:14.887569Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### K-Nearest-Neighbours after K-Means\n\nTo speed-up the process and still keep information of the whole dataset, we have first performed K-means on the images of each digit individually, and we can now perform K-means on this restricted dataset (from ~40k to 1k). Notice that the # of clusters used in the K-means step should be an hyper-parameter as well.\n\nLet's first get the dataset in a format that is accepted by the sklearn method,","metadata":{}},{"cell_type":"code","source":"y_cent = df_centroids.pop('label').to_numpy().astype(np.int)\nX_cent = df_centroids.to_numpy()\n\nstandardizer_cent = StandardScaler()\nX_cent = standardizer_cent.fit_transform(X_cent)","metadata":{"execution":{"iopub.status.busy":"2021-06-21T14:46:06.742049Z","iopub.execute_input":"2021-06-21T14:46:06.742477Z","iopub.status.idle":"2021-06-21T14:46:07.018085Z","shell.execute_reply.started":"2021-06-21T14:46:06.74244Z","shell.execute_reply":"2021-06-21T14:46:07.016855Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"and the test set,","metadata":{}},{"cell_type":"code","source":"X_centriod_test = standardizer_cent.transform(X_test)","metadata":{"execution":{"iopub.status.busy":"2021-06-21T14:46:11.592154Z","iopub.execute_input":"2021-06-21T14:46:11.592698Z","iopub.status.idle":"2021-06-21T14:46:11.770963Z","shell.execute_reply.started":"2021-06-21T14:46:11.592661Z","shell.execute_reply":"2021-06-21T14:46:11.77015Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's optimize the hyper-parameter of the model,","metadata":{}},{"cell_type":"code","source":"estimator = KNeighborsClassifier()\n\nest_params = {}\nest_params['n_neighbors'] = [1,2,3,4,5,6]\nest_params['weights'] = ['uniform','distance']\nest_params['p'] = [1,2,3]\n\ngrid = GridSearchCV(\n    estimator,\n    est_params,\n    scoring='accuracy',\n    n_jobs=-1,\n    refit=True,\n    cv=5,\n    verbose=2,\n    return_train_score=False\n)\n\ngrid.fit(X_cent,y_cent)\n\ndf_results_cent = pd.DataFrame(grid.cv_results_)\ndf_results_cent.sort_values('rank_test_score')\n\ndf_results_cent.head()","metadata":{"execution":{"iopub.status.busy":"2021-06-17T12:48:50.7396Z","iopub.execute_input":"2021-06-17T12:48:50.739951Z","iopub.status.idle":"2021-06-17T16:07:13.249976Z","shell.execute_reply.started":"2021-06-17T12:48:50.739922Z","shell.execute_reply":"2021-06-17T16:07:13.248685Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = px.line(\n    df_results_cent,\n    x='param_n_neighbors',\n    y='mean_test_score',\n    error_y='std_test_score',\n    color='param_p',\n    facet_col='param_weights'\n)\n\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2021-06-17T12:21:04.690967Z","iopub.execute_input":"2021-06-17T12:21:04.691476Z","iopub.status.idle":"2021-06-17T12:21:04.802619Z","shell.execute_reply.started":"2021-06-17T12:21:04.691445Z","shell.execute_reply":"2021-06-17T12:21:04.80154Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"opt_model = grid.best_estimator_\ny_test = opt_model.predict(X_centriod_test)\n\nsubmission['Label'] = y_test\nsubmission.to_csv('clustering_knn_results.csv')","metadata":{"execution":{"iopub.status.busy":"2021-06-17T16:14:37.581807Z","iopub.execute_input":"2021-06-17T16:14:37.582192Z","iopub.status.idle":"2021-06-17T16:18:22.826377Z","shell.execute_reply.started":"2021-06-17T16:14:37.58216Z","shell.execute_reply":"2021-06-17T16:18:22.825424Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We get an accuracy of 0.94421 from this model.\n\n## 2. PCA and K-NN\n\nHere, we try to speed up things in a different way than before. If in the previous method we reduced the size of the training matrix by reducing the number of observations, here we reduce it by performing SVD and keeping a number of components such that the explained variance is ~95%.","metadata":{"execution":{"iopub.status.busy":"2021-06-17T16:20:19.925114Z","iopub.execute_input":"2021-06-17T16:20:19.925491Z","iopub.status.idle":"2021-06-17T16:20:19.929673Z","shell.execute_reply.started":"2021-06-17T16:20:19.925459Z","shell.execute_reply":"2021-06-17T16:20:19.928486Z"}}},{"cell_type":"code","source":"y = df_train.pop('label').to_numpy()\nX = df_train.to_numpy()\n\nstandardizer_train = StandardScaler()\nX = standardizer_train.fit_transform(X)","metadata":{"execution":{"iopub.status.busy":"2021-06-18T07:55:07.993219Z","iopub.execute_input":"2021-06-18T07:55:07.99357Z","iopub.status.idle":"2021-06-18T07:55:08.918869Z","shell.execute_reply.started":"2021-06-18T07:55:07.993539Z","shell.execute_reply":"2021-06-18T07:55:08.917996Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### PCA\n\nLet's see how many component we should keep to have an explained variance of ~ 95%,","metadata":{}},{"cell_type":"code","source":"decomposition = PCA()\ndecomposition.fit(X)\n\nPCA_explained_variance = np.cumsum(decomposition.explained_variance_ratio_)\n\nplt.plot(PCA_explained_variance)\nplt.hlines(y=.95,\n           xmin=0,\n           xmax=df_train.columns.size,\n           colors='r',\n           linestyles='--'\n)\nplt.xlabel('# component in PCA')\nplt.ylabel('explained variance')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-06-18T07:55:13.031825Z","iopub.execute_input":"2021-06-18T07:55:13.032177Z","iopub.status.idle":"2021-06-18T07:55:17.998323Z","shell.execute_reply.started":"2021-06-18T07:55:13.032148Z","shell.execute_reply":"2021-06-18T07:55:17.997259Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can keep ~300. Then,","metadata":{}},{"cell_type":"code","source":"reduced_comp = 300\ndecomposition = PCA(n_components=reduced_comp)\n\nX_pca = decomposition.fit_transform(X)","metadata":{"execution":{"iopub.status.busy":"2021-06-18T07:55:17.999651Z","iopub.execute_input":"2021-06-18T07:55:17.99992Z","iopub.status.idle":"2021-06-18T07:55:25.395809Z","shell.execute_reply.started":"2021-06-18T07:55:17.999895Z","shell.execute_reply":"2021-06-18T07:55:25.394921Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's see what the first few principal components look like,","metadata":{}},{"cell_type":"code","source":"n_pca_show = 9\nsize = 28\n\nfirst_components = slice(n_pca_show)\nfig, axes = plt.subplots(n_pca_show//3,n_pca_show//3,figsize=(15,15))\n\nfor i, content in enumerate(zip(decomposition.components_[first_components],axes.flatten())):\n    pca_component, ax = content\n    ax.imshow(pca_component.reshape((size,size)))\n    ax.set_title('PCA component {}'.format(i+1))\n    ax.set_xticks([])\n    ax.set_yticks([])\n    \nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-06-18T08:58:00.612246Z","iopub.execute_input":"2021-06-18T08:58:00.61276Z","iopub.status.idle":"2021-06-18T08:58:01.281583Z","shell.execute_reply.started":"2021-06-18T08:58:00.612728Z","shell.execute_reply":"2021-06-18T08:58:01.280934Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### K-NN after PCA\n\nWe can now run K-NN on the reduced dataset, ","metadata":{}},{"cell_type":"code","source":"estimator = KNeighborsClassifier()\n\nest_params = {}\nest_params['n_neighbors'] = [1,2,3,4,5,6]\nest_params['weights'] = ['uniform','distance']\nest_params['p'] = [1,2]\n\ngrid = GridSearchCV(\n    estimator,\n    est_params,\n    scoring='accuracy',\n    n_jobs=-1,\n    refit=True,\n    cv=5,\n    verbose=2,\n    return_train_score=False\n)\n\ngrid.fit(X_pca,y)\n\ndf_results_pca = pd.DataFrame(grid.cv_results_)\ndf_results_pca.sort_values('rank_test_score')\n\ndf_results_pca.head()","metadata":{"execution":{"iopub.status.busy":"2021-06-17T16:50:54.966505Z","iopub.execute_input":"2021-06-17T16:50:54.96685Z","iopub.status.idle":"2021-06-17T16:57:53.403916Z","shell.execute_reply.started":"2021-06-17T16:50:54.966819Z","shell.execute_reply":"2021-06-17T16:57:53.400434Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = px.line(\n    df_results_pca,\n    x='param_n_neighbors',\n    y='mean_test_score',\n    error_y='std_test_score',\n    color='param_p',\n    facet_col='param_weights'\n)\n\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2021-06-17T16:57:53.405069Z","iopub.status.idle":"2021-06-17T16:57:53.405549Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's map the test set in the correct space,","metadata":{}},{"cell_type":"code","source":"X_pca_test = standardizer_train.transform(X_test)\nX_pca_test = decomposition.transform(X_pca_test)","metadata":{"execution":{"iopub.status.busy":"2021-06-18T08:00:45.801997Z","iopub.execute_input":"2021-06-18T08:00:45.802358Z","iopub.status.idle":"2021-06-18T08:00:46.247232Z","shell.execute_reply.started":"2021-06-18T08:00:45.80233Z","shell.execute_reply":"2021-06-18T08:00:46.246096Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"And use the optimal model to classify the digits,","metadata":{}},{"cell_type":"code","source":"opt_model = KNeighborsClassifier(\n    n_neighbors = 4,\n    weights = 'distance',\n    p = 2\n)\n\nopt_model.fit(X_pca,y)\n\ny_test = opt_model.predict(X_pca_test)\n\nsubmission['Label'] = y_test\nsubmission.to_csv('pca_knn_results.csv')","metadata":{"execution":{"iopub.status.busy":"2021-06-18T08:02:00.922396Z","iopub.execute_input":"2021-06-18T08:02:00.922904Z","iopub.status.idle":"2021-06-18T08:02:32.683149Z","shell.execute_reply.started":"2021-06-18T08:02:00.922869Z","shell.execute_reply":"2021-06-18T08:02:32.682189Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We get an accuracy of 0.95153 from this model!\n\n### 3. Syntetic data and tangent method\n\nWe could additionally create syntetic data to imporove the preformance of our classifier. To do so, we can implement a function that rotate, magnify, and sheer the images. Let's see how it looks,","metadata":{}},{"cell_type":"code","source":"size = 28\nimage = df_train.sample(1).to_numpy()\n\n# We can transform the image as we please, for example\nnew_image = image_transform(\n    image,\n    theta = -0.4,\n    mx = 0.8,\n    my = 1.2,\n    sheer = 0.2,\n    gamma = 3)\n\nfig,(ax1,ax2) = plt.subplots(1,2,figsize=(10,5))\n\nax1.imshow(image.reshape((size,size)),cmap='gray')\nax1.set_xticks([])\nax1.set_yticks([])\nax1.set_title('orginal image')\n\nax2.imshow(new_image.reshape((size,size)),cmap='gray')\nax2.set_xticks([])\nax2.set_yticks([])\nax2.set_title('transformed image')\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-06-21T14:19:39.441815Z","iopub.execute_input":"2021-06-21T14:19:39.442193Z","iopub.status.idle":"2021-06-21T14:19:39.620246Z","shell.execute_reply.started":"2021-06-21T14:19:39.442161Z","shell.execute_reply":"2021-06-21T14:19:39.61957Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In principle, we could now generate a pletora of new data, by perturbing the original images via rotations, magnification, and sheering. However, since we are already struggling with the number of observations we can use for out K-NN, maybe this is not such a great idea. Could be valuable for a neural network though!\n\n### Tangent method\n\nLet's try the tangent method for a single parameter (rotations).\n\nThe idea is that digits might be written with slightly different orientation, and we might account for these small-angle corrections by considering the distance between tangent lines generated by a first-order expansion of the rotation around the original image. In the following, we show what the tangent vector to the image wrt a rotation looks like, and how the rotated image compare to its first-order expansion.\n\nFor more information on this method, we refer to section 13.3.3. of 'The Element of Statistical Learning'.","metadata":{}},{"cell_type":"code","source":"# rotation parameter\ntheta = 0.15\n\n# New image obtained by rotating the input one\nrotated_image = image_transform(image,theta,mx=1,my=1,sheer=0,gamma=3)\n\n# New image obtained by first order expansion of the above rotation\nderivative_image = first_derivative_rotation(image,gamma=3)\n\ntangent_image = (image + theta * derivative_image).astype(np.int)\ntangent_image = rescale_image_values(image,tangent_image)\n\nfig, axs = plt.subplots(2,2,figsize=(12,12))\n(ax1,ax2,ax3,ax4) = axs.flatten()\n\nax1.imshow(image.reshape((size,size)),cmap='gray')\nax1.set_xticks([])\nax1.set_yticks([])\nax1.set_title('orginal image')\n\nax2.imshow(derivative_image.reshape((size,size)),cmap='gray')\nax2.set_xticks([])\nax2.set_yticks([])\nax2.set_title('derivative image')\n\nax3.imshow(rotated_image.reshape((size,size)),cmap='gray')\nax3.set_xticks([])\nax3.set_yticks([])\nax3.set_title('rotated image')\n\nax4.imshow(tangent_image.reshape((size,size)),cmap='gray')\nax4.set_xticks([])\nax4.set_yticks([])\nax4.set_title('tangent image')\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-06-21T15:06:17.687418Z","iopub.execute_input":"2021-06-21T15:06:17.687848Z","iopub.status.idle":"2021-06-21T15:06:18.051493Z","shell.execute_reply.started":"2021-06-21T15:06:17.687814Z","shell.execute_reply":"2021-06-21T15:06:18.050051Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's apply the method. We first compute the derivative of each image when rotated around its centre, for both training set and test set,","metadata":{}},{"cell_type":"code","source":"X_tg_centroid_train = format_data_tangent_method(X_cent)\nX_tg_centroid_test = format_data_tangent_method(X_centriod_test)","metadata":{"execution":{"iopub.status.busy":"2021-06-21T15:29:50.190666Z","iopub.execute_input":"2021-06-21T15:29:50.191071Z","iopub.status.idle":"2021-06-21T15:29:56.694854Z","shell.execute_reply.started":"2021-06-21T15:29:50.191031Z","shell.execute_reply":"2021-06-21T15:29:56.693934Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now, we use K-NN using the 10k centrids as training set, and as distance the one between two tangent lines passing by the corresponding images and parametrised by the rotation angle theta.\n\nWe do not optimize over the hyper-parameters since the distance we are using is substatially more computationally demanding than the previous ones. We instead select hyper-parameters that worked well for the previous examples.","metadata":{}},{"cell_type":"code","source":"model = KNeighborsClassifier(\n    n_neighbors=5,\n    weights='distance',\n    metric=optimal_line_distance,\n    n_jobs=-1,\n)\n\n# It would take to long to run this together with all the other models, so we comment this out.\n\n# model.fit(X_tg_centroid_train,y_cent)\n# y_test = model.predict(X_tg_centroid_test)\n\n# submission['Label'] = y_test\n# submission.to_csv('tg_knn_results.csv')","metadata":{"execution":{"iopub.status.busy":"2021-06-21T15:42:13.445836Z","iopub.execute_input":"2021-06-21T15:42:13.446305Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We get an accuracy of 0.94600 from this model. Although this does not beat the best method used here, we still see an improvement wrt the same method where euclidean distance was used instead of tangent distance.","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}