Свёрточная нейронная сеть с нуля. Часть 1. Свёрточный слой

В прошлой статье мы рассказывали о том, как устроена свёрточная нейронная сеть. В этой статье мы подробно разберём алгоритмы работы свёрточного слоя, напишем прямое распространение сигнала, обратное распространение ошибки с вычислением градиентов функции потерь по весам и по входам и научимся обновлять коэффициенты фильтров с помощью метода градиентного спуска.

В основном мы будем использовать код, написанный на С++, или псевдокод, однако это никак не должно помешать созданию аналогичных структур и алгоритмов на любом другом языке программирования.

Этап 0. Подготовительный

Прежде чем переходить к созданию свёрточного слоя, необходимо создать структуру для тензора. Тензор, как было сказано в прошлой статье, представляет из себя 3D массив. У тензора имеются три размерности, а значит нам потребуется структура для хранения размера тензора, мы назовём её TensorSize:

// размерность тензора
struct TensorSize {
    int depth; // глубина
    int height; // высота
    int width; // ширина
};

Создав размерность, можно создавать и сам тензор. В нём нам потребуется хранить depth·height·width вещественных чисел, которые мы будем хранить в одномерном векторе для уменьшения количества обращений к элементам по индексам. Для того, чтобы можно было работать с тензором, получать и изменять числа, необходима индексация. Индексация будет производиться по трём индексам в следующем порядке: глубина, высота, ширина. Поскольку в качестве хранения мы выбрали одномерный вектор, то потребуется применить несколько умножений и сложений.

// тензор
class Tensor {
    TensorSize size; // размерность тензора
    std::vector<double> values; // значения тензора

    int dw; // произведение глубины на ширину для индексации

    void Init(int width, int height, int depth);

public:
    Tensor(int width, int height, int depth); // создание из размеров
    Tensor(const TensorSize &size); // создание из размера

    double& operator()(int d, int i, int j); // индексация
    double operator()(int d, int i, int j) const; // индексация

    TensorSize GetSize() const; // получение размера

    friend std::ostream& operator<<(std::ostream& os, const Tensor &tensor); // вывод тензора
};

// инициализация по размерам
void Tensor::Init(int width, int height, int depth) {
    size.width = width; // запоминаем ширину
    size.height = height; // запоминаем высоту
    size.depth = depth; // запоминаем глубину

    dw = depth * width; // запоминаем произведение глубины на ширину для индексации

    values = std::vector<double>(width * height * depth, 0); // создаём вектор из width * height * depth нулей
}

// создание из размеров
Tensor::Tensor(int width, int height, int depth) {
    Init(width, height, depth);  
}

// создание из размера
Tensor::Tensor(const TensorSize &size) {
    Init(size.width, size.height, size.depth);
}

// индексация
double& Tensor::operator()(int d, int i, int j) {
    return values[i * dw + j * size.depth + d];
}

// индексация
double Tensor::operator()(int d, int i, int j) const {
    return values[i * dw + j * size.depth + d];
}

// получение размера
TensorSize Tensor::GetSize() const {
    return size;
}

// вывод тензора
std::ostream& operator<<(std::ostream& os, const Tensor &tensor) {
    for (int d = 0; d < tensor.size.depth; d++) {
        for (int i = 0; i < tensor.size.height; i++) {
            for (int j = 0; j < tensor.size.width; j++)
                os << tensor.values[i * tensor.dw + j * tensor.size.depth + d] << " ";
            
            os << std::endl;
        }

        os << std::endl;
    }

    return os;
}

Теперь, когда у нас в распоряжении есть тензор, можно приступать к созданию свёрточного слоя.

Этап 1. Свёрточный слой

