DSP на OpenCL -- строительные блоки и блочные диаграммы

Но по сути, мы тут не про язык разговариваем. У меня копится материал для изложения, который я выкладываю для обсуждения. Первая часть изложения - про объекты в протоколе и в описании модели. А вторая часть продолжение темы выкручивания алгоритмов DSP для примения в составе блочной диаграммы.

Блоки в языке OpenCL C - это такие чудные вложенные функции со шляпкой. В GCC поддерживается расширение языка Си, т.н. "nested functions", похожие сущности -- гнездящиеся функции. В языке OpenCL C 2.0+ можно на стороне устройства описывать из блоков асинхронный процесс обработки.

Достаточно крупные операции в DSP надо делать inline функциями, вернее давать возможность компилятору оптимизировать(векторизовать) большие блоки кода. Один из основных методов дающих существенный прирост производительности - векторизация циклов. Это значит что одиночные циклы, для которых число итераций определяется снаружи следует писать в виде inline функций. Мне нравится набор функций, который предлагает ARM CMSIS DSP. Но есть оговорки и вопросы оптимизации под аппаратуру. У нас есть задача переносить коэффициенты и модели из настольного компьютера в контроллер. В компьютере мы используем векторизацию float, в устройстве управления и сбора данных, сервере автоматизации используем ARM Neon, планируем использовать OpenCL. Нам нужен объект и формат данных для переноса "тензоров" и организации "потока тензоров" средствами протокола обмена, между платформами, нам нужен единый API (оптимизированный под конкретную архитектуру) для представления блочных диаграмм. Концептуально, мы говорим про внешний интерфейс для представления функционального блока FunctionalBlockObject, объекты-контейнеры данных TensorObject. А также мы будем рассматривать некоторый набор имен программных блоков, представленных в библиотеке DSP, чтобы ссылаться на их использование в блочной диаграмме.
Говоря про формат данных для внешнего представления блочных диаграмм, мы на сегодня будем ориентироваться на форматы принятые в BACnet: правила кодирования ASN.1, формат Json и CSML. Эти форматы мы выбираем прежде всего, чтобы не плодить сущности, поддержка этих форматов уже присутствует в составе нашего ПО.

Тензоры

Тензор - это массив данных определенного типа: float, double, int, char. Мы говорим про внешнее представление данных. Тнезор обладает свойством выравнивания данных - это его внутреннее свойство, определяется автоматически при создании объекта. Тензор имеет порядок укладки данных при записи и чтении, формат записи. Порядок записи включает: число каналов (feature), пакет данных (batch), пространственную размерность (spatial). Например, изображение с камеры имеет размерность 2d (spatial), организовано как поток кадров (batch), может иметь несколько каналов (feature). Звук имеет размерность =1 (spatial), организовано в поток по времени (batch), и разделено на каналы (features). Тензор - это единица обрабатываемой информации, с точки зрения функционального блока DSP или NN.
Мы будем исходить из того, что передать тензор целиком в одном пакете может оказаться не возможно. Поэтому обращение к тензору может быть как к массиву, ReadRange - запрос на чтение фрагмента массива. Мы предполагаем некоторую процедуру заполнения входных данных, объектов связанных между собой блочной диаграммой, затем запуск преобразования. Результатом преобразования может являться набор тензоров. Результат преобразования может распространяться по подписке CoV. Преобразование запускается по готовности входных данных, т.е. когда все входные данные(тензоры) записаны.
Тензоры могут быть вида TensorOutputObject, TensorInputObject. TensorOutputObject - обладает командным свойством на PresentValue. Когда записывается содержимое тензора, происходит запуск зависимостей по готовности.

Функциональные блоки

