Метод внутренней точки

Пример поиска решения. Синие линии показывают ограничения, красные точки — повторяющиеся решения.

Методы внутренней точки (также называемые барьерными методами или IPM ) — это алгоритмы для решения линейных и нелинейных выпуклой оптимизации задач . IPM сочетают в себе два преимущества ранее известных алгоритмов:

• Теоретически время их выполнения полиномиально — в отличие от симплексного метода , который в худшем случае имеет экспоненциальное время выполнения.
• На практике они работают так же быстро, как и симплексный метод — в отличие от метода эллипсоида , который в теории имеет полиномиальное время выполнения, но очень медленный на практике.

В отличие от симплекс-метода, который пересекает границу допустимой области, и метода эллипсоида, который ограничивает допустимую область снаружи , IPM достигает наилучшего решения, пересекая внутреннюю часть допустимой области — отсюда и название.

История

Метод внутренней точки был открыт советским математиком И. И. Дикиным в 1967 году. Этот метод был заново изобретен в США в середине 1980-х годов. В 1984 году Нарендра Кармаркар разработал метод линейного программирования, названный алгоритмом Кармаркара . который работает за доказуемо полиномиальное время ( $O(n^{3.5}L)$ операции над L -битными числами, где n — количество переменных и констант), а также очень эффективны на практике. Статья Кармаркара вызвала всплеск интереса к методам внутренних точек. Два года спустя Джеймс Ренегар изобрел первый метод отслеживания внутренней точки по пути с использованием времени выполнения. $O(n^{3}L)$. Позже метод был расширен от линейных до задач выпуклой оптимизации на основе самосогласованной барьерной функции , используемой для кодирования выпуклого множества .

Любую задачу выпуклой оптимизации можно преобразовать в минимизацию (или максимизацию) линейной функции на выпуклом множестве путем преобразования к форме надграфика . Идея кодирования допустимого множества с использованием барьера и разработки барьерных методов изучалась Энтони В. Фиакко, Гартом П. Маккормиком и другими в начале 1960-х годов. Эти идеи в основном развивались для общего нелинейного программирования , но позже от них отказались из-за наличия более конкурентоспособных методов для этого класса задач (например, последовательного квадратичного программирования ).

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

Класс примитивно-двойственных методов следования по внутренней точке считается наиболее успешным. Алгоритм предиктора-корректора Мехротры обеспечивает основу для большинства реализаций этого класса методов.

Определения

Нам дана выпуклая программа вида:

$${\begin{aligned}{\underset {x\in \mathbb {R} ^{n}}{\text{minimize}}}\quad &f(x)\\{\text{subject to}}\quad &x\in G.\end{aligned}}$$

где f — выпуклая функция , а G — выпуклое множество . Без ограничения общности можно предположить, что цель f является линейной функцией . Обычно выпуклое множество G представляется набором выпуклых неравенств и линейных равенств; линейные равенства можно устранить с помощью линейной алгебры, поэтому для простоты мы предполагаем, что существуют только выпуклые неравенства, а программу можно описать следующим образом, где gi _— выпуклые функции:

$${\begin{aligned}{\underset {x\in \mathbb {R} ^{n}}{\text{minimize}}}\quad &f(x)\\{\text{subject to}}\quad &g_{i}(x)\leq 0{\text{ for }}i=1,\dots ,m.\\\end{aligned}}$$

Мы предполагаем, что функции ограничений принадлежат некоторому семейству (например, квадратичных функций), так что программа может быть представлена ​​конечным вектором коэффициентов (например, коэффициентов квадратичных функций). Размерность этого вектора коэффициентов называется размером программы. Численный решатель для данного семейства программ — это алгоритм, который по вектору коэффициентов генерирует последовательность приближенных решений x _t для t = 1,2,..., используя конечное число арифметических операций. Численный решатель называется сходящимся , если для любой программы из семейства и любого положительного ε > 0 существует некоторое T (которое может зависеть от программы и от ε ) такое, что для любого t > T приближенное решение x _t является ε-приблизительным, то есть:

е ( x_t ) - е ^* ≤ е

г _я ( x_t ) ≤ ε для i в 1,..., m ,