Из прошлой статьи мы узнали, что свёрточный слой содержит в себе набор фильтров и смещений. Также для создания слоя необохдимо задать несколько гиперпараметров, а именно: fs — размер фильтров, fc — количество фильтров, S — шаг свёртки и P — дополнение нулями. Поскольку фильтры являются тензорами, а смещения всего лишь числами, то добавим фильтры в виде вектора тензоров и назовём W, а смещения добавим в виде векторов вещественных чисел и назовём их b. Поскольку эти веса являются обучаемыми, нам потребуется аналогичные наборы векторов для хранения градиентов.

Весовые коэффициенты должны заполняться случайными числами, так что добавим в класс генератор случайных чисел с нормальным распределением и метод инициализации весов. Для работы слоя также необходимо определить метод прямого распространения сигнала, метод обратного распространения ошибки и метод для обновления весовых коэффициентов:

class ConvLayer {
    std::default_random_engine generator; // генератор случайных чисел
    std::normal_distribution<double> distribution; // с нормальным распределением

    TensorSize inputSize; // размер входа
    TensorSize outputSize; // размер выхода

    std::vector<Tensor> W; // фильтры
    std::vector<double> b; // смещения

    std::vector<Tensor> dW; // градиенты фильтров
    std::vector<double> db; // градиенты смещений
    
    int P; // дополнение нулями
    int S; // шаг свёртки

    int fc; // количество фильтров
    int fs; // размер фильтров
    int fd; // глубина фильтров

    void InitWeights(); // инициализация весовых коэффициентов

public:
    ConvLayer(TensorSize size, int fc, int fs, int P, int S); // создание слоя

    Tensor Forward(const Tensor &X); // прямое распространение сигнала
    Tensor Backward(const Tensor &dout, const Tensor &X); // обратное распространение ошибки
    void UpdateWeights(double learningRate); // обновление весовых коэффициентов
};

Само создание слоя довольно простое: необходимо запомнить переданные параметры в полях класса и инициализировать весовые коэффициенты фильтров и смещений:

// создание свёрточного слоя
ConvLayer::ConvLayer(TensorSize size, int fc, int fs, int P, int S) : distribution(0.0, sqrt(2.0 / (fs*fs*size.depth))) {
    // запоминаем входной размер
    inputSize.width = size.width;
    inputSize.height = size.height;
    inputSize.depth = size.depth;

    // вычисляем выходной размер
    outputSize.width = (size.width - fs + 2 * P) / S + 1;
    outputSize.height = (size.height - fs + 2 * P) / S + 1;
    outputSize.depth = fc;

    this->P = P; // сохраняем дополнение нулями
    this->S = S; // сохраняем шаг свёртки

    this->fc = fc; // сохраняем число фильтров
    this->fs = fs; // сохраняем размер фильтров
    this->fd = size.depth; // сохраняем глубину фильтров

    // добавляем fc тензоров для весов фильтров и их градиентов
    W = std::vector<Tensor>(fc, Tensor(fs, fs, fd));
    dW = std::vector<Tensor>(fc, Tensor(fs, fs, fd));
        
    // добавляем fc нулей для весов смещения и их градиентов
    b = std::vector<double>(fc, 0);
    db = std::vector<double>(fc, 0);

    InitWeights(); // инициализируем весовые коэффициенты
}

Для заполнения весовых коэффициентов случайными числами просто проходимся по всем тензорам и присваиваем каждому элементу очередное случайное число:

// инициализация весовых коэффициентов
void ConvLayer::InitWeights() {
    // проходимся по каждому из фильтров
    for (int index = 0; index < fc; index++) {
        for (int i = 0; i < fs; i++)
            for (int j = 0; j < fs; j++)
                for (int k = 0; k < fd; k++)
                    W[index](k, i, j) = distribution(generator); // генерируем случайное число и записываем его в элемент фильтра

        b[index] = 0.01; // все смещения устанавливаем в 0.01
    }
}

Этап 2. Прямое распространение сигнала

При прямом распространении сигнала нам необходимо выполнить свёртку входного тензора с каждым из фильтров, а значит примерный алгоритм должен быть следующим:

