NDArray and Neural networks

About This Module

Prework

Prework Reading

Pre-lecture Reflections

Lecture

Learning Objectives

• Simple neural networks can be easily represented by matrix arithmentic.

• Let's see how we can build a relatively simple neural network to recognize handwritten digits.

Biological Neuron

An _artificial neuron_is loosely modeled on a biological neuron.

From Stanford's cs231n

• The dendrites carry impulses from other neurons of different distances.
• Once the collective firing rate of the impulses exceed a certain threshold, the neuron fires its own pulse through the axon to other neurons.

Artificial Neuron

The artificial neuron is loosely based on the biological one.

From Stanford's cs231n.

The artifical neuron:

• collects one or more inputs,
• each multiplied by a unique weight,
• sums the weighted inputs,
• adds a bias,
• then applies a nonlinear activation function.

Multi-Layer Perceptron (MLP) or Fully Connected Network (FCN)

From Stanford's cs231n.

• Multiple artificial neurons can act on the same inputs,. This defines a layer. We can have more than one layer until we produce one or more outputs.

• The example above shows a network with 3 inputs, two layers of neurons, each with 4 neurons, followed by one layer that produces a single value output.

This architecture can be used as a binary classifier.

FCN Detailed View

Here's a more detailed view of a fully connected network making biases explicit and with weight matrix and bias vector dimensions

We can collect all the inputs into a column vector.

$$\mathbf{x}=\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]$$

Then we can gather the weights for the first layer:

$$h_1=\text{a}\left[{\mathbf{\Omega }}_0\mathbf{x}+{\mathbf{\beta }}_0\right]=\text{a}\left[\left[\begin{array}{ccc}{\omega }_{11}& {\omega }_{12}& {\omega }_{13}\\ {\omega }_{21}& {\omega }_{22}& {\omega }_{23}\\ {\omega }_{31}& {\omega }_{32}& {\omega }_{33}\\ {\omega }_{41}& {\omega }_{42}& {\omega }_{43}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]+\left[\begin{array}{c}{\beta }_1\\ {\beta }_2\\ {\beta }_3\\ {\beta }_4\end{array}\right]\right]$$

where $a[z]$ is a nonlinear activation function, such as ReLU shown below.

The equation for which is.

