Документация Engee

NLsolve.jl

Решение нелинейных систем уравнений в Julia.

NLsolve.jl является частью семейства JuliaNLSolvers.

Нелинейные системы уравнений

Пакет NLsolve решает системы нелинейных уравнений. Формально, если F является многозначной функцией, этот пакет ищет некий вектор x, который удовлетворяет F(x)=0 с некоторой точностью.

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

Существует также идентичный API для решения фиксированных точек (принимающий в качестве входных данных функцию F(x) и решающий F(x) = x).

Если требуется решить одно уравнение, а не систему, и производительность не критична, рекомендуется использовать Roots.jl. Для решения небольшой системы уравнений с высокой производительностью используйте NonlinearSolve.jl.

Очень простой пример

Рассмотрим следующую двумерную функцию с двумя переменными.

(x, y) -> [(x+3)*(y^3-7)+18, sin(y*exp(x)-1)]

Чтобы найти нуль этой функции и отобразить его, нужно написать следующую программу.

using NLsolve

function f(x)
    [(x[1]+3)*(x[2]^3-7)+18,
    sin(x[2]*exp(x[1])-1)]
end

sol = nlsolve(f, [ 0.1, 1.2])
sol.zero

Первый аргумент для nlsolve — это решаемая функция, которая принимает вектор в качестве входных данных и возвращает невязку в виде вектора. Второй аргумент является начальным предположением алгоритма. sol.zero извлекает решение, если сходится.

Простой пример

Продолжаем работу с той же системой уравнений, но теперь используем функцию на месте и заданный пользователем якобиан для повышения производительности:

using NLsolve

function f!(F, x)
    F[1] = (x[1]+3)*(x[2]^3-7)+18
    F[2] = sin(x[2]*exp(x[1])-1)
end

function j!(J, x)
    J[1, 1] = x[2]^3-7
    J[1, 2] = 3*x[2]^2*(x[1]+3)
    u = exp(x[1])*cos(x[2]*exp(x[1])-1)
    J[2, 1] = x[2]*u
    J[2, 2] = u
end

nlsolve(f!, j!, [ 0.1; 1.2])

Во-первых, обратите внимание, что функция f! вычисляет невязки нелинейной системы и сохраняет их в предварительно выделенном векторе, переданном в качестве первого аргумента. Аналогичным образом функция j! вычисляет якобиан системы и сохраняет его в предварительно выделенной матрице, переданной в качестве первого аргумента. Невязки и якобианы могут иметь различную форму (см. ниже).

Во-вторых, функция nlsolve возвращает объект типа SolverResults. В частности, поле zero этой структуры содержит решение, если произошла сходимость. Если r является объектом типа SolverResults, converged(r) указывает, произошла ли сходимость.

Способы задания функции и ее якобиана

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

С использованием функций, изменяющих аргументы на месте

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

Далее предполагается, что вы определили функцию f!(F::AbstractVector, x::AbstractVector) или, в более общем случае, функцию f!(F::AbstractArray, x::AbstractArray), вычисляющую невязку системы в точке x и помещающую ее в аргумент F.

В свою очередь, существует три способа указания варианта вычисления якобиана:

Вычисление конечных разностей

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

nlsolve(f!, initial_x)

В качестве альтернативы можно создать объект типа OnceDifferentiable и передать его функции nlsolve, как в примере:

initial_x = ...
initial_F = similar(initial_x)
df = OnceDifferentiable(f!, initial_x, initial_F)
nlsolve(df, initial_x)

Обратите внимание, что мы передали initial_x и initial_F конструктору для df. Это не обязательно должен быть реальный начальный x и вектор невязки в x, но он используется для инициализации переменных кэша в df, поэтому их типы и измерения должны быть такими, как если бы они существовали.

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

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

nlsolve(f!, initial_x, autodiff = :forward)

Доступный якобиан

