# Neural Networks Optimized with Matrix Calculation

## Abstract


I implemented a neural network optimized with matrix calculation. Though neural networks can be implemented with many 'for' loops, which may be easy to read, the traininig time may become longer because we need at leas 3 'for' loops. When we implement such a network with 30 epochs, 50000 mnist training data, and mini_batch size 10, the network has to implement 150000 loops in tota.

This is quite expensive, so I implemented a neural network which uses two for loops and some matrix calculation. I used matrix caculation on update over mini-batch example. The update shares the same weights and biases in the network, so we can conduct the update simultaneously over the examples in the mini-batch.

## Introduction


Usually, update over mini-batch is conducted using one 'for' loop. We iterate over the examples in the mini-batch. In each loop, we calculate the gradient of the cost function with respect to weight matrix nad biases of the neurons. 

Suppose that ${C_x}_i$ is the cost calculated over a specific example $x_i$ in the mini_batch. Althoug we want to calculate the gradient of the cost $C_{all}$ over the entier training examples, it would be quite expensive computation. So we estimate the gradient with some randomely chosen examples. Then out estimate of the gradient is: 
$$
\widehat{ \nabla_{w} C_{all} } = \frac{1}{m}\sum_{i=0}^{m}\nabla_{w}C_{xi} 
$$
$$
\widehat{ \nabla_{b} C_{all} } = \frac{1}{m}\sum_{i=0}^{m}\nabla_{b}C_{xi} 
$$
where m is the mini-batch size. Each gradient of $C_x$ with respect to $\mathbf{w}$ and b can be calculated using backpropagation. Suppose that we use the corss entropy as the cost function and that add a L2 regularization term. Then the cost over the entire training examples is representes as:


$$C_{all} = \sum_{x}^{}C_{0x} + \frac{\lambda}{2n}\sum_{\mathbf{w}}^{}\mathbf{w}^2 $$
$$\text{where }C_{0x} = \sum_{j}^{}[ y_j log a_j^L + (1-y_j)log( 1 - a_j^L )] $$

$a_j^L$ is the jth activation output from the output layer. So we update $\mathbf{w}$ as :

$$\mathbf{w} \leftarrow \mathbf{w} - \eta \nabla_{w}C_{all} $$
$$= \mathbf{w} - \eta ( \nabla_{w}C_{0all} + \frac{\lambda}{n}\mathbf{w} )$$
$$= ( 1 - \frac{\eta\lambda}{n} )\mathbf{w} - \eta\nabla_{w}C_{0all}$$


Because we can not clculate $C_{0all}$, we estimate it over the examples in a mini_batch.

$$\mathbf{w} = ( 1 - \frac{\eta\lambda}{n} )\mathbf{w} - \eta\frac{1}{m}\sum_{i=0}^{m}\nabla_{w}C_{0xi} $$

and biases are also updated as:

$$\mathbf{b} \leftarrow \mathbf{b} - \frac{\eta}{m}\sum_{x}^{}\nabla_{b}C_{0xi}$$

But we cannot directly drive these 
gradient of the cost function, so we use the back propagation algorithm to get them. Let $\mathbf{z}^l$ and $\mathbf{a}^l$ be the inputs to the neurons in the lth layer and its sigmoid activation respectively. And we also define $\mathbf{\delta}^l$ as the gradient of the const function with respect to $\mathbf{z}^l$
$$\mathbf{z}^{l+1} = \mathbf{w}^{l}\mathbf{z}^{l} + \mathbf{b}^{l}$$
$$\mathbf{a}^{l} = \sigma ( \mathbf{z}^{l} )$$
$$\mathbf{\delta}^{l} = \nabla_{z^l}C_{0xi}$$

 

## Python code
The python code is composed with three separete file. 
The first code $\textbf{network_optimized.py}$ defines a class called 'Network'. This class holds basic parameters that defines the structure of the network. We give the number of neurons in the hidden layer when we instanciate it. The main function for training is 'SGD'. We give the function some hyper parameters specifying how the network learns, such as learing late $\eta$, the number of epochs, mini-bach size, and lambda $\lambda$. $\lambda$ defines how much the regularization term affects the way the networks learns. When it is large, it will be more beneficial, when we minimize the cost function, to make the weights small rather than minimize the non-regularized part of the cost function. And hence, the network will have strong generality. Conversely, when $\lambda = 0$, the network will behave as a non-regularized one and the weights may become larger as the training proceeds.

