Автоматическое дифференцирование

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

f(a : double) : double
{
 a*a + 3 * a
}

df(a : double) : double
{
 2*a + 3
}


Очевидно, что, как минимум, код f и df функционально зависим, а, возможно, мы наблюдаем дублирование функциональности. Ситуацию можно исправить, если вспомнить определение производной и, соответственно, исправить код df:

df(a : double) : double
{
 def delta = 0.1;
 (f(a+delta)-f(a))/delta
}

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

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

Мнимая единица i полагается решением уравнения x²+1 = 0. По-аналогии введем другую мнимую единицу d, которая будет решать уравнение x²=0 и при этом не равняться нулю. Затем вспомним разложение функции в ряд Тейлора:

Если в качестве h подставить нашу мнимую единицу и воспользоваться тем, что d²=0, то получим:

Таким образом, что бы узнать значение производной функции f в точке x достаточно взять коэффициент при мнимой части числа f(x+d). В современных языках не проблема определить класс, для работы с комплексными числами, так же легко можно определить класс для работы с этими расширенными числами. Ниже переведен пример на Nemerle:

 public class Dual
 {
  public this(real : double, imaginary : double)
  {
   this.real = real;
   this.imaginary = imaginary;
  }
  
  real : double;
  imaginary : double;
  
  public Real : double { get { real } }
  public Imaginary : double { get { imaginary } }
  
  public static @+(a : Dual, b : Dual) : Dual
  {
   Dual(a.real + b.real, a.imaginary + b.imaginary)
  }
  
  public static @-(a : Dual, b : Dual) : Dual
  {
   Dual(a.real - b.real, a.imaginary - b.imaginary)
  }
  
  public static @*(a : Dual, b : Dual) : Dual
  {
   Dual(a.real * b.real, a.real * b.imaginary + a.imaginary * b.real)
  }
  
  public static @/(a : Dual, b : Dual) : Dual
  {
   a*Dual(1/b.real, - b.imaginary / (b.real*b.real))
  }
  
  public static @:(k : double) : Dual
  {
   Dual(k,0)
  }
 }


Теперь можно переписать пример из начала статьи:

f(a : Dual) : Dual
{
 a*a + 3.0 * a
}
  
f(a : double) : double
{
 f(Dual(a,0)).Real
}
  
df(a : double) : double
{
 f(Dual(a,1)).Imaginary
}

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

public Sin(a : Dual) : Dual
{
 Dual(Math.Sin(a.Real), Math.Cos(a.Real)*a.Imaginary)
}

Хочу еще:

	Automatic differentiation
	Dual numbers
	Automatic Differentiation
	Automatic Differentiation, C++ Templates and Photogrammetry
	Nemerle

Хотел создать блог Computer Science и поместить этот пост туда, но на хабре, надеюсь, что временно, нельзя создавать новые блоги, поэтому помещаю в блог «Ненормальное программирование». Блог Computer Science планировался быть местом, куда можно писать пересказы, переводы, аннотации к статям попадающим под определение Computer Science, ссылками на pdf и так далее. Всем, кто хотел бы писать статьи в этот блог советую не дожидаться его открытия, а писать в «Ненормальное программирование» с тегом «to computer science».