Interior-point method: Difference between revisions

Revision as of 06:19, 21 December 2023

Latest revision as of 00:20, 20 June 2025

(25 intermediate revisions by 14 users not shown)

Line 9:

== History ==

An interior point method was discovered by Soviet mathematician I. I. Dikin in 1967.<ref>{{Cite journal |last1=Dikin |first1=I.I. |year=1967 |title=Iterative solution of problems of linear and quadratic programming. |url=https://zbmath.org/?q=an:0189.19504 |journal=Dokl. Akad. Nauk SSSR |volume=174 |issue=1 |pages=747–748|zbl=0189.19504 }}</ref> The method was reinvented in the U.S. in the mid-1980s. In 1984, [[Narendra Karmarkar]] developed a method for [[linear programming]] called [[Karmarkar's algorithm]],<ref>{{cite conference |last1=Karmarkar |first1=N. |year=1984 |title=Proceedings of the sixteenth annual ACM symposium on Theory of computing – STOC '84 |pages=302 |doi=10.1145/800057.808695 |isbn=0-89791-133-4 |archive-url=https://web.archive.org/web/20131228145520/http://retis.sssup.it/~bini/teaching/optim2010/karmarkar.pdf |archive-date=28 December 2013 |doi-access=free |chapter-url=http://retis.sssup.it/~bini/teaching/optim2010/karmarkar.pdf |chapter=A new polynomial-time algorithm for linear programming |url-status=dead}}</ref> which runs in provablyprobably polynomial time (<math>O(n^{3.5} L)</math> operations on ''L''-bit numbers, where ''n'' is the number of variables and constants), and is also very efficient in practice. Karmarkar's paper created a surge of interest in interior point methods. Two years later, [[James Renegar]] invented the first ''path-following'' interior-point method, with run-time <math>O(n^{3} L)</math>. The method was later extended from linear to convex optimization problems, based on a [[self-concordant]] [[barrier function]] used to encode the [[convex set]].<ref name=":0">{{Cite book |last=Arkadi Nemirovsky |url=https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=8c3cb6395a35cb504019f87f447d65cb6cf1cdf0 |title=Interior point polynomial-time methods in convex programming |year=2004}}</ref>

Any convex optimization problem can be transformed into minimizing (or maximizing) a [[linear function]] over a convex set by converting to the [[Epigraph (mathematics)|epigraph]] form.<ref name=":3">{{cite book |lastlast1=Boyd |firstfirst1=Stephen |title=Convex Optimization |last2=Vandenberghe |first2=Lieven |publisher=[[Cambridge University Press]] |year=2004 |isbn=978-0-521-83378-3 |___location=Cambridge |pages= |mr=2061575}}</ref>{{Rp|___location=143}} The idea of encoding the [[candidate solution|feasible set]] using a barrier and designing barrier methods was studied by Anthony V. Fiacco, Garth P. McCormick, and others in the early 1960s. These ideas were mainly developed for general [[nonlinear programming]], but they were later abandoned due to the presence of more competitive methods for this class of problems (e.g. [[sequential quadratic programming]]).

[[Yurii Nesterov]] and [[Arkadi Nemirovski]] came up with a special class of such barriers that can be used to encode any convex set. They guarantee that the number of [[iteration]]s of the algorithm is bounded by a polynomial in the dimension and accuracy of the solution.<ref>{{Cite journal |mr=2115066 |doi=10.1090/S0273-0979-04-01040-7 |title=The interior-point revolution in optimization: History, recent developments, and lasting consequences |year=2004 |last1=Wright |first1=Margaret H. |journal=Bulletin of the American Mathematical Society |volume=42 |pages=39–57|doi-access=free }}</ref><ref name=":0" />

The class of primal-dual path-following interior-point methods is considered the most successful. [[Mehrotra predictor–corrector method|Mehrotra's predictor–corrector algorithm]] provides the basis for most implementations of this class of methods.<ref>{{cite journal |last=Potra |first=Florian A. |author2=Stephen J. Wright |title=Interior-point methods |journal=Journal of Computational and Applied Mathematics |volume=124 |year=2000 |issue=1–2 |pages=281–302 |doi=10.1016/S0377-0427(00)00433-7|doi-access=free |bibcode=2000JCoAM.124..281P }}</ref>

== Definitions ==

Line 28:

\text{subject to}\quad & g_i(x) \leq 0 \text{ for } i = 1, \dots, m. \\

\end{aligned}

</math>We assume that the constraint functions belong to some family (e.g. quadratic functions), so that the program can be represented by a finite ''vector of coefficients'' (e.g. the coefficients to the quadratic functions). The dimension of this coefficient vector is called the ''size'' of the program. A ''numerical solver'' for a given family of programs is an algorithm that, given the coefficient vector, generates a sequence of approximate solutions ''x_t'' for ''t''=1,2,..., using finitely many arithmetic operations. A numerical solver is called ''convergent'' if, for any progarmprogram from the family and any positive ''ε''>0, there is some ''T'' (which may depend on the program and on ''ε'') such that, for any ''t''>''T'', the approximate solution ''x_t'' is ''ε-approximate,'' that is: <blockquote>''f''(''x'')math - f<supdisplay="block">*</sup> ≤ ''ε''

\begin{aligned}

& f(x_{t}) - f^{*} \leq \epsilon, \\

''g_i''(''x'') ≤ ''ε'' for ''i'' in 1,...,''m'',

& g_{i}(x_{t}) \leq \epsilon \quad \text{for} \quad i = 1, \dots, m, \\

& x \in G,

''x'' in ''G'',</blockquote>where f^* is the optimal solution. A solver is called ''polynomial'' if the total number of arithmetic operations in the first ''T'' steps is at most<blockquote>poly(problem-size) * log(''V''/''ε''),</blockquote>where ''V'' is some data-dependent constant, e.g., the difference between the largest and smallest value in the feasible set. In other words, ''V''/''ε'' is the "relative accuracy" of the solution - the accuracy w.r.t. the largest coefficient. log(''V''/''ε'') represents the number of "accuracy digits". Therefore, a solver is 'polynomial' if each additional digit of accuracy requires a number of operations that is polynomial in the problem size.

\end{aligned}

''x'' in ''G'',</blockquotemath>where f<supmath>f^{*}</supmath> is the optimal solution. A solver is called ''polynomial'' if the total number of arithmetic operations in the first ''T'' steps is at most<blockquote>poly(problem-size) * log(''V''/''ε''),</blockquote>where ''V'' is some data-dependent constant, e.g., the difference between the largest and smallest value in the feasible set. In other words, ''V''/''ε'' is the "relative accuracy" of the solution - the accuracy w.r.t. the largest coefficient. log(''V''/''ε'') represents the number of "accuracy digits". Therefore, a solver is 'polynomial' if each additional digit of accuracy requires a number of operations that is polynomial in the problem size.

== Types ==

Line 38 ⟶ 40:

* '''Potential reduction methods''': [[Karmarkar algorithm|Karmarkar's algorithm]] was the first one.

* '''Path-following methods''': the algorithms of [[James Renegar]]<ref name=":1">{{Cite journal |last=Renegar |first=James |date=1988-01-01 |title=A polynomial-time algorithm, based on Newton's method, for linear programming |url=https://doi.org/10.1007/BF01580724 |journal=Mathematical Programming |language=en |volume=40 |issue=1 |pages=59–93 |doi=10.1007/BF01580724 |issn=1436-4646|url-access=subscription }}</ref> and Clovis Gonzaga<ref name=":2">{{Citation |last=Gonzaga |first=Clovis C. |title=An Algorithm for Solving Linear Programming Problems in O(n3L) Operations |date=1989 |url=https://doi.org/10.1007/978-1-4613-9617-8_1 |work=Progress in Mathematical Programming: Interior-Point and Related Methods |pages=1–28 |editor-last=Megiddo |editor-first=Nimrod |access-date=2023-11-22 |place=New York, NY |publisher=Springer |language=en |doi=10.1007/978-1-4613-9617-8_1 |isbn=978-1-4613-9617-8|url-access=subscription }}</ref> were the first ones.

* '''Primal-dual methods'''.

Line 50 ⟶ 52:

* The barrier function b(x).

* A policy for determining the penalty parameters ''t_i''.

* The unconstrained-optimization solver used to solve (''P_i'') and find ''x_i'', such as [[Newton's method]]. Note that we can use each ''x_i'' as a starting-point for solving the next problem (''P_i+1'').

The main challenge in proving that the method is polytime is that, as the penalty parameter grows, the solution gets near the boundary, and the function becomes steeper. The run-time of solvers such as [[Newton's method]] becomes longer, and it is hard to prove that the total runtime is polynomial.

Line 58 ⟶ 60:

* The constraints (and the objective) are linear functions;

* The barrier function is [[Logarithmic barrier function|logarithmic]]: b(x) := - sum''_j'' log(''-g_j''(''x'')).

* The penalty parameter ''t'' is updated geometrically, that is, <math>t_{i+1} := \mu \cdot t_i</math>, where ''μ'' is a constant (they took <math>\mu = 1+0.001\cdot \sqrt{m}</math>, where ''m'' is the number of inequality constraints);

* The solver is Newton's method, and a ''single'' step of Newton is done for each single step in ''t''.

They proved that, in this case, the difference ''x_i'' - ''x''*(''t_i'') remains at most 0.01, and f(''x_i'') - f* is at most 2*''m''/''t_i''. Thus, the solution accuracy is proportional to 1/''t_i'', so to add a single accuracy-digit, it is suffiicentsufficient to multiply ''t_i'' by 2 (or any other constant factor), which requires O(sqrt(''m'')) Newton steps. Since each Newton step takes O(''m n''²) operations, the total complexity is O(''m^3/2 n''²) operations for accuracy digit.

[[Yurii Nesterov|Yuri Nesterov]] extended the idea from linear to non-linear programs. He noted that the main property of the logarithmic barrier, used in the above proofs, is that it is [[self-concordant]] with a finite barrier parameter. Therefore, many other classes of convex programs can be solved in polytime using a path-following method, if we can find a suitable self-concordant barrier function for their feasible region.<ref name=":0" />{{Rp|___location=Sec.1}}

=== Details ===

We are given a convex optimization problem (P) in "standard form":<blockquote>'''minimize ''c''^T''x'' s.t. ''x'' in ''G''''', </blockquote>where ''G'' is convex and closed. We can also assume that ''G'' is bounded (we can easily make it bounded by adding a constraint |''x''|≤''R'' for some sufficiently large ''R'').<ref name=":0" />{{Rp|___location=Sec.4}}

To use the interior-point method, we need a [[self-concordant barrier]] for ''G''. Let ''b'' be an ''M''-self-concordant barrier for ''G'', where ''M''≥1 is the self-concordance parameter. We assume that we can compute efficiently the value of ''b'', its gradient, and its [[Hessian matrix|Hessian]], for every point x in the interior of ''G''.

For every ''t''>0, we define the ''penalized objective'' '''f_t(x) := t''c''^T''x +'' b(''x'')'''. We define the path of minimizers by: '''x*(t) := arg min f_t(x)'''. We apporimateapproximate this path along an increasing sequence ''t_i''. The sequence is initialized by a certain non-trivial two-phase initialization procedure. Then, it is updated according to the following rule: <math>t_{i+1} := \mu \cdot t_i</math>.

For each ''t_i'', we find an approximate minimum of ''f_ti'', denoted by ''x_i''. The approximate minimum is chosen to satisfy the following "closeness condition" (where ''L'' is the ''path tolerance''):<blockquote><math>\sqrt{[\nabla_x f_t(x_i)]^T [\nabla_x^2 f_t(x_i)]^{-1} [\nabla_x f_t(x_i)]} \leq L</math>.</blockquote>To find ''x_i''₊₁, we start with ''x_i'' and apply the [[damped Newton method]]. We apply several steps of this method, until the above "closeness relation" is satisfied. The first point that satisfies this relation is denoted by ''x_i''₊₁.<ref name=":0" />{{Rp|___location=Sec.4}}

Line 98 ⟶ 100:

=== Practical considerations ===

The theoretic guarantees assume that the penalty parameter is increased at the rate <math>\mu = \left(1+r/\sqrt{M}\right)</math>, so the worst-case number of required Newton steps is <math>O(\sqrt{M})</math>. In theory, if ''μ'' is larger (e.g. 2 or more), then the worst-case number of required Newton steps is in <math>O(M)</math>. However, in practice, larger ''μ'' leads to a much faster convergence. These methods are called ''long-step methods''.<ref name=":0" />{{Rp|___location=Sec.4.6}} In practice, if ''μ'' is between 3 and 100, then the program converges within 20-40 Newton steps, regardless of the number of constraints (though the runtime of each Newton step of course grows with the number of constraints). The exact value of ''μ'' within this range has little effect on the performaneperformance.<ref name=":3" />{{Rp|___location=chpt.11}}

== Potential-reduction methods ==

For potential-reduction methods, the problem is presented in the ''conic form'':<ref name=":0" />{{Rp|___location=Sec.5}} <blockquote>'''minimize ''c''^T''x'' s.t. ''x'' in ''{b+L} ᚢ∩ K''''', </blockquote>where ''b'' is a vector in R^''n'', L is a [[linear subspace]] in R^''n'' (so ''b''+''L'' is an [[affine plane]]), and ''K'' is a closed pointed [[convex cone]] with a nonempty interior. Every convex program can be converted to the conic form. To use the potential-reduction method (specifically, the extension of [[Karmarkar's algorithm]] to convex programming), we need the following assumptions:<ref name=":0" />{{Rp|___location=Sec.6}}

* A. The feasible set ''{b+L} ᚢ∩ K'' is bounded, and intersects the interior of the cone ''K''.

* B. We are given in advance a strictly-feasible solution ''x''^, that is, a feasible solution in the interior of ''K''.

* C. We know in advance the optimal objective value, c*, of the problem.

* D. We are given an ''M''-logarithmically-homogeneous [[self-concordant barrier]] ''F'' for the cone ''K''.

Assumptions A, B and D are needed in most interior-point methods. Assumption C is specific to Karmarkar's approach; it can be alleviated by using a "sliding objective value". It is possible to further reduce the program to the ''Karmarkar format'':<blockquote>'''minimize ''s''^T''x'' s.t. ''x'' in ''M ᚢ∩ K'' and ''e''^T''x'' = 1''' </blockquote>where ''M'' is a [[linear subspace]] of in R^''n'', and the optimal objective value is 0.

The method is based on the following [[scalar potential]] function:<blockquote>''v''(''x'') = ''F''(''x'') + ''M'' ln (''s''^T''x'')</blockquote>where ''F'' is the ''M''-self-concordant barrier for the feasible cone. It is possible to prove that, when ''x'' is strictly feasible and ''v''(''x'') is very small (- very negative), ''x'' is approximately-optimal. The idea of the potential-reduction method is to modify ''x'' such that the potential at each iteration drops by at least a fixed constant ''X'' (specifically, ''X''=1/3-ln(4/3)). This implies that, after ''i'' iterations, the difference between objective value and the optimal objective value is at most ''V'' * exp(-''i X'' / ''M''), where ''V'' is a data-dependent constant. Therefore, the number of Newton steps required for an ''ε''-approximate solution is at most <math>O(1) \cdot M \cdot \ln\left(\frac{V}{\varepsilon} + 1\right)+1 </math>.

Line 169 ⟶ 171:

:<math>(x,\lambda) \to (x + \alpha p_x, \lambda + \alpha p_\lambda).</math>[[File:Interior_Point_Trajectory.webm|center|thumb|400x400px|Trajectory of the iterates of ''x'' by using the interior point method.]]

== Types of Convexconvex Programsprograms Solvablesolvable via Interiorinterior-Pointpoint Methodsmethods ==

Here are some special cases of convex programs that can be solved efficiently by interior-point methods.<ref name=":0" />{{Rp|___location=Sec.10}}

Line 183 ⟶ 185:

The function <math>b</math> is self-concordant with parameter ''M''=''m'' (the number of constraints). Therefore, the number of required Newton steps for the path-following method is O(''mn''²), and the total runtime complexity is O(''m''^3/2 ''n''²).{{Clarify|reason=This is the cost for an approximate solution - not an exact solution. The text does not elaborate on this.|date=November 2023}}

===[[Quadratically constrained quadratic program]]s===

Given a quadratically constrained quadratic program of the form:

<math display="block">\begin{aligned}

Line 190 ⟶ 192:

& 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}</math>

where all matrices ''A_j'' are [[Positive semidefinite matrices|positive-semidefinite programmatrices]].

We can apply path-following methods with the barrier

<math display="block">b(x) := -\sum_{j=1}^m \ln(-f_j(x)).</math> The function <math>b</math> is a self-concordant barrier with parameter ''M''=''m''. The Newton complexity is O(''(m+n)n''²), and the total runtime complexity is O(''m''^1/2 (m+n) ''n''²).

===L_p norm approximation===

Consider a problem of the form

<math display="block">\begin{aligned}

\operatorname{minimize}\quad & \sum_j |v_j - u_j^\top x|_p

\end{aligned},</math>

where each <math>u_j</math> is a vector, each <math>v_j</math> is a scalar, and <math>|\cdot|_p</math> is an [[Lp norm|L_p norm]] with <math>1< p < \infty.</math> After converting to the standard form, we can apply path-following methods with a self-concordant barrier with parameter ''M''=4''m''. The Newton complexity is O(''(m+n)n''²), and the total runtime complexity is O(''m''^1/2 (m+n) ''n''²).

===[[Geometric program]]s===

Consider the problem

<math display="block">\begin{aligned}

Line 210 ⟶ 212:

\end{aligned}</math>

There is a self-concordant barrier with parameter 2''k''+''m''. The path-following method has Newton complexity O(''mk''²+''k''³+''n''³) and total complexity O((''k+m'')^1/2[''mk''²+''k''³+''n''³]).

=== [[Semidefinite program]]s ===

Line 227 ⟶ 229:

* {{cite book |last1=Bonnans |first1=J.&nbsp;Frédéric |last2=Gilbert |first2=J.&nbsp;Charles |last3=Lemaréchal |first3=Claude |authorlink3=Claude Lemaréchal |last4=Sagastizábal |first4=Claudia&nbsp;A. |author4-link=Claudia Sagastizábal |title=Numerical optimization: Theoretical and practical aspects |url=https://www.springer.com/mathematics/applications/book/978-3-540-35445-1 |edition=Second revised ed. of translation of 1997 <!-- ''Optimisation numérique: Aspects théoriques et pratiques'' --> French |series=Universitext |publisher=Springer-Verlag |___location=Berlin |year=2006 |pages=xiv+490 |isbn=978-3-540-35445-1 |doi=10.1007/978-3-540-35447-5 |mr=2265882}}

* {{cite book |title=Numerical Optimization |first=Jorge |last=Nocedal |author2=Stephen Wright |year=1999 |publisher=Springer |___location=New York, NY |isbn=978-0-387-98793-4}}

*{{Cite book | last1=Press | first1=WH | last2=Teukolsky | first2=SA | last3=Vetterling | first3=WT | last4=Flannery | first4=BP | year=2007 | title=Numerical Recipes: The Art of Scientific Computing | edition=3rd | publisher=Cambridge University Press | ___location=New York | isbn=978-0-521-88068-8 | chapter=Section 10.11. Linear Programming: Interior-Point Methods | chapter-url=http://apps.nrbook.com/empanel/index.html#pg=537 | access-date=12 August 2011 | archive-date=11 August 2011 | archive-url=https://web.archive.org/web/20110811154417/http://apps.nrbook.com/empanel/index.html#pg=537 | url-status=dead }}

{{Optimization algorithms|convex}}