for index = 0 to fc - 1 do
    output[index] = Convolve(X, W[index]) + b[index]

Теперь необходимо реализовать саму свёртку. Нам известно, что необходимо пройтись по всем возможным элементам на входном тензоре с заданным шагом, однако давайте сначала разберёмся, как действовать, когда параметр дополнения нулями отличен от нуля.

Дополнение нулями, padding

Самым простым решением дополнения нулями было бы создание тензора с высотой и шириной на 2P больше и последующим копированием во внутреннюю часть исходного тензора. Увы, подобные копирования сильно замедлят программу, так как требуют постоянных копирований памяти. Мы, конечно, хотим написать простую свёрточную сеть, однако всё-таки постараемся описать её наиболее эффективным способом.

Поскольку входные тензоры дополняются нулями, то элементы, выходящие за границу тензора, можно просто игнорировать. Для проверки на выход за границу потребуется всего навсего один if, что не идёт ни в какое сравнение с полным копированием исходного тензора в другой, дополненный нулями. Теперь же, когда дополнение нулями нам больше не помеха, можем написать прямое распространение:

// прямое распространение
Tensor ConvLayer::Forward(const Tensor &X) {
    Tensor output(outputSize); // создаём выходной тензор

    // проходимся по каждому из фильтров
    for (int f = 0; f < fc; f++) {
        for (int y = 0; y < outputSize.height; y++) {
            for (int x = 0; x < outputSize.width; x++) {
                double sum = b[f]; // сразу прибавляем смещение

                // проходимся по фильтрам
                for (int i = 0; i < fs; i++) {
                    for (int j = 0; j < fs; j++) {
                        int i0 = S * y + i - P;
                        int j0 = S * x + j - P;

                        // поскольку вне границ входного тензора элементы нулевые, то просто игнорируем их
                        if (i0 < 0 || i0 >= inputSize.height || j0 < 0 || j0 >= inputSize.width)
                            continue;

                        // проходимся по всей глубине тензора и считаем сумму
                        for (int c = 0; c < fd; c++)
                            sum += X(c, i0, j0) * W[f](c, i, j);
                    }
                }

                output(f, y, x) = sum; // записываем результат свёртки в выходной тензор
            }
        }
    }

    return output; // возвращаем выходной тензор
}

Этап 3. Обратное распространение ошибки

Цель обратного распространения ошибки заключается в том, что от следующего слоя нам приходят градиенты некоторой функции потерь L. Исходя из этих градиентов, нам нужно расчитать градиенты по входному тензору () и градиенты весовых коэффициентов ( и ).

Это становится вполне логичным, если рассматривать нейросеть как обыкновенную функцию с огромным количеством параметров, которую мы хотим минимизировать. А чтобы найти её минимум, нам нужно найти производные по каждому из весовых коэффициентов, чтобы затем шагнуть в направлении антиградиента (направление, в котором функция уменьшается). Поскольку нейросеть состоит из большого числа блоков, сходу тяжело предрассчитать производные, а потому на помощь приходит расчёт производной сложной функции: пусть у нас есть некоторая функция y = f(x) и функция g(y), и нам требуется найти производную по x функции z = g(f(x)). Для этого сначала нужно найти производную функции g, а затем умножить её на производную функции f: z' = g'(f(x)) * f'(x).

В нейросети в качестве функций выступают слои. Каждый из слоёв умеет вычислять производную по своему аргументу (входному тензору), основываясь на градиентах следующего слоя (или градиентах функции ошибки, если слой является выходным). Некоторые слои не имеют обучаемых параметров, поэтому они вычисляют только градиенты по своему аргументу. Свёрточный слой является обучаемым, а потому помимо входных градиентов ему необходимо вычислять и градиенты по весам, чтобы затем их можно было обновить и уменьшить ошибку.

Вычисление градиента по весовым коэффициентам

