Laboratorio 09

Il Grande Refactor — Template e Matrix

Docente

Soluzione del Laboratorio 08

Box-Muller e il calcolo della scala con lo switch su enum class:

static double gauss() {
    double x1, x2;
    do { x1 = rng.ranf(); x2 = rng.ranf(); } while (x1 == 0 || x2 == 0);
    return std::cos(2 * M_PI * x2) * std::sqrt(-2.0 * std::log(x1));
}

double s = 0.0;
switch (init) {
    case InitType::HE:
        s = (distro == DistroType::UNIFORM) ? std::sqrt(6.0 / input_size)
                                            : std::sqrt(2.0 / input_size);
        break;
    case InitType::XAVIER:
        s = (distro == DistroType::UNIFORM)
                ? std::sqrt(6.0 / (input_size + output_size))
                : std::sqrt(2.0 / (input_size + output_size));
        break;
}

weights(i, j) = (distro == DistroType::UNIFORM) ? 2.0 * s * (rng.ranf() - 0.5)
                                                : s * gauss();

Chi ha fatto la verifica quantitativa ha trovato la deviazione standard dei pesi vicina a \(\sqrt{2/784} \approx 0.0505\).

Obiettivo di Oggi

Questo è il laboratorio più invasivo del corso. Il framework prende la forma definitiva che manterrà fino alla fine, e per farlo dobbiamo cambiare molte cose insieme.

Non è un capriccio: ogni singola modifica serve a rendere possibile la backpropagation efficiente che implementeremo fra due laboratori.

Warning

Fate un commit (o una copia) del progetto prima di iniziare. Oggi si rompe tanto codice.

Note

Sul tempo. È il laboratorio più denso della prima metà del corso: circa 150 righe da toccare. Non si può spezzare, perché il codice compila solo quando tutte le modifiche sono fatte. L’ordine dei Task è pensato per ridurre il tempo in cui nulla compila; i tre layer di attivazione sono quasi identici fra loro, scrivetene uno con cura e adattate gli altri due. Se non finite in aula, finite a casa: il Lab 10 riparte da qui e ha bisogno che tutto compili.

Le sei cose che cambiano

	prima (Lab 05–08)	da oggi
Tipo numerico	Eigen::MatrixXd (solo double)	Matrix<ntype> (float o double)
Interfaccia	ILayer (ilayer.hpp)	layer<ntype> (layer.hpp)
Pesi	weights + biases	Wt + dWt, senza bias
Costruttore	DenseLayer(input, output)	DenseLayer(Nnodes, Nin) — invertito
Init	enum class InitType/DistroType	int type, int distro
forward	ritorna una copia	ritorna Matrix<ntype>&

Le ultime tre meritano una spiegazione, perché sembrano passi indietro.

Perché sparisce il bias

Il framework BASIC_NEURAL++ calcola semplicemente:

\[ \text{outputs} = W_t \cdot \text{inputs} \]

Niente \(+\,b\). Non è una dimenticanza: è la scelta del libro “Make Your Own Neural Network” di Rashid, su cui il progetto è modellato.

Il motivo pratico: su MNIST con input normalizzati in \([0, 1]\) e attivazioni sigmoidi, il bias contribuisce poco all’accuratezza finale, mentre raddoppia i tensori da propagare, derivare e salvare. Togliendolo, il gradiente diventa una sola formula invece di due, e il codice della backpropagation si dimezza.

Note

Rimetterlo è un buon esercizio, e la struttura lo consente: basta aggiungere un tensore bt accanto a Wt e derivarlo. Ma il framework di riferimento — quello con cui confronterete i vostri risultati — non lo ha.

Perché si torna a int type, int distro

Gli enum class del Lab 08 erano didatticamente giusti, e il ragionamento sulla type safety resta valido. Il framework però usa due int con valori di default:

ntype randfunc(int type = 0, int distro = 1);   // type: 0=He, 1=Xavier
void init_weights(int type = 0, int distro = 1); // distro: 0=uniform, 1=gauss