SGD function calls 'update_mini_batch' function to update $\mathbf{w}$ and $b$ over an mini-batch. We use matrix calculation to derive the estimated gradient of the cost function rather than iterate over the examples in the mini-batch. The update function takes mini-batch examples minibatch_x and minibatch_y, which are matrices. The number of rows is the mini-batch size. minibatch_x and minibatch_y have 784 and 10 columns respectively.

We first feedforward through the network to calculate the $\mathbf{z}$ vectors and its activations. Note that we calculate them simultaneously over all the exmples in the mini-batch using matrix calculation. Afther the feedforward, we conduct the back-propagation. Also in the back-propagation, $\delta$ vector is expanded to a matrixs, where jth column in the matrix represents $\delta$ of the jth example in the mini-batch. $\mathbf{z}$ vecotors and its activations also have the same structure.
$$
\begin{align}
\mathbf{z} &= ( \mathbf{z_1}, \dotsc, \mathbf{z_m} ) \\
\mathbf{a} &= \sigma( \mathbf{z} )

\end{align}
$$


Calculating the gradient of the cost function with respect to $\mathbf{w}$ is a little complicated. Usually, when we do that with 'for' loop, we first calculate the gradients for each example in the mini-batch and then take its average. However, with the matrix compulation we conduct the two operations simultraneously in a single line of code.





