Давайте попробуем немного разобраться в теме импульсных нейронных сетей - ИНС (англ. spiking neural network - SNN). Напишем простую импульсную нейронную сеть, используя NumPy и Pandas, для классической задачи машинного обучения с использованием кодирования рецептивными полями

Пара слов вместо введения

Эту статью я написал для тех, кто только начинает знакомиться с импульсными нейронными сетями. Когда я самостоятельно начинал разбираться в этой теме, мне встречались примеры кода и статьи подобные этой, но с одним значимым отличием - практически всегда авторы использовали готовые решения типа PyNN или snnTorch. Мне хотелось лучше разобраться в принципах работы ИНС, поэтому я написал код, который сможет повторить и понять каждый, кто хочет разобраться в данной теме. Надеюсь эта статья найдет своего читателя и будет полезна.

Также статья доступна на английском языке тут.

С чего начать?

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

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

Первый подход к снаряду

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

Можно привести несколько примеров моделей нейронов, используемых в ИНС: одна из самых простых и популярных моделей - это, так называемая, модель LIF - нейрона (leaky integrate-and-fire neuron) или модель нейрона – порогового интегратора с утечкой, более продвинутая - модель нейрона Ижикевича, возможно самая подробная - модель нейрона Ходжкина-Хаксли - одна из самых реалистичных моделей биологического нейрона, однако редко используется для вычислительных экспериментов с ИНС (особенно большими) из-за их высокой вычислительной сложности.

Первая и возможно главная особенность ИНС - использование спайков в качестве основного средства передачи информации между нейронами и момента их генерации как единственного атрибута. Функционирование нейрона в этой парадигме ограничено только генерацией и приемом спайков. Момент генерации спайка нейроном зависит от спайков, полученных от других нейронов через специальные связи - синапсы.
Вторая особенность - это асинхронность ИНС: функционирование разных нейронов не является явно синхронизированным процессом.

Я кратко вернусь к некоторым теоретическим моментам ниже, так что не волнуйтесь, если что-то не очень понятно сейчас.

Какую задачу решаем?

Я предлагаю оригинальное решение одной из самых известных задач машинного обучения - Iris Species - с использованием импульсной нейронной сети. Я буду использовать одну из самых популярных моделей нейронов - модель LIF - нейрона. Для преобразования данных я буду использовать фазовое кодирование с использованием гауссовых рецептивных полей (англ. Gaussian receptive fields). Я не буду сильно углубляться в теоретическую часть и постараюсь передать суть подхода максимально кратко. Список используемых мною источников в конце статьи поможет более глубоко погрузиться в тему всем желающим.

Представленный код и преобразованный набор исходных данных доступны на Kaggle (код, датасет) и GitHub (код и датасет).

Наша простейшая ИНС

В общем, наша простая сеть будет иметь следующую архитектуру:

Набор данных содержит три типа цветов: Setosa, Versicolor и Virginica.
Каждый тип содержит 4 характеристики: “Sepal Length”, “Sepal Width”, “Petal Length”, “Petal Width”.
Цель классификации цветов ириса - предсказать тип цветов на основе их конкретных характеристик - уверен, большинству прекрасна знакома эта задачка.

Мы закодируем значения каждой характеристики с помощью 10 пресинаптических нейронов, всего их будет 40. В выходном слое у нас будет 3 постсинаптических нейрона для каждого типа цветка.

Ещё раз, но чуть более подробно, как это работает?

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

Основная суть модели: пресинаптические нейроны генерируют спайки (входная информация), которые стимулируют постсинаптические нейроны -> мембранный потенциал (Vm) постсинаптических нейронов повышается до порогового значения -> постсинаптические нейроны генерируют собственные спайки (выходная информация). Тип цветка, которому соответствует первый сработавший (испустивший спайк) постсинаптический нейрон на заданном входном слое (закодированном наборе из четырех значений каждой характеристики) или еще можно сказать на заданном периоде (конкретно в моем случае он равен 10 мс, но, по-сути, это сейчас не очень принципиально) является предсказанным типом цветка нашей импульсной сетью. Суть понятия период в данном контексте - это время в течение которого мы "показываем" нашей сети закодированный набор данных и в течение которого должен сработать один из постсинаптических нейронов, определив тип цветка. О том, как закодировать числовые значения входных данных в моменты генерации спайков пресинаптическими нейронами будет описано ниже.