х в G ,

где f ^* является оптимальным решением. Решатель называется полиномиальным , если общее количество арифметических операций за первые T шагов не превосходит

poly(problem-size) * log( V / ε ),

где V — некоторая константа, зависящая от данных, например, разница между наибольшим и наименьшим значением в допустимом наборе. Другими словами, V / ε — это «относительная точность» решения — точность по наибольшему коэффициенту. log( V / ε ) представляет количество «цифр точности». Следовательно, решатель является «полиномиальным», если каждая дополнительная цифра точности требует количества операций, полиномиального по размеру задачи.

Типы

Типы методов внутренней точки включают в себя:

• Возможные методы уменьшения : алгоритм Кармаркара был первым.
• Методы следования по пути : алгоритмы Джеймса Ренегара и Кловис Гонзага были первыми.
• Первично-двойственные методы .

Методы следования по пути

Идея

Учитывая программу выпуклой оптимизации (P) с ограничениями, мы можем преобразовать ее в программу без ограничений , добавив барьерную функцию . В частности, пусть b — гладкая выпуклая функция, определенная внутри допустимой области G , такая, что для любой последовательности { x _j in Interior(G)}, предел которой находится на границе G : $\lim _{j\to \infty }b(x_{j})=\infty$. Мы также предполагаем, что b невырожден, то есть: $b''(x)$ положительно определен для всех x в интерьере (G). Теперь рассмотрим семейство программ:

( P _t ) минимизировать t * f(x) + b(x)

Технически программа ограничена, поскольку определен только внутри G. b Но практически ее можно решить как программу без ограничений, поскольку любой решатель, пытающийся минимизировать функцию, не приблизится к границе, где b стремится к бесконечности. Следовательно, ( P _t ) имеет единственное решение — обозначим его x *( t ). Функция x * является непрерывной функцией t , которая называется центральным путем . Все предельные точки x * при стремлении t к бесконечности являются оптимальными решениями исходной программы (P).

Метод следования по пути — это метод отслеживания функции x * вдоль некоторой возрастающей последовательности t ₁ ,t ₂ ,..., то есть: вычисление достаточно хорошего приближения x _i к точке x *( t _i ), такая, что разность x _i - x *( t _i ) приближается к 0, когда i приближается к бесконечности; тогда последовательность x _i приближается к оптимальному решению (P). Для этого необходимо указать три вещи:

• Барьерная функция b(x).
• Политика определения штрафных параметров t _i .
• Решатель неограниченной оптимизации, используемый для решения ( P _i ) и поиска xi _метод , например, Ньютона . что мы можем использовать каждый xi _{Обратите внимание ,} в качестве отправной точки для решения следующей задачи ( P _i+1 ).

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

Отступник и Гонзага доказал, что конкретный экземпляр метода следования по пути является политаймом:

• Ограничения (и цель) являются линейными функциями;
• Барьерная функция является логарифмической : b(x) := - sum _j log( -g _j ( x )).
• Параметр штрафа t обновляется геометрически, т.е. $t_{i+1}:=\mu \cdot t_{i}$, где µ — константа (они взяли $\mu =1+0.001\cdot {\sqrt {m}}$, где m — количество ограничений-неравенств);
• Решающая программа — это метод Ньютона, и один выполняется для каждого шага t шаг Ньютона .

Они доказали, что в этом случае разница x _i - x *( t _i ) остается не более 0,01, а f( x _i ) - f* составляет не более 2* m / t _i . Таким образом, точность решения пропорциональна 1/ t _i , поэтому для добавления одной цифры точности достаточно умножить ti _. на 2 (или любой другой постоянный коэффициент), что требует O(sqrt( m )) шагов Ньютона . Поскольку каждый шаг Ньютона занимает O( mn ² ) операций, общая сложность O( m ^3/2 н ² ) операции для точности разряда.

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

Подробности

Нам дана задача выпуклой оптимизации (П) в «стандартной форме»:

минимизировать с ^Т x st x в G ,

где G выпуклая и замкнутая. Мы также можем предположить, что G ограничена (мы можем легко сделать ее ограниченной, добавив ограничение | x |≤ R для некоторого достаточно большого R ).

