Examples

Single input, single output

In this section, the goal is to use the FNN on small examples in order to have a feeling for the expected results.

Hint

On some systems, even if the matplotlib library is installed, a window manager backend is missing. It leads to an error when plt.show().

The solution is to add the pyqt6 dependency to the project.

Exercise (intermediate)

1. Create the smallest and simplest possible FNN and learn the parameters. Comment. (no hidden layer, linear activation, only two entries in the train set). Make an inference. Comment.

2. Increase the size of the train set and use non aligned points. Plot the result and comment.

3. Build a train set around a known non linear function (for example \(f(x)=x^2+5\)). Use one hidden layer and with a logistic activation. Plot and comment.

FNN code

The FNN implementation used in the exercises below is:

src

"""Forward Neural Network implementation from scratch using numpy.
Three sections:
- the vectorized activation functions and their derivative
- the Layer class to be used in the FNN class
- the FNN class"""

from dataclasses import dataclass, field
from collections.abc import Callable, Sequence

from tqdm import tqdm
import numpy as np
import numpy.typing as npt

FLOAT = np.float64
type Array2D = npt.NDArray[FLOAT]
type Vector = npt.NDArray[FLOAT]
type FuncVector = Callable[[Vector], Vector]


class IncompatibleDimensions(Exception): ...


### vectorized activation and their derivative


def relu(vector: Vector) -> Vector:
    result: Vector = np.zeros(vector.shape)
    is_pos = vector > 0.0
    result[is_pos] = vector[is_pos]
    return result


def relu_derivative(vector: Vector) -> Vector:
    result: Vector = np.zeros(vector.shape)
    is_pos = vector > 0.0
    result[is_pos] = 1.0
    return result


def logistic(vector: Vector) -> Vector:
    return 1.0 / (1 + np.exp(-vector))


def logistic_derivative(vector: Vector) -> Vector:
    v: Vector = logistic(vector)
    return v * (1.0 - v)


def tanh(vector: Vector) -> Vector:
    return np.tanh(vector)


def tanh_derivative(vector: Vector) -> Vector:
    return 1.0 - np.tanh(vector) ** 2


def linear(vector: Vector) -> Vector:
    return vector


def linear_derivative(vector: Vector) -> Vector:
    return np.ones(vector.shape)


def softmax(vector: Vector) -> Vector:
    num: Vector = np.exp(vector)
    return num / sum(num)


def softmax_derivative(vector: Vector) -> Vector:
    return vector * (1.0 - vector)


_activations_and_derivatives: dict[FuncVector, FuncVector] = {
    relu: relu_derivative,
    linear: linear_derivative,
    logistic: logistic_derivative,
    softmax: softmax_derivative,
    tanh: tanh_derivative,
}


@dataclass
class Layer:
    """A layer has a certain number of nodes on the left,
    and another on the right. Every node on the left is
    fully connected to every node on the right.

    A layer has the ownership of the nodes on the right only.

    The nodes on the right have two states: z then a
    - before activation z_i = sum_k w_ik a_k_previous_layer
    - affter activation a_i = sigma(z_i)

    The weights and biases are initialized randomly,
    potentially using a seed for repeatability."""

    nb_nodes_left: int
    nb_nodes_right: int
    seed: int | None = None

    def __post_init__(self) -> None:
        self.init_parameters_random()
        self.z: Vector = np.zeros(self.nb_nodes_right, dtype=FLOAT)
        self.a: Vector = np.zeros(self.nb_nodes_right, dtype=FLOAT)

    def init_parameters_random(self) -> None:
        rng = np.random.default_rng(self.seed)
        self.weights: Array2D = (
            rng.random((self.nb_nodes_right, self.nb_nodes_left), dtype=FLOAT) - 0.5
        )
        self.biases: Vector = rng.random(self.nb_nodes_right, dtype=FLOAT) - 0.5


