<!--
Abstract:

CS 461
Introduction to Deep Learning
Lecture 11
-->

# CS 461 - Lecture 15

## Machine Learning Principles

### The Kernel Trick

Bernhard Firner

2025-10-27

---

## Suggested Reading

* Murphy (Machine Learning: A Probabilistic Perspective)
  * Intro to chapter 14.4
* Flach
  * Chapters 1.3 and 7.2

---

## Review: High Dimensional Data

* Images are high-dimensional, structured data
  * Classifying based on individual pixels discards structure
  * Creating hidden states for a Markov model sounds intractable
* So we turn to comparison methods

---

## Similarity Scoring

* `exemplar` techniques use similarity scoring for classification
* How can we formalize similarity?
  * $k(x, x') \geq 0$ and $k(x, x')$ is real valued
  * That function is a `kernel function`
    * It would be nice if $k(x, x') = k(x', x)$ for reasons we'll see later

---

## Linear Kernel

* The matching detector, $X_i^TX_j$, is a linear kernel
  * We normalize it, turning it into the `cosine similarity`
  * $cos\theta = \frac{x\cdot y}{||x||_2\cdot ||y||_2}$
  * We tried this with Digits last time
* `cosine similarity` *is* positionally dependent, unlike Euclidean Distance
  * e.g. translating one image over the other changes their similarity

---

## String Similarity

* Define similarity as the number of shared substrings
* The matched words are:
  * 'a', 2*'the', 'White', 'Whale'
* So $K(x_1,x_2)$, the similarity of the first sentence to the second, is 5
* The similarity is also symmetric, since $K(x_1,x_2) = K(x_1,x_2)$

---

## Make a String Classifier

* Moby Dick Vs Pride and Prejudice
* Only use words greater length 5 or greater

---

## Moby Dick

<pre>
Almost simultaneously, with a mighty volition of ungraduated, instantaneous
swiftness, the White Whale darted through the weltering sea.

In this attitude the White Whale now shook the slight cedar as a mildly cruel
cat her mouse.
</pre>

Becomes

<pre>
Almost simultaneously mighty volition ungraduated instantaneous
swiftness White Whale darted through weltering

attitude White Whale shook slight cedar mildly cruel mouse
</pre>

---

## Pride and Prejudice

<pre>
The evening altogether passed off pleasantly to the whole family.

Elizabeth blushed and blushed again with shame and vexation.
</pre>

Becomes

<pre>
evening altogether passed pleasantly whole family

Elizabeth blushed blushed again shame vexation
</pre>

---

## Code

```python
import argparse
import nltk
import numpy as np
import sys

def filter_words(text, min_length):
    sentences = nltk.sent_tokenize(text)
    min_words = []
    for sentence in sentences:
        words = nltk.word_tokenize(sentence)
        filtered_words = [word for word in words if len(word) >= min_length]
        min_words.append(set(filtered_words))
    return min_words

def sentence_similarity(a_words, b_words):
    return len(a_words.intersection(b_words))

if __name__ == "__main__":
    parser = argparse.ArgumentParser()
    parser.add_argument(
        "--class_a",
        required=True,
        help="Text file with class A data")
    parser.add_argument(
        "--class_b",
        required=True,
        help="Text file with class B data")
    parser.add_argument(
        "--test_strings",
        required=True,
        help="Text file with test strings")
    parser.add_argument(
        "--min_length",
        default=5,
        type=int,
        help="The minimum word length to use.")

args = parser.parse_args()

#########
    # Note: You need the nltk parser
    # Open a python terminal and run:
    # import nltk
    # nltk.download('punkt_tab')

try:
        with open(args.class_a, 'r') as file:
            text_a = file.read()
        with open(args.class_b, 'r') as file:
            text_b = file.read()
        with open(args.test_strings, 'r') as file:
            test_text = file.read()
    except IOError as e:
        print(f"Error: Could not open or read file at path.")
        sys.exit(1)
    except Exception as e:
        print(f"An unexpected error occurred during parsing")
        sys.exit(1)

# Parse all of the files into sentences
    sentences_a = filter_words(text_a, args.min_length)
    sentences_b = filter_words(text_b, args.min_length)
    test_sentences = filter_words(test_text, args.min_length)
    original_test_sentences = nltk.sent_tokenize(test_text)

for idx, t in enumerate(test_sentences):
        max_similarity = 0
        classified = None
        closest_sentence = None
        for a in sentences_a:
            sim = sentence_similarity(t, a)
            if sim > max_similarity:
                classified = 'A'
                closest_sentence = a
                max_similarity = sim
        for b in sentences_b:
            sim = sentence_similarity(t, b)
            if sim > max_similarity:
                classified = 'B'
                closest_sentence = b
                max_similarity = sim

if closest_sentence is None:
            print(f"unmatched: {original_test_sentences[idx]}")
        else:
            print(f"{classified}: {original_test_sentences[idx]}; closest sentence is {' '.join(closest_sentence)}")
```

