Метод Ньютона

Нахождение корня функции

Пусть f ⁣:R→Rf \colon \RR \to \RR. Мы хотим найти корень ff, то есть такую точку x∗∈Rx^* \in \RR, что

f(x∗)=0f(x^*) = 0

Скалярная функция многих переменных

Нахождение корня скалярной функции многих переменных

Пусть f ⁣:Rn→Rf \colon \RR^n \to \RR. Мы хотим найти корень ff, то есть такой вектор x∗∈Rnx^* \in \RR^n, что

f(x∗)=0f(x^*) = 0

Будем опять находить решение итеративно. Пусть xkx_k — какое-то приближённое решение задачи, и мы хотим улучшить это решение. То есть мы хотим найти такое Δx\Delta x, чтобы f(xk+Δx)=0f(x_k + \Delta x) = 0.

По формуле Тейлора

f(xk+Δx)≈f(xk)+∇f(xk)T⋅Δxf(x_k + \Delta x) \approx f(x_k) + \nabla f (x_k) ^\T \cdot \Delta x

Мы очень хотим выбрать такой шаг Δx\Delta x, чтобы f(xk+Δx)=0f(x_k + \Delta x) = 0. Используем предыдущую формулу

f(xk)+∇f(xk)T⋅Δx≈0f(x_k) + \nabla f (x_k) ^\T \cdot \Delta x \approx 0

Теперь нам нужно выбрать вектор Δx\Delta x, по которому нужно двигаться. Учитывая, что градиент — направление наибыстрейшего подъема, кажется естественным выбрать вектор Δx\Delta x так, чтобы он был коллинеарен градиенту

Δx=−α⋅∇f(xk)\Delta x = - \alpha \cdot \nabla f (x_k)

Осталось решить уравнение f(xk)+∇f(xk)T⋅(−α⋅∇f(xk))=0f(x_k) + \nabla f (x_k) ^\T \cdot (- \alpha \cdot \nabla f (x_k)) = 0. Получаем

α=f(xk)∇f(xk)⋅∇f(xk)T=f(xk)∥∇f(xk)∥22\alpha = \frac{f(x_k)}{\nabla f (x_k) \cdot \nabla f (x_k) ^\T} = \frac{f(x_k)}{ \|\nabla f (x_k)\|_2^2}

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

xk+1=xk−f(xk)∥∇f(xk)∥22⋅∇f(xk)x_{k+1} = x_k - \frac{f(x_k)}{\| \nabla f (x_k) \|_2^2} \cdot \nabla f (x_k)

Напишем алгоритм в псевдокоде. Для этого нам нужно уметь вычислять функцию ff и её градиент ∇f\nabla f.

function f(vector x) -> real: ... function gradient(vector[real] x) -> vector[real]: ... function solve_by_newton(f, gradient) -> real: select x: real while f(x) >= EPSILON: x = f(x) / (norm(gradient(x)) ** 2) * gradient(x) return x

Далёкие от корня значения

Наш метод будет работать только если начальное значение x0x_0 уже близко к корню. Если это не так, то у нас точно возникнет большое количество проблем.

Вместо поиска корня ff можно найти минимум функции 12f(x)2\frac{1}{2} f(x)^2. В принципе, для этого можно использовать и другие алгоритмы, но мы посмотрим на алгоритм Ньютона.

∇(12⋅f2)=f⋅∇fиH⁡(12⋅f2)=∇f⋅∇fT+f⋅H⁡f\nabla \left( \frac{1}{2} \cdot f^2 \right) = f \cdot \nabla f \quad\text{и}\quad \hess \left( \frac{1}{2} \cdot f^2 \right) = \nabla f \cdot \nabla f ^\T + f \cdot \hess f

Многомерные функции

Пусть F ⁣:Rn→RnF \colon \RR^n \to \RR^n — многомерная функция. Мы хотим научиться решать уравнение F(x)=0F(x) = \0.

Действовать придётся итеративно, начиная с какого-то приближения x0∈Rnx_0 \in \RR^n.

Вспомним формулу для линеаризации, где h∈Rnh \in \RR^n — какой-то небольшой шаг:

F(x+h)=F(x)+J⁡F(x)⋅h+o(∥h∥)F(x + h) = F(x) + \jacobi F (x) \cdot h + o \bigl( \| h \| \bigr)

Попробуем выбрать hh так, чтобы F(x+h)=0F(x+h) = \0. Используем приближение

0=F(x+h)≈F(x)+J⁡F(x)⋅h\0 = F(x+h) \approx F(x) + \jacobi F (x) \cdot h

У нас получилась система линейных уравнений относительно hh. Решаем, находя обратную матрицу

J⁡F(x)⋅h=−F(x)  ⟹  h=−(J⁡F(x))−1⋅F(x)\jacobi F (x) \cdot h = - F(x) \implies h = - \bigl( \jacobi F (x) \bigr)^{-1} \cdot F(x)

Получается, для решения уравнения F(x)=0F(x) = \0, нужно выполнять итеративный процесс, начиная с точки x0x_0

xj+1=xj−(J⁡F(xj))−1⋅F(xj)x_{j+1} = x_j - \bigl( \jacobi F (x_j) \bigr)^{-1} \cdot F(x_j)