$$a[z]=\text{ReLU}[z]=\{\begin{array}{ll}0& z<0\\ z& z\ge 0\end{array}$$

or in pseudocode simply:

x = max(0, x)

The nonlinear activation functions are crucial, because without them, the neural network can be simplified down to a single linear system of equations, which is severely limited in its representational power.

The subsequent layers can be defined similarly interms of the activations of the previous activation functions.

$$\begin{array}{rl}{\mathbf{h}}_2& =\text{a}\left[{\mathbf{\Omega }}_1{\mathbf{h}}_1+\mathbf{\beta }\right]\\ {\mathbf{h}}_3& =\text{a}\left[{\mathbf{\Omega }}_2{\mathbf{h}}_2+\mathbf{\beta }\right]\\ {\mathbf{h}}_4& ={\mathbf{\Omega }}_3{\mathbf{h}}_3+\mathbf{\beta }\\ \mathbf{y}& =\text{softmax}\left[{\mathbf{h}}_4\right]\end{array}$$

where the last function is a softmax, which limits the output to between 0 and 1, and all the outputs sum up to 1.

$$\mathrm{softmax}(y_i)=\frac{e^{y_i}}{{\sum }_{j=1}^K{e}^{y_j}}$$

Here's a visualization.

From TDS Post.

Much of neural network training and evaluation comes down to vector operations in which ndarray excels.

Training a Neural Network

• We won't go into all the details of training a network, but we'll show an example.

• The steps involve producing the outpur from some input -- Forward Pass

• Then updating the weights and biases of the network, taking the partial derivatives of the entire model with respect to each parameter and then updating each parameter proportional to its partial derivative -- Backward Pass

Simple weight update example

We'll build our understanding of the gradient by considering derivatives of single variable functions.

Let's start with a quadratic function

$$f(x)=3x^2-4x+5.$$

The derivative with respect to $x$ of our example is

$$f^{\prime }(x)=6x-4.$$

This is the asymptotic slope at any particular value of $x$.

Question: Which way does $x$ need to move at $x=-2$ to decrease $f(x)$? Positive or negative direction?

Question: What about at $x=2$? Positive or negative direction?

An easy way to update this is:

$$x\leftarrow x-\alpha \frac{\partial f}{\partial x}$$

where $\alpha$ is some small number (e.g. 0.01) which we call the learning rate.

Minimizing the Network Loss

To train our network we want to minize some function that computes the difference between our target values $t_i$ and what the network produces $y_i$ for each input $x_i$.

An appropriate loss function for classification models is the categorical cross-entropy loss.

The following image from the Medium Post illustrates the full pipeline:

except in the figure $s$ is the model output and $y$ are the targets.

We won't go into details but a nice output is that

$$\frac{\partial {\ell }_i}{\partial {\mathbf{h}}_4}=\mathbf{y}-\mathbf{t}$$

Forward propagation

• Weights would be represented by an Array2<f32> with dimensions: (# neurons $\times$ # inputs)

  • For example if you have 100 input neurons and 50 2nd layer neurons then you need a $50\times 100$ matrix of weights

• Bias vectors would be represented by an Array2<f32> with dimensions (# neurons $\times$ 1)

Backwards propagation

• As mentioned, in practice the calculation of partial derivatives must be calculated starting with the end of the pipeline and then working backwards.

• The chain rule for differentiation in calculus gives a modular, scalable way of doing this calculation, again starting from the end and working backwards.

Just a reminder, we have

$$\begin{array}{rl}{\mathbf{f}}_0& ={\mathbf{\Omega }}_0\mathbf{x}+{\mathbf{\beta }}_0\\ {\mathbf{h}}_1& =\text{a}\left[{\mathbf{f}}_0\right]\\ {\mathbf{f}}_1& ={\mathbf{\Omega }}_1{\mathbf{h}}_1+\mathbf{\beta }\\ {\mathbf{h}}_2& =\text{a}\left[{\mathbf{f}}_1\right]\\ {\mathbf{f}}_2& ={\mathbf{\Omega }}_2{\mathbf{h}}_2+\mathbf{\beta }\\ {\mathbf{h}}_3& =\text{a}\left[{\mathbf{f}}_2\right]\\ {\mathbf{h}}_4& ={\mathbf{\Omega }}_3{\mathbf{h}}_3+\mathbf{\beta }\\ \mathbf{y}& =\text{softmax}\left[{\mathbf{h}}_4\right]\end{array}$$

We already did

$$\begin{array}{rl}\frac{\partial {\ell }_i}{\partial {\mathbf{h}}_4}& =\mathbf{y}-\mathbf{t}\\ \frac{\partial {\ell }_i}{\partial {\mathbf{\Omega }}_3}& =\frac{\partial {\ell }_i}{\partial {\mathbf{h}}_4}\frac{\partial {\mathbf{h}}_4}{\partial {\mathbf{\Omega }}_3}=\frac{\partial {\ell }_i}{\partial {\mathbf{h}}_4}{\mathbf{h}}_3^T\end{array}$$

And the process continues.

Example Implementaiton

Let's look at an example implementation.

We'll start with some utility functions.

#![allow(unused)]
fn main() {
:dep ndarray = { version = "^0.15" }
:dep rand = { version = "0.8" }

use ndarray::prelude::*;
use rand::Rng;

// Helper function to populate arrays with random values
// We could use ndarray-rand crate to do this, but it doesn't seem to work in a Jupyter notebook
fn populate_array(arr: &mut Array2<f32>, m: usize, n: usize) {
    let mut rng = rand::thread_rng();
    for i in 0..m {
        for j in 0..n {
            arr[(i, j)] = rng.gen_range(-1.0..1.0);
        }
    }
}

// ReLU activation function
fn relu(x: &Array2<f32>) -> Array2<f32> {
    x.mapv(|x| if x > 0.0 { x } else { 0.0 })
}

// Derivative of ReLU
fn relu_derivative(x: &Array2<f32>) -> Array2<f32> {
    x.mapv(|x| if x > 0.0 { 1.0 } else { 0.0 })
}

// Softmax function
fn softmax(x: &Array2<f32>) -> Array2<f32> {
    let exp_x = x.mapv(|x| x.exp());
    let sum_exp_x = exp_x.sum();
    exp_x / sum_exp_x
}

}

Then we'll define a NeuralNetwork struct and associated methods.

#![allow(unused)]

fn main() {
struct NeuralNetwork {
    input_size: usize,
    output_size: usize,
    weights1: Array2<f32>,
    biases1: Array2<f32>,
    weights2: Array2<f32>,
    biases2: Array2<f32>,
    learning_rate: f32,
}

impl NeuralNetwork {
    /// Create a shallow neural network with one hidden layer
    /// 
    /// # Arguments
    /// 
    /// * `input_size` - The number of input neurons
    /// * `hidden_size` - The number of neurons in the hidden layer
    /// * `output_size` - The number of output neurons
    /// * `learning_rate` - The learning rate for the neural network
    fn new(input_size: usize, hidden_size: usize, output_size: usize, learning_rate: f32) -> Self {
        let mut weights1 = Array2::zeros((hidden_size, input_size));
        let mut weights2 = Array2::zeros((output_size, hidden_size));

        // Initialize weights with random values
        populate_array(&mut weights1, hidden_size, input_size);
        populate_array(&mut weights2, output_size, hidden_size);

        let biases1 = Array2::zeros((hidden_size, 1));
        let biases2 = Array2::zeros((output_size, 1));

        NeuralNetwork {
            input_size,
            output_size,
            weights1,
            biases1,
            weights2,
            biases2,
            learning_rate,
        }
    }

    fn forward(&self, input: &Array2<f32>) -> (Array2<f32>, Array2<f32>, Array2<f32>) {
        // First layer
        let pre_activation1 = self.weights1.dot(input) + &self.biases1;
        let hidden = relu(&pre_activation1);

        // Output layer
        let pre_activation2 = self.weights2.dot(&hidden) + &self.biases2;
        let output = softmax(&pre_activation2);

        (hidden, pre_activation2, output)
    }

    fn backward(
        &mut self,
        input: &Array2<f32>,
        hidden: &Array2<f32>,
        pre_activation2: &Array2<f32>,
        output: &Array2<f32>,
        target: &Array2<f32>,
    ) {
        let batch_size = input.shape()[1] as f32;

        // Calculate gradients for output layer
        let output_error = output - target;
        
        // Gradients for weights2 and biases2
        let d_weights2 = output_error.dot(&hidden.t()) / batch_size;
        let d_biases2 = &output_error.sum_axis(Axis(1)).insert_axis(Axis(1)) / batch_size;

        // Backpropagate error to hidden layer
        let hidden_error = self.weights2.t().dot(&output_error);
        let hidden_gradient = &hidden_error * &relu_derivative(hidden);

        // Gradients for weights1 and biases1
        let d_weights1 = hidden_gradient.dot(&input.t()) / batch_size;
        let d_biases1 = &hidden_gradient.sum_axis(Axis(1)).insert_axis(Axis(1)) / batch_size;

        // Update weights and biases using gradient descent
        self.weights2 = &self.weights2 - &(d_weights2 * self.learning_rate);
        self.biases2 = &self.biases2 - &(d_biases2 * self.learning_rate);
        self.weights1 = &self.weights1 - &(d_weights1 * self.learning_rate);
        self.biases1 = &self.biases1 - &(d_biases1 * self.learning_rate);
    }

    fn train(&mut self, input: &Array2<f32>, target: &Array2<f32>) -> f32 {
        // Forward pass
        let (hidden, pre_activation2, output) = self.forward(input);

        // Calculate loss (cross-entropy)
        let epsilon = 1e-15;
        let loss = -target * &output.mapv(|x| (x + epsilon).ln());
        let loss = loss.sum() / (input.shape()[1] as f32);

        // Backward pass
        self.backward(input, &hidden, &pre_activation2, &output, target);

        loss
    }
}

}

fn main() {
    // Create a neural network
    let mut nn = NeuralNetwork::new(6, 6, 4, 0.01);

    // Create sample input
    let mut input = Array2::zeros((nn.input_size, 1));
    populate_array(&mut input, nn.input_size, 1);

    // Create sample target
    let mut target = Array2::zeros((nn.output_size, 1));
    populate_array(&mut target, nn.output_size, 1);
    
    // Calculate initial cross entropy loss before training
    let (_, _, initial_output) = nn.forward(&input);
    let epsilon = 1e-15;
    let initial_loss = -&target * &initial_output.mapv(|x| (x + epsilon).ln());
    let initial_loss = initial_loss.sum() / (input.shape()[1] as f32);
    println!("Initial loss before training: {}", initial_loss);

    
    // Train for one iteration
    let loss = nn.train(&input, &target);
    //println!("Loss: {}", loss);

    // Calculate loss after training
    let (_, _, final_output) = nn.forward(&input);
    let epsilon = 1e-15;
    let final_loss = -&target * &final_output.mapv(|x| (x + epsilon).ln());
    let final_loss = final_loss.sum() / (input.shape()[1] as f32);
    println!("loss after training: {}", final_loss);

}

We can now run an iteration of a forward pass and backward pass.

You can repeatedly run this cell to see the loss decreasing.

#![allow(unused)]
fn main() {
main()
}

Initial loss before training: 1.5886599
loss after training: 1.5812341





()

Has anyone implemented neural networks in Rust?

See for example candle. Perhaps a pun on "Torch" from "PyTorch".

But the point isn't to use an existing package but learn how to build a simple one from basic linear algebra principles.

All of this and more...

• Quick start: https://github.com/rust-ndarray/ndarray/blob/master/README-quick-start.md
• Numpy equivalence: https://docs.rs/ndarray/latest/ndarray/doc/ndarray_for_numpy_users/index.html
• Linear Algebra: https://github.com/rust-ndarray/ndarray-linalg/blob/master/README.md

Technical Coding Challenge

Coding Challenge

Coding Challenge Review