---

## Not Great Results

* Only took 4 chapters from each book, then attempted to classify some sentences from the last chapter
* Failures:

<pre>
P&P: I thank you again and again, for not going to the Lakes.
MD: Your idea of the ponies is delightful.
</pre>

* The matching sentences tended to be large, overbearing sentences spanning many lines
  * Similar to digits problems!

---

## MNIST DIGITS

* Hand written numeric digits converted to 28x28 black and white images
  * Half drawn by U.S. Census Bureau employees, other half by high school students
* 60,000 training images
* 10,000 validation
* Equally distributed over classes 0-9

---

## Matching Detector

* Just take the cosine similarity of everything in the test set compared to everything in the training set
  * Sounds reasonable; 8's should line up better with other 8's than with 7's
* Just perform $X_i^TX_j$ for I in the training data and J in the testing data

---

## The Code

```python

import argparse
import gzip
import os
import numpy as np
from PIL import Image

def load_mnist_ubyte(image_path):
    """
    Loads MNIST images from the raw ubyte files.

Args:
        image_path (str): Path to the image file (e.g., 'train-images-idx3-ubyte.gz').

Returns:
        images (numpy.ndarray)
    """
    with gzip.open(image_path, 'rb') as f:
        # Read the header: magic number (4 bytes) + num images (4 bytes) +
        # num rows (4 bytes) + num cols (4 bytes) = 16 bytes.

# Read the entire file content into a buffer
        image_data = f.read()

# The image data starts at byte 16.
        images = np.frombuffer(image_data, dtype=np.uint8, offset=16)

# We need the dimensions to reshape. We can extract them from the header bytes,
        # which are big-endian ('>'). We use struct.unpack if we were being strict,
        # but here we'll assume the standard MNIST format and calculate the dimensions
        # for a clean numpy approach.

# The number of images is in the 5th to 8th byte (4 bytes)
        num_images = np.frombuffer(image_data, dtype='>i4', offset=4, count=1)[0]
        # Rows and columns are 28x28 for MNIST, stored in bytes 9-12 and 13-16.
        # num_rows = np.frombuffer(image_data, dtype='>i4', offset=8, count=1)[0]
        # num_cols = np.frombuffer(image_data, dtype='>i4', offset=12, count=1)[0]
        num_rows = 28
        num_cols = 28

# Reshape the 1D array into a 3D array (num_images, rows, columns)
        images = images.reshape(num_images, num_rows, num_cols)
    return images

def load_mnist_labels(label_path):
    """
    Loads MNIST labels from the raw ubyte files.

Args:
        label_path (str): Path to the label file (e.g., 'train-labels-idx1-ubyte.gz').

Returns:
        labels (numpy.ndarray)
    """
    with gzip.open(label_path, 'rb') as f:
        # Read the header: magic number (4 bytes) + num items (4 bytes) = 8 bytes.
        # Skip these 8 bytes.

# Read the entire file content into a buffer
        label_data = f.read()

# The label data starts at byte 8. The data type is unsigned byte ('B' or np.uint8).
        # Labels are a 1D vector.
        labels = np.frombuffer(label_data, dtype=np.uint8, offset=8)

return labels

args = parser.parse_args()

X_train = load_mnist_ubyte(args.train)/255
    X_test = load_mnist_ubyte(args.test)/255

for j in range(args.matches):
        max_similarity = 0
        idx = -1
        for i in range(60000):
            similarity = np.sum(X_train[i].T@X_test[j])/(np.linalg.norm(X_train[i]) * np.linalg.norm(X_test[j]))
            if similarity > max_similarity:
                max_similarity = similarity
                idx = i

print(f"{j} similarity {max_similarity}")
        Image.fromarray((255*X_train[idx]).astype(np.uint8)).save(f"mnist_match_{j}_train.png")
        Image.fromarray((255*X_test[j]).astype(np.uint8)).save(f"mnist_match_{j}_test.png")
```