La ragione è che init_weights è dichiarata virtuale nella classe base layer<ntype>, dove un enum definito in parametricLayer.hpp non sarebbe visibile. Mantenere gli int tiene l’interfaccia base indipendente dai layer concreti.

È un compromesso reale fra type safety e disaccoppiamento: sapere che esiste, e perché, vale più che avere sempre la risposta bella.

L’ordine invertito degli argomenti

Questo è l’errore che vi costerà più tempo se non ci fate attenzione adesso.

// PRIMA (Lab 05–08):  (quanti input, quanti neuroni)
DenseLayer(784, 128);      // 784 input -> 128 neuroni

// DA OGGI:            (quanti neuroni, quanti input)
DenseLayer<float>(128, 784);   // 784 input -> 128 neuroni

Il perché è coerente: la matrice Wt ha forma [Nnodes × Nin], quindi gli argomenti seguono l’ordine delle sue dimensioni, come in Eigen.

Important

Wt è la \(\mathbf{W}\) delle lezioni, non la sua trasposta. (Lo schema con gli indici è nella Lezione 2.) Nelle Lezioni 2 e 3 la matrice dei pesi ha una riga per ogni neurone del layer e una colonna per ogni ingresso, e il forward è \(\mathbf{z} = \mathbf{W}\mathbf{x}\): esattamente la forma [Nnodes × Nin] di Wt e l’istruzione outputs = Wt * inputs. Sugli indici: Wt(i, j) ha i = neurone di questo layer (la riga, da 0 a Nnodes-1) e j = ingresso (la colonna, da 0 a Nin-1), come si vede nel doppio ciclo di init_weights. La “t” nel nome viene dal libro di Rashid, Make Your Own Neural Network, il testo di riferimento del corso: lì i pesi si scrivono \(w_{ij}\) con \(i\) = neurone del layer \(L\) (partenza) e \(j\) = neurone del layer \(L+1\) (arrivo), e la matrice che entra nel prodotto \(\mathbf{X} = \mathbf{W}\cdot\mathbf{I}\) è la trasposta di quella tabella. Il framework chiama quella matrice Wt per non nasconderlo; le lezioni la chiamano direttamente \(\mathbf{W}\). Nel Lab 13 ritroverete le formule della Lezione 3 tali e quali: dWt = delta * inputs.transpose() e delta_prev = Wt.transpose() * delta.

Per una rete 784 → 200 → 10 si scrive:

net.add(new DenseLayer<float>(200, 784));   // primo layer
net.add(new SigmoidLayer<float>());
net.add(new DenseLayer<float>(10, 200));    // secondo layer
net.add(new SigmoidLayer<float>());

Notate come i numeri si “incastrano”: 784 → (200, 784) → (10, 200) → 10. Se sbagliate l’ordine, forward() lancerà un std::runtime_error sulle dimensioni — il controllo c’è proprio per questo.

Task 1: eigen_wrap.hpp

Il file che rende possibile tutto il resto: un solo alias.

#ifndef __WRAP_EIGEN__
#define __WRAP_EIGEN__
#include <Eigen/Dense>

template <typename ntype>
using Matrix = Eigen::Matrix<ntype, Eigen::Dynamic, Eigen::Dynamic>;

#endif

Da qui in avanti Matrix<double> è esattamente Eigen::MatrixXd, e Matrix<float> è Eigen::MatrixXf. Una sola scrittura per entrambi.

Task 2: layer.hpp — la nuova interfaccia

Sostituisce ilayer.hpp. Notate i quattro membri protected: sono lo stato che ogni layer possiede.

#ifndef __LAYER__
#define __LAYER__
#include "./eigen_wrap.hpp"
#include <string>
#include <iostream>

template<typename ntype>
class layer {
protected:
    bool is_training = true;
    Matrix<ntype> outputs;                       // il risultato del forward
    Matrix<ntype> delta_l;                       // il gradiente del backward
    const Matrix<ntype>* inputs_ptr = nullptr;   // PUNTATORE all'input, non copia

public:
    virtual std::string get_type() const = 0;