Функциональный блок содержит ссылку на исполняемый код(kernel) - DeviceObjectReference, список аргументов, список фиксированный идентификаторов объектов ARRAY_OF(DeviceObjectReference) - зависимости. В процессе исполнения, блок оформляет подписку на соответствующие объекты списка зависимостей. Функциональный блок имеет статус(состояние,контекст) - это некоторая структура данных, в которой могут хранится промежуточные величины. Контекст сохраняется и может быть перезаписан в режиме out-of-service.
Для блока мы определяем фазы: Инициализация, обновление контекста, обработка, пауза/завершение/сериализация контекста. Блочная диаграмма это не сам конвейер, это пакетное задание, которое может выполняться на конвейере с использованием данных. Конвейер это очередь исполнения по готовности.
Функциональный блок может обладать отладочной информацией (опция), для отладки используются штампы времени: queued, submit, start, end.

Программы

Программа - это текст или бинарник. Со стороны хост-приложения возможно определить из бинарника список блоков и список агрументов(входных/выходных портов). Каждый порт описывается 'форматом записи', однако формат записи из программы не определить.

Строительные блоки

Список функций API

dot(a,b)

операция свертки вектора sum_{0}^{N-1}(a[i]·b[i])

conv(a,b,n)

операция частичной свертки, выполняется с разворотом вектора sum_{0}^{N-1}(a[i]·b[n-i-1]). В это определение попадает корреляционная функция

float cascade_biquad_df2T(z,x)

вычисление значения каскада y = z+b0·x

float cascade_biquad_df2T_update(S,x)

обновление регистров каскада (S)

float fir(S,x)

FIR фильтр. Интерполяция.

pid(S,x)

PID регулятор.

lms_update(S,x)

LMS подстройка параметров адаптивного фильтра

IIR фильтры обрабатывают batch составляющую, время.
Для задач мультирейт используются фильтры LPF_decimate (по времени).
Аудио процессор может обрабатывать полосы и сшивать их на выходе. Для этого используются фильтры Allpass. В каких-то случаях могут использоваться фильтры iir_lattice.
Пид регулятор - это фильтр типа biquad. LMS фильтр - подстройка.

PID регулятор

 

Передаточная функция:
H(s) = ( K_P + K_I·1/s + K_D·s ) Замена: 1/s = T/2·(1+z^-1)/(1-z^-1), s = 1/T·(1-z^-1)

Передаточная функция в дискретном пространстве:
H(z) = KP + KI·T/2·(1+z-1)/(1-z-1) + KD/T·(1-z-1)
или 
H(z) = KP + KI·T·1/(1-z-1) + KD/T·(1-z-1)
Ki := KI·T
Kd := KD/T
Может быть представлена биквадратной секцией

void pid_init(S, x){
  S->b0 =  Kp + Ki + Kd;

  S->b1 = -Kp -2Kd
  S->b2 =  Kd;

  S->z1 = 0;
  S->z2 = 0;
}
// первый проход
float pid(S, x){
  return b0*x+z1;
}
// второй проход, обновление регистров
void pid_update(S, x){
  z1 =(b1-b0)*x + (z2 - z1);
  z2 = b2*x;
}

При разработке PID регулятора мы также используем два прохода: первый проход встроен в поток обработки, второй - поток обновления регистров и состояния фильтров. Пока что у нас получается что абсолютно все фильтры представлены выражением: y = g*x+Z, включая и PID регулятор. Далее мы хотим добавить свойство автоматической подстройки каэффициентов FIR фильтров и IIR фильтров. На втром цикле обновления регистров, мы добавляем возможность подстройки коэффциентов по методу LMS (). Есть известное преобразование FIR фильтра в LMS - получается адаптивный фильтр. Мы создаем сеть, в которой каждый функциональный блок должен обладать адаптивностью. Точнее можно сказать, что в процессе работы возникают периоды обучения - подстройки коэффициентов модели. Таким образом, фаза UPDATE может сопровождаться "обучением", если известна величина ошибки выходного сигнала. Для каждого слоя сети, мы можем определить этот признак "обучения".

LMS фильтры

Фильтры адаптивные (LMS) мы будем изготавливать из фильтров типа IIR и FIR. Метод LMS можно применять для обучения (подстройки весовых коэффициентов) нейронной сети. Мы вводим единую методику обучения/распознавания для нейронных сетей (NN) и для структур цифровой обработки сигналов (DSP). Нам известен способ изготовления адаптивного цифрового фильтра из FIR фильтра эту методику мы хотим расширить на остальные задачи.