---

## The results

* Works well with typeset fonts, but not with human writing
* Small rotations of a letter make it difficult to match, thick digits still match too much
* So we want a more robust approach

</div>
<div class="col">

</div>
</div>

---

## Trainable Kernel Methods

* One problem with similarity on its own is that we cannot tell which features are more important than others
* The other problem is that "large" samples are too similar to everything else

---

## Types of Kernels

* That was a `bag of words` or `bag of features` kernel
  * If we restricted comparisons to be words of length $k$, that is a $k$-spectrum kernel
* Not useful yet
* Let's see how we can train something better using kernels

---

## Perceptron

* Recall that the perception was a line that estimates a class using the sign function
  * $\hat{y} = sign(bias + \sum weight_i x$)
* Every time the perceptron is wrong, we adjust the bias and weights in the "correct" direction by a learning rate, $\lambda$
  * $weight_i \leftarrow weight_i + \lambda y x_i$
* Predictions are
  * $\hat{y} = sign(bias + WX)$

---

## Perceptron Example

* $Y = {1, -1, 1, -1}$
* $X = {(1,2), (1,1), (2,4), (3,2)}$

---

## Perception Example

* Initialize perceptron parameters

<p><math display="block" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mtable><mtr><mtd columnalign="center" style="text-align: center"><mi>b</mi><mi>i</mi><mi>a</mi><mi>s</mi><mo>←</mo><mn>0</mn></mtd></mtr><mtr><mtd columnalign="center" style="text-align: center"><msub><mi>w</mi><mn>1</mn></msub><mo>←</mo><mn>0</mn></mtd></mtr><mtr><mtd columnalign="center" style="text-align: center"><msub><mi>w</mi><mn>2</mn></msub><mo>←</mo><mn>0</mn></mtd></mtr><mtr><mtd columnalign="center" style="text-align: center"><mi>λ</mi><mo>←</mo><mn>0.25</mn></mtd></mtr><mtr><mtd columnalign="center" style="text-align: center"></mtd></mtr></mtable><annotation encoding="application/x-tex">\begin{gather*}
bias \leftarrow 0\\
w_1 \leftarrow 0\\
w_2 \leftarrow 0\\
\lambda \leftarrow 0.25\\
\end{gather*}</annotation></semantics></math></p>

---

## Perceptron Example

* Iterate, updating when incorrect