Чтобы использовать метод внутренней точки, нам нужен барьер для G. самосогласованный Пусть b — M -самосогласованный барьер для G , где M ≥1 — параметр самосогласования. что можем эффективно вычислить значение b , его градиент и гессиан для каждой точки x внутри G. Мы предполагаем ,

Для каждого t >0 мы определяем штрафную цель f _t (x) := c ^Т х + б( х ) . Мы определяем путь минимизаторов следующим образом: x*(t) := arg min f _t (x) . Аппроксимируем этот путь вдоль возрастающей последовательности t _i . Последовательность инициализируется определенной нетривиальной двухфазной процедурой инициализации. Затем он обновляется по следующему правилу: $t_{i+1}:=\mu \cdot t_{i}$.

Для каждого t _i находим приблизительный минимум f _ti , обозначаемый xi _мы . Приблизительный минимум выбирается так, чтобы удовлетворять следующему «условию близости» (где L — допуск пути ):

${\sqrt {[\nabla _{x}f_{t}(x_{i})]^{T}[\nabla _{x}^{2}f_{t}(x_{i})]^{-1}[\nabla _{x}f_{t}(x_{i})]}}\leq L$.

Чтобы найти x _{i +1} , мы начинаем с x _i и применяем демпфированный метод Ньютона . Мы применяем несколько шагов этого метода, пока вышеуказанное «отношение близости» не будет удовлетворено. Первая точка, удовлетворяющая этому соотношению, обозначается x _{i +1} .

и

Скорость сходимости метода определяется следующей формулой для каждого i :

$c^{T}x_{i}-c^{*}\leq {\frac {2M}{t_{0}}}\mu ^{-i}$

принимая $\mu =\left(1+r/{\sqrt {M}}\right)$, количество шагов Ньютона, необходимое для перехода от , _{фиксированное} не более _{чем xi к xi +1} число, которое зависит только от r и L . В частности, общее количество шагов Ньютона, необходимое для нахождения ε -приближенного решения (т. е. нахождения x в G такого, что c ^Т x - c* ≤ ε ) не более:

$O(1)\cdot {\sqrt {M}}\cdot \ln \left({\frac {M}{t_{0}\varepsilon }}+1\right)$

где постоянный множитель O(1) зависит только r и L. от Число шагов Ньютона, необходимое для двухэтапной процедуры инициализации, не превышает:

$O(1)\cdot {\sqrt {M}}\cdot \ln \left({\frac {M}{1-\pi _{x_{f}^{*}}({\bar {x}})}}+1\right)+O(1)\cdot {\sqrt {M}}\cdot \ln \left({\frac {M{\text{Var}}_{G}(c)}{\epsilon }}+1\right)$

где постоянный множитель O(1) зависит только от r и L , а ${\text{Var}}_{G}(c):=\max _{x\in G}c^{T}x-\min _{x\in G}c^{T}x$, и ${\bar {x}}$ является некоторой точкой внутри G . В целом, общая ньютоновская сложность поиска ε -приближенного решения не превышает

$O(1)\cdot {\sqrt {M}}\cdot \ln \left({\frac {V}{\varepsilon }}+1\right)$, где V — некоторая константа, зависящая от задачи: $V={\frac {{\text{Var}}_{G}(c)}{1-\pi _{x_{f}^{*}({\bar {x}})}}}$.

Каждый шаг Ньютона занимает O( n ³ ) арифметические операции.

фазы

Чтобы инициализировать методы следования по пути, нам нужна точка внутри допустимой G. области Другими словами: если G определяется неравенствами gi _x ( ) ≤ 0, то нам нужен некоторый x , для которого gi ₍ ) x < 0 для всех i из 1,..., m . Если у нас нет такой точки, нам нужно найти ее, используя так называемый метод фазы I. Простой метод фазы I заключается в решении следующей выпуклой программы:

$${\begin{aligned}{\text{minimize}}\quad &s\\{\text{subject to}}\quad &g_{i}(x)\leq s{\text{ for }}i=1,\dots ,m\end{aligned}}$$

Обозначим оптимальное решение через x*, s *.