Для некоторого преобразования c передаточной функцией представленной полиномом H(z) = h₀+h₁z^-1+h₂z^-2+...+h_Nz^-N , можно подстраивать коэффициенты, если известна ошибка выходного значения. Ошибка может быть пересчитана в поправку каждого конкретного коэффициента.
Ошибку на выходе определим, как квадрат отклонения.
E = (d[n] - y[n])² = (d[n] - sum(b[i]z^-i))²
изменение коэффицента для фильтра можно определить как градиент. Используем метод "градиентного спуска"
w[n+1]=w[n]-μ·∂E/∂w
Здесь μ - это величина шага подстройки. Для стабильного проявления методики, бывает полезно ограничить ошибку в каких то пределах, используя функцию clamp(). Берем частную производную, не сильно увлекаясь частностями. w_i[n+1] = w_i[n] - 2μ·e[n]·∂e[n]/∂w_i
∂E/∂w_i = -2(e[n])(x[n-i])

b_i[n+1]=b_i[n]+2μ (e[n])(x[n-i]).
Это все! Каждый коэффицент подстраиваем исходя из ошибки на выходе. IIR фильтр - это комбинация двух фильтров FIR. Для подстройки коэффициентов A(z) в IIR-фильтре с передаточной функцией H(z)=B(z)/A(z) используем аналогичный вывод:
a_i[n+1]=a_i[n]+2μ (e[n])(y[n-i]).