<p><math display="block" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mtable><mtr><mtd columnalign="center" style="text-align: center"><mi>e</mi><mi>r</mi><mi>r</mi><mi>o</mi><mi>r</mi><mi>s</mi><mo>←</mo><mi>∞</mi></mtd></mtr><mtr><mtd columnalign="center" style="text-align: center"><mi>w</mi><mi>h</mi><mi>i</mi><mi>l</mi><mi>e</mi><mrow><mo stretchy="true" form="prefix">(</mo><mi>e</mi><mi>r</mi><mi>r</mi><mi>o</mi><mi>r</mi><mi>s</mi><mo>></mo><mn>0</mn><mo stretchy="true" form="postfix">)</mo></mrow><mo>:</mo></mtd></mtr><mtr><mtd columnalign="center" style="text-align: center"><mspace width="2.0em"></mspace><mi>f</mi><mi>o</mi><mi>r</mi><mrow><mo stretchy="true" form="prefix">(</mo><mi>i</mi><mo>∈</mo><mrow><mo stretchy="true" form="prefix">|</mo><mi>Y</mi><mo stretchy="true" form="postfix">|</mo></mrow><mo stretchy="true" form="postfix">)</mo></mrow><mo>:</mo></mtd></mtr><mtr><mtd columnalign="center" style="text-align: center"><mspace width="2.0em"></mspace><mspace width="2.0em"></mspace><msub><mover><mi>y</mi><mo accent="true">̂</mo></mover><mi>i</mi></msub><mo>=</mo><mn>0</mn></mtd></mtr><mtr><mtd columnalign="center" style="text-align: center"><mspace width="2.0em"></mspace><mspace width="2.0em"></mspace><mi>i</mi><mi>f</mi><mrow><mo stretchy="true" form="prefix">(</mo><msub><mover><mi>y</mi><mo accent="true">̂</mo></mover><mi>i</mi></msub><mo>≠</mo><msub><mi>y</mi><mi>i</mi></msub><mo stretchy="true" form="postfix">)</mo></mrow><mo>:</mo></mtd></mtr><mtr><mtd columnalign="center" style="text-align: center"><mspace width="2.0em"></mspace><mspace width="2.0em"></mspace><mspace width="2.0em"></mspace><mi>b</mi><mi>i</mi><mi>a</mi><mi>s</mi><mo>←</mo><mi>b</mi><mi>i</mi><mi>a</mi><mi>s</mi><mo>+</mo><mi>λ</mi><msub><mi>y</mi><mi>i</mi></msub></mtd></mtr><mtr><mtd columnalign="center" style="text-align: center"><mspace width="2.0em"></mspace><mspace width="2.0em"></mspace><mspace width="2.0em"></mspace><msub><mi>w</mi><mn>1</mn></msub><mo>←</mo><msub><mi>w</mi><mn>1</mn></msub><mo>+</mo><mi>λ</mi><msub><mi>y</mi><mi>i</mi></msub><msub><mi>x</mi><mrow><mi>i</mi><mo>,</mo><mn>1</mn></mrow></msub></mtd></mtr><mtr><mtd columnalign="center" style="text-align: center"><mspace width="2.0em"></mspace><mspace width="2.0em"></mspace><mspace width="2.0em"></mspace><msub><mi>w</mi><mn>2</mn></msub><mo>←</mo><msub><mi>w</mi><mn>2</mn></msub><mo>+</mo><mi>λ</mi><msub><mi>y</mi><mi>i</mi></msub><msub><mi>x</mi><mrow><mi>i</mi><mo>,</mo><mn>2</mn></mrow></msub></mtd></mtr><mtr><mtd columnalign="center" style="text-align: center"></mtd></mtr></mtable><annotation encoding="application/x-tex">\begin{gather*}
errors \gets \infty\\
while(errors > 0):\\
\qquad for(i \in |Y|):\\
\qquad\qquad \hat{y}_i = 0\\
\qquad\qquad if(\hat{y}_i \neq y_i):\\
\qquad\qquad\qquad bias \gets bias + \lambda y_i\\
\qquad\qquad\qquad w_1 \gets w_1 + \lambda y_i x_{i,1}\\
\qquad\qquad\qquad w_2 \gets w_2 + \lambda y_i x_{i,2}\\
\end{gather*}</annotation></semantics></math></p>

---

## Perceptron Iteration 1

* At $i=1$ (point 1,2, class 1)
  * $bias, w_1, w_2 \gets 0.25, 0.25, 0.5$
* At $i=2$, (point 1,1, class -1)
  * $bias, w_1, w_2 \gets 0, 0, 0.25$
* At $i=3$, (point 2,4, class 1) no error
* At $i=4$, (point 3,2, class -1)
  * $bias, w_1, w_2 \gets -0.25, -0.75, -0.25$

---

## Perceptron Iteration 2

* $bias, w_1, w_2 = -0.25, -0.75, -0.25$
* At $i=1$ (point 1,2, class 1)
  * $bias, w_1, w_2 \gets 0, -0.5, 0.25$
* At $i=2$, (point 1,1, class -1) (no error)
* At $i=3$, (point 2,4, class 1)
  * $bias, w_1, w_2 \gets 0.25, 0, 1.25$
* At $i=4$, (point 3,2, class -1)
  * $bias, w_1, w_2 \gets 0, -0.75, 0.75$

---

## Perceptron Iteration 3

* $bias, w_1, w_2 = 0, -0.75, 0.75$
* At $i=1$ (point 1,2, class 1) no error
* At $i=2$, (point 1,1, class -1)
  * $bias, w_1, w_2 \gets -0.25, -1, 0.5$