Если дополнительно к f!(F::AbstractArray, x::AbstractArray) у вас есть функция j!(J::AbstractArray, x::AbstractArray) для вычисления якобиана системы, то синтаксис будет таким, как в примере выше:

nlsolve(f!, j!, initial_x)

Опять же, можно указать две функции f!(F::AbstractArray, x::AbstractArray) и j!(J::AbstractArray, x::AbstractArray), которые работают с произвольными массивами x.

Обратите внимание, что не следует считать, что якобиан J, переданный в j!, инициализирован нулевой матрицей. Необходимо задать все элементы матрицы в функции j!.

В качестве альтернативы можно создать объект типа OnceDifferentiable и передать его функции nlsolve, как в примере:

df = OnceDifferentiable(f!, j!, initial_x, initial_F)
nlsolve(df, initial_x)

Оптимизация одновременных невязок и якобиана

Если помимо f! и j! у вас есть функция fj!(F::AbstractArray, J::AbstractArray, x::AbstractArray), которая одновременно вычисляет невязку и якобиан, вы можете использовать следующий синтаксис.

df = OnceDifferentiable(f!, j!, fj!, initial_x, initial_F)
nlsolve(df, initial_x)

Если функция fj! использует некоторую оптимизацию, благодаря которой она обходится дешевле, чем последовательный вызов f! и j!, этот синтаксис, возможно, улучшит производительность.

Указание только функции fj!

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

function myfun!(F, J, x)
    ## начало общих вычислений
    ## ...
    ## конец общих вычислений
    if !(F == nothing)
        ## выполняются изменяемые вычисления, характерные для f!
    end
    if !(J == nothing)
        ## выполняются изменяемые вычисления, характерные для j!
    end
end

и выполните решение:

nlsolve(only_fj!(myfun), initial_x)

Это позволит функции nlsolve вычислять F(x) и J(x) вместе, но при этом эффективно вычислять F(x) или J(x) по отдельности.

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

Здесь предполагается, что у вас есть функция f(x::AbstractArray), которая возвращает только что выделенный вектор, содержащий невязки. Просто передайте его функции nlsolve, и она автоматически определит, определен ли f для одного или двух аргументов:

nlsolve(f, initial_x)

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

С помощью ключевого слова autodiff для вычисления якобиана в этом случае можно использовать как нахождение конечных разностей, так и автоматическое дифференцирование.

Если наряду с f(x::AbstractArray) имеется функция j(x::AbstractArray), возвращающая только что выделенную матрицу, содержащую якобиан, мы снова просто передаем их функции nlsolve:

nlsolve(f, j, initial_x)

Если наряду с f и j имеется функция fj, возвращающая кортеж только что выделенного вектора невязок и только что выделенную матрицу, содержащую якобиан, подход останется прежним:

nlsolve(f, j, fj, initial_x)

С использованием функций, принимающих несколько скалярных аргументов

Если у вас есть функция f(x::Float64, y::Float64, ...), которая принимает нужную точку в виде нескольких скаляров и возвращает вектор или кортеж, содержащий невязки, можно использовать вспомогательную функцию n_ary. Полный синтаксис выглядит следующим образом.

nlsolve(n_ary(f), initial_x)

Для вычисления якобиана используется метод нахождения конечных разностей.

Если якобиан разрежен

Если якобиан функции разрежен, можно настроить подпрограммы для работы с разреженными матрицами, а не полными, чтобы увеличить производительность в больших системах. Это означает, что необходимо обязательно указать соответствующий тип якобиана, чтобы решатель знал, что передавать функции j!.

df = OnceDifferentiable(f!, j!, x0, F0, J0)
nlsolve(df, initial_x)

Можно передать конструктору необязательную третью функцию fj!, как для случая с полным якобианом.

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

fill!(a, 0)
dropzeros!(a) ## если вы также хотите удалить шаблон разреженности

Тонкие настройки

