8.2. Дуальные числа

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

Определение 8.3 (Множество дуальных чисел)

Дуальные числа \(\dualset\) расширяют[1] действительные числа \(\real\) числом \(\dual\), для которого определяются следующие правила

\[\dual{} \neq 0, \quad \dual{} \times 0 = 0, \quad \text{но} \quad \dual{} \times \dual{} = \dual{}^2 = 0.\]

Например, в дуальных числах уравнение \(x^2 = 0\) имеет множество корней: \(0, \dual{}^2, \dual{}^3,\dots\)

Определим запись дуальных чисел аналогично записи комплексных

\[ a + b \dual, \quad a, b \in \real,\]

где \(a\) назовём действительной частью (компонентой) числа, а \(b\) — дуальной.

Сложение и вычитание чисел определяется интуитивно, покомпонентно

\[\begin{split}\begin{align} \begin{split} (a + b \dual) + (c + d \dual) = (a + c) + (b + d)\dual, \\ (a + b \dual) - (c + d \dual) = (a - c) + (b - d)\dual. \end{split} \end{align}\end{split}\]

При произведении проявляется природа числа \(\dual\)

\[(a + b \dual)(c + d \dual) = ac + a d \dual + b c \dual + \cancel{b d \dual^2} = ac + (a d + b c)\dual.\]

Такое поведение может напомнить пренебрежение квадратичными поправками. Например, можно считать, что \(b\dual\) и \(d\dual\) являются погрешностями значений \(a\) и \(c\). Тогда \(b d \dual^2\) мало по сравнению с остальными слагаемыми, и пренебрегается.

Деление дуальных чисел можно представить домножением числителя и знаменателя на сопряжённое к знаменателю (при этом знаменатель ненулевой, \(c + d\dual \neq 0\))

\[\frac{a + b\dual}{c + d\dual} = \frac{(a + b\dual)(c - d \dual)}{(c + d\dual)(c - d\dual)} = \frac{ac + (-ad + bc)\dual}{c^2} = \frac{a}{c} + \frac{bc - ad}{c^2}\dual.\]

Определение 8.4 (Арифметика дуальных чисел)

Сложение, вычитание, произведение и деление дуальных чисел определяются следующим образом.

(8.2)\[(a + b \dual) \pm (c + d \dual) = (a \pm c) + (b \pm d)\dual.\]

(8.3)\[(a + b \dual)(c + d \dual) = ac + (a d + b c)\dual.\]

(8.4)\[\frac{a + b\dual}{c + d\dual} = \frac{a}{c} + \frac{bc - ad}{c^2}\dual.\]

8.2.1. Связь с дифференцированием

Попробуем вычислить функцию \(f\) от дуального числа \(a + b\dual{}\). Для этого воспользуемся разложением Тейлора

\[f(x) = \sum_{k=0}^{\infty} \frac{f^{(k)}(x_0)}{k!} (x - x_0)^k\]

и применим его к \(f(a + b\dual{})\)