* At $i=3$, (point 2,4, class 1)
  * $bias, w_1, w_2 \gets 0, -0.5, 1.5$
* At $i=4$, (point 3,2, class -1)
  * $bias, w_1, w_2 \gets -0.25, -1.25, 1.0$

---

## Perceptron Iteration 4

* $bias, w_1, w_2 = -0.25, -1.25, 1.0$
* At $i=1$ (point 1,2, class 1) no error
* At $i=2$, (point 1,1, class -1) no error
* At $i=3$, (point 2,4, class 1) no error
* At $i=4$, (point 3,2, class -1) no error

---

## Perceptron Result

---

## What is a perceptron?

* $bias, w_1, w_2 = -0.25, -1.25, 1.0$
  * What is that?
  * It's a linear combination of the training inputs
    * Weighted by the number of times each point was wrong, their class, and $\lambda$
* $B = \lambda(3 - 4)$
* $W = \lambda(2X_1 - 2X_2 + 2X_3 - 3X_4) = (-5\lambda, 4\lambda)$

---

## Instance weights, not features

* We can think of the perceptron as a linear combination of X
  * $\alpha_i$ is the weight of the $X_i$
  * $\alpha = [2, 2, 2, 3]$ in our example
* This is called the `dual form`
* The bias becomes a prior, useful when we cannot normalize inputs (as in strings)

---

## New Perceptron Prediction

* Given input, t, and training set, X, with classes y
* $\hat{y} = sign(bias + \sum_{i=1}^n \alpha_i y_i X_i \cdot t)$
* What is $X_i \cdot t$?
  * It's a kernel, but a linear one
  * Could we use other kernels?
    * Yes! And some are fantastic!

---

## Cool Things about Kernels

* Imagine two roughly circularly grouped classes 
  * They are not linearly separable

---

## Making things separable

* If we look at $(x^2, y^2)$ instead, things are linearly separable

---

## Transformed Space

* To solve with linear regression or a perceptron, we would need to preprocess the data
* Is there any way to represent out weight vector, $w$, in the original space?
  * Yes! With a kernel!
* Let's use $k(a, b) = (a \cdot b)^2 = (a_x b_x + a_y b_y)^2$

---

## Linear Classifier

* A linear classifier will draw a line between the means of class 1 and class 2
  * $sign(w*x) = class$, where $w = \mu_1 - \mu_2$
* In x,y space we can see that the mean of class 1 and 2 are the same
  * So they aren't separable
* How about in transformed space?

---

## Choosing some Exemplars

* The mean of class 1 is $n = (0, 0)$, so let's have that represent the class
* Represent class 2 with four points, $P = (-1,-1), (-1,1), (1,1), (1,-1)$
* The mean of class 2 is: $\frac{1}{4}(p_1 + p_2 + p_3 + p_4)$, which is the same
* $w = \mu_1 - \mu_2$ hasn't changed

---

## Rewriting W

* Replace $w$ with the exemplars
* $\frac{1}{4}p_1 \cdot x + \frac{1}{4}p_2 \cdot x +\frac{1}{4}p_3 \cdot x +\frac{1}{4}p_4 \cdot x - n\cdot x$
* Now use the kernel
* $\frac{1}{4}k(p_1,x) + \frac{1}{4}k(p_2,x) +\frac{1}{4}k(p_3,x) +\frac{1}{4}k(p_4,x) - k(n,x) = class$

---

## Expanding

* $n = (0, 0)$
* $P = (-1,-1), (-1,1), (1,1), (1,-1)$
* $\frac{1}{4}k(p_1,x) + \frac{1}{4}k(p_2,x) +\frac{1}{4}k(p_3,x) +\frac{1}{4}k(p_4,x) - k(n,x)$
* $\frac{1}{4}[(-x-y)^2 + (-x+y)^2 + (x+y)^2 + (x-y)^2] - 0$
* $\frac{1}{4}[(x^2+y^2+2xy) + (x^2+y^2-2xy) + (x^2+y^2+2xy) + (x^2+y^2-2xy)]$
* $x^2+y^2$

---

## Result

* The decision boundary is the unit circle, $x^2+y^2$
  * Our regression parameters would be: $bias=-1, w=(1,1)$
  * $sign(bias + k(w,x))$ will correctly classify $x$