На данный момент доступно три алгоритма. Выбрать нужный можно путем задания необязательного аргумента method функции nlsolve. Алгоритмом по умолчанию является метод доверительной области.

Метод доверительной области

Это хорошо известный метод решения, основанный на квадратичной аппроксимации задачи наименьших квадратов, которая считается действительной в компактной области с центром вокруг текущей итерации.

Этот метод выбирается следующим образом: method = :trust_region.

Метод принимает следующие пользовательские параметры.

  • factor: определяет размер начальной доверительной области. Этот размер устанавливается как произведение коэффициента и евклидовой нормы initial_x, если она ненулевая, или в противном случае как сам коэффициент. Значение по умолчанию: 1.0.

  • autoscale: если true, переменные будут автоматически перемасштабированы. Коэффициенты масштабирования представляют собой нормы столбцов якобиана. Значение по умолчанию: true.

Метод Ньютона с линейным поиском

Это классический алгоритм Ньютона с дополнительным линейным поиском.

Этот метод выбирается следующим образом: method = :newton.

Он принимает пользовательский параметр linesearch, который должен быть равен функции, вычисляющей линейный поиск. В настоящее время источником доступных значений является пакет LineSearches. По умолчанию линейный поиск не выполняется. Примечание. Предполагается, что переданная функция линейного поиска по крайней мере обновит вектор решения и вычислит функцию в новой точке.

Ускорение Андерсона

Этот метод выбирается следующим образом: method = :anderson.

Он также известен как DIIS или смешение Пулая и основан на ускорении итерации с фиксированной точкой xₙ₊₁ = xₙ + beta*f(xₙ), где по умолчанию beta=1. Он не использует информацию якобиана или линейный поиск, но имеет историю, размер которой контролируется параметром m: m=0 соответствует простой итерации с фиксированной точкой, а более высокие значения используют больший размер истории для ускорения итераций. Большие значения m обычно увеличивают скорость сходимости, но повышают требования к хранению и вычислениям и могут привести к нестабильности. Этот метод полезен для ускорения итерации с фиксированной точкой xₙ₊₁ = g(xₙ) (в этом случае используйте этот решатель при f(x) = g(x) - x).

Справочные материалы: H. Walker, P. Ni, Anderson acceleration for fixed-point iterations, SIAM Journal on Numerical Analysis, 2011

Общие параметры

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

  • xtol: разность норм x между двумя последовательными итерациями, при которой объявляется сходимость. Значение по умолчанию: 0.0.

  • ftol: бесконечная норма невязок, при которой объявляется сходимость. Значение по умолчанию: 1e-8.

  • iterations: максимальное количество итераций. Значение по умолчанию: 1_000.

  • store_trace: следует ли хранить трассировку состояния алгоритма

    1. Значение по умолчанию: false.

  • show_trace: следует ли отображать трассировку состояния алгоритма оптимизации в STDOUT. Значение по умолчанию: false.

  • extended_trace: следует ли добавлять дополнительные внутренние компоненты алгоритма в трассировку состояния. Значение по умолчанию: false.

Фиксированные точки

Существует оболочка fixedpoint() вокруг nlsolve(), которая сопоставляет входную функцию F(x) с G(x) = F(x) - x, и то же самое для версии на месте. Это позволяет решать задачи с фиксированной точкой, например такие, которые часто встречаются в вычислительной экономике. Далее приведены некоторые примечания.

  • По умолчанию используется метод :anderson с m = 5. Собственная итерация в стиле Пикара может быть выполнена, если задать m=0, но это не рекомендуется делать для сжатий, константы Липшица которых близки к 1. Однако если сходимость завершается сбоем, можно подумать о ее снижении.

  • Поддерживается автоматическое дифференцирование, например fixedpoint(f!, init_x; method = :newton, autodiff = :forward).

  • Допуски и границы итераций могут быть заданы точно так же, как в nlsolve(), поскольку эта функция является оболочкой, например fixedpoint(f, init_x; iterations = 500, ...).