{Тут в изложении присутсвует некоторое искусство, с которым надо поразбираться. Работать будет в любом случае, однако математически правильный вариант брать частную производную не от уравнения в конечных разностях. Варианты, которые можно получить исследуя конечные разности и интегралы от дельта-функций - это линии задержки разного вида и фильтры на поток ошибок. К сожалению, быстро ответ не нашел. Нашел, что этим когда-то занимался I.-D. Landau, }
{1976 Feintuch’s algorithm}
{A parallel-type adaptive algorithm is proposed which utilizes the fast least-squares method. Its computational complexity is much less than that of I. Landau's (1978) method, which is based on hyperstability theory Hyperstability requires much computation because it involves matrix operations. The proposed method has nothing to do with hyperstability. It is also shown theoretically and using computer simulation, that the convergence performance of the proposed method is better than that of the series-parallel-type fast least-squares method.} Линия задержки будет разная для разных форм представления фильтров. Глядя на презентацию, мне почему-то кажется, что тут может быть ошибка в интерпретации, что-то не то с частной производной. Должна быть простая форма. Мне видится простая форма при вторичном использовании регистров z. Предположительно выглядит так, доказать надо на практике:

void lms_dt2t_update(S, x, ue){
  b0+= ue*x;
  b += ue*x;// вектор
  a += ue*y;// вектор
}
void lms_update(S, x, ue){
  b0+= ue*x;
  b += ue*z;// вектор
  a += ue*z;// вектор
}

Если мы применяем LMS к фильтру типа PID должны быть веса разные для коэффициентов. Коэффциенты A не подстравиваются.

DSP на OpenCL - Изоморфные преобразования алгоритмов

Хочу поделиться знаниями, вернее акцентировать внимание на известном факте.. Каково место цифровой обработки сигналов в коде. Хотим затолкать цифровую обработку в поток OpenCL, выполнять обработку на потоке данных. При этом не дает покоя одна мысль. Мы хотим получить наивысшую скорость работы и малые задержки. Поэтому мы каждый алгоритм DSP хотим разделить на две части. В первой части рассчитываем результат за минимальное число операций, во второй обновляем значения переменных. Я начну писать, потом видно будет.
В разделе DSP есть два вида фильтров: IIR- и FIR- фильтры. Прежде всего, оговорюсь все это относится к анализу LTI-систем, для которых можно описать передаточную функцию Y(s)=H(s)X(s) в непрерывном виде и далее перейти к дискретному виду передаточной функции H(z) путем замены переменной s←2/T(z-1)/(z+1). Передаточную функцию можно представить в виде отношения полиномов B(z)/A(z). Т.е.в виде дроби (дробно-рациональной функции) - сверху и снизу дроби полином z. Причем буква Z^-1 на самом деле представляет из себя оператор задержки в дискретном по времени процессе. Как вообще люди в процессе своего образования доходят до этого - загадка. Думаю все сознательные электронщики и программисты должны это знать, но полагаю мало, кто знает. Наше образование не было полным в том виде, в каком оно может служить инструментом программиста или электронщика. Статья для избранных. На самом деле я уже начинаю сомневаться, что вообще этому учат. Когда-то учили анализу систем, учили передаточным функциям. Предлагаю, сейчас найти учебник на английском (по теме Linear time-invariant system и Digital Signal Processing). Для работы, чтобы получить готовый набор коэффициентов использую GNU Octave или MATLAB с модулем Signal, он умеет готовить матрицы коэффициентов. Использую конспекты лекций по курсу DSP & Digital Filters, потому что мне это близко и понятно, именно в таком виде. Автор курса ссылается на учебник.
[#] "Digital Signal Processing" Sanjit K. Mitra, 4th edition, McGraw Hill, 2011. ISBN 0071289461
Для стационарных систем (LTI) передаточная функция — это дробно-рациональная функция комплексной переменной s. Тернистый путь познания лежит через преобразования Фурье и Лапласа. Однако более короткий путь можно найти через дифференциальные уравнения в конечных разностях, только этот путь неведом выпускникам со специальностью кибернетика, так учат физиков. У меня в образовании был курс ТФКП(теория функции комплексного переменного) и курс математической физики (который символизирует наезд дельта функцией на уравнения), от которого ничего тут использовать не буду.
1) Систему можно представить в виде дробно рациональной функции N-ого порядка в странном пространстве наполненном операторами задержки и умножителями, -- это понятно программисту. Причем порядок числителя не выше чем у знаменателя. 2) Безумную дробь можно разложить на полюса и нули(полюса - это решения полинома в знаменателе, нули - решение полинома в числителе), путем факторизации полинома(еще одно умное слово - разложить на скобки, - так говорят в школе). Длинную дробно рациональную функцию большого порядка можно представить в виде суммы или произведения этих самых дробно-рациональных функций второго порядка. (точка!! - никакой больше теории, только картинки).
Иногда я чувствую, что являюсь хранителем некоего сокровенного невербального знания :). Могу доказать теорему не произнося слова, а рисуя только квадратики, треугольнички и стрелочки. Передаточные функции в дискретном пространстве можно строить в виде графов и представлять математические преобразования над функциям, как тождественные преобразования графов. Граф в точности отражает архитектуру программной или аппаратной реализации алгоритма цифровой обработки. На картинке буквой Z обозначаем оператор задержки, фактически это будет переменная (буфер, в аппаратуре - D-триггер), в которой хранится значение с предыдущего шага обработки.

Рис.1 Биквадратный фильтр. Прямое представление 1-типа DF1 (слева) и прямое представление 2-типа DF2 (справа)

Рисунок #1 получен путем прямой записи выражения B(z)X = (b₀+b₁·z^-1+b₂·z^-2)X и затем выражения A(z)Y = (1-a₁·z^-1-a₂·z^-2)Y, где z^-1 - оператор задержки на такт.
В первом случае (DF1)сложность, 5 умножений и 4 сложения: 5m+4a. Во втором (DF2) сложность та же, но для получения хранения промежуточных значений достаточно всего двух регитров Z.

Рис.2 Биквадратный фильтр. Прямое представление 2-типа DF2 (слева) и прямое представление второго типа транспонированное DF2T (справа)

Метод Транспонирования цифрового фильтра:
1) Меняем направление обхода, все стрелочки и треугольники разворачиваем. 
2) Слева будет выход, справа вход. Меняем местами буквы X и Y.
3) Разворачиваем буферы Z.
4) Где были узлы ставим "плюсик", где плюсики - ставим узел.