@dataclass
class FNN:
    """A FNN has a sequence of layers. It owns all the hidden layers
    and the output layer, not the input layer. However, to build
    the first layer, the number of nodes of the input layer must be known.

    Every layer has an activation, and its derivative can be retrieved in
    the activation_derivatives dict.

    The three main methods are:
    - infer: forward execution of the FNN from input to output
      with the current weights and biases
    - fit: learn all the weights and biases on a train set
    - mean_accuracy: assess the accuracy of the estimations
    """

    nb_nodes_per_layer_including_input: Sequence[int]
    activations: Sequence[FuncVector]
    seed: int | None = None

    activation_derivatives: dict[FuncVector, FuncVector] = field(
        default_factory=_activations_and_derivatives.copy
    )

    def __post_init__(self) -> None:
        """creates a sequence of layers."""
        self.layers: list[Layer] = []
        self.nb_layers: int = len(self.nb_nodes_per_layer_including_input) - 1
        if self.nb_layers != len(self.activations):
            raise IncompatibleDimensions(
                "The number of activation functions is not compatible with the number of nodes."
            )

        for ind in range(self.nb_layers):
            layer = Layer(
                self.nb_nodes_per_layer_including_input[ind],
                self.nb_nodes_per_layer_including_input[ind + 1],
                self.seed + ind if self.seed is not None else None,
                # to avoid using the same sequence of random numbers
                # for every layer which could introduce unwanted symetries
            )
            self.layers.append(layer)

    def infer(self, input_: Vector) -> Vector:
        """forward sweep across the FNN,
        z_right = w . a_left
        a_right = sigma(z_right)

        changes the attributes z and a of every layer."""
        if self.bad_input_or_target_size(
            self.nb_nodes_per_layer_including_input[0], [input_]
        ):
            raise IncompatibleDimensions("The size of the input is incompatible.")

        left: Vector = input_
        for ind_layer, layer in enumerate(self.layers):
            layer.z = layer.weights.dot(left) + layer.biases
            layer.a = self.activations[ind_layer](layer.z)
            left = layer.a

        return left

    def init_parameters_gradient_zero(self) -> tuple[list[Array2D], list[Vector]]:
        """initialize the arrays for the weight and biases gradients to zero
        with the correct dimension to be able to compute the mean gradients
        accross the train set."""
        biases_all: list[Vector] = []
        weights_all: list[Array2D] = []
        # init to zero for cumsum over training examples
        for ind_layer in range(0, self.nb_layers):
            layer = self.layers[ind_layer]
            weights_all.append(
                np.zeros((layer.nb_nodes_right, layer.nb_nodes_left), dtype=FLOAT)
            )
            biases_all.append(np.zeros(layer.nb_nodes_right, dtype=FLOAT))
        return weights_all, biases_all

    def get_left_a_or_input(self, ind_layer: int, input_: Vector) -> Vector:
        """helper function, the FNN does not own the input,
        but uses it when left nodes of the first layer are required for an operation"""
        return input_ if ind_layer == 0 else self.layers[ind_layer - 1].a

    def get_layer_activation_activation_derivative(
        self, ind_layer: int
    ) -> tuple[Layer, FuncVector, FuncVector]:
        """helper function, code factorisation"""
        layer: Layer = self.layers[ind_layer]
        activation: FuncVector = self.activations[ind_layer]
        activation_derivative: FuncVector = self.activation_derivatives[activation]
        return layer, activation, activation_derivative

    def parameter_gradient(
        self, inputs: Sequence[Vector], targets: Sequence[Vector]
    ) -> tuple[list[Array2D], list[Vector]]:
        """The loss gradient with respect to the parameters: either weights or biases.

        dJ/dz of the last layer is first calculated.
        The detail of this computation depends on the loss function itself: usually
        either multi-class or regression. It is performed in the last_layer_dJdz method.

        The dJdz of the last layer allows to obtain the dJdb and dJdw of the last layer:
        dJ/dw_jk = dJ/dz_j a_k  with index k on the left and j on the right
        (which can be obtained by broadcast
                    a_0 ... a_k ...
        dJ/dz_0
        ...
        dJ/dz_j
        ...
        )

        Then, using back-propagation, every layer can be obtained from one before last to first:
        dJ/dw_kl = [ sum_j [w_jk s'(z_k) ] dJ/dz_j ] a_l  with index k on the left and j on the right

        c_jk = w_jk s'(z_k) simple broadcasting

        sum_j c_jk  dJ/dz_j is almost a dot product, requires a transpose -> d = c^T . dJ/dz

        dJ/dw_kl = c_k a_l can be obtained by broadcast similar as the one for the last layer

        At the end, the sum of the gradients for every case of the train set is divided by
        the number of cases of the train set to get the average gradient.
        """
        # first store the sum for all training example of gradients
        # then at the end, the mean gradient
        gradients_weights: list[Array2D]
        gradients_biases: list[Vector]
        gradients_weights, gradients_biases = self.init_parameters_gradient_zero()

        for input_, target in zip(inputs, targets):
            # update all a and z for given training example
            self.infer(input_)

            ## last layer
            ind_layer = self.nb_layers - 1
            layer, activation, activation_derivative = (
                self.get_layer_activation_activation_derivative(ind_layer)
            )
            left: Vector = self.get_left_a_or_input(ind_layer, input_)

            dJdz: Vector = self.last_layer_dJdz(
                layer, activation, activation_derivative, target
            )

            gradients_biases[ind_layer] += dJdz
            dJdw: Array2D = dJdz[:, np.newaxis] * left  # broadcast
            gradients_weights[ind_layer] += dJdw

            dJdz_prev = dJdz

            ## layers from one before last included, to first
            for ind_layer in range(self.nb_layers - 2, -1, -1):
                layer, _, activation_derivative = (
                    self.get_layer_activation_activation_derivative(ind_layer)
                )
                left: Vector = self.get_left_a_or_input(ind_layer, input_)
                dJdz = (
                    self.layers[ind_layer + 1].weights
                    * activation_derivative(layer.z)  # broadcast
                ).T.dot(dJdz_prev)

                gradients_biases[ind_layer] += dJdz
                dJdw = dJdz[:, np.newaxis] * left  # broadcast
                gradients_weights[ind_layer] += dJdw

                dJdz_prev = dJdz

        # from sum to average
        nb_training = len(inputs)
        for ind_layer in range(0, self.nb_layers):
            gradients_weights[ind_layer] /= nb_training
            gradients_biases[ind_layer] /= nb_training

        return gradients_weights, gradients_biases

    @staticmethod
    def last_layer_dJdz(
        layer: Layer,
        activation: FuncVector,
        activation_derivative: FuncVector,
        target: Vector,
    ) -> Vector:
        """Depending on the last layer activation, one assumes the loss is:
        + multi-class cross entropy for softmax
        + mean square error for linear"""
        error: Vector = layer.a - target
        if activation == softmax:  # multi-class cross entropy gradient
            dJdz: Vector = error * target
        elif activation == linear:  # regression mean square gradient
            dJdz: Vector = error
        else:  # regression mean square gradient, other activations
            dJdz: Vector = error * activation_derivative(layer.z)
        return dJdz

    def fit(
        self,
        inputs: Sequence[Vector],
        targets: Sequence[Vector],
        learning_rate: float = 0.1,
        nb_iter_gradient: int = 1_000,
        regularization: float = 0.0,
        rel_convergence: float = 1e-3,
        display_interval: int | None = 100,
    ) -> None:
        """Gradient descent algorithm on both the weights and biases, no stopping criteria,
        no stochastic pooling.

        The learning_rate is for the gradient from the train set.

        The regularization is the lambda_r of the basic regularization:
        J_r = lambda_r/2 (sum_all w^2 + sum_all b^2)
        """
        if self.bad_input_or_target_size(
            self.nb_nodes_per_layer_including_input[0], inputs
        ):
            raise IncompatibleDimensions("The size of the input is incompatible.")
        if self.bad_input_or_target_size(
            self.nb_nodes_per_layer_including_input[-1], targets
        ):
            raise IncompatibleDimensions("The size of the target is incompatible.")

        first_gradient_norm = 0.0

        for ind_iter in tqdm(range(nb_iter_gradient), desc="fnn fit"):
            gradients_weights, gradients_biases = self.parameter_gradient(
                inputs, targets
            )
            gradient_norm = 0.0
            for ind_layer in range(self.nb_layers):
                self.apply_regularization(regularization)
                self.layers[ind_layer].weights -= (
                    learning_rate * gradients_weights[ind_layer]
                )
                self.layers[ind_layer].biases -= (
                    learning_rate * gradients_biases[ind_layer]
                )
                gradient_norm += np.linalg.norm(
                    gradients_weights[ind_layer]
                    + regularization * self.layers[ind_layer].weights,
                    "fro",
                )

            if ind_iter == 0:
                first_gradient_norm = gradient_norm

            if display_interval is not None and ind_iter % display_interval == 0:
                tqdm.write(str(gradient_norm / first_gradient_norm))

            if gradient_norm / first_gradient_norm < rel_convergence:
                break

    def apply_regularization(self, regularization: float):
        """The gradient descent for the regularization used is:
        w_new = w_current - lambda_r w_current = (1-lambda_r) w_current
        w_new *= 1-lambda_r"""
        for ind_layer in range(self.nb_layers):
            self.layers[ind_layer].weights *= 1.0 - regularization
            self.layers[ind_layer].biases *= 1.0 - regularization

    @staticmethod
    def bad_input_or_target_size(expected_size: int, vectors: Sequence[Vector]) -> bool:
        """check that every vector in a sequence is of the expected size"""
        return any([vector.size != expected_size for vector in vectors])

    def jacobian(self, input_: Vector) -> Array2D:
        self.infer(input_)
        jacobian: Array2D = np.eye(self.nb_nodes_per_layer_including_input[0])
        for ind_layer in range(0, self.nb_layers):
            layer, _, activation_derivative = (
                self.get_layer_activation_activation_derivative(ind_layer)
            )
            jacobian = (
                layer.weights.dot(jacobian).T
                * activation_derivative(layer.z)  # broadcast
            ).T

        return jacobian

    def mean_accuracy(
        self,
        inputs: list[Vector],
        targets: list[Vector],
        metric_single_example: Callable[[Vector, Vector], float | bool],
    ) -> float:
        """with a given metric callable on every pair (prediction, target)
        which gives either a measure of their distance or a measure of exact match,
        this method computes a mean accuracy of the FNN over a train or test set."""
        residual: float = 0.0
        for input_, target in zip(inputs, targets):
            prediction: Vector = self.infer(input_)
            residual += metric_single_example(prediction, target)
        return residual / len(inputs)