* This `kernel trick` allows a linear classifier, such as our perceptron, to have a curved decision boundary
  * Really, we are projecting from a linear space onto a curved space

---

## Kernel Trick + Perceptron

* Given input, t, and training set, X, with classes y
* $\hat{y} = sign(bias + \sum_{i=1}^n \alpha_i y_i k(X_i, t))$

---

## Gram Matrix

* During training, we do not need to recalculate the kernels repeatedly
* Instead, populate a matrix, `G`, where $G_{ij} = k(x_i, x_j)$
  * G is NxN, where N is the number of training samples
  * For sample $i$, the class prediction is 
    * $\hat{y} = sign(bias + \alpha y G(N, i))$
    * $G(N, i)$ is the i-th column of the matrix

---

## Kernel Perceptron

```python
"""
A perceptron class.
It will accept a kernel function to use when processing inputs.
"""

import numpy as np

class KernelPerceptron:
    """
    Perceptron Classifier using the Kernel Trick.
    """
    def __init__(self, kernel_func):
        """
        Arguments:
            kernel_func (function): The kernel function K(x_i, x_j) to use.
        """
        self.kernel_func = kernel_func
        self.alphas = None      # Dual coefficients (alpha_i)
        self.bias = 0           # Bias term (b)
        self.X_train = None     # Store training data (required for prediction)
        self.y_train = None     # Store training labels (required for prediction)

def fit(self, X, y, max_epochs, random_seed=None):
        """
        Trains the Kernel Perceptron.

Arguments:
            X (numpy.ndarray): Training data matrix (N_samples, N_features).
            y (numpy.ndarray): Target labels, either -1 or 1
            max_epochs  (int): Number of epochs, unless convergence is faster
            random_see  (int): Seed for shuffling the data during each epoch.
        """
        if not np.all(np.in1d(y, [-1, 1])):
            raise ValueError("Labels must be in {-1, 1}.")

N = X.shape[0]
        # We need to remember this for classification and for training
        self.X_train = X
        self.y_train = y
        # These are the counts of how wrong we are
        self.alphas = np.zeros(N)
        self.bias = 0

# Calculate the Gram Matrix K, which holds K(x_i, x_j) for all pairs
        # This pre-calculation speeds up the training process.
        # It also means that we no longer require X
        self.gram_matrix = self.precalculate_gram_matrix(X, X)

# Randomize the order of presentations to maybe avoid loops
        rng = np.random.default_rng(random_seed)
        indices = np.arange(N)

for epoch in range(self.max_epochs):
            n_updates = 0
            # Randomize the presentation order with each epoch
            rng.shuffle(indices)

for i in indices:
                # 1. Prediction using the dual form and Gram matrix:
                # Decision Function: f(x_i) = sum(alpha_j * y_j * K(x_j, x_i)) + b
                # In matrix form: alpha * y * K[:, i] + b
                # The term K[:, i] is the i-th column of the Gram matrix
                decision_value = np.sum(self.alphas * self.y_train * self.gram_matrix[:, i]) + self.bias
                
                prediction = np.sign(decision_value)
                
                # We could predict 0 (for example, on the first step)
                # Count it as an error
                if prediction == 0:
                    prediction = -1 * self.y_train[i]

# 2. Update rule:
                if prediction != self.y_train[i]:
                    # Increment the mistakes made on this exemplar
                    self.alphas[i] += 1
                    # Also increment the bias
                    self.bias += self.y_train[i]
                    n_updates += 1
            print(f"Epoch {epoch} had {n_updates}/N updates")
            
            # Stop if the perceptron has converged
            if n_updates == 0:
                print(f"Converged at epoch {epoch + 1}.")
                break
        else:
            print(f"Finished training after max_epochs ({self.max_epochs}).")

def precalculate_gram_matrix(self, X1, X2):
        """Calculates the matrix of pairwise kernel evaluations K(X1[i], X2[j])."""
        n1 = X1.shape[0]
        n2 = X2.shape[0]
        K = np.zeros((n1, n2))
        
        for i in range(n1):
            for j in range(n2):
                K[i, j] = self.kernel_func(X1[i], X2[j])
        return K

def predict(self, X):
        """
        Predicts labels for the test data X.

Arguments:
            X (numpy.ndarray): Test data matrix (M_samples, N_features).
        Returns:
            Predicted labels (list[M_samples])
        """
        # Calculate the kernel values between training and test data: K(X_train, X_test)
        kernel_test = self.precalculate_gram_matrix(self.X_train, X)
        
        # Decision function: f(x) = sum(alpha_i * y_i * K(x_i, x)) + b
        # In matrix form: (alpha * y) @ K.T + b
        weighted_alphas = self.alphas * self.y_train
        decision_values = weighted_alphas @ kernel_test + self.bias

# Return the sign of the decision function
        return np.sign(decision_values)
```