Вычисление градиента по весовым коэффициентам можно разбить на две части: по весам фильтров и по коэффициентам смещения. Проще всего вычисляются градиенты по весам смещения: это всего лишь сумма всех градиентов, которые относятся к данному весу смещения:

for (int f = 0; f < fc; f++)
    for (int y = 0; y < size.height; y++)
        for (int x = 0; x < size.width; x++)
            db[f] += deltas(f, y, x);

Для вычисления же градиентов по весам фильтров, необходимо выполнить слегка модифицированную свёртку входного тензора и тензора градиентов:

for (int f = 0; f < fc; f++) {
    for (int y = 0; y < size.height; y++) {
        for (int x = 0; x < size.width; x++) {
            double delta = deltas(f, y, x); // запоминаем значение градиента

            for (int i = 0; i < fs; i++) {
                for (int j = 0; j < fs; j++) {
                    int i0 = i + y - P;
                    int j0 = j + x - P;

                    // игнорируем выходящие за границы элементы
                    if (i0 < 0 || i0 >= inputSize.height || j0 < 0 || j0 >= inputSize.width)
                        continue;

                    // наращиваем градиент фильтра
                    for (int c = 0; c < fd; c++)
                        dW[f](c, i, j) += delta * X(c, i0, j0);
                }
            }
        }
    }
}

Вычисление градиента по входному тензору

Градиент по входному тензору считается обыкновенной свёрткой тензора дельт и весов с небольшим отличием: для вычисления фильтры поворачиваются на 180 градусов, а дополнение нулями заменяется на величину fs - 1 - P. Данную операцию называют транспонированной свёрткой или кросскорелляцией.

Tensor dX(inputSize);
int pad = fs - 1 - P;

for (int y = 0; y < inputSize.height; y++) {
    for (int x = 0; x < inputSize.width; x++) {
        for (int c = 0; c < fd; c++) {
            double sum = 0; // сумма для градиента

            // идём по всем весовым коэффициентам фильтров
            for (int i = 0; i < fs; i++) {
                for (int j = 0; j < fs; j++) {
                    int i0 = y + i - pad;
                    int j0 = x + j - pad;

                    // игнорируем выходящие за границы элементы
                    if (i0 < 0 || i0 >= size.height || j0 < 0 || j0 >= size.width)
                        continue;

                    // суммируем по всем фильтрам
                    for (int f = 0; f < fc; f++)
                        sum += W[f](c, fs - 1 - i, fs - 1 - j) * deltas(f, i0, j0); // добавляем произведение повёрнутых фильтров на дельты
                }
            }

            dX(c, y, x) = sum; // записываем результат в тензор градиента
        }
    }
}

Внимательный читаель мог заметить, что в заголовке метода использовался параметр dout, а в приведённом выше коде используется параметр deltas, а также неизвестный размер size. Действительно, deltas и dout не совсем одно и то же. Если шаг свёртки равен 1, то они совпадают, при другом шаге size и deltas вычисляются следующим образом:

TensorSize size; // размер дельт

// расчитываем размер для дельт
size.height = S * (outputSize.height - 1) + 1;
size.width = S * (outputSize.width - 1) + 1;
size.depth = outputSize.depth;

Tensor deltas(size); // создаём тензор для дельт

// расчитываем значения дельт
for (int d = 0; d < size.depth; d++)
    for (int i = 0; i < outputSize.height; i++)
        for (int j = 0; j < outputSize.width; j++)
            deltas(d, i * S, j * S) = dout(d, i, j);

Итоговый код для расчёта градиентов будет таким:

// обратное распространение
Tensor ConvLayer::Backward(const Tensor &dout, const Tensor &X) {
    TensorSize size; // размер дельт

    // расчитываем размер для дельт
    size.height = S * (outputSize.height - 1) + 1;
    size.width = S * (outputSize.width - 1) + 1;
    size.depth = outputSize.depth;

    Tensor deltas(size); // создаём тензор для дельт

    // расчитываем значения дельт
    for (int d = 0; d < size.depth; d++)
        for (int i = 0; i < outputSize.height; i++)
            for (int j = 0; j < outputSize.width; j++)
                deltas(d, i * S, j * S) = dout(d, i, j);

    // расчитываем градиенты весов фильтров и смещений
    for (int f = 0; f < fc; f++) {
        for (int y = 0; y < size.height; y++) {
            for (int x = 0; x < size.width; x++) {
                double delta = deltas(f, y, x); // запоминаем значение градиента

                for (int i = 0; i < fs; i++) {
                    for (int j = 0; j < fs; j++) {
                        int i0 = i + y - P;
                        int j0 = j + x - P;
                        
                        // игнорируем выходящие за границы элементы
                        if (i0 < 0 || i0 >= inputSize.height || j0 < 0 || j0 >= inputSize.width)
                            continue;

                        // наращиваем градиент фильтра
                        for (int c = 0; c < fd; c++)
                            dW[f](c, i, j) += delta * X(c, i0, j0);
                    }
                }

                db[f] += delta; // наращиваем градиент смещения
            }
        }
    }

    int pad = fs - 1 - P; // заменяем величину дополнения
    Tensor dX(inputSize); // создаём тензор градиентов по входу

    // расчитываем значения градиента
    for (int y = 0; y < inputSize.height; y++) {
        for (int x = 0; x < inputSize.width; x++) {
            for (int c = 0; c < fd; c++) {
                double sum = 0; // сумма для градиента

                // идём по всем весовым коэффициентам фильтров
                for (int i = 0; i < fs; i++) {
                    for (int j = 0; j < fs; j++) {
                        int i0 = y + i - pad;
                        int j0 = x + j - pad;

                        // игнорируем выходящие за границы элементы
                        if (i0 < 0 || i0 >= size.height || j0 < 0 || j0 >= size.width)
                            continue;

                        // суммируем по всем фильтрам
                        for (int f = 0; f < fc; f++)
                            sum += W[f](c, fs - 1 - i, fs - 1 - j) * deltas(f, i0, j0); // добавляем произведение повёрнутых фильтров на дельты
                    }
                }

                dX(c, y, x) = sum; // записываем результат в тензор градиента
            }
        }
    }

    return dX; // возвращаем тензор градиентов
}

Этап 4. Обновление весовых коэффициентов

После того, как градиенты вычислены, необходимо обновить весовые коэффициенты, чтобы уменьшить ошибку на текущем обучающем примере. Для этого из каждого веса вычитается значение градиента, умноженное на скорость обучения — вещественное число, гиперпараметр метода градиентного спуска. Данный алгоритм является наиболее простым и довольно медленным с точки зрения сходимости методом оптимизации, однако ничего не мешает реализовать более быстрые алгоритмы, как, например, Adam, Adagrad, RMSprop и другие. В одной из статей мы обязательно расскажем, как их реализовать и добавить к создаваемой нами сети, а пока посмотрим на код обновления параметров:

// обновление весовых коэффициентов
void ConvLayer::UpdateWeights(double learningRate) {
    for (int index = 0; index < fc; index++) {
        for (int i = 0; i < fs; i++) {
            for (int j = 0; j < fs; j++) {
                for (int d = 0; d < fd; d++) {
                    W[index](d, i, j) -= learningRate * dW[index](d, i, j); // вычитаем градиент, умноженный на скорость обучения
                    dW[index](d, i, j) = 0; // обнуляем градиент фильтра
                }
            }
        }

        b[index] -= learningRate * db[index]; // вычитаем градиент, умноженный на скорость обучения
        db[index] = 0; // обнуляем градиент веса смещения
    }
}

Этап 5. Проверка