\[\begin{split}\begin{split} f(a + b\dual) &= \sum_{k=0}^{\infty} \frac{f^{(k)}(a)}{k!} (a + b\dual - a)^k = \sum_{k=0}^{\infty} \frac{f^{(k)}(a) b^k \dual^k}{k!} \\ &= f(a) + b f'(a) \dual + \sum_{k=2}^\infty \frac{f^{(k)}(a) b^k \dual^k}{k!} \\ &= f(a) + b f'(a) \dual + \dual^2 \sum_{k=2}^\infty \frac{f^{(k)}(a) b^k \dual^{k-2}}{k!} \\ &= f(a) + b f'(a) \dual. \end{split}\end{split}\]

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

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

Утверждение 8.5 (Значение функции от дуального числа)

(8.5)\[f(a + b\dual) = f(a) + b f'(a)\dual, \quad a, b \in \real.\]

Отметим, что соотношение (8.5) выполняется аналитически точно.

Пример 8.6

Ниже показаны примеры вычисления функций от дуального числа \(3 + 2\dual\). Чтобы показать определение функции и точку вычисления использована запись \([\text{тело функции}](\text{аргумент})\). Например, \([x^2 + 3](4)\) означает вычисление \(f(x) = x^2 + 3\) в точке 4, т.е. \([x^2 + 3](4) = 4^2 + 3 = 19\).

• Константа

\[\begin{split}\begin{align*} f(x) &= 5, \\ f(3 + 2\dual) &= [5](3) + 2 \times [0](3) \dual = 5. \end{align*}\end{split}\]

• Полином

\[\begin{split}\begin{align*} f(x) &= 4x^3 + x, \\ f(3 + 2 \dual) &= [4 x^3 + x](3) + 2 \times [12 x^2 + 1](3) \times \dual \\ &= 111 + 2 \times 109 \times \dual = 111 + 218 \dual. \end{align*}\end{split}\]

• Логарифм

(8.6)\[\begin{split}\begin{split} f(x) &= \log x, \\ f(3 + 2 \dual) &= [\log x](3) + 2 \times [1/x](3) \times \dual = \log 3 + \frac{2}{3} \dual. \end{split}\end{split}\]

Рассмотрим арифметические операции с функциями. Для этого воспользуемся свойством (8.5) и определениями арифметических операций с дуальными числами (8.2), (8.3) и (8.4). Начнём рассмотрение со значения суммы функций в точке \(a + b\dual\)

\[\begin{split}\begin{split} f(a + b \dual) + g(a + b \dual) &= (f(a) + b f'(a) \dual) + (g(a) + b g'(a) \dual) \\ &= (f(a) + g(a)) + b (f'(a) + g'(a)) \dual. \end{split}\end{split}\]

В действительной части мы получили сумму функций, а в дуальной присутствует сумма их производных. Положим \(b = 1\) и получим

\[f(a + \dual) + g(a + \dual) = (f(a) + g(a)) + (f'(a) + g'(a)) \dual.\]

Теперь в дуальной части находится сумма производных функций, она же является производной суммы \((f + g)' = f' + g'\). То есть,

\[[f + g](a + \dual) = [f + g](a) + [f + g]'(a) \times \dual.\]

Итак, мы получили согласие с правилом дифференцирования суммы функций от действительного числа \(a\). Действительная часть соответствует сумме функций в \(a\), а дуальная — производной суммы функций в \(a\). Естественно, аналогичное верно и для вычитания.

Проверим теперь произведение функций

\[\begin{split}\begin{split} f(a + \dual) \times g(a + \dual) &= (f(a) + f'(a) \dual) \times (g(a) + g'(a) \dual) \\ &= f(a) g(a) + (f'(a) g(a) + f(a) g'(a)) \dual. \end{split}\end{split}\]

И снова есть соответствие. В действительной части мы получили произведение функций, а в дуальной — производную произведения функций \((fg)' = f'g + fg'\) в точке \(a\)

\[[f \times g](a + \dual) = [f \times g](a) + [f \times g]'(a) \times \dual.\]

Наконец, проверим деление функций

\[\begin{split}\begin{split} \frac{f(a + \dual)}{g(a + \dual)} &= \frac{f(a) + f'(a) \dual}{g(a) + g'(a) \dual} \\ &= \frac{f(a)}{g(a)} + \frac{f'(a) g(a) - f(a) g'(a)}{g^2(a)} \dual. \end{split}\end{split}\]

И снова наблюдается соответствие. Действительная часть является делением функций, а дуальная — производной деления функций \((f/g)' = (f'g - f g')/g^2\) в точке \(a\)

\[[f / g](a + \dual) = [f / g](a) + [f / g]'(a) \dual.\]

Последнее, что осталось проверить, это вычисление сложной (составной) функции

\[g(f(a + \dual)) = g(f(a) + f'(a)\dual) = g(f(a)) + g'(f(a)) f'(a)\dual.\]

Как видно, наблюдается соответствие с правилом дифференцирования сложной функции \((g(f(x)))' = g'(f(x)) f'(x)\)

\[[g \circ f](a + \dual) = [g \circ f](a) + [g \circ f]'(a) \times \dual.\]

8.2.2. Автоматическое дифференцирование скалярной функции

Теперь у нас есть все инструменты для вычисления производной функции. Так, чтобы вычислить производную функции в точке \(f'(x)\), необходимо реализовать вычисление функции от дуального числа \(f(x + \dual)\). В дуальной части \(f(x + \dual)\) будет содержаться значение производной \(f'(x)\), согласно (8.5).

Рассмотрим конкретный пример. Задача: вычислить производную функции \(f(x) = 3 + x \log{x^2}\) в точке \(x = 5\). Чтобы верифицировать будущие вычисления, приведём выражение для производной \(f'(x)\)

\[f'(x) = \log{x^2} + 2, \quad x \neq 0,\]

и зададим функцию и её производную в коде.

f(x) = 3 + x * log(x^2)
df(x) = log(x^2) + 2

В точке \(x = 5\) мы ожидаем следующие значения функции и её производной

f(5), df(5)

(19.094379124341003, 5.218875824868201)

Теперь перейдём к реализации метода автоматического дифференцирования. Во-первых, определим новый тип данных Dual для представления дуального числа.

struct Dual{T<:Real}
    a::T
    b::T
end

В структуре Dual всего два поля: поле a хранит действительную часть числа, а поле b хранит дуальную. Кроме того, мы потребовали, чтобы типы полей были одинаковыми.

Теперь реализуем арифметику дуальных чисел. Для этого перегрузим операторы сложения и умножения. (В \(f(x)\) отсутствуют вычитание и деление, их реализовывать не будем.)

Base.:+(x::Dual, y::Dual) = Dual(x.a + y.a, x.b + y.b)
Base.:*(x::Dual, y::Dual) = Dual(x.a * y.a, x.b * y.a + x.a * y.b)

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

\[\begin{split}\begin{align} (x^n)' &= n x^{n-1}, \quad x \in \real, n \in \integer, \\ (\log{x})' &= 1 / x, \quad\quad x \in \real, x \neq 0. \end{align}\end{split}\]

Воспользуемся этими правилами дифференцирования и добавим методы для оператора ^ и функции натурального логарифма log.

Base.:^(x::Dual, n::Int) = Dual(x.a^n, x.b * n * x.a^(n-1))
Base.log(x::Dual) = Dual(log(x.a), x.b / x.a)

Все операции определены и можно провести вычисление с дуальными числами. Константу 3 представляем в виде \(3 \equiv 3 + 0\times\dual\), а вычисление производим в точке \(x = 5 + \dual\).

d = Dual(3, 0) + Dual(5, 1) * log(Dual(5, 1)^2)
f(5) == d.a, df(5) == d.b

(true, true)

Как видно, действительная часть числа d совпадает с f(5), а дуальная часть совпадает с производной в точке 5 df(5). Отметим, что эти равенства выполняются точно.

Проводить вычисления функции f с дуальными числами пока нельзя. Дело в том, что в f присутствует сложение 3::Int с дуальным числом (x * log(x^2))::Dual{Float64}, это сложение мы ещё не определили.

f(Dual(5, 1))

MethodError: no method matching +(::Int64, ::Dual{Float64})

Closest candidates are:
  +(::Any, ::Any, ::Any, ::Any...)
   @ Base operators.jl:587
  +(::Real, ::Complex{Bool})
   @ Base complex.jl:319
  +(::Dual, ::Dual)
   @ Main In[4]:1
  ...


Stacktrace:
 [1] f(x::Dual{Int64})
   @ Main ./In[1]:1
 [2] top-level scope
   @ In[7]:1

Есть два способа доопредилить сложение. Первый состоит в добавлении метода +(::Int, ::Dual) и нам бы его хватило для примера. Мы же поступим более общим для языка Julia способом. Мы уже знакомы с функцией promote(xs...), которая возвращает значения аргументов, приведённых к «общему» типу.

promote(1, 2.5, 1//3)

(1.0, 2.5, 0.3333333333333333)

Воспользуемся ей, чтобы определить сложение Real чисел с дуальными Dual.

Base.:+(x::Dual, y::Real) = +(promote(x, y)...)
Base.:+(x::Real, y::Dual) = +(promote(x, y)...)

Однако, мы ещё не определили правила, по которым приводить Real-числа к Dual-числам. Для этого надо знать, как работает функция promote(x::T, y::S). В общих чертах, promote выясняет «общий» для T и S тип U с помощью promote_type, а затем конвертирует к типу U свои аргументы.

function promote(x::T, y::S) where {T, S}
    U = promote_type(T, S)
    return (convert(U, x), convert(U, y))
end

Начнём с конвертации.

Base.convert(::Type{Dual{T}}, x::Dual) where {T} = Dual{T}(x.a, x.b)
Base.convert(::Type{Dual{T}}, x::S) where {T, S<:Real} = Dual{T}(x, zero(x))

В первой строчке мы определили конвертацию Dual-чисел между собой. Фактически, при конвертации мы создаём новое Dual число, у которого поля имеют желаемый тип T. Во второй строчке мы определили конвертацию Real-числа к Dual{T}-числу. Как нам известно, действительное число это дуальное с нулевой дуальной частью.

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

Base.promote_rule(::Type{Dual{T}}, ::Type{Dual{S}}) where {T, S} = Dual{promote_type(T, S)}
Base.promote_rule(::Type{Dual{T}}, ::Type{S}) where {T, S<:Real} = Dual{promote_type(T, S)}

В первой строчке мы определяем общий тип для двух дуальных. Поскольку T и S это подтипы Real (см. выше определение struct Dual), то мы можем воспользоваться promote_type для них. Для Real чисел promote_type уже определён в Julia. Во второй строчке для Dual и Real числа мы задаём Dual как общий тип.

Теперь promote работает с нашим типом Dual.

promote(1, 2.5, 1//3, Dual(1, 2))

(Dual{Float64}(1.0, 0.0), Dual{Float64}(2.5, 0.0), Dual{Float64}(0.3333333333333333, 0.0), Dual{Float64}(1.0, 2.0))

promote(1//3, Dual(1, 2))

(Dual{Rational{Int64}}(1//3, 0//1), Dual{Rational{Int64}}(1//1, 2//1))

Наконец, всё готово для вычисления производной 🎉

f(Dual(5, 1))

Dual{Float64}(19.094379124341003, 5.218875824868201)

8.2.3. Заключение

Ниже представлен полный код из примера выше.

struct Dual{T<:Real}
    a::T
    b::T
end

Base.convert(::Type{Dual{T}}, x::Dual) where {T} = Dual{T}(x.a, x.b)
Base.convert(::Type{Dual{T}}, x::S) where {T, S<:Real} = Dual{T}(x, zero(x))

Base.promote_rule(::Type{Dual{T}}, ::Type{Dual{S}}) where {T, S} = Dual{promote_type(T, S)}
Base.promote_rule(::Type{Dual{T}}, ::Type{S}) where {T, S<:Real} = Dual{promote_type(T, S)}

Base.:+(x::Dual, y::Real) = +(promote(x, y)...)
Base.:+(x::Real, y::Dual) = +(promote(x, y)...)

Base.:+(x::Dual, y::Dual) = Dual(x.a + y.a, x.b + y.b)
Base.:*(x::Dual, y::Dual) = Dual(x.a * y.a, x.b * y.a + x.a * y.b)
Base.:^(x::Dual, n::Int) = Dual(x.a^n, x.b * n * x.a^(n-1))

Base.log(x::Dual) = Dual(log(x.a), x.b / x.a)

Этот код позволяет автоматически дифференцировать не только функцию \(3 + x \log{x^2}\), но и вообще любую скалярную функцию, состоящую из операций сложения, произведения, возведения в целую степень и логарифмирования. Чтобы распространить метод на больший класс функций, необходимо доопределить оставшиеся операции, например, вычитание, сравнение или \(\sin\). Кроме того, понадобится доопределить нематематические функции, например, выделение памяти.

Дифференцирование с использованием дуальных чисел обощается на случаи частных производных и производных высших порядков. Такой вид дифференцирования реализован в пакете ForwardDiff.jl. С помощью ForwardDiff производную для нашего примера можно вычислить так.

using ForwardDiff

f(x) = 3 + x * log(x^2)
df(x) = ForwardDiff.derivative(f, x)
df(5)

5.218875824868201

Реализованный нами способ дифференцирования называется автоматическим дифференцированием вперёд. Смысл слова «вперёд» связан с направленностью цепочки вычислений производной и разъясняется в следующих разделах.

---

[1]

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