Цель можно считать достигнутой, мы развернули фильтр таким образом (DF2T), что на первом проходе можно получить в одно действие результат: Y = b₀·X+z₀ От входа X до выхода Y одно суммирование и одно умножение. Сложность первого прохода: 1m+1a. Этот проход совершается в той же самой математике, что и процесс свертки слоя нейронной сети.
К чему все это?. А, к нейролизации цифровой обработки сигналов в дискурсе для непосвещенных. Все что ниже первой линии будем считать частью процесса обучения, т.е. расчитывать и обновлять значения Z будем на втором проходе сети.
Какие у нас открываются возможности при использовании параллельной архитектуры. Правильно, фильтры могли бы считать паралельно, а не последовательно. Привычное разложение передаточной функции высокого порядка - в каскад (произведение) фильтров второго порядка. Каскад - это последовательная обработка фильтров, предполагает загрузку одного ядра/потока. Возможно представление в виде суммы, можно разложить в сумму фильтров первого и второго порядка. Cмотрю на структуру каскада фильтров. Каскад расчитывается, как sum(z_1k+b_0kX). Переобозначим переменную, S=sum(z_1k) -- сумму можно расчитывать на втором проходе.

Рис.3 Каскад биквадратных фильтров.

Получаем финально такой результат: каскад фильтров можно представить на первом проходе Y =Z+ gX. Весь фильтр в одно действие!

// Алгоритм 
float cascaded_iir_df2t(float x){
// первый проход
   y = g*x+s;
   return y;
}
float cascaded_iir_df2t_update(float x){
// второй проход, вычисляем s
   s=0;
   x=g*x;
   for (i=0; i<len; i++){
     y = x+z1[i];
     z1[i] = b1[i]*x - a1[i]*y + z2[i];
     z2[i] = b2[i]*x - a2[i]*y;
     s+=z1[i];
     x = y;
   }
   return s;
}

Преобразования фильтров

Рис.4 Тождественные преобразования. a) замена переменных, b) объединение идентичных буферов, c) свойство линейности преобразования.

Предлагаю для практики проделать преобразования фильтров. На рисунке #5 привожу последовательное, шаг за шагом, преобразование фильтра Allpass первого порядка из одной формы в другую. Считается, что каждый студент прослушавший курс лекций должен понимать, как проделать подобный переход.

Рис.5 Преобразования фильтра Allpass

На рисунке #6 представлена одиночная операция транспонирования над одним каскадом фильтра FIR.

Рис.6 Преобразования фильтра FIR. Транспонирование каскада.

Векторизация фильтров

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

Рис.7 Прямая и транспонированная форма FIR-фильтра

Ниже пример реализации FIR фильтра. Также разделен на два прохода.

// FIR фильтр с явной векторизацией
float  fir_dft(float x){
  y = g*x+z.s0;
  return y;
}
float4 fir_dft_update(float4 z, float x){
  Z0=(0.0f);
  z = shuffle2(z,Z0, (int4){1,2,3,4});
  z =z + H*x;
  return z;
}

В данном варианте расчета, для обновления регистров фильтра нужно всего две процессорных инструкции shuffle и fma.

Функция активации в слое нейронной сети

В статье дано описание функций активации. Выбран класс гладких функций удобный для реализации на CPU с поддержкой инструкций AVX512 и на GPU Intel Gen8-Gen11. При работе с GPU мы ориентируемся на язык OpenCL C. Дается некоторое количество утверждений, основанных на предположениях автора. Мы бездоказательно утверждаем, что использование класса гладких Smoothie-функций дает достаточную точность, для получения адекватных результатов работы сети.