Давайте рассмотрим небольшой пример визуализации одного периода обучения/тестирования, соответствующего одному набору входных данных (в случае набора данных Iris - одному из 150). Посмотрим на зависимость мембранного потенциала постсинаптических нейронов от входных спайков пресинаптических нейронов с учетом их весов. Если мембранный потенциал достигает порогового значения, постсинаптический нейрон генерирует спайк:

Теперь еще раз попробуем осмыслить все это, но с новыми подробностями, а затем начнем писать код.

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

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

Отлично! Мы рассмотрели все основы. Теперь перейдем к первой части, которая ответит на вопрос: как преобразовать числовые значения исходного набора данных в моменты спайков пресинаптических нейронов, а если точнее, то в задержки при генерации спайков.

Кодирование с использованием гауссовых рецептивных полей

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

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
import warnings
warnings.filterwarnings("ignore")

Загрузите исходный набор данных из любого удобного источника, по аналогии с ноутбуком из Kaggle, можно сделать это так:

URL = '/kaggle/input/iris/Iris.csv'
df = pd.read_csv(URL)
df = df.iloc[:,1:]

Посмотрим на первые несколько строк нашего датасета:

df.head()

Сделаем копию датасета, оставив только количественную часть детасета:

df_ = df.drop(columns=['Species']).copy()

Постройте гистограмму данных:

df_.plot.hist(alpha = 0.4, figsize = (12, 4))
plt.legend(title = "Dataset cilumns:" ,bbox_to_anchor = (1.0, 0.6),
                                                   loc = 'upper left')
plt.title('Iris dataset', fontsize = 20)
plt.xlabel('Input value', fontsize = 15)
plt.show()

Напомню еще раз, что мы закодируем каждую из четырех характеристик (“Sepal Length”, “Sepal Width”, “Petal Length”, “Petal Width”) с помощью 10 пресинаптических нейронов, всего будет 40 пресинаптических нейронов. Давайте напишем функцию, которая генерирует 10 гауссиан для каждой входной характеристики так, чтобы:

1. средние значения каждой гауссианы равномерно распределялись между крайними значениями диапазона, включая границы (минимальное и максимальное значение) для каждой характеристики;

2. высота каждой гауссианы равна 1, что является максимальным значением "возбуждения" пресинаптического нейрона, на основе которого мы затем рассчитаем задержку генерации спайка для каждого пресинаптического нейрона;

3. из-за отсутствия явной концентрации в исходных данных для каждой характеристики ширина каждой гауссины должна быть одинаковой.

def Gaus_neuron(df, n, step, s):

    neurons_list = list()
    x_axis_list = list()
    t = 0

    for col in df.columns:

        vol = df[col].values
        min_ = np.min(vol)
        max_ = np.max(vol)
        x_axis = np.arange(min_, max_, step)
        x_axis[0] = min_
        x_axis[-1] = max_
        x_axis_list.append(np.round(x_axis, 10))
        neurons = np.zeros((n, len(x_axis)))

        for i in range(n):

            loc = (max_ - min_) * (i /(n-1)) + min_
            neurons[i] = norm.pdf(x_axis, loc, s[t])
            neurons[i] = neurons[i] / np.max(neurons[i])

        neurons_list.append(neurons)
        t += 1

    return neurons_list, x_axis_list

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

sigm = [0.1, 0.1, 0.2, 0.1]
d = Gaus_neuron(df_, 10, 0.001, sigm)

Теперь визуализируем наши гауссианы для каждой характеристики:

fig, (ax1, ax2, ax3, ax4) = plt.subplots(4)

fig.set_figheight(8)
fig.set_figwidth(10)

k = 0

for ax in [ax1, ax2, ax3, ax4]:

    ax.set(ylabel = f'{df_.columns[k]} \n\n "Возбуждение" \nнейрона')

    for i in range(len(d[0][k])):

        ax.plot(d[1][k], d[0][k][i], label = i + 1)

    k+=1

plt.legend(title = "Номер пресинаптического нейрона \n         каждой характеристики:" ,bbox_to_anchor = (1.05, 3.25), loc = 'upper left')
plt.suptitle(' \n\n  Гауссианы для кодирования каждой характеристики', fontsize = 15)
ax.set_xlabel(' Диапозон входных данных\n каждой характеристики', fontsize = 12, labelpad = 15)


plt.show()

Подробно рассмотрим логику кодирования, используя в качестве примера первые пять точек данных характеристики «SepalWidthCm»: мы нарисуем пунктирные вертикальные сегменты, соответствующие первым пяти значениям характеристики «SepalWidthCm» и найдем их пересечения с гауссианами пресинаптических нейронов. Эти точки будут отмечены красным:

x_input = 5
fig, ax = plt.subplots(1)

fig.set_figheight(5)
fig.set_figwidth(15)

ax.set(ylabel = df_.columns[1])

for i in range(len(d[0][1])):
    ax.plot(d[1][1], d[0][1][i])

for n in range(x_input):

    plt.plot(np.tile(df_['SepalWidthCm'][n], (10,1)), 
         d[0][1][np.tile(d[1][1] == df_['SepalWidthCm'][n], (10,1))], 'ro', markersize=4)

    plt.vlines(x = df_['SepalWidthCm'][n], ymin = - 0.1, ymax = 1.1, 
               colors = 'purple', ls = '--', lw = 1, label = df_['SepalWidthCm'][n])

    plt.text(df_['SepalWidthCm'][n] * 0.997, 1.12, n + 1, size = 10)


plt.legend(title = "Первые пять значений:", bbox_to_anchor = (1.0, 0.7), loc = 'upper left')

plt.suptitle('Гауссианы для кодирования. \n \
                Подробное описание идеи на примере первых пяти значений "SepalWidthCm"',
            fontsize = 15)

ax.set_xlabel('Входное значение X ∈ [x_min, x_max] характеристики', fontsize = 12, labelpad = 15)
ax.set_ylabel('"Возбуждение" нейрона ∈ [0,1]', fontsize = 12, labelpad = 15)

plt.show()

Теперь выведим числовые значения точек пересечения каждого вертикального сегмента с каждой гауссианой:

np.set_printoptions(formatter={'float_kind':'{:f}'.format})
five_x = np.zeros((5, 10))

for n in range(x_input):
    five_x[n,:] = d[0][1][np.tile(d[1][1] == df_['SepalWidthCm'][n], (10,1))]

five_x

Найдем задержки каждого пресинаптического нейрона как 1 — (возбуждение пресинаптического нейрона) при условии, что возбуждение нейрона больше 0.1, в противном случае считаем пресинаптический нейрон невозбужденным на данной итерации:

five_x = np.where(five_x > 0.1, 1 - five_x, np.nan)
five_x[five_x == 0] = 0.0001
five_x

Теперь визуализируем процесс возникновения спайка с учетом задержек. Черные точки на графике – моменты генерации спайка пресинаптическими нейронами:

fig, ax = plt.subplots(5, figsize=(10, 8))

for i in range(5):
    ax[i].scatter(x = five_x[i], y = np.arange(1, 10 + 1), s = 10, color = 'black')
    ax[i].hlines(xmin = 0, xmax=1, y=np.arange(1, 11, 1),
               colors = 'purple', ls = '--', lw = 0.25)
    ax[i].yaxis.set_ticks(np.arange(0, 11, 1))
    ax[i].set_ylabel(f'x{i+1} = {df_.iloc[i,1]}\n (период{i+1}) \n\n № \nпресинаптического\nнейрона', fontsize = 7)
    ax[i].set_xlim(0, 1)
    ax[i].set_ylim(0, 10 * 1.05)
    ax[i].tick_params(labelsize = 7)

ax[i].set_xlabel('Задержка пресинаптического нейрона')
plt.suptitle(' \n\nВходные данные после кодирования \nс помощью гауссовых рецептивных полей', fontsize = 12)
plt.show()

Теперь проделаем это для всех значений всех характеристик:

def Lat_Spike(df, d, n):

    for i in range(len(df.columns)):

        k = len(df.iloc[:, i])
        st1 = np.tile(d[1][i], (k, 1))
        st2 = df.iloc[:, i].values.reshape(-1, 1)
        ind = (st1 == st2)
        exc = np.tile(d[0][i], (k, 1)).reshape(k, n, len(d[0][i][0]))[
            np.repeat(ind, n, axis=0).reshape(k, n, len(ind[0]))].reshape(k, n)
        lat_neuron = np.transpose(np.where(exc > 0.1, 1 - exc, np.nan))

        if i == 0:
            lat_neuron_total = lat_neuron
        else:
            lat_neuron_total = np.concatenate((lat_neuron_total, lat_neuron), axis = 0)

    lat_neuron_total[lat_neuron_total == 0] = 0.0001

    return lat_neuron_total

fin = Lat_Spike(df_, d, 10)

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

Отлично, теперь представим результаты в виде финального фрейма данных.
По строкам — пресинаптические нейроны, по столбцам — номер набора входных данных, значения внутри - задержки:

Final_df = pd.DataFrame(fin)
Final_df