Примечание. Если вы предоставляете собственную производную, убедитесь, что она соответствующим образом преобразована (например, в настоящее время мы сопоставляем f -> f - x, но ждем стабилизации API, прежде чем сопоставлять J -> J - I, поэтому вам придется сделать это самостоятельно).

Смешанные задачи комплементарности

Если задана многомерная функция f и два вектора a и b, решением смешанной задачи дополнительности (MCP) является вектор x такой, что для каждого индекса i выполняется одно из следующих условий:

  • либо f_i(x) = 0 и a_i <= x_i <= b_i

  • или f_i(x) > 0 и x_i = a_i

  • или f_i(x) < 0 и x_i = b_i

Вектор a может содержать элементы, равные -Inf, а вектор b может содержать элементы, равные Inf. В частном случае, когда все элементы a равны -Inf, а все элементы b равны Inf, задача MCP в точности эквивалентна многомерной задаче нахождения корней, описанной выше.

Пакет решает задачи MCP, переформулируя их как решение системы нелинейных уравнений (как описано в статье Миранды (Miranda) и Факлера (Fackler) от 2002 г., хотя NLsolve использует соглашение о знаках, противоположное их).

Функция mcpsolve решает задачи MCP. Она принимает те же аргументы, что и nlsolve, за исключением того, что векторы a и b должны сразу же следовать за аргументами, соответствующими f (и, возможно, ее производной). Также есть дополнительный необязательный аргумент reformulation, который может принимать два значения.

  • reformulation = :smooth: использовать корректную переформулировку задачи с помощью функции Фишера. Это значение используется по умолчанию, так как оно более надежно для сложных задач.

  • reformulation = :minmax: использовать минимально-максимальную переформулировку задачи. Это быстрее, чем сглаженное приближение, поскольку используется меньше алгебры, но менее надежно, поскольку переформулированная задача имеет перегибы.

Вот полный пример.

using NLsolve

function f!(F, x)
    F[1]=3*x[1]^2+2*x[1]*x[2]+2*x[2]^2+x[3]+3*x[4]-6
    F[2]=2*x[1]^2+x[1]+x[2]^2+3*x[3]+2*x[4]-2
    F[3]=3*x[1]^2+x[1]*x[2]+2*x[2]^2+2*x[3]+3*x[4]-1
    F[4]=x[1]^2+3*x[2]^2+2*x[3]+3*x[4]-3
end

r = mcpsolve(f!, [0., 0., 0., 0.], [Inf, Inf, Inf, Inf],
             [1.25, 0., 0., 0.5], reformulation = :smooth, autodiff = :forward)

Вот решение:

julia> r.zero
4-element Array{Float64,1}:
  1.22474
  0.0
 -1.378e-19
  0.5

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

julia> F = similar(r.zero)

julia> f!(F, r.zero)

julia> F
4-element Array{Float64,1}:
 -1.26298e-9
  3.22474
  5.0
  3.62723e-11

Список задач

Связанные пакеты

  • JuMP.jl также может решать нелинейные уравнения. Просто переформулируйте вашу задачу как задачу оптимизации с нелинейными ограничениями: используйте набор уравнений в качестве ограничений и введите 1.0 в качестве целевой функции. В настоящее время JuMP поддерживает ряд решателей с открытым исходным кодом и коммерческих решателей.

  • Complementarity.jl использует мощный язык моделирования JuMP.jl для решения задач комплементарности. Он поддерживает два решателя: PATHSolver.jl и NLsolve.jl.

Справочные материалы

Nocedal, Jorge and Wright, Stephen J. (2006): Numerical Optimization, second edition, Springer

MINPACK от Jorge More', Burt Garbow, and Ken Hillstrom at Argonne National Laboratory

Miranda, Mario J. and Fackler, Paul L. (2002): Applied Computational Economics and Finance, MIT Press