---

## Problems

* It's not fast
  * One could even say it's incredibly slow
  * N operations per prediction, N predictions per epoch
  * Each operation is also scaled by the number of features
  * n=60k, features=784
  * Into the $10^{12}$ range per epoch

---

## More problems

* Only guessing level accuracy
  * May not be converging with current features
  * Meaning it's not linearly separable with the square kernel
* Even if it converges, the decision boundary could be "just barely" good enough for the current data
  * e.g. no margin between the boundary and the training points

---

## Some Improvements Possible

* We could be computing things on the GPU
* And we could keep a moving average of the weight to improve stability
* And we could use a different kernel

---

## Speedup on GPU

* We can just throw everything onto a GPU
  * Can only train with half the samples using fp16 on 8GB memory
  * Fine

---

## GPU Version

```python
import torch
import numpy as np

# --- 1. Device Configuration ---
# Set the device to CUDA/GPU if available, otherwise use CPU.
# This ensures all tensors are created on the correct hardware.
DEVICE = torch.device("cuda" if torch.cuda.is_available() else "cpu")
print(f"Using device: {DEVICE}")

class KernelPerceptron():
    """
    A Kernel Perceptron implementation using PyTorch for GPU acceleration.

The model parameters (alpha) and data are stored on the selected device.
    Uses the Gaussian (RBF) Kernel for non-linear separation.
    """
    def __init__(self, gamma: float = 0.1, degree: int = 2, offset: float = 0.0):
        """
        Initializes the Kernel Perceptron.

Arguments:
            gamma: The kernel parameter (often 1 / (2*sigma^2)) for the RBF kernel.
            degree: The power (d) of the polynomial kernel.
            offset: The constant offset (c) of the polynomial kernel.
        """
        super().__init__()
        self.gamma = gamma
        self.degree = degree
        self.offset = offset
        self.alpha = None      # The weight vector (dual variables)
        self.X_train = None    # Training data (support vectors)
        self.y_train = None    # Training labels
        self.kernel_matrix = self._poly_kernel_matrix
        #self.kernel_matrix = self._rbf_kernel_matrix
        self.learning_rate = 0.01

def _poly_kernel_matrix(self, X1, X2):
        """
        Calculates the Polynomial Kernel matrix K(X1, X2) using vectorized
        PyTorch operations.

K_ij = (x1_i.T @ x2_j + offset)^degree
        """
        # 1. Calculate the dot product matrix (X1 @ X2.T)
        # This computes x1_i.T @ x2_j for all pairs
        dot_product = torch.mm(X1, X2.T)

# 2. Apply the polynomial formula element-wise
        # K = (dot_product + offset)^degree
        K = (dot_product + self.offset).pow(self.degree)

return K

def _rbf_kernel_matrix(self, X1, X2):
        """
        Calculates the Gaussian (RBF) Kernel matrix K(X1, X2) using vectorized
        PyTorch operations.

K_ij = exp(-gamma * ||x1_i - x2_j||^2)
        """
        # Calculate squared Euclidean distances (||x1_i - x2_j||^2) using
        # the identity: ||a - b||^2 = ||a||^2 - 2*a.T*b + ||b||^2

# 1. Calculate ||X1||^2 and ||X2||^2 (sum of squares for each row)
        X1_sq_norm = torch.sum(X1**2, dim=1, keepdim=True)
        X2_sq_norm = torch.sum(X2**2, dim=1, keepdim=True)
        
        # 2. Calculate -2 * X1 @ X2.T (Dot product)
        dot_product = torch.mm(X1, X2.T)

# 3. Calculate the squared distance matrix (D^2)
        # Note: X1_sq_norm is N x 1, X2_sq_norm is M x 1.
        # X1_sq_norm + X2_sq_norm.T results in an N x M matrix due to broadcasting.
        sq_distances = X1_sq_norm - 2 * dot_product + X2_sq_norm.T
        
        # 4. Apply the RBF formula
        K = torch.exp(-self.gamma * sq_distances)
        return K

def fit(self, X: np.ndarray, y: np.ndarray, max_epochs: int = 10):
        """
        Trains the Kernel Perceptron using the iterative Perceptron rule.
        
        Arguments:
            X (numpy.ndarray): Training data matrix (N_samples, N_features).
            y (numpy.ndarray): Target labels, either -1 or 1. If 0 and 1, 0 will be converted to -1.
            max_epochs  (int): Number of epochs, unless convergence is faster
        """
        # --- Data Preparation ---
        # Convert numpy arrays to PyTorch tensors and move to the device
        X_tensor = torch.tensor(X, dtype=torch.float16, device=DEVICE)
        # Ensure labels are +1 or -1 (Perceptron standard)
        y_tensor = torch.tensor(y, dtype=torch.float16, device=DEVICE)
        y_tensor[y_tensor == 0] = -1

N_samples = X_tensor.shape[0]
        
        self.X_train = X_tensor
        self.y_train = y_tensor
        
        # Initialize dual variables (alpha) on the device
        self.alpha = torch.zeros(N_samples, dtype=torch.float16, device=DEVICE)
        self.bias = 0

# --- Gram Matrix Calculation ---
        # Calculate the kernel matrix K = K(X_train, X_train) once
        print(f"Calculating {N_samples}x{N_samples} Gram Matrix on {DEVICE}...")
        K = self.kernel_matrix(self.X_train, self.X_train)
        print("Gram Matrix calculation complete. Starting training...")

# --- Training Loop ---
        for epoch in range(max_epochs):
            n_misclassified = 0
            
            # The decision function is f(x_i) = sum_j(alpha_j * y_j * K_ji)
            # This is computed efficiently as: f = K @ (alpha * y)
            f = K @ (self.alpha * self.y_train)

# Perceptron update rule check: y_i * f(x_i) <= 0
            misclassified_mask = (self.y_train * f <= 0)
            
            # Update alpha.
            # We want to add y_i to alpha_i only where misclassified_mask is True
            self.alpha[misclassified_mask] += 1
            
            n_misclassified = torch.sum(misclassified_mask).item()
            
            # Break condition
            if n_misclassified == 0:
                print(f"Epoch {epoch+1}/{max_epochs}: Converged!")
                break
            
            print(f"Epoch {epoch+1}/{max_epochs}: Misclassified {n_misclassified}/{N_samples}")

def predict(self, X: np.ndarray) -> np.ndarray:
        """
        Predicts labels for a new dataset X.

Arguments:
            X: New data features (N_test_samples, N_features).

Returns:
            Predicted labels (+1 or -1).
        """
        if self.alpha is None:
            raise RuntimeError("Model has not been trained. Call .fit() first.")

# Convert test data to tensor and move to the device
        X_test_tensor = torch.tensor(X, dtype=torch.float16, device=DEVICE)

# 1. Calculate the test kernel matrix K_test = K(X_test, X_train)
        # K_test has shape (N_test, N_train)
        K_test = self.kernel_matrix(X_test_tensor, self.X_train)

# 2. Calculate the decision function f(x)
        # f(x) = K_test @ (alpha * y_train)
        # K_test (N_test x N_train) * (alpha * y_train) (N_train x 1) -> N_test x 1
        decision_values = K_test @ (self.alpha * self.y_train)

# 3. Apply sign function: sign(f(x))
        predictions = torch.sign(decision_values)
        
        # Convert back to numpy array and ensure 0 is returned for -1 for consistency
        # if the user expects binary (0, 1) labels
        predictions_np = predictions.cpu().numpy()
        predictions_np[predictions_np == -1] = 0
        
        return predictions_np
```

---

## Training Code

```python
import argparse
import gzip
import os
import numpy as np
from PIL import Image
#import perceptron_funcs
import gpu_perceptron_funcs as perceptron_funcs