def mean_square_loss_single(prediction: Vector, target: Vector) -> float:
    return float(np.linalg.norm(prediction - target) ** 2 / target.size)


def is_correct_class_single(prediction: Vector, target: Vector) -> bool:
    """is the most probable estimated class the correct one?"""
    return np.argmax(prediction) == np.argmax(target)

Correction

1. Careful, even if each input and target are of length 1, they still must be Vector.

   Two entries (meaning two inputs and two targets), means that the parameters represent exactly the slope and intercept of the straight line going the two points.

   This could have been written as a test.

src

def R_to_R_linear_two_points(
    learning_rate: float = 0.1,
    nb_iter_gradient: int = 1_000,
    regularization: float = 0.0,
    rel_convergence: float = 1e-6,
    display_interval: int | None = 100,
) -> None:
    fnn = FNN([1, 1], [linear])
    inputs: list[Vector] = [np.array([0.0]), np.array([1.0])]
    targets: list[Vector] = [np.array([1.0]), np.array([3.0])]
    fnn.fit(
        inputs,
        targets,
        learning_rate=learning_rate,
        nb_iter_gradient=nb_iter_gradient,
        regularization=regularization,
        rel_convergence=rel_convergence,
        display_interval=display_interval,
    )
    middle: Vector = np.array([0.5])
    inferred: Vector = fnn.infer(middle)
    print("R to R, no hidden layer, linear activation")
    print(f"inputs: {inputs}")
    print(f"targets: {targets}")
    print(f"weights: {fnn.layers[0].weights}")
    print(f"biases: {fnn.layers[0].biases}")
    print(f"inference:  in={middle}, out={inferred}")
    print("-" * 30)