    virtual void save_config(std::ostream& os) const {}
    virtual void save_weights(std::ostream& os) const {}
    virtual void load_weights(std::istream& is) {}

    virtual ~layer() = default;

    virtual Matrix<ntype>& forward(const Matrix<ntype>& inp) = 0;

    // Per ora un corpo di default: la backpropagation arriva nel Lab 13
    virtual Matrix<ntype>& backward(const Matrix<ntype>& out) { return delta_l; }

    virtual void train() { is_training = true; }
    virtual void eval()  { is_training = false; }

    virtual void init_weights(int type = 0, int distro = 1) {}
};
#endif

clone() non c’è più: la copia polimorfica verrà risolta diversamente, salvando e ricaricando il modello da file (Lab 11).

I tre metodi save_config, save_weights, load_weights sono ganci per quella serializzazione: oggi hanno il corpo vuoto. Prendono std::ostream& e std::istream&, non std::ofstream&/std::ifstream&: sono le classi base degli stream (anche std::cout è un ostream, std::cin un istream), così lo stesso codice potrà scrivere su file o a schermo senza cambiare una riga.

Perché forward ritorna un riferimento

Prima ogni layer restituiva una matrice per valore: una rete di 4 layer faceva 4 copie di matrici a ogni immagine, per 60.000 immagini per 20 epoche.

Ora ogni layer scrive nel proprio membro outputs e ne restituisce un riferimento:

Matrix<ntype>& forward(const Matrix<ntype>& inp) override {
    inputs_ptr = &inp;                    // memorizzo DOVE è l'input
    outputs = Wt * (*inputs_ptr);         // scrivo nel mio buffer
    return outputs;                       // restituisco un riferimento ad esso
}

E inputs_ptr è un puntatore: il layer non copia l’input, si ricorda soltanto dove sta. Servirà nel backward, dove la formula del gradiente ha bisogno dell’input di quel layer.

Warning

Il rovescio della medaglia: inputs_ptr punta a memoria che appartiene a qualcun altro. Se quella matrice viene distrutta prima del backward(), il puntatore diventa dangling (Lezione 8). Nel nostro flusso forward → backward questo non accade, ma è il motivo per cui non si può chiamare backward() senza aver prima chiamato forward().

Task 3: parametricLayer.hpp

Due classi in un file: la base che raggruppa i layer con pesi, e il DenseLayer.

#include "./layer.hpp"
#include "./randnumgen.hpp"
#include <cmath>

template <typename ntype>
class ParametricLayer : public layer<ntype> {
protected:
    ntype gauss(void) {                       // Box-Muller, come nel Lab 08
        double x1, x2;
        do { x1 = rng.ranf(); x2 = rng.ranf(); } while (x1 == 0 || x2 == 0);
        return std::cos(2 * M_PI * x2) * std::sqrt(-2.0 * std::log(x1));
    }
    // Contratto per i layer parametrici: chi ha pesi deve saper riceverli
    virtual void set_weights(const Matrix<ntype>& wm) = 0;

public:
    virtual ~ParametricLayer() = default;
};

Important

set_weights è puro virtuale: qualunque layer parametrico è obbligato a implementarlo. Serve a imporre i pesi dall’esterno senza distruggere e ricreare l’oggetto — lo useremo per caricare un modello da file (Lab 11) e vi servirà nel gradient check (Lab 13).

template <typename ntype>
class DenseLayer : public ParametricLayer<ntype> {
    using bc = ParametricLayer<ntype>;
    using bc::gauss;
    using bc::outputs, bc::delta_l, bc::inputs_ptr;   // <- IMPORTANTE

    ntype randfunc(int type = 0, int distro = 1) {
        double s;
        if (type == 0)   // He (Kaiming), per ReLU
            s = (distro == 0) ? sqrt(6.0 / Nin) : sqrt(2.0 / Nin);
        else             // Xavier (Glorot), per tanh e sigmoide
            s = (distro == 0) ? sqrt(6.0 / (Nnodes + Nin))
                              : sqrt(2.0 / (Nnodes + Nin));
        if (distro == 0) return 2.0 * s * (rng.ranf() - 0.5);
        else             return s * gauss();
    }

public:
    Matrix<ntype> Wt, dWt;      // pesi e gradienti
    int Nnodes, Nin;