• Если s *<0, то мы знаем, что x* является внутренней точкой исходной задачи и можем перейти к «фазе II», которая решает исходную задачу.
• Если s *>0, то мы знаем, что исходная программа неосуществима — допустимая область пуста.
• Если s *=0 и оно достигается некоторым решением x*, то задача разрешима, но не имеет внутренней точки; если оно не достигнуто, то задача невыполнима.

Для этой программы легко получить внутреннюю точку: мы можем произвольно взять x = 0 и принять s за любое число, большее max( f ₁ (0),..., f _m (0)). Поэтому ее можно решить методами внутренних точек. Однако время выполнения пропорционально log(1/ s *). Когда s* приближается к 0, становится все труднее и труднее найти точное решение проблемы фазы I и, следовательно, труднее решить, выполнима ли исходная задача.

соображения

Теоретические гарантии предполагают, что параметр штрафа увеличивается со скоростью $\mu =\left(1+r/{\sqrt {M}}\right)$, поэтому наихудшее количество требуемых шагов Ньютона равно $O({\sqrt {M}})$. Теоретически, если µ больше (например, 2 или более), то наихудшее количество требуемых шагов Ньютона находится в $O(M)$. Однако на практике большее значение ц приводит к гораздо более быстрой сходимости. Эти методы называются методами длинного шага . На практике, если ц находится между 3 и 100, то программа сходится за 20–40 шагов Ньютона, независимо от количества ограничений (хотя время выполнения каждого шага Ньютона, конечно, растет с количеством ограничений). Точное значение μ в этом диапазоне мало влияет на производительность.

потенциала

Для методов потенциальной редукции задача представляется в конической форме :

минимизировать с ^Т x st x в {b+L} ᚢ K ,

где b — вектор из R ^н , L — линейное подпространство в R ^н (поэтому b + L — аффинная плоскость ), а K — замкнутый заостренный выпуклый конус с непустой внутренностью. Любую выпуклую программу можно преобразовать к конической форме. Чтобы использовать метод редукции потенциала (в частности, расширение алгоритма Кармаркара до выпуклого программирования), нам необходимы следующие предположения:

• A. Допустимое множество {b+L} ᚢ K ограничено и пересекает внутреннюю часть конуса K .
• Б. Нам заранее дано строго допустимое решение х т. е. допустимое решение внутри К. ^ ,
• C. Мы заранее знаем оптимальное целевое значение c* задачи.
• D. Дан M -логарифмически-однородный самосогласованный барьер F для конуса K .

Допущения A, B и D необходимы в большинстве методов внутренней точки. Предположение C специфично для подхода Кармаркара; его можно смягчить, используя «скользящее объективное значение». Возможно дальнейшее сведение программы к формату Кармаркара :

минимизировать s ^Т х ст х в М ᚢ К и е ^Т х = 1

где M — линейное подпространство в R ^н , а оптимальное целевое значение равно 0.В основе метода лежит следующая скалярная потенциальная функция:

v ( Икс ) знак равно F ( Икс ) + M пер ( s ^Т х )

где F — M -самосогласованный барьер для допустимого конуса. Можно доказать, что, когда x строго осуществима и v ( x ) очень мало (- очень отрицательно), x приблизительно оптимален. Идея метода уменьшения потенциала состоит в том, чтобы изменить x так, чтобы потенциал на каждой итерации уменьшался как минимум на фиксированную константу X (в частности, X = 1/3-ln(4/3)). Это означает, что после i итераций разница между целевым значением и оптимальным целевым значением составляет не более V * exp(- i X / M ), где V — константа, зависящая от данных. Следовательно, число шагов Ньютона, необходимое для получения ε -приближенного решения, не превосходит $O(1)\cdot M\cdot \ln \left({\frac {V}{\varepsilon }}+1\right)+1$.

Обратите внимание, что в методах следования по пути выражение имеет вид ${\sqrt {M}}$ а не M , что теоретически лучше. Но на практике метод Кармаркара позволяет сделать гораздо большие шаги к цели, поэтому он может сойтись гораздо быстрее, чем теоретические гарантии.

Первично-двойственные методы

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

