SML_Assignment4

Task 1: Neural Networks

1a) Multi-layer Perceptron

We use mean cross entropy loss, because it is good for multi-category problems;
We use ReLu as activation function due to excellent adaptability to multi-class classification. We have also tested sigmoid, but there are problems with vanishing gradients.
The input layer has 784 neurons, because the count of pixel is (28 imes 28);
The output layer has 10 neurons, because we need to classify 10 categories.
The selection strategy of learning rate is constantly changing in the process of network training. At the beginning, the parameters are relatively random, so we should choose a relatively large learning rate, so that loss will fall faster. When we train for a period of time, the parameter update should have a smaller range of change, and then we choose a smaller learning rate

Full code is below:

import numpy as np
from matplotlib import pyplot as plt


def initialize(hidden_dim, output_dim):
    # retrieve mnist data
    X_train = np.loadtxt('./dataSets/mnist_small_train_in.txt', delimiter=',')
    X_test = np.loadtxt('./dataSets/mnist_small_test_in.txt', delimiter=',')
    y_train = np.loadtxt('./dataSets/mnist_small_train_out.txt', delimiter=',', dtype="int32")
    y_test = np.loadtxt('./dataSets/mnist_small_test_out.txt', delimiter=',', dtype="int32")

    # normalize data
    X_train = (X_train - np.mean(X_train)) / np.std(X_train)
    X_test = (X_test - np.mean(X_test)) / np.std(X_test)

    # one hot encode labels
    temp_y_train = np.eye(10)[y_train]
    temp_y_test = np.eye(10)[y_test]
    y_train = temp_y_train.reshape(list(y_train.shape) + [10])
    y_test = temp_y_test.reshape(list(y_test.shape) + [10])

    # weights for single hidden layer
    W1 = np.random.randn(hidden_dim, X_train.shape[1]) * 0.01
    b1 = np.zeros((hidden_dim,))
    W2 = np.random.randn(output_dim, hidden_dim) * 0.01
    b2 = np.zeros((output_dim,))

    parameters = [W1, b1, W2, b2]

    # return to column vectors
    return parameters, X_train.T, X_test.T, y_train.T, y_test.T


def sigmoid(x):
    return 1.0 / (1.0 + np.exp(-x))


def dsigmoid(x):
    # input x is already sigmoid, no need to recompute
    return x * (1.0 - x)


def ReLU(x):
    return x * (x > 0)


def dReLU(x):
    return 1. * (x > 0)

# Cross entropy loss
def closs(pred, y):
    return np.squeeze(-np.sum(np.multiply(np.log(pred), y)) / len(y))

def dcloss(pred, y):
    return pred - y

# Mean squared error loss
def loss(pred, y):
    return np.sum(.5 * np.sum((pred - y) ** 2, axis=0), axis=0) / y.shape[1]

def dloss(pred, y):
    return (pred - y) / y.shape[1]


class NeuralNet(object):

    def __init__(self, hidden_dim, output_dim):
        # size of layers
        self.hidden_dim_1 = hidden_dim
        self.output_dim = output_dim

        # weights and data
        parameters, self.x_train, self.x_test, self.y_train, self.y_test = initialize(hidden_dim, output_dim)
        self.W1, self.b1, self.W2, self.b2 = parameters

        # activations
        batch_size = self.x_train.shape[0]

        self.ai = np.ones((self.x_train.shape[1], batch_size))
        self.ah1 = np.ones((self.hidden_dim_1, batch_size))
        self.ao = np.ones((self.output_dim, batch_size))

        # classification output for transformed OHE
        self.classification = np.ones(self.ao.shape)

        # container for loss progress
        self.loss = None
        self.test_error = None

    def forward_pass(self, x):
        # input activations
        self.ai = x

        # hidden_1 activations
        self.ah1 = sigmoid(self.W1 @ self.ai + self.b1[:, np.newaxis])

        # output activations
        self.ao = sigmoid(self.W2 @ self.ah1 + self.b2[:, np.newaxis])

        # transform to OHE for classification
        self.classification = (self.ao == self.ao.max(axis=0, keepdims=0)).astype(float)

    def backward_pass(self, target):
        # calculate error for output
        out_error = dloss(self.ao, target)
        out_delta = out_error * dsigmoid(self.ao)

        # calculate error for hidden_1
        hidden_1_error = self.W2.T @ out_delta
        hidden_1_delta = hidden_1_error * dReLU(self.ah1)

        # derivative for W2/b2 (hidden_1 --> out)
        w2_deriv = out_delta @ self.ah1.T
        b2_deriv = np.sum(out_delta, axis=1)

        # derivative for W1/b1 (input --> hidden_1)
        w1_deriv = hidden_1_delta @ self.ai.T
        b1_deriv = np.sum(hidden_1_delta, axis=1)

        return [w1_deriv, b1_deriv, w2_deriv, b2_deriv]

    def train(self, epochs, lr, batch_size=128):

        self.loss = np.zeros((epochs,))
        self.test_error = np.zeros((epochs,))

        for epoch in range(epochs):

            indices = np.arange(self.x_train.shape[1])
            np.random.shuffle(indices)

            for i in range(0, self.x_train.shape[1] - batch_size + 1, batch_size):
                excerpt = indices[i:i + batch_size]

                batch_x, batch_y = self.x_train[:, excerpt], self.y_train[:, excerpt]

                # compute output of forward pass
                self.forward_pass(batch_x)

                # back prop error
                w1_deriv, b1_deriv, w2_deriv, b2_deriv = self.backward_pass(batch_y)

                # adjust weights with simple SGD
                self.W1 -= lr * w1_deriv
                self.b1 -= lr * b1_deriv
                self.W2 -= lr * w2_deriv
                self.b2 -= lr * b2_deriv

            # compute error
            self.forward_pass(self.x_test)
            error = loss(self.ao, self.y_test)
            self.loss[epoch] = error

            t = 0
            for i, pred in enumerate(self.classification.T):
                if np.argmax(self.y_test.T[i]) == np.argmax(pred):
                    t += 1

            acc = t / self.classification.shape[1]
            self.test_error[epoch] = 1 - acc

            # print progress
            if (epoch + 1) % 50 == 0:
                print("Epoch {}/{} -> loss: {:.5f} -> val_acc: {:.5f}".format(epoch + 1, epochs, error, acc))

    def predict(self, x):
        self.forward_pass(x)
        return self.classification

    def visualize_prediction(self, expected, predicted):
        # plot loss progress
        plt.figure(2)
        plt.plot(self.loss)
        plt.xlabel("Epochs")
        plt.ylabel("Loss")
        plt.show()

        # plot loss progress
        plt.figure(3)
        plt.plot(self.test_error)
        plt.xlabel("Epochs")
        plt.ylabel("Test Error")
        plt.show()