Визуализируем полученные задержки для первых трех периодов обучения. Считаем один период обучения равным 10 мс, поэтому масштабируем задержку каждого пресинаптического нейрона (суммарный временной интервал по трем периодам равен 30 мс):

fig, ax = plt.subplots(1, figsize=(10, 6))
h = 3

for i in range(h):
    ax.scatter(x = (i+Final_df.iloc[:,i].values)*10, y = np.arange(1, 41), s = 6, color = 'black')

    plt.vlines(x = (i) * 10, ymin = 0, ymax = 40, 
               colors = 'purple', ls = '--', lw = 0.5)
    ax.tick_params(labelsize = 7)

ax.yaxis.set_ticks(np.arange(1, 41, 1))
ax.xaxis.set_ticks(np.arange(0, (h+1)*10, 10))
ax.set_xlabel('Время (мс)')
ax.set_ylabel('№ пресинаптического нейрона')
plt.suptitle(' \n\nМоменты спайков пресинаптических нейронов втечение первых 30 мс', fontsize = 10)
plt.gca().invert_yaxis()
plt.show()

Отлично! Мы успешно закодировали наш набор данных!
Теперь перейдем к обучению нашего LIF-нейрона!

Модель LIF-нейрона

Часть первая. Размер подвыборки 60: первые 20 значений каждого типа цветка.

Данные в исходном наборе распределены последовательно по типам ирисов: всего 150 наборов, из которых первые 50 записей принадлежат к Iris-setosa, 50–100 — Iris-versicolor, 100–150 — Iris-virginica. Я написал функцию, которая выберет нужный нам набор данных из портфеля, в котором будет равное количество экземпляров каждого типа:

def model_data(ind, ind_type, lat_ne, start, end):
    
    train_stack = np.vstack((lat_ne[ind_type[ind, 0] + start:ind_type[ind, 0] + end],
                            lat_ne[ind_type[ind, 1] + start:ind_type[ind, 1] + end],
                            lat_ne[ind_type[ind, 2] + start:ind_type[ind, 2] + end]))
    train_stack = np.where(train_stack > 0, train_stack, 0)
    
    return train_stack

На первом этапе нам необходимо каким-то образом увеличить веса синапсов активных пресинаптических нейронов, сформировав набор весов для каждого постсинаптического нейрона для дальнейшей корректировки с помощью STDP (Spike-time-dependable Plasticity) на втором этапе.

Как мы можем это сделать?
Вариантов может быть много, я буду это делать следующим образом:

1. Все веса изначально равны 0.1;

2. Для каждого постсинаптического нейрона в период обучения, соответствующий его типу цветка, увеличиваем веса всех активных пресинаптических нейронов на одинаковую для всех константу, скорректированную в зависимости от задержки конкретного сработавшего пресинаптического нейрона: меньше задержка -> пресинаптический нейрон сгенерировал спайк раньше -> получает больший вес. Общая формула выглядит следующим образом: delta_weight = +Const * (1 — latency);

3. Для каждого постсинаптического нейрона в период обучения, не соответствующий его типу цветка, мы накладываем штраф на все активные пресинаптические нейроны, имеющие вес, отличный от значения по умолчанию (вес > 0.1). То есть мы «наказываем» пресинаптические нейроны, которые активируются для различных типов цветов, уменьшая их вес -> уменьшая их вклад в стимуляцию мембранного потенциала постсинаптических нейронов других типов цветов. Аналогично пункту 2 размер штрафа равен той же константе с поправкой на задержку активного нейрона в данный период обучения. Общая формула выглядит следующим образом: delta_weight = -Const * (1 — latency);

4. «Штраф» из пункта 3 равномерно распределяется среди неактивных в этот период пресинаптических нейронов, имеющих вес, отличный от значения по умолчанию (вес > 0.1). Общая формула выглядит следующим образом: delta_weight = -Const * (1 — latency) / N, где latency — задержка сработавшего пресинаптического нейрона, N — количество «молчащих» пресинаптических нейронов с весом больше, чем по умолчанию;

5. Скорректированные веса не могут быть меньше значения по умолчанию: если вес уменьшается, то 0.1 является нижним пределом. У нас нет «тормозящих» пресинаптических нейронов с отрицательными весами связей;

6. В каждой итерации (каждом периоде) «молчащие» нейроны с дефолтными весами никак не корректируются.

Фу! Если я вас немного запутал, картинка ниже поможет разобраться:

Напишем код:

lat_ne = np.transpose(Final_df.values)
ind_type = np.array(([0, 50, 100], [50, 100, 0], [100, 0, 50]))
list_weight = np.zeros((3,40))

for ind in range(3):
    
    train_stack = model_data(ind, ind_type, lat_ne, 0, 20)
    tr_ar = np.where(np.transpose(train_stack) > 0, 2 * (1 - np.transpose(train_stack)), 0)
    tr_ar[:, 20:] = tr_ar[:, 20:] * (-1)
    tr_ar = pd.DataFrame(tr_ar)
    tr_ar[20] = tr_ar.iloc[:,:20].sum(axis = 1) + 0.1
    tst_ar = np.float64(np.transpose(np.array(tr_ar.iloc[:,20:])))
    
    for i in range(1, len(tst_ar)):
        
        tst_ar[0][((np.round(tst_ar[0], 4) > 0.1) & (tst_ar[i] == 0))] += - np.float64(
            np.sum(tst_ar[i][np.round(tst_ar[0], 4) > 0.1]) / len(tst_ar[0][((
                np.round(tst_ar[0], 4) > 0.1) & (tst_ar[i] == 0))]))
        tst_ar[0][np.round(tst_ar[0], 4) > 0.1] += tst_ar[i][np.round(tst_ar[0], 4) > 0.1]
        tst_ar[0][tst_ar[0] < 0.1] = 0.1
        
    list_weight[ind, :] = tst_ar[0]

list_weight

Мы получили наш первый набор весов! Давайте напишем функцию для нашего нейрона, которая учитывает синаптические веса и время моментов генерации спайков пресинаптическими нейронами. Затем посмотрим, как изменяется мембранный потенциал каждого постсинаптического нейрона.

Но сначала давайте введем формулу для изменения мембранного потенциала постсинаптического нейрона. Если говорить просто: спайк с учетом веса синапса приходит на постсинаптический нейрон -> мембранный потенциал постсинаптического нейрона увеличивается, если спайки не приходят, потенциал падает до своего минимального уровня [Vmin] (в нашем случае это 0) - эта логика видна из следующей формулы:

- эта формула является слегка измененной и упрощенной версией формулы из [12], но вся логика классического LIF-нейрона здесь сохранена. В моей формуле отсутствует потенциал покоя (Vmin), я предполагаю его равным 0.
Параметры tau, dt, Vmin могут быть подобраны для приближения модели к более близкой к биологической модели нейрона. В данном случае, у меня такой задачи нет, поэтому я выбрал параметры такими, чтобы результаты были максимально понятными, а модель простая. Теперь давайте реализуем эту логику в коде:

def LIF_SNN(n, l, data, weight, v_spike):
    
    V_min = 0
    V_spike = v_spike
    r = 5
    tau = 2.5
    dt = 0.01
    t_max = 10
    time_stamps = t_max / dt
    time_relax = 10
    v = np.zeros((n, l, int(time_stamps)))
    t_post = np.zeros((n, l))
    t_post_ = np.zeros((n, int(l / 3)))
    v[:, :, 0] = V_min
    
    for n in range(n):
        for u in range(l):
            
            t = 0
            f0 = (np.round(data[u][np.newaxis].T, 3) * 1000).astype(int)
            f1 = np.tile(np.arange(1000), (40, 1))
            f2 = np.where(((f1 == f0) & (f0 > 0)), 1, 0)
            f2 = f2 * weight[n][np.newaxis].T
            spike_list = np.sum(f2.copy(), axis = 0)

            for step in range(int(time_stamps) - 1):
                if v[n, u, step] > V_spike:
                    t_post[n, u] = step
                    v[n, u, step] = 0
                    t = time_relax / dt
                elif t > 0:
                    v[n, u, step] = 0
                    t = t - 1

                v[n, u, step + 1] = v[n, u, step] + dt / tau * (-v[n, u, step] + r * spike_list[step])
        t_post_[n, :] = t_post[n, n * int(l / 3):n * int(l / 3) + int(l / 3)]
    
    return v, t_post_, t_post

Функция визуализации моментов спайков постсинаптическими нейронами:

def spike_plot(spike_times, one_per, n, cur_type):
    
    fig, (ax1, ax2, ax3) = plt.subplots(3, figsize = (25, 10))#, dpi = 70)
    
    if one_per:
        k, t, a  = 1, n, 0
        cur = cur_type
    else:
        k, t, a = len(spike_times[0]), 0, 1
        cur = 1
        
    spike_times[spike_times == 0] = np.nan
    di = {0: 'blue', 1: 'red', 2: 'black'}
    di_t = {0: 'Iris-setosa', 1: 'Iris-versicolor', 2: 'Iris-virginica'}
    p = 0
    
    for ax in [ax1, ax2, ax3]:
        for i in range(k * t, k + t):
            ax.vlines(x = spike_times[p, i] / 100 + i * a * 10, ymin = 0.0, ymax = 1.1, 
                       colors = di[p], ls = '-', lw = 3)
            ax.set_ylabel(f'Нейрон {p + 1} \n {di_t[p]}', fontsize = 15)
            
        if one_per:
            ax.axvspan(0, int(k * 10), color = di[cur - 1], alpha = 0.05, label = di_t[cur - 1])
            ax.margins(0)
        else:
            ax.axvspan(0, int(k * 10 / 3), color = di[0], alpha = 0.05, label = di_t[0])
            ax.axvspan(int(k * 10 / 3), int(k * 10 * 2 / 3), color = di[1], alpha = 0.05, label = di_t[1])
            ax.axvspan(int(k * 10 * 2 / 3), int(k * 10 * 3 / 3), color = di[2], alpha = 0.05, label = di_t[2])
            ax.set_xlim(0, k * 10)
            ax.margins(0)
            
        p += 1
        
    
    if one_per:
        plt.suptitle(f' \n\n Моменты генерации постсинаптическими нейронами спайков в течение периода {n}', fontsize = 20)
        plt.legend(title = "    Часть датасета,\n соответствующая типу:" ,bbox_to_anchor = (1, 1.9), loc = 'upper left',
               fontsize = 15, title_fontsize = 15)
    else:
        plt.suptitle(f' \n\n Моменты генерации постсинаптическими нейронами спайков на используемой части датасета', fontsize = 20)
        plt.legend(title = "    Часть датасета,\n соответствующая типу:" ,bbox_to_anchor = (1, 2.1), loc = 'upper left',
               fontsize = 15, title_fontsize = 15)
    
    plt.xlabel('Время (мс)', fontsize = 15)
    plt.show()

Функция визуализации уровня мембранного потенциала каждого постсинаптического нейрона:

def v_plot(v):
    
    fig, (ax1, ax2, ax3) = plt.subplots(3, figsize = (25, 10))#, dpi = 70)
    k = len(v[0,:,:])
    di = {0: 'blue', 1: 'red', 2: 'black'}
    di_t = {0: 'Iris-setosa', 1: 'Iris-versicolor', 2: 'Iris-virginica'}
    p = 0
    
    for ax in [ax1, ax2, ax3]:
        for i in range(k):
            ax.plot(np.arange(i * 10, (i + 1) * 10, 0.01), v[p, i, :], di[p], linewidth = 1)
            ax.set_ylabel(f' Нейрон {p + 1} \n {di_t[p]} \nV (mV)', fontsize = 15)

        ax.axvspan(0, int(k * 10 / 3), color = di[0], alpha = 0.05, label = di_t[0])
        ax.axvspan(int(k * 10 / 3), int(k * 10 * 2 / 3), color = di[1], alpha = 0.05, label = di_t[1])
        ax.axvspan(int(k * 10 * 2 / 3), int(k * 10 * 3 / 3), color = di[2], alpha = 0.05, label = di_t[2])
        ax.margins(0)

        p += 1
    
    plt.legend(title = "    Часть датасета,\n соответствующая типу:" ,bbox_to_anchor = (1, 2), loc = 'upper left', fontsize = 15, title_fontsize = 15)
    plt.xlabel('Время (мс)', fontsize = 15)
    plt.suptitle(' \n Активность постсинаптических нейронов на используемой части датасета \n (Мембранный потенциал)', fontsize = 20)

Функция точности. Главный принцип: если несколько постсинаптических нейронов генерируют спайки в течение одного периода, считается, что постсинаптический нейрон, который первым сгенерировал спайк, был единственно активным:

def accuracy_snn(spike_time, start, end, df, ind_type, ind):
    
    type_dict = {'Iris-setosa': 1, 'Iris-versicolor': 2, 'Iris-virginica': 3}
    target_type_total = np.array(df.replace({'Species': type_dict}).iloc[:, - 1])
    target_type = np.vstack((target_type_total[ind_type[ind, 0] + start:ind_type[ind, 0] + end],
                            target_type_total[ind_type[ind, 1] + start:ind_type[ind, 1] + end],
                            target_type_total[ind_type[ind, 2] + start:ind_type[ind, 2] + end])).flatten()
    
    spike_time_ = np.where(spike_time > 0, np.array(([1], [2], [3])), np.nan)
    final_test = np.full([len(spike_time[0])], np.nan).astype(int)
    for i in range(len(spike_time[0])):
        try:
            final_test[i] = spike_time_[:, i][spike_time[:, i] == np.min(spike_time[:, i][spike_time[:, i] > 0])][0]
        except:
            final_test[i] = 0
    
    ac = np.sum(np.where(final_test == target_type, 1, 0)) / len(target_type)

    return final_test, target_type, print('accur.:', np.round(ac * 100, 2), '%')

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

train_stack = model_data(0, ind_type, lat_ne, 0, 20)
res = LIF_SNN(3, 60, train_stack, list_weight, 100)
v = res[0]

v_plot(v)

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

res = LIF_SNN(3, 60, train_stack, list_weight, 0.25)
spike_time = res[2]
spike_plot(spike_time, False, False, False)
accuracy_snn(spike_time, 0, 20, df, ind_type, 0)[2]

Точность 93.33% действительно неплохая! Давайте рассмотрим несколько периодов, когда один из постсинаптических нейронов имеет "ложную" активацию. Попытаемся понять, что происходит и как это влияет на точность.

Посмотрим на последний ложный спайк первого постсинаптического нейрона, который возникает в 46 периоде первой части обучения:

spike_plot(spike_time, True, 46, 3)

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

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

spike_plot(spike_time, True, 24, 2)

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

В целом дела идут неплохо, так что плавно переходим ко второй части обучения.

Часть вторая. Размер подвыборки 60: вторые 20 значений для каждого типа цветка.

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

train_stack = model_data(0, ind_type, lat_ne, 20, 40)
res = LIF_SNN(3, 60, train_stack, list_weight, 100)
v = res[0]

v_plot(v)
res = LIF_SNN(3, 60, train_stack, list_weight, 0.25)
spike_time = res[2]
spike_plot(spike_time, False, False, False)
accuracy_snn(spike_time, 20, 40, df, ind_type, 0)[2]

В целом профиль аналогичен профилю первой обучающей части. Немного снизилась точность до 91.67% — ожидаемый результат.

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

Чтобы применить STDP, мы будем записывать время моментов спайков каждого постсинаптического нейрона при пороговом значении 0.25, получая значения t_post. У нас уже есть значения t_pre (задержки), которые мы рассчитали ранее. Мы рассчитаем изменение веса каждого пресинаптического нейрона для каждого постсинаптического нейрона, используя следующие формулы:

В зависимости от значения задержки, веса пресинаптических нейронов будут скорректированы в сторону увеличения или уменьшения. Константы А+ и А- подбираются индивидуально для каждой задачи: если обучающая выборка большая, можно выбрать их маленькими значениями и наоборот для малой обучающей выборки (как у нас сейчас). A+ и A- обычно выражаются друг через друга:

res = LIF_SNN(3, 60, train_stack, list_weight, 0.25)
t_post = res[1]
A_p = 0.8
A_m = A_p * 1.1

for n in range(3):
    for u in range(20):
        
        t1 = np.round(train_stack[u + 10 * n] * 1000)
        t2 = t1.copy()
        
        t2[((t1 <= t_post[n, u]) & (t1 > 0))] = A_p * np.exp((t1[((t1 <= t_post[n, u]) & (t1 > 0))] - t_post[n, u]) / 1000)
        t2[((t1 > t_post[n, u]) & (t1 > 0))] = - A_m * np.exp((t_post[n, u] - t1[((t1 > t_post[n, u]) & (t1 > 0))]) / 1000)
        
        list_weight[n, :] += t2
        
list_weight[list_weight < 0] = 0
list_weight

Мы скорректировали веса, давайте посмотрим, как теперь изменилась точность модели на втором наборе обучающих экземпляров:

res = LIF_SNN(3, 60, train_stack, list_weight, 0.25)
spike_time = res[2]
spike_plot(spike_time, False, False, False)
accuracy_snn(spike_time, 20, 40, df, ind_type, 0)[2]

Отлично! Точность увеличилась до прежнего уровня 93.33%. Теперь проверим точность классификации на всем обучающем наборе (все первые 40 экземпляров каждого типа), используя эти веса:

Напишем код:

train_stack = model_data(0, ind_type, lat_ne, 0, 40)
res = LIF_SNN(3, 120, train_stack, list_weight, 100)
v = res[0]