$${\begin{aligned}\operatorname {minimize} \quad &f(x)\\{\text{subject to}}\quad &x\in \mathbb {R} ^{n},\\&c_{i}(x)\geq 0{\text{ for }}i=1,\ldots ,m,\\{\text{where}}\quad &f:\mathbb {R} ^{n}\to \mathbb {R} ,\ c_{i}:\mathbb {R} ^{n}\to \mathbb {R} .\end{aligned}}\quad (1)$$

Эта задача оптимизации, ограниченная неравенством, решается путем преобразования ее в неограниченную целевую функцию, минимум которой мы надеемся эффективно найти.В частности, логарифмическая барьерная функция , связанная с (1), равна

$$B(x,\mu )=f(x)-\mu \sum _{i=1}^{m}\log(c_{i}(x)).\quad (2)$$

Здесь $\mu$ — небольшая положительная скалярная величина, иногда называемая «параметром барьера». Как $\mu$ сходится к нулю минимум $B(x,\mu )$ должно сходиться к решению (1).

Градиент дифференцируемой функции $h:\mathbb {R} ^{n}\to \mathbb {R}$ обозначается $\nabla h$.Градиент барьерной функции равен

$$\nabla B(x,\mu )=\nabla f(x)-\mu \sum _{i=1}^{m}{\frac {1}{c_{i}(x)}}\nabla c_{i}(x).\quad (3)$$

В дополнение к исходной («основной») переменной $x$ мы вводим множителе Лагранжа , основанную на двойную переменную $\lambda \in \mathbb {R} ^{m}$

$$c_{i}(x)\lambda _{i}=\mu ,\quad \forall i=1,\ldots ,m.\quad (4)$$

Уравнение (4) иногда называют условием «возмущенной дополнительности» из-за его сходства с «дополнительной нежесткостью» в условиях ККТ .

Мы стараемся найти такие $(x_{\mu },\lambda _{\mu })$ для которого градиент барьерной функции равен нулю.

Замена $1/c_{i}(x)=\lambda _{i}/\mu$ из (4) в (3) получаем уравнение для градиента:

$$\nabla B(x_{\mu },\lambda _{\mu })=\nabla f(x_{\mu })-J(x_{\mu })^{T}\lambda _{\mu }=0,\quad (5)$$

где матрица $J$ является якобианом ограничений $c(x)$.

Интуиция, лежащая в основе (5), заключается в том, что градиент $f(x)$ должно лежать в подпространстве, охватываемом градиентами ограничений. «Нарушенная дополнительность» с малыми $\mu$ (4) можно понимать как условие, что решение должно либо лежать вблизи границы $c_{i}(x)=0$, или что проекция градиента $\nabla f$ на компоненте ограничения $c_{i}(x)$ нормальный должен быть почти нулевым.

Позволять $(p_{x},p_{\lambda })$ быть направлением поиска для итеративного обновления $(x,\lambda )$.Применяя метод Ньютона к (4) и (5), получаем уравнение для $(p_{x},p_{\lambda })$:

$${\begin{pmatrix}H(x,\lambda )&-J(x)^{T}\\\operatorname {diag} (\lambda )J(x)&\operatorname {diag} (c(x))\end{pmatrix}}{\begin{pmatrix}p_{x}\\p_{\lambda }\end{pmatrix}}={\begin{pmatrix}-\nabla f(x)+J(x)^{T}\lambda \\\mu 1-\operatorname {diag} (c(x))\lambda \end{pmatrix}},$$

где $H$ представляет собой Гессе матрицу $B(x,\mu )$, $\operatorname {diag} (\lambda )$ представляет собой матрицу диагональную $\lambda$, и $\operatorname {diag} (c(x))$ диагональная матрица $c(x)$.

Ввиду (1), (4) условие

$\lambda \geq 0$

должны соблюдаться на каждом этапе. Это можно сделать, выбрав соответствующий $\alpha$:

$(x,\lambda )\to (x+\alpha p_{x},\lambda +\alpha p_{\lambda }).$

Продолжительность: 11 секунд. 0:11

Траектория итераций x с использованием метода внутренней точки.

Типы выпуклых программ, решаемых