R to R, no hidden layer, linear activation
inputs: [array([0.]), array([1.])]
targets: [array([1.]), array([3.])]
weights: [[1.99999114]]
biases: [1.00000547]
inference:  in=[0.5], out=[2.00000105]

2. Using a larger train set and building it around a knwown slope and intercept adding a noise, one can check that the obtained model is in fact the least square best stright line.

src

def R_to_R_linear_with_noise(
    noise_amplitude: float = 0.4,
    seed: int = 0,
    path: str | None = None,
    learning_rate: float = 0.1,
    nb_iter_gradient: int = 1_000,
    regularization: float = 0.0,
    rel_convergence: float = 1e-12,
    display_interval: int | None = 100,
) -> None:
    fnn = FNN([1, 1], [linear])
    left_bound, right_bound = 0.0, 5.0
    nb_points: int = 11
    abscissas = np.linspace(
        left_bound, right_bound, num=nb_points, endpoint=True, dtype=FLOAT
    )
    inputs: list[Vector] = [np.array([abscissa]) for abscissa in abscissas]
    slope: float = 2.0
    intercept: float = 1.0
    rng = np.random.default_rng(seed=seed)

    def generate_y(abscissa):
        return (
            intercept
            + slope * abscissa
            + rng.uniform(-noise_amplitude, noise_amplitude)
        )

    targets: list[Vector] = [np.array([generate_y(abscissa)]) for abscissa in abscissas]
    fnn.fit(
        inputs,
        targets,
        learning_rate=learning_rate,
        nb_iter_gradient=nb_iter_gradient,
        regularization=regularization,
        rel_convergence=rel_convergence,
        display_interval=display_interval,
    )

    matrix_for_least_squares = np.vstack([abscissas, np.ones(nb_points)]).T
    least_squares_slope, least_sqares_intercept = np.linalg.lstsq(
        matrix_for_least_squares, np.array(targets).squeeze()
    )[0]

    print("R to R, no hidden layer, linear activation")
    print(f"base slope: {slope}")
    print(f"base intercept: {intercept}")
    print(f"inputs: {inputs}")
    print(f"targets: {targets}")
    print(f"weights: {fnn.layers[0].weights}")
    print(f"biases: {fnn.layers[0].biases}")
    print(f"least squares slope: {least_squares_slope}")
    print(f"least squares intercept: {least_sqares_intercept}")
    print("-" * 30)

    abscissas_plot = np.linspace(
        left_bound, right_bound, num=21, endpoint=True, dtype=FLOAT
    )
    ordinates_plot = []
    for abscissa in abscissas_plot:
        infered_ordinate = fnn.infer(np.array([abscissa]))
        ordinates_plot.append(infered_ordinate[0])

    fig, ax = plt.subplots()
    for abscissa, target in zip(inputs, targets):
        ax.plot(abscissa[0], target[0], "ok")
    ax.plot(abscissas_plot, ordinates_plot, "-+r")

    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.set_title("strictly linear model")

    if path is not None:
        fig.savefig(path)