Применяют два варианта функции активации в промежуточных слоях, остальные сводятся к производным. Мы хотим обобщить понятие о функции активации. В обобщенном виде к функции предъявляется ряд требований, таких как гладкость, монотонность и уровни насыщения. Функция задана на области определения и насыщается при определенных значениях. Важен внешний вид функции и важно, чтобы производная функции работала, давала соответствующий градиент в правильном направлении. Величина градиента не сильно важна и сказывается только на времени сходимости при обратном распространении ошибки.
Вариант 1) Сигмоида. Функция типа ступенька с размытыми краями. В обобщенном виде сигмоида может быть заменена на полное отсутствие сигмоиды лишь бы насыщалась на уровне 0 и 1. Т.е. вполне реально активацию можно заменить на операцию умножения с насыщением. Производная сигмоиды: s'(x) = s(x)(1-s(x)) выражается через значение самой функции, т.е. не надо знать ничего, кроме значения на выходе.
Вариант 2) Тангенс гиперболический. Это такая функция, которая монотонно возрастает имеет один перегиб. Имеет область значений [-1,+1].Все что реально нужно знать про эту функцию - она монотонная, нечетная и производная этой функции: (1+t(x))(1-t(x)). Тангенс гиперболический можно выразить через сигмоиду и наоборот. По сути это одна и та же функция.
Надо ли говорить, почему функция сигмоида и тангенс гиперболический не подходит под описание?! Потому что мы говорим про насыщение, а насыщение у сигмоиды не наблюдается. Это значит что нейрон будет "обучаться"(подстраивать коэффициенты) только в области, где производная не равна нулю.
Давайте предположим, что традиция использовать экспоненту в функции обусловлена исключительно требованием гладкости и дифференцируемости функции, которая в свою очередь продиктована методом вычисления градиента. Давайте предположим, что нам удаться сбалансировать систему таким образом, чтобы она не уходила в насыщение в процессе обучения.

tanh(x) = (1-e-2x)/(1+e-2x)
sigm(x) = 1/(1+e-x)
relu(x) = ln(1+ex) -- SmoothReLU, softplus, выпрямитель с гладким переходом.

Следует отметить, что производная от SmoothReLU дает сигмоиду. Функцию tanh можно выразить через сигмоиду и наоборот:

tanh(x) = (1-e-2x)/(1+e-2x) = 2sigm(2x) - 1
sigm(x) = 0.5*(tanh(x/2)+1)

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

d t(x)/d(x) = (1+t(x))(1-t(x))
d s(x)/d(x) =    s(x)·(1-s(x))

-- Производная имеет смысл градиента.

Кусочно-полиномиальные приближения функции

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

// Кусочно-полиномиальные приближения функций
float activation(float x)
{
const int N=8;
   x = clamp(X_MIN, X_MAX,x);
// округляет в меньшую сторону до целого числа.
   idx = floor(x)-X_MIN;
   t  = x-floor(x);// дробная часть
   v0 = Coeff0[idx];
   v1 = Coeff1[idx];
   v2 = Coeff2[idx];
// квадратная интерполяция 
// v(t) = t*(t*v2+v1)+v0
   v = fma(t, v2,v1);
   v = fma(t, v, v0);
   return v;
}

Примерно в таком виде Intel описывает некоторую оптимизацию в собрании сочинений Intel Architectures Optimization Manual, глава 7.5.2 Fused post GEMM, Sigmoid Approximation with Minimax Polynomials. Которая оптимизация сводится к интерполяции полиномами второй степени. Этот вариант запасной. Он работает, не имеет противопоказаний, это доказано на опыте.
Можно определить функцию активации следующим образом.
На диапазоне [0+eps,1-eps] - происходит линейный рост с производной =1.
(0-eps, 0+eps) и (1-eps, 1+eps) -- сглаженный переход к насыщению - квадратным полиномом.
Назовем функцию, которая осуществляет плавный переход к насыщению Smoothie-выпрямитель. В масштабе 1.0 Smoothie-выпрямитель мы определим на интервале (-1.0,1.0) со значениями производной на конце интервала 0 и 1, соответственно.