    DenseLayer(int Nn, int Ni, int type = 0, int distro = 1)
        : Nnodes(Nn), Nin(Ni) {
        resize_weights(Nn, Ni);
        if (type >= 0) init_weights(type, distro);
    }

    // Ridimensiona pesi e gradienti: estratto dal costruttore perché
    // serve anche quando si cambia la forma del layer a posteriori
    void resize_weights(int Nn, int Ni) {
        Nnodes = Nn;
        Nin = Ni;
        Wt.resize(Nn, Ni);
        dWt.resize(Nn, Ni);
    }

    // L'implementazione del contratto di ParametricLayer
    void set_weights(const Matrix<ntype>& wm) override { Wt = wm; }

    void init_weights(int type = 0, int distro = 1) override {
        for (int i = 0; i < Nnodes; ++i)
            for (int j = 0; j < Nin; ++j)
                Wt(i, j) = randfunc(type, distro);
    }

    std::string get_type() const override { return "DenseLayer"; }

    Matrix<ntype>& forward(const Matrix<ntype>& inp) override {
        if (inp.rows() != Nin)
            throw std::runtime_error("[DenseLayer] Input size (" +
                std::to_string(inp.rows()) + ") non corrisponde a " +
                std::to_string(Nin));
        inputs_ptr = &inp;
        outputs = Wt * (*inputs_ptr);
        return outputs;
    }
};

La trappola dei template: using e this->

Questa riga non è decorativa:

using bc::outputs, bc::delta_l, bc::inputs_ptr;

In C++, quando una classe template deriva da un’altra classe template, il compilatore non cerca automaticamente i membri della base. Scrivere outputs = ... dà un errore di variabile non dichiarata, anche se il membro esiste. Due soluzioni equivalenti:

using bc::outputs;      // 1. la dichiari una volta e la usi normalmente
outputs = Wt * inp;

this->outputs = Wt * inp;   // 2. oppure la qualifichi ogni volta con this->

Il framework usa la prima nel DenseLayer e la seconda nelle attivazioni: vedrete entrambe.

Task 4: activation_layers.hpp

template <typename ntype>
class ActivationLayer : public layer<ntype> {
protected:
    Matrix<ntype> outputs, delta_l;
    int total_elements = 0;

    // Rialloca solo se la forma è cambiata: evita allocazioni inutili
    void match_shape(const Matrix<ntype>& inp) {
        if (total_elements != inp.size()) {
            outputs.resize(inp.size(), 1);
            delta_l.resize(inp.size(), 1);
            total_elements = inp.size();
        }
    }
};

template<typename ntype>
class SigmoidLayer : public ActivationLayer<ntype> {
    using bc = ActivationLayer<ntype>;
    using bc::delta_l, bc::outputs;
public:
    std::string get_type() const override { return "SigmoidLayer"; }

    Matrix<ntype>& forward(const Matrix<ntype>& inp) override {
        this->match_shape(inp);
        outputs = static_cast<ntype>(1.0)
                / (static_cast<ntype>(1.0) + ((-inp.array()).exp()));
        return outputs;
    }
};

template<typename ntype>
class tanhLayer : public ActivationLayer<ntype> {
    using bc = ActivationLayer<ntype>;
    using bc::delta_l, bc::outputs;
public:
    std::string get_type() const override { return "tanhLayer"; }

    Matrix<ntype>& forward(const Matrix<ntype>& inp) override {
        this->match_shape(inp);
        outputs = inp.array().tanh();
        return outputs;
    }
};

Note

Il nome è tanhLayer con la t minuscola, non TanhLayer. È il nome nel framework, e conta: nel Lab 11 la stringa get_type() verrà usata per ricostruire il modello da file.

Notate gli static_cast<ntype>(1.0): con ntype = float un 1.0 nudo è un double e forzerebbe Eigen a promuovere tutto il calcolo a doppia precisione, vanificando il risparmio di memoria.