base slope: 2.0
base intercept: 1.0
inputs: [array([0.]), array([0.5]), array([1.]), array([1.5]), array([2.]), array([2.5]), array([3.]), array([3.5]), array([4.]), array([4.5]), array([5.])]
targets: [array([1.10956935]), array([1.81582937]), array([2.63277882]), array([3.61322211]), array([5.25061619]), array([6.33020446]), array([7.08530862]), array([8.18359725]), array([9.03489999]), array([10.34805794]), array([11.25268284])]
weights: [[2.09138705]]
biases: [0.83123845]
least squares slope: 2.091387054024167
least squares intercept: 0.8312384509871413

FNN fit and inference: linear plus noise

3. With the proposed non-linear function and some noise, one hidden layer with two nodes with logistic activation.

src

def R_to_R_nonlinear_with_noise(
    nb_nodes_per_layer_including_input: Sequence[int],
    activations: Sequence[FuncVector],
    func: Callable,
    noise_amplitude: float = 3.0,
    seed: int = 0,
    path: str | None = None,
    learning_rate: float = 0.01,
    nb_iter_gradient: int = 100_000,
    regularization: float = 0.0,
    rel_convergence: float = 1e-3,
    display_interval: int | None = 100,
):
    fnn = FNN(nb_nodes_per_layer_including_input, activations)
    left_bound, right_bound = 0.0, 5.0
    nb_points: int = 11
    abscissas = np.linspace(
        left_bound, right_bound, num=nb_points, endpoint=True, dtype=FLOAT
    )
    inputs: list[Vector] = [np.array([abscissa]) for abscissa in abscissas]

    rng = np.random.default_rng(seed=seed)

    def generate_y(abscissa):
        return func(abscissa) + rng.uniform(-noise_amplitude, noise_amplitude)

    targets: list[Vector] = [np.array([generate_y(abscissa)]) for abscissa in abscissas]

    fnn.fit(
        inputs,
        targets,
        learning_rate=learning_rate,
        nb_iter_gradient=nb_iter_gradient,
        regularization=regularization,
        rel_convergence=rel_convergence,
        display_interval=display_interval,
    )

    left_extra, right_extra = -2.0, 3.0
    abscissas_plot = np.linspace(
        left_bound + left_extra,
        right_bound + right_extra,
        num=201,
        endpoint=True,
        dtype=FLOAT,
    )
    ordinates_plot = []
    for abscissa in abscissas_plot:
        infered_ordinate = fnn.infer(np.array([abscissa]))
        ordinates_plot.append(infered_ordinate[0])

    # breakpoint()

    fig, ax = plt.subplots()
    for abscissa, target in zip(inputs, targets):
        ax.plot(abscissa[0], target[0], "ok")
    ax.plot(abscissas_plot, ordinates_plot, "-+r")

    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.set_title(
        (
            f"with hidden layer ({fnn.nb_nodes_per_layer_including_input} nodes) "
            f"with {fnn.activations[0].__name__} activation model"
        )
    )
    if path is not None:
        fig.savefig(path)

