7.2. Метод Бройдена#
Метод Ньютона послужил основой для многих других методов. Недостатком метода является требовательное по времени вычисление якобиана на каждой итерации. Для некоторых приложений это неприемлимо. Вместо прямого вычисления якобиана, квазиньютоновские методы заменяют якобиан некоторым приближением.
7.2.1. Приближение якобиана конечной разностью#
В 6.2.2 Методе секущих мы фактически сталкивались с приближением производной функции в точке
Численный метод, приближающий производную «из определения» (
Подобным образом можно поступить для приближения якобиана. В общем случае для
где
Якобиан состоит из частных производных, таким образом, можно выписать его приближение через конечные разности
Значение
Приближение якобиана конечной разностью
"""
jacobianfd(f, x[; y, δ])
Вычисляет якобиан функции `f` в точке `x` через конечную разность в точках `x` и `x + δ
I[:, j]`, где `δ::Number` - скаляр. Опционально можно подать `y == f(x)`.
"""
function jacobianfd(f, x; y=f(x), δ=sqrt(eps())*max(norm(x), 1))
m, n = size(y, 1), size(x, 1)
J = zeros(m, n)
x = float(copy(x))
for j in 1:n
x[j] += δ
J[:, j] .= (f(x) .- y) ./ δ
x[j] -= δ
end
return J
end
Сравним на небольшом примере точность приближения якобиана. В качестве функции и точки приближения возьмём
function f(x)
x₁, x₂, x₃ = x
return [
x₂*exp(x₁) + sin(x₁),
cos(x₂) + x₁*x₃,
log(x₃) + x₂^2,
]
end
function J(x)
x₁, x₂, x₃ = x
return [
x₂*exp(x₁)+cos(x₁) exp(x₁) 0 ;
x₃ -sin(x₂) x₁ ;
0 2x₂ 1/x₃ ;
]
end
Вычислим якобиан напрямую Jexact
, и через конечные разности Jfd
x = [2/3, -1.0, sqrt(2)]
Jexact, Jfd = J(x), jacobianfd(f, x)
Jexact |> display
Jfd |> display
@show norm(Jexact - Jfd);
3×3 Matrix{Float64}:
-1.16185 1.94773 0.0
1.41421 0.841471 0.666667
0.0 -2.0 0.707107
3×3 Matrix{Float64}:
-1.16185 1.94773 0.0
1.41421 0.841471 0.666667
0.0 -2.0 0.707107
norm(Jexact - Jfd) = 4.4920040357185245e-8
Отличие в данной точке получилось небольшим. В качестве меры невязки здесь приведена норма Фробениуса для матриц
7.2.2. Метод Бройдена#
Обозначим приближение якобиана в точке
На
здесь обозначены
Квазиньютоновское условие на новое приближение якобиана
Выбор
такой, что условие (7.5) удовлетворится.
Самыми ресурсоемкими операциями в правиле Бройдена является матричное умножение (
Пусть дана функция
В качестве начального приближения якобиана можно использовать как его конечную разность, так и настоящий якобиан.
7.2.3. Реализация#
Метод Бройдена
"""
broydensys(f, x, J[; maxiter, xtol, ftol])
Решает нелинейную систему уравнений `f`(x) = 0 методом Бройдена.
Требует начального приближения корня `x` уравнения и якобиана `J` в этой точке.
Выполняет итерации, пока норма решения `> xtol` или норма функции `> ftol`.
В случае превышения числа итераций `maxiter` вызывает ошибку.
"""
function broydensys(f, x, J; maxiter=50, xtol=1e-6, ftol=1e-6)
δx = float(similar(x))
yp, yn = similar.((δx, δx))
x = float(copy(x))
B = float(copy(J))
yn .= f(x)
for i in 1:maxiter
yp .= yn
δx .= .- (B \ yp)
x .+= δx
yn .= f(x)
norm(δx) < xtol && return x
norm(yn) < ftol && return x
g = B * δx
B .+= (1 / dot(δx, δx)) .* (yn .- yp .- g) .* δx'
end
error("Превышено число итераций.")
end
В данной реализации для простоты изложения явно не использовано QR-разложение, имеющее преимещуство перед LU-разложением при обновлении B
.
Рассмотрим пример использования метода Бройдена. Возьмём функцию и якобиан из Демонстрации 7.3, где использовался метод Ньютона
function f(x)
x₁, x₂ = x
return [
x₁^2 - 2x₂^2 - x₁*x₂ + 2x₁ - x₂ + 1,
2x₁^2 - x₂^2 + x₁*x₂ + 3x₂ - 5,
]
end
function J(x)
x₁, x₂ = x
return [
2x₁-x₂+2 -4x₂-x₁-1;
4x₁+x₂ -2x₂+x₁+3
]
end
Сравним решения методом Бройдена, где в качестве
x = [10.0, 10.0]
rjac = broydensys(f, x, J(x))
rfdjac = broydensys(f, x, jacobianfd(f, x))
@show rjac
@show rfdjac
@show f(rjac);
rjac = [0.9999999480032711, 1.0000000136175216]
rfdjac = [0.999999948003263, 1.0000000136175238]
f(rjac) = [-2.376953129878956e-7, -2.327485972841714e-7]
Поищем также другие корни, не используя подлинный якобиан
x = [-10.0, 10.0]
root = broydensys(f, x, jacobianfd(f, x))
root, f(root)
([-1.4999998826044518, 0.5000004907722365], [-9.122522022231294e-7, -4.0028955261561805e-7])
x = [10.0, -10.0]
root = broydensys(f, x, jacobianfd(f, x))
root, f(root)
([-1.499999989517552, 0.5000000339688456], [-6.667694307793681e-8, -4.066904146782235e-8])
x = [-10.0, -10.0]
root = broydensys(f, x, jacobianfd(f, x))
root, f(root)
([-1.666666678734911, -0.3333330286952676], [6.213441938740516e-7, 6.937537442297526e-7])
7.2.4. Другие версии метода#
Изложенная выше версия метода (7.7) требует решения линейной системы (7.4). Однако, она сочетается с линейным поиском, обеспечивающим глобальную сходимость метода, т.е. сходимостью к корню из произвольного начального приближения
Если же линейный поиск не используется, распространены две версии алгоритма, использующие вместо матрицы
Самим Бройденом была предложена формула для
здесь
Вторая версия метода на основе обратной матрицы
которая минимизирует норму Фробениуса уже между обратными приближениями якобиана
Эта версия распространилась под названием «плохого метода Бройдена» (bad Broyden’s method). Хотя, слова «плохой» и «хороший» в названии методов имеют нейтральную коннотацию.