nn = NeuralNet(hidden_dim=784, output_dim=10)
nn.train(epochs=500, lr=.05)

y_hat = nn.predict(nn.x_test)
t = 0
for i, pred in enumerate(y_hat.T):
    if np.argmax(nn.y_test.T[i]) == np.argmax(pred):
        t += 1

print("Accuracy:", t / y_hat.shape[1])
nn.visualize_prediction(nn.y_test.T, y_hat.T)

Show how the misclassification error (in percent) on the testing set evolves during learning:

Task 3: Gaussian Processes

3a) GP Regression

Full code is below:

import numpy as np
from matplotlib import pyplot as plt

def target(x):
    s = np.sin(x)
    return s + np.square(s)


def kernel(x, z):
    return np.exp(-(x - z) ** 2)


def compute_c(x, noise):
    return kernel(x, x) + noise


def predict(x, X, Y, C_inv, noise):
    k = np.array([kernel(x, val) for val in X])
    mean = k.T.dot(C_inv).dot(Y)
    std = np.sqrt(compute_c(x, noise) - k.T.dot(C_inv).dot(k))
    return mean, std


def plot_gp(X, mean, std, iteration):
    y1 = mean - 2 * std
    y2 = mean + 2 * std

    plt.plot(X, target(X), "-", color="red", label="$sin(x) + sin^2(x)$")
    plt.plot(X, mean, color="black", label="$mu$")
    plt.fill_between(X, y1, y2, facecolor='blue', interpolate=True, alpha=.5, label="$2sigma$")
    plt.title("$mu$ and $sigma$ for iteration={}".format(iteration))
    plt.xlabel("x")
    plt.ylabel("y")


def gpr():
    noise = 0.001
    step_size = 0.005
    X = np.arange(0, 2 * np.pi + step_size, step_size)
    iterations = 15

    std = np.array([np.sqrt(compute_c(x, noise)) for x in X])
    j = np.argmax(std)

    Xn = np.array([X[j]])
    Yn = np.array([target(Xn)])

    C = np.array(compute_c(X[j], noise)).reshape((1, 1))
    C_inv = np.linalg.solve(C, np.identity(C.shape[0]))

    for iteration in range(0, iterations):
        mean = np.zeros(X.shape[0])
        std = np.zeros(X.shape[0])

        for i, x in enumerate(X):
            mean[i], std[i] = predict(x, Xn, Yn, C_inv, noise)

        if iteration + 1 in [1, 2, 5, 10, 15]:
            plot_gp(X, mean, std, iteration + 1)
            plt.plot(Xn, Yn, "o", c="blue", label="sampled")
            plt.legend()
            plt.show()

        j = np.argmax(std)
        Xn = np.append(Xn, X[j])
        Yn = np.append(Yn, target(Xn[-1]))

        # update C matrix
        x_new = Xn[-1]
        k = np.array([kernel(x_new, val) for val in Xn])
        c = kernel(x_new, x_new) + noise

        dim = C.shape[0]
        C_new = np.empty((dim + 1, dim + 1))
        C_new[:-1, :-1] = C
        C_new[-1, -1] = c
        C_new[:, -1] = k
        C_new[-1:] = k.T

        C = C_new
        C_inv = np.linalg.solve(C, np.identity(C.shape[0]))

gpr()