Task 5: sequential.hpp

Sequential deriva da layer<ntype>: una rete è un layer, e potrà quindi essere inserita in un’altra rete.

template <typename ntype>
class Sequential : public layer<ntype> {
protected:
    layer<ntype>** layers;      // ancora un array crudo: cambierà nel Lab 10
    int capacity;
    int size;

public:
    Sequential(int max_layers = 10) : capacity(max_layers), size(0) {
        layers = new layer<ntype>*[capacity];
        for (int i = 0; i < capacity; ++i) layers[i] = nullptr;
    }

    ~Sequential() override {
        for (int i = 0; i < size; ++i) delete layers[i];
        delete[] layers;
    }

    std::string get_type() const override { return "Sequential"; }

    void add(layer<ntype>* l) {
        if (size >= capacity) throw std::out_of_range("Network capacity esaurita!");
        layers[size++] = l;
    }

    Matrix<ntype>& forward(const Matrix<ntype>& input) override {
        if (size == 0) throw std::runtime_error("Rete vuota!");

        Matrix<ntype>* current_out = &layers[0]->forward(input);
        for (int i = 1; i < size; ++i)
            current_out = &layers[i]->forward(*current_out);
        return *current_out;
    }
};

Il forward ora lavora con puntatori a matrici: &layers[0]->forward(...) prende l’indirizzo del buffer outputs del layer, e lo passa al successivo. Zero copie da un capo all’altro della rete.

Task 6: Il test che giustifica tutto

Il senso dei template si vede qui: la stessa rete, due precisioni.

#include "eigen_wrap.hpp"
#include "sequential.hpp"

int main() {
    rng.seed(42);

    // ATTENZIONE all'ordine: (Nnodes, Nin)
    Sequential<float> net32(4);
    net32.add(new DenseLayer<float>(128, 784));
    net32.add(new SigmoidLayer<float>());
    net32.add(new DenseLayer<float>(10, 128));
    net32.add(new tanhLayer<float>());

    Matrix<float> X32 = Matrix<float>::Random(784, 1);
    std::cout << "Float:\n" << net32.forward(X32).transpose() << std::endl;

    Sequential<double> net64(4);
    net64.add(new DenseLayer<double>(128, 784));
    net64.add(new SigmoidLayer<double>());
    net64.add(new DenseLayer<double>(10, 128));
    net64.add(new tanhLayer<double>());

    Matrix<double> X64 = Matrix<double>::Random(784, 1);
    std::cout << "Double:\n" << net64.forward(X64).transpose() << std::endl;
}

Verifica: entrambe compilano ed eseguono, e i due output sono vicini ma non identici — è l’errore di arrotondamento della singola precisione.

Confrontate anche la memoria occupata: Matrix<float> usa metà dei byte di Matrix<double>. Su una rete con 1500 neuroni nel primo layer sono milioni di parametri, e la differenza si sente.

Task da Completare

• Creare eigen_wrap.hpp con l’alias Matrix<ntype>
• Creare layer.hpp con is_training, outputs, delta_l, inputs_ptr
• Eliminare ilayer.hpp, dense_layer.*, activation_layer.*, sequential_network.*
• Creare parametricLayer.hpp con ParametricLayer (che ha gauss() e set_weights puro virtuale) e DenseLayer
• Implementare resize_weights() e set_weights() in DenseLayer
• Rinominare i pesi in Wt, aggiungere dWt, eliminare i bias
• Passare a int type, int distro e al costruttore (Nnodes, Nin)
• Creare activation_layers.hpp con match_shape, SigmoidLayer, tanhLayer
• Riscrivere sequential.hpp come Sequential<ntype> : layer<ntype>
• Far ritornare a forward() un Matrix<ntype>& e usare inputs_ptr
• Verificare che Sequential<float> e Sequential<double> compilino entrambi
• Trappola: provate a scrivere DenseLayer<float>(784, 128) e leggete l’eccezione

La soluzione di riferimento verrà discussa all’inizio della prossima lezione.