Квадратичная интерполяция функции Smoothie-выпрямитель
F (x)  =  a2 x^2 + a1 x + a0
F'(x)  = 2a2 x   + a1
// граничные значения функции
F ( 1) =  a2 + a1 + a0 = 1
F (-1) =  a2 - a1 + a0 = 0
// производная непрерывная
F'( 1) = 2a2 + a1      = 1
F'(-1) =-2a2 + a1      = 0

a0 = 1/4, a1 = 1/2, a2 = 1/4, 

На области определения (-1,1)
SmoothieRect (x) = (x+1)(x+1)·0.25
SmoothieRect'(x) = (x+1)·0.5

В общем виде этот пример можно рассматривать, как метод вычисления коэффициентов полинома на интервале. Мы задаем производную и значение функции в каждой строчке таблицы. Таким образом рассчитываем коэффициенты интерполяции полиномом второй или третьей степени.
Позволю себе небольшое отступление для обоснования расчета коэффициентов табличного метода. Для каждого интервала нам нужно рассчитать коэффициенты полинома. Для этого нам нужно решить систему линейных уравнений полученных подстановкой граничных условий на интервале (0,1). Кто не помнит, надо посмотреть метод Гаусса-Жордана для решения систем линейных уравнений.

// Решение системы линейных уравнений.
Запишем функцию и его производную. Надо найти коэффициенты ai
F (x)  =  a3 x3 + a2 x2 + a1 x + a0
F'(x)  = 3a3 x2 +2a2 x  + a1
Шаг 1. Строим матрицу из множителей коэффициентов полинома ai, 
при заданном x и заданном значении функции bi:
  F(0) = (0, 0, 0, 1 | b0)  
 dF(0) = (0, 0, 1, 0 | b1)  
  F(1) = (1, 1, 1, 1 | b2)  
 dF(1) = (3, 2, 1, 0 | b3)  
Шаг 2. Вычитаем так чтобы получилась диагональная матрица.
b3 = b3-3b2
  F(0) = (0, 0, 0, 1 | b0)  
 dF(0) = (0, 0, 1, 0 | b1)  
  F(1) = (1, 1, 1, 1 | b2)  
 dF(1) = (0,-1,-2,-3 | b3-3b2)  
b2 = -2b2+b3
  F(0) = (0, 0, 0, 1 |  b0)  
 dF(0) = (0, 0, 1, 0 |  b1)  
  F(1) = (1, 0,-1,-2 |-2b2+ b3)  
 dF(1) = (0,-1,-2,-3 |  b3-3b2) 
b3 = b3-3b2+2b1+3b0
b2 = 2b0+b1-2b2+b3
          a3 a2 a1 a0
  F(0) = (0, 0, 0, 1 |  b0 )  
 dF(0) = (0, 0, 1, 0 |  b1 )  
  F(1) = (1, 0, 0, 0 | 2b0+ b1-2b2+b3 )  
 dF(1) = (0,-1, 0, 0 | 3b0+2b1-3b2+b3 ) 

Итого получаем: 
a0 = b0;
a1 = b1;
a3 = 2b0+ b1-2b2+b3
a2 =-3b0-2b1+3b2-b3

Далее мы рассмотрим вариант функции tanh-like активации заданный на интервале (-1,1) полиномом третьей степени.

Гладкие функции активации

Обращаю внимание на существование иных функций, smoothie-функций, которые подпадают под определение функции активации. Например, кубическая интерполяция Эрмита, функция smoothstep, S(x)=3x² - 2x³. Эта функция введена в язык OpenCL и достойна применения. Функция имеет хорошую производную S'(x) = 6(x - x²) = 6x(1-x). Даже производная похожа на сигмойду, sigmoid-like function.
Продемонстрируем вывод полиномиальной функции активации, представляем функцию активации полином третьей степени. Приравниваем производную вне области определения (-1, 1) нулю и приравниваем значение вне области определения -1 и 1 соответственно, чтобы получить tanh-like функцию.
Вычисляем коэффициенты полинома для граничных условий:

Для полинома и его производной
T (x) = a3 x3 + a2 x2 + a1 x + a0
T'(x) =3a3 x2 +2a2 x  + a1
Составляем систему уравнений для граничных условий
T ( 1) =   a3 +  a2 + a1 + a0 = 1
T (-1) = - a3 +  a2 - a1 + a0 =-1 
T'( 1) =  3a3 + 2a2 + a1      = 0
T'(-1) = -3a3 + 2a2 - a1      = 0
Из выражения для производной получаем:
a1 = -3a3, a2 = 0.
Подставляем в выражение для функции:
T ( 1) = -2a3 + a0 = 1
T (-1) =  2a3 + a0 =-1
Получаем набор коэффициентов:
a0 = 0, a1 = 3/2, a2 = 0, a3 = -1/2
Функция Tanh-like на множестве (-1,1)
T (x)= (3x - x3)1/2 = 2*Smoothstep(x + 1/2)-1
T'(x)= (1  - x2)3/2 = 3/2 (1-x)(1+x)

Интерполяция Sigmoid-like на (-1,1).
S (x)= (T(x)+1)/2 = Smoothstep(x + 1/2)
S'(x)=  T'(x)/2

На языке OpenCL C кубическую интерполяцию функции активации Sigmoid и Tanh можно выразить с использованием базовой функции smoothstep, которая может быть реализована аппаратно или хорошо оптимизирована для векторных типов. В спецификации Khronos указано, что функции реализуются с использованием инструкций fma и mad.

T(x) = 2*smoothstep(-1.0, 1.0, x+0.5)-1.0;
S(x) = smoothstep(-1.0, 1.0, x+0.5);

float smooth_tanh(float x){
  x = clamp(x, -1.0f, 1.0f);
  return (3.0f - x*x)*x*0.5f;
}
float smooth_sigm(float x){
  x = clamp(x, -1.0f, 1.0f);
  return (3.0f - x*x)*x*0.25f+0.5f;
}
float smooth_step(float edge0, float edge1, float x){
  float t;
  t = clamp((x - edge0)/(edge1 - edge0), 0, 1);
  return t * t * (3 - 2 * t );
}

Коль скоро затронули тему Smoothie-функций, можем сюда же определить смузи-функцию выпрямителя, SmoothReLU. Производная от SmoothReLU должна давать функцию SmoothStep, а производная от ReLU - функцию Step, пороговую активацию. Строго говоря, SmoothReLU должна представлять собой полином четвертой степени.

Тестирование

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

// расcчет коэффициентов линейной интерполяции
. . .

Заключение

Я хочу поставить току в обсуждении. Самая простая функция активации - это ее отсуствие. Это и есть самая эффективная реалиазция.

float activation_s(float x){
  return clamp(-1.0f, 1.0f, x);
}

float activation_u(float x){
  return clamp(0.0f, 1.0f, x);
}

float activation_relu(float x){
  return fmax(0.0f, x);
}

Все что нам нужно - ограничивать значение на выходе.
Кстати, чтоб развеять мистический туман вокруг встроенной функции clamp, скажу, что ее можно реализовать двумя функциями: fmax и fmin, которые в свою очередь представлены соответствующими векторными инструкциями.

// Эмуляция функции clamp
float clamp(float x, float minval, float maxval) { 
  return fmin(fmax(x, minval), maxval); 
}

Такой вид функции активации позволяет легко перейти к реализации алгоритмов в целых числах. Активация в целых числах - это просто насыщение.
Пример такой операции оптимизированной под AVX-512:

activation_u8(a,b,c,d){
// добавление константы относится к свертке.
// a = _mm512_adds_epi32(a,a0);
// b = _mm512_adds_epi32(b,b0);
// c = _mm512_adds_epi32(c,c0);
// d = _mm512_adds_epi32(d,d0);
 v0= _mm512_packus_epi32(a, b);
 v1= _mm512_packus_epi32(c, d);
 return _mm512_packus_epi16(v0, v1);
}

activation_i16(a,b){
// a = _mm512_adds_epi32(a,a0);
// b = _mm512_adds_epi32(b,b0);
 return _mm512_packs_epi32(a, b);
}