v_plot(v)
res = LIF_SNN(3, 120, train_stack, list_weight, 0.25)
spike_time = res[2]
spike_plot(spike_time, False, False, False)
accuracy_snn(spike_time, 0, 40, df, ind_type, 0)[2]

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

Теперь запустим сеть на тестовом наборе (последние 10 экземпляров каждого типа):

Напишем код:

train_stack = model_data(0, ind_type, lat_ne, 40, 50)
res = LIF_SNN(3, 30, train_stack, list_weight, 100)
v = res[0]
res = LIF_SNN(3, 30, train_stack, list_weight, 0.25)
spike_time = res[2]

v_plot(v)
spike_plot(spike_time, False, False, False)
accuracy_snn(spike_time, 40, 50, df, ind_type, 0)[2]

Точность 100%!

Кажется, есть одно ложное срабатывание второго постсинаптического нейрона в периоде 27, давайте рассмотрим его более внимательно:

spike_plot(spike_time, True, 27, 3)

Картина вполне ожидаемая, в этом случае из-за того, что третий постсинаптический нейрон сработал раньше "ложного" второго (8.50 мс против 8.66 мс), считается, что именно он сработал в текущем периоде и ошибка не происходит.

Поздравляю! Мы написали простую импульсную нейронную сеть с нуля и научились кодировать данные, используя рецептивные поля с помощью только NumPy и Pandas!)

Источники мудрости:

References:

[1] Alexander Sboev, Danila Vlasov, Roman Rybka, Alexey Serenko, "Spiking neural network reinforcement learning method based on Abstract temporal coding and STDP", Procedia Computer Science Volume 145, 2018, Pages 458-463

[2] Stefan Schliebs, Nikola Kasabov, "Evolving spiking neural networks: A Survey", Article in Evolving Systems, June 2013 DOI: 10.1007/s12530-013-9074-9

[3] Sander M. Bohte, Joost N. Kok, Han La Poutre, "Error-backpropagation in temporally encoded networks of spiking neurons", Neurocomputing 48 (2002) 17–37

[4] S. M. Bohte, H. La Poutre and J. N. Kok, "Unsupervised clustering with spiking neurons by sparse temporal coding and multilayer RBF networks" in IEEE Transactions on Neural Networks, vol. 13, no. 2, pp. 426-435, March 2002, doi: 10.1109/72.991428

[5] M. Kiselev, "Spiking neural networks - Information Representation, Learning, Memory" (manuscript)

[6] Eugene M. Izhikevich, "Simple Model of Spiking Neurons", IEEE TRANSACTIONS ON NEURAL NETWORKS, VOL. 14, NO. 6, NOVEMBER 2003

[7] Dmitry Ivanov, Aleksandr Chezhegov, Mikhail Kiselev, Andrey Grunin, Denis Larionov, "Neuromorphic artificial intelligence systems", Front. Neurosci., 14 September 2022, Sec. Neuromorphic Engineering, Volume 16 - 2022

[8] Laurent U Perrinet, Arnaud Delorme, Manuel Samuelides, Simon Jonathan Thorpe, "Networks of Integrate-and-Fire Neuron using Rank Order Coding A: How to Implement Spike Time Dependent Hebbian Plasticity", June 2001, Neurocomputing 38-40:817-822, DOI:10.1016/S0925-2312(01)00460-X

[9] Senglan Li, Qiang Yu, "New Efficient Multi-Spike Learning for Fast Processing and Robust Learning", April 2020Proceedings of the AAAI Conference on Artificial Intelligence 34(04):4650-4657, DOI:10.1609/aaai.v34i04.5896

[10] Gütig R, Sompolinsky H., "The tempotron: a neuron that learns spike timing-based decisions", Nat Neurosci. 2006 Mar;9(3):420-8. doi: 10.1038/nn1643. Epub 2006 Feb 12. PMID: 16474393.

[11] Baktash Babadi, L. F. Abbott, "Stability and Competition in Multi-spike Models of Spike-Timing Dependent Plasticity", Published: March 3, 2016, https://doi.org/10.1371/journal.pcbi.1004750

[12] Neuromatch Academy: Computational Neuroscience

[13] Zexiang Yi, Jing Lian, Qidong Liu, Hegui Zhu, Dong Liang, Jizhao Liu, "Learning rules in spiking neural networks: A survey", Neurocomputing, February 2023 DOI: 10.1016/j.neucom.2023.02.026

Надеюсь, вам понравилось читать мою статью. Если вам интересен подобный контент, вы можете подписаться на меня здесь или Linkedin