После написания всех методов хотелось бы удостовериться, что они работают корректно. Логичнее всего было бы написать небольшую процедуру для тестирования слоя, однако на данный момент у нас нет возможности задавать свои собственные параметры для весовых коэффициентов, а проверять работу слоя на случайных числах крайне неприятно. Поэтому мы создадим два метода для задания значений весовых коэффициентов фильтров и смещений соответственно:

// установка веса фильтра по индексу
void ConvLayer::SetWeight(int index, int i, int j, int k, double weight) {
    W[index](i, j, k) = weight;
}

// установка веса смещения по индексу
void ConvLayer::SetBias(int index, double bias) {
    b[index] = bias;
}

После этого напишем небольшую тестирующую функцию, которая проверит значения после применения прямого распространения и обратного распространения. Мы не будем приводить её код ввиду простоты и большого объёма проверок, но покажем значения, на которых будет производиться проверка.

Для начала создадим свёрточный слой, который будет принимать на вход тензор размером 5x5x3 (три канала 5x5), содержать 2 фильтра размером 3x3, шаг смещения 2 и единичное дополнение нулями. На выходе мы ожидаем тензор размером (5-3+2*1)/2+1 x (5-3+2*1)/2+1 x 2 = 3x3x2. Возьмём следующие фильтры и входной тензор для проверки:

      |-1 1 1 |  |  0 -1  1 |  |  1 -1  0 |
W[0]: |-1 1 1 |  | -1  0 -1 |  | -1  0 -1 |
      | 0 0 1 |  |  1  0  0 |  |  0  1 -1 |

      | 0  0 -1 |  |  1 -1  0 |  |  1  0  1 |
W[1]: | 1  0  0 |  | -1  1  0 |  |  1 -1  1 |
      | 0 -1  0 |  | -1  1  1 |  | -1  0  0 |

b[0]: 1
b[1]: 0

   | 1 2 0 1 0 |  | 1 2 2 1 2 |  | 0 2 0 1 1 |
   | 2 0 0 0 1 |  | 0 2 2 0 2 |  | 2 0 2 1 1 |
X: | 1 2 2 0 2 |  | 1 2 2 1 1 |  | 2 0 1 2 1 |
   | 2 2 2 0 1 |  | 2 2 0 1 0 |  | 0 1 0 1 2 |
   | 2 0 1 0 1 |  | 2 2 1 0 0 |  | 1 0 1 2 1 |        

Если всё сделано правильно, то на выходе должен получиться следующий тензор:

    | 2 -3 1 |  | 3 5  2 |
Y:  | 5 -7 2 |  | 1 0 -3 |
    | 5  0 0 |  |-2 4  3 |

Пусть теперь слой принимает тензоры размером 4x4x1 и содержит 1 фильтр 3х3 с единичным шагом без дополнения нулями.

      | 1 4 1 |
W[0]: | 1 4 3 |
      | 3 3 1 |

b[0]: 0

   | 4 5 8 7 |
X: | 1 8 8 8 |
   | 3 6 6 4 |
   | 6 5 7 8 |

Результат прямого распространения должен быть таким:

Y: | 122 148 |
   | 126 134 |

Запустим же теперь обратное распространение с такими дельтами:

dout: | 2 1 |
      | 4 4 |

В результате должны получить такой тензор:

    |  2  9  6  1 |
dX: |  6 29 30  7 |
    | 10 29 33 13 |
    | 12 24 16  4 |

Интерактивный пример

Чтобы операция свёртки стала ещё более понятной, мы добавили интерактивный пример ниже. В нём вы можете задать гиперпараметры слоя и увидеть анимацию прохода фильтра по тензорам.

Итоги

В результате мы получили класс для работы с тензорами, и, что самое главное, описали класс для работы свёрточного слоя с различными гиперпараметрами, а также проверили корректность работы методов прямого и обратного распространения. В следующей части мы расскажем вам о том, как создать слои подвыборки (пулинга) различного типа (максимума, среднего или суммы).

Следующая часть: Свёрточная нейронная сеть с нуля. Часть 2. Слой подвыборки