Автоматическое дифференцирование «на пальцах»

в 5:19, , рубрики: Блог компании Intel, математика, ненормальное программирование, производная, численные методы, метки: , ,

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

Простые соображения

Итак, пусть у нас есть какая-нибудь функция f(x). Пусть мы в некоторой точке x_0 знаем ее значение f(x_0)=f_0, и знаем значение ее производной f'(x_0)=f_0'. Мы знаем только только два числа: f_0 и f_0', больше ничего знать нам не надо — ни выражения для f(x), ни даже значения x_0. Рассмотрим функцию g(x)=2f(x) и зададимся вопросом: чему равно ее значение и значение ее производной в точке x_0? Очевидно: g(x_0)=2f_0 и g'(x_0)=2f_0'. Обратите внимание, что в правой части здесь стоят только те числа, которые мы знаем.

Вопрос чуть сложнее: рассмотрим функцию h(x)=f^2(x), и зададимся про нее тем же вопросом. Несложно видеть, что и в этом случае мы легко находим h(x_0)=f_0^2, и h'(x_0)=2f_0f_0'. Аналогично, для l(x)=cos(f(x)) имеем l(x_0)=cos(f_0) и l'(x_0)=-sin(f_0)f_0'. И так далее: для любой элементарной функции, примененой к f(x), мы можем вычислить значение и производную, используя только два числа: f_0 и f_0'.

Далее, пусть есть еще какая-нибудь функция m(x), и про нее мы тоже знаем только два числа: m(x_0)=m_0 и m'(x_0)=m_0'. Тогда для функции p(x)=f(x)+m(x) мы столь же легко можем вычислить и значение, и производную в той же точке; и аналогично для любой бинарной операции. Например, для умножения имеем: если q(x)=f(x)cdot m(x), то q(x_0)=f_0m_0 и q'(x_0)=f_0m_0'+f_0'm_0.

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

Код

Легко пишется класс, который реализует такую арифметику:

class Derivable {
        double val, deriv;

    public:
        Derivable(double _val, double _deriv) : val(_val), deriv(_deriv) {}

        double getValue() {return val;}
        double getDerivative() {return deriv;}

        Derivable operator+(Derivable f) { 
           return Derivable(val + f.val, deriv + f.deriv);
        }
        Derivable operator-(Derivable f) { 
           return Derivable(val - f.val, deriv - f.deriv);
        }
        Derivable operator*(Derivable f) { 
           return Derivable(val * f.val, deriv * f.val + val * f.deriv);
        }
        Derivable operator/(Derivable f) { 
           return Derivable(val / f.val, (deriv * f.val - val * f.deriv) / f.val / f.val);
        }

        friend Derivable cos(Derivable f);
};

Derivable cos(Derivable f) {
    return Derivable(cos(f.val), -sin(f.val)*f.deriv);
}

Теперь, если у нас есть код, вычисляющий некоторую функцию, то достаточно просто заменить везде double на Derivable — и получится код, вычисляющий ту же функцию вместе с ее производной. Правда, конечно, возникнет вопрос: с чего начать? Мы пока умеем по Derivable получать новые Derivable, но откуда взять самые первые Derivable? На самом деле, все понятно. В выражения для нашей функции входят, во-первых, различные константы, не зависящие от x, и, во-вторых, собственно сама x. Любую константу c, не зависящую от x, надо, конечно, заменить на Derivable(c, 0); а вхождения собственно переменной x — на Derivable(x0, 1). (Здесь для понятности x0 обозначает то значение x, для которого мы вычисляем функцию. В реальной программе, скорее всего, соответствующая переменная будет называться тоже x).

Вот пример вычисления функции ax^3-cos(x/2) вместе с ее производной:

Derivable f(double x, double a) {
    Derivable xd(x, 1);
    Derivable ad(a, 0);
    Derivable two(2, 0);
    return ad*xd*xd*xd - cos(xd/two);
}

Естественно, чтобы не городить такой код, проще добавить в наш класс неявное преобразование типов:

	Derivable(double c): val(c), deriv(0) {}

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

	static Derivable IndependendVariable(double x) { return Derivable(x,1); }

после чего код функции f становится еще проще:

Derivable f(double x, double a) {
    Derivable xd = Derivable::IndependendVariable(x);
    return xd*xd*xd*a - cos(xd/2);
}

Вот пример кода, решающего уравнение ax^3=cos(x/2) методом Ньютона (здесь я еще, естественно, операторы сделал глобальными, чтобы работало a*xd; выше они члены класса только для того, чтобы не загромождать код). Если вы захотите изменить решаемое уравнение, вам надо поменять код функции f и всё; производная будет автоматически вычисляться верно. (Конечно, в данном примере, может быть, проще было бы посчитать производную руками, потому что функция простая, — но как только уравнения у вас становятся более сложными, возможность не думать о коде вычисления производной оказывается очень удобной.)

Все это будет работать, даже если у нас в коде есть не очень сложные ветвления, циклы, вызовы других функций и т.д. — главное везде все промежуточные переменные (и в том числе результаты промежуточных функций) сделать Derivable, остальной код даже менять не потребуется!

Конечно, можно придумать и такую функцию, для которой такой подход работать не будет — вдруг вам придет в голову вычислять функцию f методом Монте-Карло? — но все равно область применения автоматического дифференцирования весьма широка. (А потом, если уж вы и правда вычисляете свою функцию методом Монте-Карло, то как вы вообще представляете себе ее производную?)

Обобщения

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

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

В-третьих, есть еще один интересный подход, который я краем глаза видел в одном из open-source кодов для автоматического дифференцирования (см. по ссылкам из википедии). Можно сделать класс Derivable шаблонным, принимающим в качестве параметра шаблона тип данных для значений функции и производной (т.е. чтобы надо было писать Derivable<double> и т.п.; запись Derivable<T> будет соответствовать функции mathbb{R}to T). Тогда если написать Derivable<Derivable<double> >, то не получатся ли вторые производные автоматически? Эдакое применение автоматического дифференцирования к автоматическому дифференцированию. Правда, при таким подходе делается лишняя работа: если расписать, например, какой получится оператор умножения, то видно, что всё получится правильно, но первая производная будет вычисляться дважды. Кстати, при правильной инициализации начальных переменных объекты типа Derivable<Derivable<double> >, видимо, можно применять и для вычисления производных по нескольким независимым переменным.

Другие реализации

Отмечу, что описанный выше подход с перегрузкой операторов не является единственным возможным; даже различают «forward» и «reverse» реализации (подробнее см. в википедии и по ссылкам оттуда). В частности, по-видимому, многие коды, на которые приведены ссылки в википедии, используют несколько более общие подходы — но я их смотрел только очень поверхностно.

Ссылки

http://habrahabr.ru/post/63055/
http://en.wikipedia.org/wiki/Automatic_differentiation, со ссылками на множество реализаций
http://www.coin-or.org/CppAD/Doc/mul_level.htm (про конструкцию вида Derivable<Derivable<double> >)

Автор: pkalinin

Источник

Поделиться

* - обязательные к заполнению поля