Метод Ньютона
Метод Ньютона
Метод Ньютона является одним из самых распространенных методов решения систем нелинейных алгебраических уравнений. Для начала мы рассмотрим классический вывод метода Ньютона через разложение в ряд Тейлора, после чего взглянем на метод Ньютона как на частный случай метода простой итерации. Как и раньше, мы в деталях рассмотрим одномерный случай, после чего обобщим результаты для многомерного случая.
Пусть нам дана функция $f(x) \in C^2[a; b]$ и точка $x^{(0)} \in [a; b]$ является приближением к корню $x^*$ функции $f(x)$, т.е. $f(x^*) = 0$. Разложим функцию $f(x)$ в ряд Тейлора в точке $x^{(0)}$:
\begin{equation}
f(x) = f(x^{(0)}) + f'(x^{(0)}) (x - x^{(0)}) + \frac{f''(\xi)}{2}(x - x^{(0)})^2,
\end{equation}
где $\xi \in (x; x^{(0)})$. Метод Ньютона основан на предположении, что дистанция $|x^* - x^{(0)}|$ от приближения до корня функции достаточно мала. Тогда, вычислив ряд в точке $x^*$, мы можем отбросить квадратичный член:
\begin{equation}
f(x^*) \approx f(x^{(0)}) + f'(x^{(0)}) (x^* - x^{(0)}).
\end{equation}
Так как $f(x^*) = 0$, мы получаем оценку для $x^*$, т.е. более точное приближение:
\begin{equation}
\begin{split}
&0 \approx f(x^{(0)}) + f'(x^{(0)}) (x^* - x^{(0)}) \\
\Longrightarrow \quad &x^* \approx x^{(0)} - \frac{f(x^{(0)})}{f'(x^{(0)})}.
\end{split}
\end{equation}
Таким образом, можно сформулировать следующий итерационный метод, называемый
методом Ньютона:
\begin{equation}
x^{(k)} = x^{(k-1)} - \frac{f(x^{(k-1)})}{f'(x^{(k-1)})}.
\end{equation}
Несложно заметить, что полученный итерационный метод является формой метода простой итерации:
\begin{equation}
\label{eq:newton_method_as_fixed_point_iter}
x^{(k)} = g(x^{(k-1)}),
\end{equation}
где $g(x) = x - \frac{f(x)}{f'(x)}$. Сравнение полученной функции $g(x)$ с общим выражением $g(x) = x - \Omega(x) f(x)$ показывает, что метода Ньютона выбирает в качестве функции $\Omega(x) = -\frac{1}{f'(x)}$.
Отметим две важных особенности метода Ньютона:
- начальное приближение $x^{(0)}$ должно быть достаточно близким к корню $x^*$;
- градиент $f'(x^{(k)})$ должен быть отличен от нуля для любого $k$. Более того, метод Ньютона тем эффективнее, чем дальше $f'(x^*)$ от нуля.
Первое условие накладывает серьезные ограничения на сходимость метода Ньютона для произвольных $x^{(0)}$. Следующая теорема доказывает существование некоторой окрестности $x^*$, внутри которой сходимость метода Ньютона является безусловной.
Пусть $f(x) \in C^2[a; b]$ и существует такое $x^* \in (a; b)$, что $f(x^*) = 0$ и $f'(x^*) \neq 0$. Тогда существует такое $\delta > 0$, что последовательность $\{x^{(k)}\}_{k=0}^{\infty}$, генерируемая метода Ньютона, сходится к $x^*$ для любого $x^{(0)} \in [x^* - \delta; x^* + \delta]$.
Рассмотрим метода Ньютона как метод простой итерации в соответствии с \eqref{eq:newton_method_as_fixed_point_iter}. Тогда, следуя теореме \ref{th:conv_condition_fixed_point_iter}, необходимо найти интервал $D = [x^* - \delta; x^* + \delta]$, который функция $g(x)$ отображает в себя, и в котором $|g'(x)| \leq \gamma$, где $\gamma \in (0; 1)$.
Так как $f'$ является непрерывной и $f'(x^*) \neq 0$, существует такое $\delta_1 > 0$, что $f'(x) \neq 0$ в замкнутой $\delta_1$-окрестности точки $x^*$, т.е. для $x \in [x^* - \delta_1; x^* + \delta_1] \subset [a; b]$. Тогда непрерывная производная $g'(x) \in C^1[a; b]$ также будет существовать в этом интервале:
\begin{equation}
\begin{split}
g'(x) &= 1 - \frac{\left(f'(x)\right)^2 - f(x) f''(x)}{\left(f'(x)\right)^2} \\
&= \frac{f(x) f''(x)}{\left(f'(x)\right)^2}.
\end{split}
\end{equation}
Найдем значение $g'(x)$ в точке $x^*$:
\begin{equation}
\begin{split}
g'(x^*) &= \frac{f(x^*) f''(x^*)}{\left(f'(x^*)\right)^2} \\
&= 0.
\end{split}
\end{equation}
Тогда существует такое $\delta \in (0; \delta_1)$, что $|g'(x)| \leq \gamma$ для $\gamma \in (0; 1)$ в замкнутой $\delta$-окрестности точки $x^*$, т.е. для $x \in [x^* - \delta; x^* + \delta] \subset [a; b]$.
Теперь докажем, что полученный интервал $D = [x^* - \delta; x^* + \delta]$ отображается функцией $g(x)$ сам в себя. По теореме Лагранжа о среднем значении для $x \in D$ существует такое $\xi \in (x; x^*)$ или $\xi \in (x^*; x)$, что:
\begin{equation}
\begin{split}
&|g'(\xi)| = \frac{|g(x) - g(x^*)|}{|x - x^*|} = \frac{|g(x) - x^*|}{|x - x^*|} \\
\Longrightarrow \quad &|g(x) - x^*| = |g'(\xi)| \cdot |x - x^*| \leq \gamma |x - x^*| < |x - x^*|.
\end{split}
\end{equation}
Так как $|x - x^*| \leq \delta$, мы получаем $|g(x) - x^*| \leq \delta$, из чего следует, что $g(x)$ отображает $D$ в себя. Тогда по теореме о сходимости метода простой итерации последовательность $\{x^{(k)}\}_{k=0}^{\infty}$ сходится к $x^*$ для любого $x^{(0)} \in [x^* - \delta; x^* + \delta]$.
Прежде чем перейти к рассмотрению скорости сходимости метода Ньютона, дадим строгое определение этому термину.
Пусть последовательность $\{x^{(k)}\}_{k=0}^{\infty}$ сходится к $x^*$. Тогда если существуют такие $\lambda, \alpha \in \mathbb{R}$, что:
\begin{equation}
\lim_{k \to \infty} \frac{|x^{(k+1)} - x^*|}{|x^{(k)} - x^*|^\alpha} = \lambda,
\end{equation}
то метод, генерирующий данную последовательность, обладает сходимостью степени $\alpha$. При $\alpha = 1, \lambda \in (0; 1)$ сходимость называется линейной, при $\alpha = 1, \lambda = 0$ сверхлинейной и при $\alpha = 2$ квадратичной.
Следующая теорема доказывает, что сходящийся метод простой итерации в общем случае сходится только линейно.
Пусть $D = [a; b]$, $g \in C(D)$, $g(D) \subset D$ и существует производная $g'(x)$ для $x \in (a; b)$ с таким $\gamma \in (0; 1)$, что $|g'(x)| \leq \gamma$ для любого $x \in (a; b)$. Если $g'(x^*) \neq 0$, то для любого $x^{(0)} \in D$ последовательность $\{x^{(k)}\}_{k=0}^{\infty}$, сгенерированная итерацией $x^{(k)} = g(x^{(k-1)})$, сходится линейно к единственной неподвижной точке $x^* \in [a; b]$.
По теореме о сходимости метода простой итерации данная последовательность действительно сходится к $x^*$, так что необходимо доказать только факт линейной сходимости. В силу непрерывности $g'(x)$ по теореме Лагранжа о среднем значении для любого $k > 0$ существует такое $\xi \in (x^{(k)}; x^*)$ или $\xi_k \in (x^*; x^{(k)})$, что:
\begin{equation}
\begin{split}
g'(\xi_k) &= \frac{g(x^{(k)}) - g(x^*)}{x^{(k)} - x^*} \\
&= \frac{x^{(k + 1)} - x^*}{x^{(k)} - x^*}.
\end{split}
\end{equation}
Так как $\{x^{(k)}\}_{k=0}^{\infty}$ сходится к $x^*$, к $x^*$ сходится и $\{\xi^{(k)}\}_{k=0}^{\infty}$. Более того, благодаря непрерывности $g'(x)$, мы имеем
\begin{equation}
\begin{split}
&\lim_{k \to \infty} g'(\xi_k) = g'(x^*) \\
\Longrightarrow \quad &\lim_{k \to \infty} \frac{x^{(k + 1)} - x^*}{x^{(k)} - x^*} = g'(x^*) \\
\Longrightarrow \quad &\lim_{k \to \infty} \frac{|x^{(k + 1)} - x^*|}{|x^{(k)} - x^*|} = |g'(x^*)|.
\end{split}
\end{equation}
Тогда, так как $0 < |g'(x^*)| < 1$, последовательность $\{x^{(k)}\}_{k=0}^{\infty}$ сходится к $x^*$ линейно.
Рассмотренная теорема предполагает, что сходимость быстрее линейной возможна только при $g'(x^*) = 0$. Действительно, как доказывает следующая теорема, сходимость в таком случае становится как минимум квадратичной.
Пусть $x^* \in (a; b)$ является неподвижной точкой функции $g(x)$, т.е. $x^* = g(x^*)$. Пусть также $g'(x^*) = 0$ и $|g''(x)| < M$ для $x \in (a; b)$. Тогда существует такое $\delta > 0$, что последовательность $\{x^{(k)}\}_{k=0}^{\infty}$, генерируемая итерацией $x^{(k)} = g(x^{(k-1)})$, $k > 1$, сходится к $x^*$ в области $[x^* - \delta; x^* + \delta]$ как минимум квадратически. Более того, для асимптотически больших $k$, верно следующее неравенство:
\begin{equation}
|x^{(k+1)} - x^*| < \frac{M}{2} |x^{(k)} - x^*|^2.
\end{equation}
По аналогии с доказательством теоремы о сходимости метода Ньютона выберем такое $\delta > 0$ и $\gamma \in (0; 1)$, что $|g'(x)| \leq \gamma$ для $x \in [x^* - \delta; x^* + \delta] \subset [a; b]$. В том же доказательстве мы установили, что последовательность $\{x^{(k)}\}_{k=0}^{\infty}$ для $x^{(0)} \in [x^* - \delta; x^* + \delta]$ будем так же содержаться в $[x^* - \delta; x^* + \delta]$. Разложим $g(x)$ в ряд Тейлора в точке $x^*$:
\begin{equation}
g(x) = g(x^*) + g'(x^*) (x - x^*) + \frac{g''(\xi)}{2} (x - x^*)^2
\end{equation}
где $x \in [x^* - \delta; x^* + \delta]$ и, следовательно, $\xi \in [x^*; x]$ или $\xi \in [x; x^*]$. Используя установки теоремы, а именно $g(x^*) = x^*$ и $g'(x^*) = 0$, мы имеем:
\begin{equation}
g(x) = x^* + \frac{g''(\xi)}{2} (x - x^*)^2.
\end{equation}
Тогда для $x = x^{(k)}$ мы получаем:
\begin{equation}
\begin{split}
&g(x^{(k)}) = x^* + \frac{g''(\xi_k)}{2} (x^{(k)} - x^*)^2 \\
\Longrightarrow \quad &x^{(k+1)} - x^* = \frac{g''(\xi_k)}{2} (x^{(k)} - x^*)^2.
\end{split}
\end{equation}
Так как $|g'(x)| \leq \gamma$ для $x \in [x^* - \delta; x^* + \delta]$, последовательность $\{x^{(k)}\}_{k=0}^{\infty}$ сходится к $x^*$ $\Longrightarrow$ $\{\xi^{(k)}\}_{k=0}^{\infty}$ сходится к $x^*$ $\Longrightarrow$ $\lim_{k \to \infty} g''(\xi_k) = g''(x^*)$. Тогда имеем:
\begin{equation}
\begin{split}
&\lim_{k \to \infty} g''(\xi_k) = g''(x^*) \\
\Longrightarrow \quad &\lim_{k \to \infty} \frac{x^{(k+1)} - x^*}{(x^{(k)} - x^*)^2} = \frac{g''(x^*)}{2} \\
\Longrightarrow \quad &\lim_{k \to \infty} \frac{|x^{(k+1)} - x^*|}{|x^{(k)} - x^*|^2} = \frac{|g''(x^*)|}{2},
\end{split}
\end{equation}
что доказывает квадратичную сходимость. Так как $|g''(x^*)| < M$, мы также получаем следующую оценку:
\begin{equation}
|x^{(k+1)} - x^*| < \frac{M}{2} |x^{(k)} - x^*|^2.
\end{equation}
Эта теорема дает нам возможность построить такую функцию $g(x)$ в соответствии с определением $g(x) = x - \Omega(x) f(x)$ для нахождения корня функции $f(x)$, что метод простой итерации будет сходиться квадратично. Как мы убедились, для этого необходимо, чтобы $g'(x^*) = 0$. Тогда для одномерного случая имеем:
\begin{equation}
\begin{split}
&g'(x^*) = 0 \\
\Longrightarrow \quad &1 - \Omega'(x^*) f(x^*) - \Omega(x^*) f'(x^*) = 0 \\
\Longrightarrow \quad &1 - \Omega(x^*) f'(x^*) = 0 \\
\Longrightarrow \quad &\Omega(x^*) = \frac{1}{f'(x^*)}.
\end{split}
\end{equation}
Тогда, выбрав $\Omega(x) = \frac{1}{f'(x)}$, мы получаем уже рассмотренный метод Ньютона:
\begin{equation}
x^{(k)} = g(x^{(k)}) = x^{(k-1)} - \frac{f(x^{(k-1)})}{f'(x^{(k-1)})}.
\end{equation}
Таким образом, метод Ньютона является ``оптимальным'' методом простой итерации и имеет квадратичную сходимость.
Рассмотрев одномерный случай, мы можем обобщить результаты на многомерный случай. Мы сразу обратимся к методу Ньютона в форме метода простой итерации и сформулируем без доказательства теорему о методе простой итерации с квадратичной сходимостью.
Пусть $\bm{x^*}$ является неподвижной точкой функции $\bm{g}(\bm{x})$, т.е. $\bm{g}(\bm{x^*}) = \bm{x}^*$, и пусть функция $\bm{g}(\bm{x})$ удовлетворяет следующим условиям:
- существует такое $\delta$, что производные $\dfrac{\partial g_i}{\partial x_j}$ непрерывны в $D = \{\bm{x} \quad | \quad ||\bm{x} - \bm{x^*}|| < \delta \}$;
- вторые производные $\dfrac{\partial^2 g_i}{\partial x_j \partial x_k}$ непрерывны и ограничены, т.е. $\left|\dfrac{\partial^2 g_i}{\partial x_j \partial x_k}\right| \leq M$, в $D = \{\bm{x} \quad | \quad ||\bm{x} - \bm{x^*}|| < \delta \}$;
- $\dfrac{\partial g_i(\bm{x^*})}{\partial x_j} = 0$,
для $i = 1, \dots, n$, $j = 1, \dots, n$ и $k = 1, \dots, n$. Тогда существует такое $\hat{\delta} \leq \delta$, что последовательность $\{\bm{x}^{(k)}\}_{k=0}^{\infty}$, генерируемая итерацией $\bm{x}^{(k)} = \bm{g}(\bm{x}^{(k-1)})$, сходится квадратично к $\bm{x^*}$ для любого $\bm{x}^{(0)}$, удовлетворяющего $||\bm{x}^{(0)} - \bm{x^*}|| < \hat{\delta}$, так, что верно следующее неравенство:
\begin{equation}
||\bm{x}^{(k)} - \bm{x^*}||_{\infty} \leq \frac{n^2 M}{2} ||\bm{x}^{(k-1)} - \bm{x^*}||_{\infty}^2, \quad k \geq 1.
\end{equation}
Для вывода метода Ньютона для многомерного случая рассмотрим формулу $\bm{g}(\bm{x}) = \bm{x} - \bm{\Omega}(\bm{x}) \bm{f}(\bm{x})$ покоординатно:
\begin{equation}
g_i(\bm{x}) = x_i - \sum_{j=1}^n \omega_{ij}(\bm{x}) f_j(\bm{x}).
\end{equation}
По условиям теоремы выше для построения метода простой итерации с квадратичной сходимостью нам необходимо обнулить производные функций $g_i(\bm{x})$ в точке $\bm{x^*}$. Найдем для начала соответствующие производные:
\begin{equation}
\frac{\partial g_i (\bm{x})}{\partial x_j}
=
\begin{cases}
\displaystyle
1 - \sum_{j=1}^n \frac{\partial \omega_{ij}(\bm{x})}{\partial x_j} f_j(\bm{x}) - \sum_{j=1}^n \omega_{ij}(\bm{x}) \frac{\partial f_j(\bm{x})}{\partial x_j}, &\quad i = j, \\
\displaystyle
-\sum_{j=1}^n \frac{\partial \omega_{ij}(\bm{x})}{\partial x_j} f_j(\bm{x}) - \sum_{j=1}^n \omega_{ij}(\bm{x}) \frac{\partial f_j(\bm{x})}{\partial x_j}, &\quad i \neq j. \\
\end{cases}
\end{equation}
Тогда обнуление производных дает уравнения:
\begin{align}
\sum_{j=1}^n \omega_{ij}(\bm{x^*}) \frac{\partial f_j(\bm{x^*})}{\partial x_j} = 1, &\quad i = j, \\
\sum_{j=1}^n \omega_{ij}(\bm{x^*}) \frac{\partial f_j(\bm{x^*})}{\partial x_j} = 0, &\quad i \neq j.
\end{align}
где мы использовали факт $f_j(\bm{x^*}) = 0$. Легко убедиться, что в матричном виде это эквивалентно следующему выражению:
\begin{equation}
\Omega(\bm{x^*}) \bm{J}(\bm{x^*}) = \bm{E}
\end{equation}
где матрица $\bm{J}$ называется
матрицей Якоби и имеет следующую форму:
\begin{equation}
\bm{J}(\bm{x})
=
\begin{bmatrix}
\frac{\partial f_1(\bm{x})}{\partial x_1} & \dots & \frac{\partial f_1(\bm{x})}{\partial x_n} \\
\vdots & \ddots & \vdots \\
\frac{\partial f_n(\bm{x})}{\partial x_1} & \dots & \frac{\partial f_n(\bm{x})}{\partial x_n}
\end{bmatrix}.
\end{equation}
Таким образом, для квадратичной сходимости в качестве $\Omega(\bm{x})$ мы можем выбрать обратную матрицу Якоби:
\begin{equation}
\Omega(\bm{x}) = \bm{J}^{-1}(\bm{x}),
\end{equation}
что дает формулировку
метода Ньютона для систем нелинейных алгебраических уравнений:
\begin{equation}
\bm{x}^{(k)} = \bm{g}(\bm{x}^{(k)}) = \bm{x}^{(k-1)} - \bm{J}^{-1}(\bm{x}^{(k-1)}) \bm{f}(\bm{x}^{(k-1)}).
\end{equation}
На практике обратная от матрицы Якоби не вычисляется, вместо чего используется следующая двухшаговая процедура:
\begin{align}
\bm{J}(\bm{x}^{(k-1)}) \bm{y}^{(k-1)} = \bm{f}(\bm{x}^{(k-1)}), \\
\bm{x}^{(k)} = \bm{x}^{(k-1)} - \bm{y}^{(k-1)},
\end{align}
где $\bm{y}^{(k-1)}$ находится через решение первого уравнения.
Заметим, что в такой форме метод Ньютона требует на каждой итерации $n$ вычислений значений векторной функции $\bm{f}(\bm{x}^{(k-1)})$, $n^2$ вычислений значений матрицы Якоби $\bm{J}(\bm{x}^{(k-1)})$ и $O(n^3)$ операций для нахождения $\bm{y}^{(k-1)}$. Подобная алгоритмическая сложность не позволяет использовать метод Ньютона в его стандартной форме для систем большой размерности.
Квазиньютоновские методы
Квазиньютоновские методы представляют собой класс методов, уменьшающих алгоритмическую сложность метода Ньютона за счет аппроксимации матрицы Якоби тем или иных способом. Уменьшение алгоритмической сложности, однако, приводит к замене квадратичной сходимости на сверхлинейную. Более того, так как аппроксимация матрицы Якоби в сущности является заменой точных производных скалярных функций на аппроксимации путем численного дифференцирования, квазиньютоновские методы являются вычислительно неустойчивыми и могут накапливать погрешность округления с каждой итерацией.
В этом разделе мы рассмотрим один из квазиньютоновских методов, а именно метод Бройдена, который обобщает метод хорд на многомерный случай. Предположим, что $\bm{x}^{(0)}$ является начальным приближением и следующее приближение $\bm{x}^{(1)}$ было рассчитано по методу Ньютона. Это также означает, что мы уже рассчитали матрицу Якоби в точке $\bm{x}^{(0)}$, т.е. $\bm{J}(\bm{x}^{(0)})$. Для нахождения аппроксимации матрицы Якоби в точке $\bm{x}^{(1)}$ нам необходимо сначала взглянуть на одномерный случай. Рассмотрим формулу для численного дифференцирования первого порядка:
\begin{equation}
f'(x^{(1)}) \approx \frac{f(x^{(1)}) - f(x^{(0)})}{x^{(1)} - x^{(0)}}.
\end{equation}
Обобщить это выражение на многомерный случай можно, если переписать его в следующей форме:
\begin{equation}
f'(x^{(1)}) \left( x^{(1)} - x^{(0)} \right) \approx f(x^{(1)}) - f(x^{(0)}),
\end{equation}
из чего следует:
\begin{equation}
\bm{J}(\bm{x}^{(1)}) \left( \bm{x}^{(1)} - \bm{x}^{(0)} \right) \approx \bm{f}(\bm{x}^{(1)}) - \bm{f}(\bm{x}^{(0)}).
\end{equation}
Введем обозначение для аппроксимации матрицы Якоби: $\bm{A}_1 \approx \bm{J}(\bm{x}^{(1)})$. Тогда мы получаем следующее равенство:
\begin{equation}
\bm{A}_1 \left( \bm{x}^{(1)} - \bm{x}^{(0)} \right) = \bm{f}(\bm{x}^{(1)}) - \bm{f}(\bm{x}^{(0)}).
\end{equation}
Это выражение не позволяет однозначно определить матрицу $\bm{A}_1$, так как оно определяет поведение матрицы лишь в направлении одного вектора $\bm{x}^{(1)} - \bm{x}^{(0)}$ в $n$-мерном пространстве. Бройден логично предположил, что действие матрицы $\bm{A}_1$ на все вектора $\bm{z}$, ортогональные вектору $\bm{x}^{(1)} - \bm{x}^{(0)}$, должны быть аналогичны действию на них матрицы $\bm{J}(\bm{x}^{(0)})$:
\begin{equation}
\bm{A}_1 \bm{z} = \bm{J}(\bm{x}^{(0)}) \bm{z}, \quad \langle \bm{x}^{(1)} - \bm{x}^{(0)}, \bm{z} \rangle = 0.
\end{equation}
Эти два условия однозначно задают матрицу $\bm{A}_1$, выражение для которой будет тогда иметь вид:
\begin{equation}
\bm{A}_1 = \bm{J}(\bm{x}^{(0)}) + \frac{\left[ \bm{f}(\bm{x}^{(1)}) - \bm{f}(\bm{x}^{(0)}) - \bm{J}(\bm{x}^{(0)}) (\bm{x}^{(1)} - \bm{x}^{(0)}) \right] \left[\bm{x}^{(1)} - \bm{x}^{(0)} \right]^T}{\langle \bm{x}^{(1)} - \bm{x}^{(0)}, \bm{x}^{(1)} - \bm{x}^{(0)} \rangle}.
\end{equation}
Действительно, легко убедиться, что такая матрица удовлетворяет заданным условиям. Например, для вектора $\bm{z}$, удовлетворяющего $\langle \bm{x}^{(1)} - \bm{x}^{(0)}, \bm{z} \rangle = 0$, мы имеем:
\begin{equation}
\begin{split}
\bm{A}_1 \bm{z} &= \bm{J}(\bm{x}^{(0)}) \bm{z} + \frac{\left[ \bm{f}(\bm{x}^{(1)}) - \bm{f}(\bm{x}^{(0)}) - \bm{J}(\bm{x}^{(0)}) (\bm{x}^{(1)} - \bm{x}^{(0)}) \right] \langle \bm{x}^{(1)} - \bm{x}^{(0)}, \bm{z} \rangle}{\langle \bm{x}^{(1)} - \bm{x}^{(0)}, \bm{x}^{(1)} - \bm{x}^{(0)} \rangle} \\
&= \bm{J}(\bm{x}^{(0)}) \bm{z},
\end{split}
\end{equation}
в то время как для $\bm{x}^{(1)} - \bm{x}^{(0)}$ перемножение дает:
\begin{equation}
\begin{split}
\bm{A}_1 \left( \bm{x}^{(1)} - \bm{x}^{(0)} \right) &= \bm{J}(\bm{x}^{(0)}) \left( \bm{x}^{(1)} - \bm{x}^{(0)} \right) + \\
&\quad \quad + \frac{\left[ \bm{f}(\bm{x}^{(1)}) - \bm{f}(\bm{x}^{(0)}) - \bm{J}(\bm{x}^{(0)}) (\bm{x}^{(1)} - \bm{x}^{(0)}) \right] \langle \bm{x}^{(1)} - \bm{x}^{(0)}, \bm{x}^{(1)} - \bm{x}^{(0)} \rangle}{\langle \bm{x}^{(1)} - \bm{x}^{(0)}, \bm{x}^{(1)} - \bm{x}^{(0)} \rangle} \\
&= \bm{J}(\bm{x}^{(0)}) \left( \bm{x}^{(1)} - \bm{x}^{(0)} \right) + \left[ \bm{f}(\bm{x}^{(1)}) - \bm{f}(\bm{x}^{(0)}) - \bm{J}(\bm{x}^{(0)}) (\bm{x}^{(1)} - \bm{x}^{(0)}) \right] \\
&= \bm{f}(\bm{x}^{(1)}) - \bm{f}(\bm{x}^{(0)}).
\end{split}
\end{equation}
Используя подобную аппроксимацию для матрицы Якоби, модифицированный метод Ньютона имеет вид:
\begin{align}
\bm{A}_k &= \bm{A}_{k-1} + \frac{\left( \bm{y}_k - \bm{A}_{k-1} \bm{s}_k \right) \bm{s}_k^T}{||\bm{s}_k||_2^2}, \\
\bm{A}_k \bm{w}^{(k)} &= \bm{f}(\bm{x}^{(k)}), \\
\bm{x}^{(k+1)} &= \bm{x}^{(k)} - \bm{w}^{(k)},
\end{align}
где $\bm{y}_k = \bm{f}(\bm{x}^{(k)}) - \bm{f}(\bm{x}^{(k-1)})$, $\bm{s}_k = \bm{x}^{(k)} - \bm{x}^{(k-1)}$. Даже после такого упрощения нам все еще необходимо $O(n^3)$ операций для нахождения $\bm{w}^{(k)}$, так что полученный метод в его текущей форме не имеет преимуществ перед методом Ньютона. Его можно улучшить, если сформулировать рекуррентное соотношение для обратной матрицы $\bm{A}_k^{-1}$. В этом нам поможет
формула Шермана–Моррисона, которая гласит, что для невырожденной матрицы $\bm{A}$ и векторов $\bm{x}, \bm{y}$ верно следующее выражение:
\begin{equation}
\left( \bm{A} + \bm{x} \bm{y}^T \right)^{-1} = \bm{A}^{-1} - \frac{\bm{A}^{-1} \bm{x} \bm{y}^T \bm{A}^{-1}}{\sigma},
\end{equation}
где $\sigma = 1 + \bm{y}^T \bm{A}^{-1} \bm{x}$. Убедимся, что формула действительно является тождественно верной. Домножение слева на $\bm{A} + \bm{x} \bm{y}^T$ дает:
\begin{equation}
\begin{split}
&\bm{E} = \left( \bm{A} + \bm{x} \bm{y}^T \right) \bm{A}^{-1} - \frac{\left( \bm{A} + \bm{x} \bm{y}^T \right) \bm{A}^{-1} \bm{x} \bm{y}^T \bm{A}^{-1}}{\sigma} \\
\Longrightarrow \quad &\bm{x} \bm{y}^T \bm{A}^{-1} - \frac{\bm{x} \bm{y}^T \bm{A}^{-1} + \bm{x} \bm{y}^T \bm{A}^{-1} \bm{x} \bm{y}^T \bm{A}^{-1}}{\sigma} = \bm{O}.
\end{split}
\end{equation}
Тогда домножение справа на $\bm{x}$ и подстановка $\bm{y}^T \bm{A}^{-1} \bm{x} = \sigma - 1$ дает:
\begin{equation}
\begin{split}
&\bm{x} \left( \bm{y}^T \bm{A}^{-1} \bm{x} - \frac{\bm{y}^T \bm{A}^{-1} \bm{x} + \bm{y}^T \bm{A}^{-1} \bm{x} \bm{y}^T \bm{A}^{-1} \bm{x}}{\sigma} \right) = \bm{0}, \\
\Longrightarrow \quad &\bm{x} \left( \sigma - 1 - \frac{\sigma - 1 + (\sigma - 1)^2}{\sigma} \right) = \bm{0},
\end{split}
\end{equation}
где, как несложно убедиться, выражение в скобках тождественно обращается в ноль.
Найдем выражение для $\bm{A}_k^{-1}$, используя формулу Шермана–Моррисона:
\begin{equation}
\begin{split}
\bm{A}_k^{-1} &= \left( \bm{A}_{k-1} + \frac{\left( \bm{y}_k - \bm{A}_{k-1} \bm{s}_k \right) \bm{s}_k^T}{||\bm{s}_k||_2^2} \right)^{-1} \\
&= \bm{A}_{k-1}^{-1} - \frac{\bm{A}_{k-1}^{-1} \left( \bm{y}_k - \bm{A}_{k-1} \bm{s}_k \right) \bm{s}_k^T \bm{A}_{k-1}^{-1}}{||\bm{s}_k||_2^2 \left( 1 + \bm{s}_k^T \bm{A}_{k-1}^{-1} \dfrac{\bm{y}_k - \bm{A}_{k-1} \bm{s}_k}{||\bm{s}_k||_2^2} \right)} \\
&= \bm{A}_{k-1}^{-1} - \frac{\bm{A}_{k-1}^{-1} \left( \bm{y}_k - \bm{A}_{k-1} \bm{s}_k \right) \bm{s}_k^T \bm{A}_{k-1}^{-1}}{||\bm{s}_k||_2^2 + \bm{s}_k^T \bm{A}_{k-1}^{-1} (\bm{y}_k - \bm{A}_{k-1} \bm{s}_k)} \\
&= \bm{A}_{k-1}^{-1} - \frac{\left( \bm{A}_{k-1}^{-1} \bm{y}_k - \bm{s}_k \right) \bm{s}_k^T \bm{A}_{k-1}^{-1}}{\bm{s}_k^T \bm{A}_{k-1}^{-1} \bm{y}_k},
\end{split}
\end{equation}
что позволяет рассчитать $\bm{A}_k^{-1}$, используя $\bm{A}_{k-1}^{-1}$, за $O(n^2)$ операций. Таким образом, мы получаем следующие рекуррентные соотношения для
метода Бройдена:
\begin{align}
\bm{A}_k^{-1} &= \bm{A}_{k-1}^{-1} - \frac{\left( \bm{A}_{k-1}^{-1} \bm{y}_k - \bm{s}_k \right) \bm{s}_k^T \bm{A}_{k-1}^{-1}}{\bm{s}_k^T \bm{A}_{k-1}^{-1} \bm{y}_k}, \\
\bm{x}^{(k+1)} &= \bm{x}^{(k)} - \bm{A}_k^{-1} \bm{f}(\bm{x}^{(k)}).
\end{align}
Метод градиентного спуска
Метод градиентного спуска, используемый для нахождения минимума нелинейной целевой функции, может быть использован и для нахождения корней нелинейной системы алгебраических уравнений так же, как и в случае со СЛАУ (см.
вывод метода градиентого спуска в разделе). Он обладает линейной сходимостью, вследствие чего редко используется для нахождения окончательного значения корня. Чаще его роль заключается в том, чтобы найти подходящее начальное приближение для метода Ньютона. Подобная страгегия работает благодаря тому, что метод градиентного спуска накладывает гораздо меньше ограничений на бассейн притяжения (т.е. область, которой должно принадлежать начальное приближение).
Так как в случае нелинейной системы вида $\bm{f}(\bm{x}) = \bm{0}$ вектор невязки формально равен $\bm{f}(\bm{x})$, корень нелинейной системы алгебраических уравнений совпадает со следующим аргументом минимизации:
\begin{equation}
\argmin_{\bm{x}} g(\bm{x}) = \argmin_{\bm{x}} \langle \bm{f}(\bm{x}), \bm{f}(\bm{x}) \rangle.
\end{equation}
Как мы уже обсуждали в выводе
метода сопряженных градиентов, в подобных задачах минимизации необходимо найти шаг поиска $t$ и направление поиска $\bm{v}$ для формирования итерации вида:
\begin{equation}
\bm{x}^{(k)} = \bm{x}^{(k-1)} + t^{(k)} \bm{v}^{(k)}.
\end{equation}
Оптимальным направлением поиска является направление, обратное вектору градиента функции $g(\bm{x})$, также называемое направлением наискорейшего спуска:
\begin{equation}
\begin{split}
\bm{v}^{(opt)} &= -\nabla g(\bm{x}) \\
&=
-\left[ \frac{\partial g}{\partial x_1}, \frac{\partial g}{\partial x_2}, \dots, \frac{\partial g}{\partial x_n} \right]^T \\
&=
-\left[ \sum_{i = 1}^n 2 f_i(\bm{x}) \frac{\partial f_i(\bm{x})}{\partial x_1},
\sum_{i = 1}^n 2 f_i(\bm{x}) \frac{\partial f_i(\bm{x})}{\partial x_2},
\dots,
\sum_{i = 1}^n 2 f_i(\bm{x}) \frac{\partial f_i(\bm{x})}{\partial x_n}
\right]^T \\
&= -2 \bm{J}^T(\bm{x}) \bm{f}(\bm{x}).
\end{split}
\end{equation}
Выбрав в качестве направления поиска направление наискорейшего спуска, мы получаем
метод градиентного спуска:
\begin{align}
\bm{z}^{(k)} = \bm{J}^T(\bm{x}^{(k-1)}) \bm{f}(\bm{x}^{(k-1)}),
\bm{x}^{(k)} = \bm{x}^{(k-1)} - t^{(k)} \frac{\bm{z}^{(k)}}{||\bm{z}^{(k)}||_2}.
\end{align}
Нахождение оптимального значения для $t^{(k)}$ не всегда является тривиальной задачей, так как это требует дифференцирования функции $h(t) = g(\bm{x}^{(k-1)} - t \frac{\bm{z}^{(k)}}{||\bm{z}^{(k)}||_2})$ относительно $t$. Для упрощения задачи, мы саппроксимируем функцию $h(t)$ квадратичным полиномом $P_2(t)$ и найдем ее минимум. Для этого нам необходимо найти три значения $t_1^{(k)}, t_2^{(k)}, t_3^{(k)}$. Мы сразу можем выбрать $t_1^{(k)} = 0$, так как оно является минимальным значением $t$ и гарантированно не является оптимальным. Затем мы выбираем такое $t_3^{(k)}$, что $h(t_1^{(k)}) > h(t_3^{(k)})$, что можно сделать, начав подбирать с $t_3^{(k)} = 1$, и после делить $t_3^{(k)}$ пополам или удваивать до тех пор, пока не выполнится указанное условие. Когда $t_3^{(k)}$, мы автоматически можем найти $t_2^{(k)} = \frac{t_3^{(k)}}{2}$. Тогда, в соответствии с интерполяцией Лагранжа, мы получаем следующее выражение для $P_2(t)$:
\begin{multline}
P_2(t) = h(t_1^{(k)}) \frac{(t - t_2^{(k)})(t - t_3^{(k)})}{(t_1^{(k)} - t_2^{(k)})(t_1^{(k)} - t_3^{(k)})} + h(t_2^{(k)}) \frac{(t - t_1^{(k)})(t - t_3^{(k)})}{(t_2^{(k)} - t_1^{(k)})(t_2^{(k)} - t_3^{(k)})} + \\ + h(t_3^{(k)}) \frac{(t - t_1^{(k)})(t - t_2^{(k)})}{(t_3^{(k)} - t_1^{(k)})(t_3^{(k)} - t_2^{(k)})}.
\end{multline}
Производная $P_2(t)$ имеет вид:
\begin{multline}
\frac{d P_2}{dt} = h(t_1^{(k)}) \frac{2t - t_2^{(k)} - t_3^{(k)}}{(t_1^{(k)} - t_2^{(k)})(t_1^{(k)} - t_3^{(k)})} + h(t_2^{(k)}) \frac{2t - t_1^{(k)} - t_3^{(k)}}{(t_2^{(k)} - t_1^{(k)})(t_2^{(k)} - t_3^{(k)})} + \\ + h(t_3^{(k)}) \frac{2t - t_1^{(k)} - t_2^{(k)}}{(t_3^{(k)} - t_1^{(k)})(t_3^{(k)} - t_2^{(k)})}.
\end{multline}
Для упрощения записи введем следующие обозначения:
\begin{align}
a^{(k)} &= \frac{h(t_1^{(k)})}{(t_1^{(k)} - t_2^{(k)})(t_1^{(k)} - t_3^{(k)})}, \\
b^{(k)} &= \frac{h(t_2^{(k)})}{(t_2^{(k)} - t_1^{(k)})(t_2^{(k)} - t_3^{(k)})}, \\
c^{(k)} &= \frac{h(t_3^{(k)})}{(t_3^{(k)} - t_1^{(k)})(t_3^{(k)} - t_2^{(k)})}.
\end{align}
Тогда путем обнуления производной мы находим следующее квазиоптимальное значение для $t^{(k)}$:
\begin{equation}
t^{(k)} = \frac{a^{(k)}\left(t_2^{(k)} + t_3^{(k)}\right) + b^{(k)}\left(t_1^{(k)} + t_3^{(k)}\right) + c^{(k)}\left(t_1^{(k)} + t_2^{(k)}\right)}{2\left(a^{(k)} + b^{(k)} + c^{(k)}\right)}.
\end{equation}