```python

# Last Updated: March 6, 2018
# Note:
# the vectorized version of the network2. Only one hidden layer is allowed,
# and we use only cross entropy cost.


import random
import sys
import numpy as np

class Network ( object ) :
    def __init__ ( self, num_hidden_neuron ) :
        # the number of neurons in the hidded layer
        self.num_hidden_neuron = num_hidden_neuron
        # the definistions of the biases vecotors and weights matrices.
        self.biases_hidlayer = np.random.randn ( num_hidden_neuron )
        self.biases_outlayer = np.random.randn ( 10 )
        self.weights_hidlayer = \
          np.random.randn ( num_hidden_neuron, 784 ) / 10.
        self.weights_outlayer = \
          np.random.randn ( 10, num_hidden_neuron ) / num_hidden_neuron


    def SGD ( self,
              # when whole_training_data is the whole mnist training data,
              # whole_training_data [ 0 ] <-shape( 50000, 784 )
              # whole_training_data [ 1 ] <-shape( 50000, 10 )
              whole_training_data,
              epochs,
              mini_batch_size, # 10 in demo.py
              eta, # 0.3 in demp.py
              lmbda = 0.0, # 3.0 in demo.py
              # evaluation data is either validation data or test data.
              evaluation_data = None,
              evaluate_on_training_data = False,
              evaluate_on_evaluation_data = False ) :

        if evaluation_data : n_data = len ( evaluation_data [ 0 ] ) # 10000
        whole_training_data_size = len ( whole_training_data [ 0 ] ) # 50000
        training_cost, training_accuracy = [], []
        evaluation_cost, evaluation_accuracy = [], []
        index = np.arange ( whole_training_data_size )
        for epc in xrange ( epochs ) :
            print "Epoch %d: Start" % epc
            np.random.shuffle ( index )
            whole_training_x = whole_training_data [ 0 ] [ index, : ]
            whole_training_y = whole_training_data [ 1 ] [ index, : ]
            for j in xrange ( 0, whole_training_data_size, mini_batch_size  ) :
                minibatch_x = \
                  whole_training_x [ j : j + mini_batch_size, : ]
                minibatch_y = \
                  whole_training_y [ j : j + mini_batch_size, : ]
                self.update_mini_batch ( minibatch_x, minibatch_y, eta,
                                    lmbda, whole_training_data_size )

            print "Epoch %d: training complete" % epc
            if evaluate_on_training_data :
                accuracy, cost = self.evaluate ( whole_training_data [ 0 ],
                                                 whole_training_data [ 1 ],
                                                 lmbda )
                training_cost.append ( cost )
                training_accuracy.append ( accuracy )
            if evaluate_on_evaluation_data and evaluation_data :
                accuracy, cost = self.evaluate ( evaluation_data [ 0 ],
                                                 evaluation_data [ 1 ],
                                                 lmbda )
                evaluation_cost.append ( cost )
                evaluation_accuracy.append ( accuracy )

        return training_cost, training_accuracy,\
               evaluation_cost, evaluation_accuracy

    def update_mini_batch ( self, minibatch_x, minibatch_y, eta, lmbda,
                            whole_training_data_size ) :
        mini_batch_size = minibatch_x.shape [ 0 ]
        ## feedforward
        a_inlayer = minibatch_x.T
        z_hidlayer = np.dot ( self.weights_hidlayer, a_inlayer ) + \
                     self.biases_hidlayer [ :, None ]
        a_hidlayer = sigmoid ( z_hidlayer )
        z_outlayer = np.dot ( self.weights_outlayer, a_hidlayer ) + \
                     self.biases_outlayer [ :, None ]
        a_outlayer = sigmoid ( z_outlayer )

        ## back propagate
        delta_outlayer = a_outlayer - minibatch_y.T
        biases_outlayer_grad = delta_outlayer.mean ( axis = 1 )
        weights_outlayer_grad = \
          np.dot ( delta_outlayer, a_hidlayer.T ) / mini_batch_size
        delta_hidlayer = \
          sigmoid_prime ( z_hidlayer ) * \
          np.dot ( self.weights_outlayer.T, delta_outlayer )
        biases_hidlayer_grad = delta_hidlayer.mean ( axis = 1 )
        weights_hidlayer_grad = \
          np.dot ( delta_hidlayer, a_inlayer.T ) / mini_batch_size

        ## update biases
        self.biases_outlayer = \
          self.biases_outlayer - eta * biases_outlayer_grad
        self.biaess_hidlayer = \
          self.biases_hidlayer - eta * biases_hidlayer_grad
        # update weights
        self.weights_outlayer = \
          ( 1. - eta * lmbda / whole_training_data_size ) * \
          self.weights_outlayer - eta * weights_outlayer_grad
        self.weights_hidlayer = \
          ( 1. - eta * lmbda / whole_training_data_size ) * \
          self.weights_hidlayer - eta * weights_hidlayer_grad

    def evaluate ( self, data_x, data_y, lmbda ) :
        '''
        calcuate accuracy and cost on the given data. the format of data_x
        and data_y is the same as the training_data. The cost function is
        the coross entropy cost
        '''
        num_data = data_x.shape [ 0 ]

        ## feedforward
        a_inlayer = data_x.T
        z_hidlayer = np.dot ( self.weights_hidlayer, a_inlayer ) + \
                     self.biases_hidlayer [ :, None ]
        a_hidlayer = sigmoid ( z_hidlayer )
        z_outlayer = np.dot ( self.weights_outlayer, a_hidlayer ) + \
                     self.biases_outlayer [ :, None ]
        a_outlayer = sigmoid ( z_outlayer )

        ## calcualte accuracy
        network_output = a_outlayer.argmax ( axis = 0 )
        desired_output = data_y.argmax ( axis = 1 )
        succeeds = ( network_output == desired_output ).sum ()
        accuracy = ( succeeds * 100.0 ) / num_data

        # calculate cost
        cost_C0 = np.nan_to_num \
          ( - data_y.T * np.log ( a_outlayer ) - \
          ( 1 - data_y.T ) * np.log ( 1 - a_outlayer ) ).sum () / num_data
        cost_regularization_term = \
          0.5 * lmbda / num_data *\
          ( np.linalg.norm ( self.weights_hidlayer ) ** 2 +\
            np.linalg.norm ( self.weights_outlayer ) ** 2 )
        cost = cost_C0 + cost_regularization_term

        return accuracy, cost




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

def sigmoid_prime ( x ) :
    # return sigmoid ( x ) * ( 1 - sigmoid ( x ) )
    return 1./ ( 2. + np.exp ( x ) + np.exp ( -x ) )

```

41ff19585acd62646342229f31ce8e680fcec662
ipynb
Jupyter Notebook
Untitled.ipynb
YoshioYamauchi/neuralnet_with_matrix
e39f7a8fc93d9dd595301bcfe90bae6bec1010da
MIT
Untitled.ipynb
YoshioYamauchi/neuralnet_with_matrix
e39f7a8fc93d9dd595301bcfe90bae6bec1010da
MIT
Untitled.ipynb
YoshioYamauchi/neuralnet_with_matrix
e39f7a8fc93d9dd595301bcfe90bae6bec1010da
MIT
Qwen/Qwen-72B
1. YES
2. YES