Вот некоторые частные случаи выпуклых программ, которые можно эффективно решить методами внутренних точек.

Линейные программы

Рассмотрим линейную программу вида:

$${\begin{aligned}\operatorname {minimize} \quad &c^{\top }x\\{\text{subject to}}\quad &Ax\leq b.\end{aligned}}.$$

Мы можем применить методы следования по пути с барьером

$$b(x):=-\sum _{j=1}^{m}\ln(b_{j}-a_{j}^{T}x).$$

Функция $b$ согласуется с параметром M = m (количеством ограничений). Следовательно, число необходимых шагов Ньютона для метода следования по пути равно O( mn ² ), а общая сложность времени выполнения равна O( m ^3/2 н ² ).

программы с

Дана квадратичная программа с квадратичными ограничениями вида:

$${\begin{aligned}\operatorname {minimize} \quad &d^{\top }x\\{\text{subject to}}\quad &f_{j}(x):=x^{\top }A_{j}x+b_{j}^{\top }x+c_{j}\leq 0\quad {\text{ for all }}j=1,\dots ,m,\end{aligned}}$$

где все матрицы Aj _{являются} матрицами положительно-полуопределенными .Мы можем применить методы следования по пути с барьером

$$b(x):=-\sum _{j=1}^{m}\ln(-f_{j}(x)).$$

Функция $b$ является самосогласованным барьером с параметром M = m . Сложность Ньютона равна O( (m+n)n ² ), а общая сложность времени выполнения равна O( m ^1/2 (м+п) н ² ).

L _p нормальная аппроксимация

Рассмотрим задачу вида

$${\begin{aligned}\operatorname {minimize} \quad &\sum _{j}|v_{j}-u_{j}^{\top }x|_{p}\end{aligned}},$$

где каждый $u_{j}$ представляет собой вектор, каждый $v_{j}$ является скаляром, и $|\cdot |_{p}$ является Lp _​ нормой с $1<p<\infty .$ После преобразования к стандартной форме мы можем применить методы следования по пути с самосогласованным барьером с параметром M =4 m . Сложность Ньютона равна O( (m+n)n ² ), а общая сложность времени выполнения равна O( m ^1/2 (м+п) н ² ).

Геометрические программы

Рассмотрите проблему

$${\begin{aligned}\operatorname {minimize} \quad &f_{0}(x):=\sum _{i=1}^{k}c_{i0}\exp(a_{i}^{\top }x)\\{\text{subject to}}\quad &f_{j}(x):=\sum _{i=1}^{k}c_{ij}\exp(a_{i}^{\top }x)\leq d_{j}\quad {\text{ for all }}j=1,\dots ,m.\end{aligned}}$$

Существует самосогласованный барьер с параметром 2 k + m . Метод следования по пути имеет сложность Ньютона O( mk ² + к ³ + н ³ ) и общая сложность O(( k+m ) ^1/2 [ МК ² + к ³ + н ³ ]).

Semidefinite programs

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

См. также

• Аффинное масштабирование
• Расширенный метод Лагранжа
• Алгоритм Шамболь-Пока
• Условия Каруша – Куна – Такера
• Метод штрафа

Ссылки

• Боннан, Ж. Фредерик; Гилберт, Дж. Чарльз; Лемарешаль, Клод ; Сагастисабал, Клаудия А. (2006). Численная оптимизация: Теоретические и практические аспекты . Universitext (Второе исправленное издание перевода французского издания 1997 г.). Берлин: Springer-Verlag. стр. xiv+490. дои : 10.1007/978-3-540-35447-5 . ISBN  978-3-540-35445-1 . МР   2265882 .
• Носедаль, Хорхе; Стивен Райт (1999). Численная оптимизация . Нью-Йорк, штат Нью-Йорк: Спрингер. ISBN  978-0-387-98793-4 .
• Пресс, WH; Теукольский, С.А.; Феттерлинг, WT; Фланнери, BP (2007). «Раздел 10.11. Линейное программирование: методы внутренних точек» . Численные рецепты: искусство научных вычислений (3-е изд.). Нью-Йорк: Издательство Кембриджского университета. ISBN  978-0-521-88068-8 .