Note:

• the two added logistic function are clearly visible

• extrapolation outside the range of training data will not give very satisfying predictions when one keeps in mind the non-linear function which has been used to generate this example

FNN fit and inference: non-linear plus noise, 1 hidden layer with 2 nodes

If the number of nodes in the hidden layer is increased:

FNN fit and inference: non-linear plus noise, 1 hidden layer with 6 nodes

Note that even if the inference sticks much closely to the training data, it does not mean that this inference is “better” than the previous one. If one recalls the function with has been used to generate the example, it means that the extra nodes here have mostly learned the “noise”. It will not generalize well one a test set with a new random noise.

One may achieve similar result with other activation function (below \(\text{atanh}\)) and with deeper and narrower network:

FNN fit and inference: non-linear plus noise, several hidden layers with 3 nodes

Deeper nets are usually harder to train.

Double input, double output

The goal is now to illustrate the FNN for an higher dimension problem: \(\mathbb{R}^2 \to \mathbb{R}^2\). The following figure shows in black a vector field:

\[\begin{split}\begin{cases} r &= \sqrt{x^2 + y^2} \\ \theta &= \arctan \left( \cfrac{y}{x} \right) \\ u(r, \theta) &= r * \cos(\theta^2) \\ v(r, \theta) &= r * \sin(\theta) \end{cases}\end{split}\]

In red, the inferred values at the same locations fit pretty well. The network has two 8 nodes hidden layer.

FNN fit and inference: 2d non-linear, several hidden layers