🎓How I Study AIHISA
📖Read
📄Papers📰Blogs🎬Courses
💡Learn
🛤️Paths📚Topics💡Concepts🎴Shorts
🎯Practice
⏱️Coach🧩Problems🧠Thinking🎯Prompts🧠Review
SearchSettings
How I Study AI - Learn AI Papers & Lectures the Easy Way
⚙️AlgorithmAdvanced

Interior Point Methods

Key Points

  • •
    Interior point methods solve constrained optimization by replacing hard constraints with a smooth barrier that becomes infinite at the boundary, keeping iterates strictly inside the feasible region.
  • •
    The logarithmic barrier transforms problems like min f(x) subject to gi​(x) ≤ 0 into a sequence of unconstrained problems min t f(x) + φ(x), where φ(x) = -∑ ln(-gi​(x)).
  • •
    As the barrier parameter t increases, solutions of the barrier problem trace the central path and approach the true constrained optimum.
  • •
    Each barrier subproblem is typically solved efficiently by Newton’s method with backtracking line search and a feasibility-preserving step to stay inside constraints.
  • •
    Primal-dual interior point methods solve KKT conditions directly, often using a Schur complement system, and enjoy strong polynomial-time properties for convex problems like LP and QP.
  • •
    Per-iteration cost is dominated by solving linear systems; for dense problems this is about O(n3) (barrier Newton) or O(m3 + m2 n) (primal-dual with m constraints, n variables).
  • •
    These methods are robust for large, well-scaled convex problems but require a strictly feasible start (or a Phase I) and careful numerical linear algebra.
  • •
    Common pitfalls include starting infeasible, taking steps that hit the boundary, poor scaling, and stopping too early before complementarity is small.

Prerequisites

  • →Linear algebra (vectors, matrices, factorizations) — Interior point methods are dominated by linear system solves (Gaussian elimination, Cholesky) and require understanding of matrix operations and conditioning.
  • →Multivariable calculus (gradients, Hessians) — Barrier Newton steps need accurate gradients and Hessians; understanding curvature is essential to implement and analyze Newton’s method.
  • →Convex optimization basics — Convergence guarantees rely on convexity, self-concordance, and KKT conditions; recognizing when assumptions hold is crucial.
  • →Newton’s method and line search — The core inner loop uses Newton directions with backtracking and feasibility-preserving step sizes.
  • →KKT conditions and duality — Primal-dual methods directly solve KKT systems; stopping criteria and duality gap use these concepts.
  • →Numerical stability and scaling — Good scaling and stable factorizations are vital to avoid breakdowns in linear solves and to ensure reliable progress.

Detailed Explanation

Tap terms for definitions

01Overview

Interior point methods are a family of algorithms for solving constrained optimization problems—especially convex ones such as linear programming (LP) and quadratic programming (QP). The core idea is to avoid touching the boundary of the feasible set by adding a barrier term to the objective that "blows up" near constraint boundaries. Instead of optimizing over a jagged feasible polytope or curved region directly, the algorithm optimizes a smooth surrogate objective inside the region. By gradually reducing the barrier’s influence (equivalently, increasing a parameter t that emphasizes the original objective), the solutions to the surrogate problems approach the true constrained optimum. Two widely used variants are the barrier (or path-following) method and the primal-dual method. The barrier method solves a sequence of unconstrained problems with Newton’s method and a line search that preserves feasibility. The primal-dual method solves the Karush–Kuhn–Tucker (KKT) conditions directly by Newton’s method, updating both primal variables and dual variables (including slack variables) to maintain complementarity. Both enjoy strong theoretical guarantees for convex problems and are the workhorses of modern LP/QP solvers. In practice, interior point methods require accurate gradients/Hessians (or structured linear algebra), good scaling, and a feasible starting point (or an auxiliary Phase I procedure). Their main computational cost lies in forming and solving linear systems, which can be accelerated using sparse/dense factorizations and exploiting problem structure.

02Intuition & Analogies

Imagine walking through a long hallway whose walls are the constraints. A naive strategy would be to walk while occasionally bouncing off the walls when you hit them—that’s like dealing with constraints directly and can be awkward. Interior point methods instead put a soft but very strong cushion on the walls. As you approach a wall, the cushion’s pressure grows dramatically, pushing you back toward the center. You never hit the wall; you glide smoothly through the middle of the hallway, heading toward the exit that minimizes your travel cost. The cushion is the barrier function. For inequality constraints like Ax ≤ b or g_i(x) ≤ 0, a logarithmic barrier -ln(b_i - a_i^T x) or -ln(-g_i(x)) becomes extremely large as you approach the boundary (where b_i - a_i^T x → 0 or g_i(x) → 0). When you add up these cushions for all constraints and also include the original cost, you get a smooth landscape with a deep valley somewhere inside the hallway. By slowly turning down the cushion’s thickness (increasing t), the valley shifts closer to the true exit while still discouraging you from scraping the walls. Newton’s method acts like a smart guide that looks at both the slope (gradient) and the local curvature (Hessian) of this landscape to decide fast, curved steps toward the valley’s bottom. A line search ensures each step stays within the hallway and makes progress. In the primal-dual version, you carry two maps at once: one for the primal hallway (x) and one for the dual pressures on the walls (s, y). Keeping the product of position and pressure moderate (complementarity) ensures you move along the central path—roughly the centerline of the hallway—toward the true exit efficiently and stably.

03Formal Definition

Consider a convex optimization problem: minimize f(x) subject to gi​(x) ≤ 0 for i=1,...,m and Aeq​ x=b. The logarithmic barrier method replaces the inequalities by a barrier term φ(x) = -∑i=1m​ ln(-gi​(x)), defined only when gi​(x) < 0. For a parameter t>0, we solve the barrier subproblem: minimize Ft​(x) = t f(x) + φ(x) subject to Aeq​ x=b and gi​(x) < 0. As t → ∞, the solution x^*(t) to this subproblem converges to a solution of the original constrained problem. When gi​ are affine (typical in LP/QP), φ is self-concordant and Newton’s method applied to Ft​ enjoys global convergence with iteration bounds tied to the barrier parameter ν (for the log barrier, ν = m). The central path is the trajectory x^*(t) satisfying the optimality conditions of the barrier subproblems; it approaches the KKT point of the original problem as t grows. Primal-dual interior point methods introduce dual variables y (for equalities) and s≥0 (for inequalities) and solve the KKT system with complementarity X S e = μ e, where μ → 0 along iterations. Each iteration computes a damped Newton step for the residual equations, often via a Schur complement system A D AT Δy=rhs with D diagonal from current (x, s). A fraction-to-the-boundary step ensures x and s remain positive, maintaining feasibility of the interior.

04When to Use

  • Linear Programming (LP): For large, dense or moderately sparse LPs, interior point methods scale well and provide reliable high-accuracy solutions.
  • Quadratic Programming (QP): Convex QPs with linear constraints benefit from the same machinery, with only modest changes in linear algebra (Hessian from the quadratic term).
  • Convex problems with smooth inequality constraints: If g_i are smooth (and preferably convex/affine), the logarithmic barrier is natural and Newton’s method is effective.
  • When you can exploit structure: Block structure, sparsity, or low-rank updates in A or the Hessian can make interior point methods extremely fast via specialized factorizations.
  • High-accuracy solutions: If you need small duality gap and certified optimality via KKT residuals, primal-dual interior point methods provide principled stopping criteria. Avoid using interior point methods when constraints are nonconvex without additional safeguards, or when you lack a strictly feasible starting point and cannot run a reliable Phase I. For very small LPs where basis information is valuable, simplex may be simpler and faster in practice.

⚠️Common Mistakes

  • Starting infeasible or on the boundary: The log barrier is undefined on the boundary. Use a strictly feasible start or a Phase I procedure (e.g., find the analytic center of {x | Ax < b, x > 0}).
  • Overshooting the boundary in line search: Always compute a fraction-to-the-boundary step size so that x + α Δx stays strictly inside (and similarly for slack variables). Backtracking alone is insufficient without feasibility checks.
  • Increasing the barrier parameter too aggressively: If t grows too fast, Newton steps become poorly conditioned and iterations may stall. Use moderate growth (e.g., t ← μ t with μ in [8, 20]) and solve each subproblem to sufficient accuracy.
  • Ignoring scaling/conditioning: Poorly scaled A or wildly different variable magnitudes make Hessians ill-conditioned, harming Newton steps. Pre-scale variables/constraints and prefer numerically stable factorizations (Cholesky/LDL^T with pivoting).
  • Weak stopping criteria: Stop only when primal residuals, dual residuals, and complementarity (or m/t in barrier form) are all small. Checking just the gradient norm of F_t can be misleading.
  • Mishandling Hessians for nonlinear g_i: For general g_i, the barrier Hessian includes both ∑ (∇g_i ∇g_i^T)/g_i^2 and -∑ ∇^2 g_i / g_i. Forgetting the second term breaks Newton’s method. For affine g_i this term vanishes.

Key Formulas

Primal Problem

xmin​f(x)s.t.gi​(x)≤0 (i=1..m), Aeq​x=b

Explanation: This is a generic constrained optimization problem with inequality and equality constraints. Interior point methods target this structure, especially when f and gi​ are convex.

Logarithmic Barrier

ϕ(x)=−i=1∑m​ln(−gi​(x))

Explanation: The barrier is finite only when all inequalities are strictly satisfied (gi​(x) < 0) and goes to infinity as any constraint approaches its boundary. It discourages iterates from touching the boundary.

Barrier Subproblem

xmin​Ft​(x)=tf(x)+ϕ(x)s.t.Aeq​x=b, gi​(x)<0

Explanation: For a given t>0, we solve this smooth problem. As t increases, its minimizer x^*(t) approaches the constrained solution along the central path.

Barrier Derivatives

∇ϕ(x)=−i=1∑m​gi​(x)∇gi​(x)​,∇2ϕ(x)=i=1∑m​gi​(x)2∇gi​(x)∇gi​(x)⊤​−i=1∑m​gi​(x)∇2gi​(x)​

Explanation: These are the gradient and Hessian of the log barrier. For affine constraints gi​(x) = a_i⊤x - bi​, the second-derivative term vanishes, simplifying computation.

Central Path Limit

x∗(t)=argxmin​Ft​(x),andx∗(t)→x⋆ as t→∞

Explanation: The minimizers of the barrier problems trace the central path and converge to an optimal solution x^⋆ of the original constrained problem.

Duality Gap for Log Barrier

gap(t)≤tm​

Explanation: For the log barrier with m inequalities, the suboptimality of x^*(t) is at most m/t. This provides a principled stopping rule for increasing t.

Primal-Dual Newton System (LP)

​0AS​A⊤00​I0X​​​ΔxΔyΔs​​=−​rd​rp​rc​−σμe​​

Explanation: This linear system arises from Newton’s method on the KKT conditions for LP. X and S are diagonal matrices of primal variables and slacks, σ is the centering parameter, and μ = (x⊤s)/n.

Schur Complement (LP)

(ADA⊤)Δy=−rp​+A(S−1rc​−Drd​),D=XS−1

Explanation: Eliminating Δx and Δs yields a symmetric positive definite system in Δy only. This is typically solved by Cholesky factorization for efficiency and stability.

Fraction-to-the-Boundary Step

αmax​=min{Δxj​<0min​−Δxj​xj​​, Δsj​<0min​−Δsj​sj​​},  α=ταmax​, 0<τ<1

Explanation: To maintain strict positivity of x and s, take only a fraction τ of the maximal feasible step that keeps all variables positive.

Newton Decrement

λ(x)=∇Ft​(x)⊤(∇2Ft​(x))−1∇Ft​(x)​

Explanation: This quantity measures how close x is to the minimizer of Ft​. A small value indicates that a Newton step will yield little improvement, suggesting convergence.

Complexity Analysis

For the pure barrier (path-following) method with affine inequalities gi​(x) = a_i⊤x - bi​, each Newton step requires building the gradient and Hessian and solving a linear system. Building the Hessian involves summing m rank-1 outer products ai​ a_i⊤ scaled by (bi​ - a_i⊤x)^{-2}, which costs O(m n2) for dense A (m constraints, n variables). Adding the diagonal contribution from positivity constraints costs O(n). Solving the resulting dense n×n Newton system by Gaussian elimination or Cholesky factorization costs O(n3). Thus, per iteration the dominant cost is O(n3 + m n2). The number of damped Newton iterations per barrier subproblem is typically small (dozens at most), and with a self-concordant barrier the iteration bound is O(m​ log(1/ε)) to reach ε−accuracy for each t. For primal-dual interior point methods applied to LP in equality form (Ax=b, x≥0), each iteration solves the Schur complement system (A D A⊤) Δy=rhs, where D=X S−1 is diagonal. Forming M=A D A⊤ costs O(m2 n) for dense matrices (first compute A D in O(m n), then multiply by A⊤ in O(m2 n)). Solving the m×m SPD system by Cholesky costs O(m3). Back-substitution to obtain Δx and Δs is O(n). Hence, the per-iteration cost is O(m3 + m2 n + n). When m ≪ n (common in many LPs), this is significantly cheaper than solving an n×n system directly. Practical implementations exploit sparsity and structure, reducing these bounds dramatically with sparse Cholesky or LDLT factorizations. Overall iteration counts are modest: primal-dual methods usually converge in 20–60 iterations to high accuracy on large problems. The total running time is then dominated by the linear algebra in each iteration. Memory use is dominated by storing A (O(m n)) and factorization data (e.g., O(m2) for dense Cholesky), plus O(n) vectors for x, s, and O(m) for y.

Code Examples

Logarithmic barrier method for LP: min c^T x s.t. A x ≤ b and x > 0
1#include <bits/stdc++.h>
2using namespace std;
3
4struct BarrierLPSolver {
5    // Solve: minimize c^T x subject to A x <= b and x > 0 using a log-barrier with Newton's method.
6    vector<vector<double>> A; // m x n
7    vector<double> b;         // m
8    vector<double> c;         // n
9    int m, n;
10
11    BarrierLPSolver(const vector<vector<double>>& A_, const vector<double>& b_, const vector<double>& c_)
12        : A(A_), b(b_), c(c_) {
13        m = (int)A.size();
14        n = (int)A[0].size();
15    }
16
17    // Basic linear algebra helpers
18    static double dot(const vector<double>& x, const vector<double>& y) {
19        double s = 0.0; for (size_t i = 0; i < x.size(); ++i) s += x[i]*y[i]; return s;
20    }
21    static double norm2(const vector<double>& x) { return sqrt(dot(x, x)); }
22    static vector<double> add(const vector<double>& x, const vector<double>& y, double a=1.0) {
23        vector<double> r(x.size()); for (size_t i=0;i<x.size();++i) r[i] = x[i] + a*y[i]; return r;
24    }
25    static vector<double> scale(const vector<double>& x, double a) {
26        vector<double> r(x.size()); for (size_t i=0;i<x.size();++i) r[i] = a*x[i]; return r;
27    }
28
29    vector<double> matVec(const vector<vector<double>>& M, const vector<double>& x) const {
30        vector<double> r(M.size(), 0.0);
31        for (int i=0;i<(int)M.size();++i) {
32            double s=0.0; for (int j=0;j<(int)x.size();++j) s += M[i][j]*x[j]; r[i]=s;
33        }
34        return r;
35    }
36
37    // Solve linear system H dx = rhs by Gaussian elimination with partial pivoting (dense, O(n^3))
38    bool solveLinear(vector<vector<double>> H, vector<double>& rhs, vector<double>& dx) const {
39        int n = (int)H.size();
40        // Augment
41        for (int i=0;i<n;++i) H[i].push_back(rhs[i]);
42        for (int col=0; col<n; ++col) {
43            // Pivot
44            int piv = col; double best = fabs(H[col][col]);
45            for (int r=col+1;r<n;++r) if (fabs(H[r][col]) > best) { best=fabs(H[r][col]); piv=r; }
46            if (best < 1e-14) return false; // singular/ill-conditioned
47            if (piv != col) swap(H[piv], H[col]);
48            // Normalize
49            double diag = H[col][col];
50            for (int j=col;j<=n;++j) H[col][j] /= diag;
51            // Eliminate below
52            for (int r=col+1;r<n;++r) {
53                double f = H[r][col];
54                for (int j=col;j<=n;++j) H[r][j] -= f * H[col][j];
55            }
56        }
57        // Back substitution
58        dx.assign(n, 0.0);
59        for (int i=n-1;i>=0;--i) {
60            double s = H[i][n];
61            for (int j=i+1;j<n;++j) s -= H[i][j]*dx[j];
62            dx[i] = s; // since diag normalized to 1
63        }
64        return true;
65    }
66
67    // Compute objective F_t(x) = t*c^T x - sum ln(b_i - a_i^T x) - sum ln x_j
68    double barrierObjective(const vector<double>& x, double t) const {
69        double val = t * dot(c, x);
70        // Inequality slacks s_i = b_i - a_i^T x must be positive
71        vector<double> Ax = matVec(A, x);
72        for (int i=0;i<m;++i) {
73            double s = b[i] - Ax[i];
74            if (s <= 0) return numeric_limits<double>::infinity();
75            val -= log(s);
76        }
77        for (int j=0;j<n;++j) {
78            if (x[j] <= 0) return numeric_limits<double>::infinity();
79            val -= log(x[j]);
80        }
81        return val;
82    }
83
84    // Compute gradient and Hessian of F_t at x
85    void gradHess(const vector<double>& x, double t, vector<double>& grad, vector<vector<double>>& H) const {
86        grad.assign(n, 0.0);
87        H.assign(n, vector<double>(n, 0.0));
88        // t * c
89        for (int j=0;j<n;++j) grad[j] += t * c[j];
90        vector<double> Ax = matVec(A, x);
91        // Sum over inequality constraints: -ln(b_i - a_i^T x)
92        for (int i=0;i<m;++i) {
93            double s = b[i] - Ax[i]; // slack > 0
94            double inv = 1.0 / s;
95            // gradient += a_i / s
96            for (int j=0;j<n;++j) grad[j] += A[i][j] * inv;
97            // Hessian += a_i a_i^T / s^2
98            double inv2 = inv * inv;
99            for (int r=0;r<n;++r) {
100                double ar = A[i][r];
101                if (ar == 0.0) continue;
102                for (int c_=0;c_<n;++c_) H[r][c_] += ar * A[i][c_] * inv2;
103            }
104        }
105        // Positivity barrier -sum ln x_j: grad -= 1/x_j, Hessian += diag(1/x_j^2)
106        for (int j=0;j<n;++j) {
107            grad[j] -= 1.0 / x[j];
108            H[j][j] += 1.0 / (x[j]*x[j]);
109        }
110    }
111
112    // Compute maximal feasible step (fraction-to-the-boundary)
113    double maxFeasibleStep(const vector<double>& x, const vector<double>& dx, double tau=0.99) const {
114        double alpha = 1.0;
115        // Maintain x + alpha*dx > 0
116        for (int j=0;j<n;++j) if (dx[j] < 0) alpha = min(alpha, -tau * x[j] / dx[j]);
117        // Maintain b_i - a_i^T (x + alpha dx) > 0 => s_i + alpha ds_i > 0
118        vector<double> Ax = matVec(A, x);
119        // ds_i = -a_i^T dx
120        vector<double> Adx = matVec(A, dx);
121        for (int i=0;i<m;++i) {
122            double s = b[i] - Ax[i];
123            double ds = -Adx[i];
124            if (ds < 0) alpha = min(alpha, -0.99 * s / ds);
125        }
126        alpha = max(0.0, min(1.0, alpha));
127        return alpha;
128    }
129
130    // Solve the barrier problem for a given t by damped Newton
131    bool solveBarrierSubproblem(vector<double>& x, double t, int max_newton=50, double tol=1e-8) {
132        for (int it=0; it<max_newton; ++it) {
133            vector<double> g; vector<vector<double>> H;
134            gradHess(x, t, g, H);
135            double gnorm = norm2(g);
136            if (gnorm < tol) return true; // first-order optimality for subproblem
137            vector<double> dx;
138            vector<double> rhs = g; for (double &v: rhs) v = -v;
139            if (!solveLinear(H, rhs, dx)) return false;
140            // Feasibility-preserving step size
141            double alpha = maxFeasibleStep(x, dx, 0.99);
142            // Backtracking on barrier objective (Armijo)
143            double f0 = barrierObjective(x, t);
144            double dec = 1e-4 * dot(g, dx);
145            while (true) {
146                vector<double> xnew = add(x, dx, alpha);
147                double f1 = barrierObjective(xnew, t);
148                if (isfinite(f1) && f1 <= f0 + dec * alpha) { x = xnew; break; }
149                alpha *= 0.5; if (alpha < 1e-12) return false; // step too small
150            }
151        }
152        return true;
153    }
154
155    // Main barrier loop
156    vector<double> solve(vector<double> x0, double t0=1.0, double mu=15.0, double tol=1e-6) {
157        vector<double> x = x0;
158        double t = t0;
159        while (true) {
160            // Stop when m/t <= tol (approximate duality gap criterion)
161            if ((double)m / t <= tol) break;
162            bool ok = solveBarrierSubproblem(x, t, 60, 1e-9);
163            if (!ok) {
164                cerr << "Newton failed; consider better scaling or initialization.\n";
165                break;
166            }
167            t *= mu; // increase emphasis on original objective
168        }
169        return x;
170    }
171};
172
173int main() {
174    // Example: min c^T x subject to A x <= b, x > 0
175    // c = [1, 1], constraints: x1 + x2 <= 1.8, x1 <= 1.2, x2 <= 1.2
176    vector<vector<double>> A = {
177        {1.0, 1.0},
178        {1.0, 0.0},
179        {0.0, 1.0}
180    };
181    vector<double> b = {1.8, 1.2, 1.2};
182    vector<double> c = {1.0, 1.0};
183
184    BarrierLPSolver solver(A, b, c);
185    vector<double> x0 = {0.5, 0.5}; // strictly feasible start
186    vector<double> x = solver.solve(x0, 1.0, 12.0, 1e-8);
187
188    cout.setf(std::ios::fixed); cout << setprecision(6);
189    cout << "Solution x*: [" << x[0] << ", " << x[1] << "]\n";
190    cout << "Objective c^T x*: " << (x[0]*c[0] + x[1]*c[1]) << "\n";
191    return 0;
192}
193

This implements a classical log-barrier method for an LP with inequalities Ax ≤ b and positivity x > 0. The barrier subproblem minimizes F_t(x) = t c^T x − ∑ ln(b_i − a_i^T x) − ∑ ln x_j. Each inner iteration computes the Newton direction by forming the gradient and Hessian and solving H Δx = −∇F_t. A fraction-to-the-boundary step ensures feasibility, while backtracking (Armijo) guarantees sufficient decrease. The outer loop increases t until the duality gap m/t is small. The example solves a simple 2D LP from a strictly feasible start.

Time: Per Newton step: O(n^3 + m n^2) for dense matrices (Hessian build + linear solve). Outer barrier iterations are O(log(1/ε)) with a moderate multiplier; total iterations are typically small.Space: O(n^2 + m n) to store dense Hessian and constraints (plus O(n) vectors).
Primal-dual interior point method (Mehrotra-style) for small LP in equality form
1#include <bits/stdc++.h>
2using namespace std;
3
4// Solve LP: minimize c^T x subject to A x = b, x >= 0.
5// Simple primal-dual IPM with Mehrotra predictor-corrector for small dense problems.
6struct PrimalDualLP {
7    vector<vector<double>> A; // m x n
8    vector<double> b;         // m
9    vector<double> c;         // n
10    int m, n;
11
12    PrimalDualLP(const vector<vector<double>>& A_, const vector<double>& b_, const vector<double>& c_)
13        : A(A_), b(b_), c(c_) { m = (int)A.size(); n = (int)A[0].size(); }
14
15    // Basic LA
16    static double dot(const vector<double>& x, const vector<double>& y) { double s=0; for(size_t i=0;i<x.size();++i) s+=x[i]*y[i]; return s; }
17    static double norm2(const vector<double>& x) { return sqrt(dot(x,x)); }
18    static vector<double> add(const vector<double>& x, const vector<double>& y, double a=1.0) { vector<double> r(x.size()); for(size_t i=0;i<x.size();++i) r[i]=x[i]+a*y[i]; return r; }
19    static vector<double> hadamard(const vector<double>& x, const vector<double>& y) { vector<double> r(x.size()); for(size_t i=0;i<x.size();++i) r[i]=x[i]*y[i]; return r; }
20
21    vector<double> ATy(const vector<double>& y) const {
22        vector<double> r(n, 0.0);
23        for (int j=0;j<n;++j) {
24            double s=0.0; for (int i=0;i<m;++i) s += A[i][j]*y[i]; r[j]=s;
25        }
26        return r;
27    }
28    vector<double> Ax(const vector<double>& x) const {
29        vector<double> r(m, 0.0);
30        for (int i=0;i<m;++i) {
31            double s=0.0; for (int j=0;j<n;++j) s += A[i][j]*x[j]; r[i]=s;
32        }
33        return r;
34    }
35
36    // Cholesky solve for SPD matrix M (dense). Returns false if fails.
37    bool choleskySolve(vector<vector<double>> M, vector<double>& rhs, vector<double>& y) const {
38        int m = (int)M.size();
39        // Decompose M = L L^T
40        for (int k=0;k<m;++k) {
41            if (M[k][k] <= 0) return false; // not SPD (basic check)
42            M[k][k] = sqrt(M[k][k]);
43            for (int i=k+1;i<m;++i) M[i][k] /= M[k][k];
44            for (int j=k+1;j<m;++j)
45                for (int i=j;i<m;++i)
46                    M[i][j] -= M[i][k]*M[j][k];
47        }
48        // Forward solve L z = rhs
49        vector<double> z(m);
50        for (int i=0;i<m;++i) {
51            double s = rhs[i];
52            for (int k=0;k<i;++k) s -= M[i][k]*z[k];
53            z[i] = s / M[i][i];
54        }
55        // Backward solve L^T y = z
56        y.assign(m, 0.0);
57        for (int i=m-1;i>=0;--i) {
58            double s = z[i];
59            for (int k=i+1;k<m;++k) s -= M[k][i]*y[k];
60            y[i] = s / M[i][i];
61        }
62        return true;
63    }
64
65    // Build Schur complement M = A * D * A^T and rhs = -r_p + A*(S^{-1} r_c - D r_d)
66    bool schurSolve(const vector<double>& x, const vector<double>& s,
67                    const vector<double>& r_p, const vector<double>& r_d, const vector<double>& r_c,
68                    vector<double>& dy, vector<double>& dx, vector<double>& ds) const {
69        vector<double> D(n), Sinv_rc(n), D_rd(n);
70        for (int j=0;j<n;++j) {
71            D[j] = x[j] / s[j];
72            Sinv_rc[j] = r_c[j] / s[j];
73            D_rd[j] = D[j] * r_d[j];
74        }
75        // M = A D A^T
76        vector<vector<double>> M(m, vector<double>(m, 0.0));
77        // First compute T = A * diag(D)
78        vector<vector<double>> T(m, vector<double>(n, 0.0));
79        for (int i=0;i<m;++i) for (int j=0;j<n;++j) T[i][j] = A[i][j] * D[j];
80        // M = T * A^T
81        for (int i=0;i<m;++i) for (int k=0;k<m;++k) {
82            double ssum = 0.0; for (int j=0;j<n;++j) ssum += T[i][j] * A[k][j]; M[i][k] = ssum;
83        }
84        // rhs = -r_p + A*(S^{-1} r_c - D r_d)
85        vector<double> v(n);
86        for (int j=0;j<n;++j) v[j] = Sinv_rc[j] - D_rd[j];
87        vector<double> Av = Ax(v);
88        vector<double> rhs(m);
89        for (int i=0;i<m;++i) rhs[i] = -r_p[i] + Av[i];
90        if (!choleskySolve(M, rhs, dy)) return false;
91        // Recover ds = -r_d - A^T dy
92        vector<double> ATdy = ATy(dy);
93        ds.resize(n);
94        for (int j=0;j<n;++j) ds[j] = -r_d[j] - ATdy[j];
95        // Recover dx = -S^{-1} (X ds + r_c)
96        dx.resize(n);
97        for (int j=0;j<n;++j) dx[j] = -(x[j]*ds[j] + r_c[j]) / s[j];
98        return true;
99    }
100
101    struct Result { vector<double> x, y, s; int iters; };
102
103    Result solve(vector<double> x, vector<double> y, vector<double> s, int max_iter=100, double tol=1e-8) {
104        // Ensure strict positivity
105        for (double &xi : x) xi = max(xi, 1e-2);
106        for (double &si : s) si = max(si, 1e-2);
107
108        for (int it=0; it<max_iter; ++it) {
109            // Residuals
110            vector<double> rp = add(Ax(x), b, -1.0); // A x - b
111            vector<double> rATy = ATy(y);
112            vector<double> rd(n); for (int j=0;j<n;++j) rd[j] = rATy[j] + s[j] - c[j];
113            vector<double> XS(n); for (int j=0;j<n;++j) XS[j] = x[j]*s[j];
114            double mu = dot(XS, vector<double>(n, 1.0)) / n; // x^T s / n
115
116            double rp_n = norm2(rp), rd_n = norm2(rd);
117            if (rp_n < tol && rd_n < tol && mu < tol) {
118                return {x, y, s, it+1};
119            }
120
121            // Predictor (affine-scaling): sigma = 0
122            vector<double> dy_aff, dx_aff, ds_aff;
123            if (!schurSolve(x, s, rp, rd, XS, dy_aff, dx_aff, ds_aff)) {
124                cerr << "Schur solve failed (predictor).\n"; break; }
125            // Step lengths for affine step
126            double alpha_aff_pri = 1.0, alpha_aff_dual = 1.0;
127            for (int j=0;j<n;++j) if (dx_aff[j] < 0) alpha_aff_pri = min(alpha_aff_pri, -x[j]/dx_aff[j]);
128            for (int j=0;j<n;++j) if (ds_aff[j] < 0) alpha_aff_dual = min(alpha_aff_dual, -s[j]/ds_aff[j]);
129            alpha_aff_pri = min(1.0, 0.99*alpha_aff_pri);
130            alpha_aff_dual = min(1.0, 0.99*alpha_aff_dual);
131            // mu_aff
132            double mu_aff = 0.0;
133            for (int j=0;j<n;++j) {
134                double xa = x[j] + alpha_aff_pri * dx_aff[j];
135                double sa = s[j] + alpha_aff_dual * ds_aff[j];
136                mu_aff += xa * sa;
137            }
138            mu_aff /= n;
139            double sigma = pow(mu_aff / mu, 3.0);
140
141            // Corrector: r_c = X S e + dx_aff ⊙ ds_aff - sigma * mu * e
142            vector<double> rc(n);
143            for (int j=0;j<n;++j) rc[j] = XS[j] + dx_aff[j]*ds_aff[j] - sigma * mu;
144
145            vector<double> dy, dx, ds;
146            if (!schurSolve(x, s, rp, rd, rc, dy, dx, ds)) {
147                cerr << "Schur solve failed (corrector).\n"; break; }
148
149            // Step sizes (fraction-to-the-boundary)
150            double alpha_pri = 1.0, alpha_dual = 1.0;
151            for (int j=0;j<n;++j) if (dx[j] < 0) alpha_pri = min(alpha_pri, -x[j]/dx[j]);
152            for (int j=0;j<n;++j) if (ds[j] < 0) alpha_dual = min(alpha_dual, -s[j]/ds[j]);
153            alpha_pri = min(1.0, 0.99*alpha_pri);
154            alpha_dual = min(1.0, 0.99*alpha_dual);
155
156            // Update
157            for (int j=0;j<n;++j) x[j] += alpha_pri * dx[j];
158            for (int i=0;i<m;++i) y[i] += alpha_dual * dy[i];
159            for (int j=0;j<n;++j) s[j] += alpha_dual * ds[j];
160        }
161        return {x, y, s, max_iter};
162    }
163};
164
165int main(){
166    // Example LP: minimize [1,2]^T x subject to [1 1] x = 1, x >= 0. Optimum is x = [1, 0].
167    vector<vector<double>> A = { {1.0, 1.0} }; // m=1, n=2
168    vector<double> b = {1.0};
169    vector<double> c = {1.0, 2.0};
170
171    PrimalDualLP solver(A, b, c);
172    vector<double> x0 = {0.5, 0.5}; // feasible and positive
173    vector<double> y0 = {0.0};
174    vector<double> s0 = {1.0, 2.0}; // c - A^T y0 = c; strictly positive
175
176    auto res = solver.solve(x0, y0, s0, 100, 1e-10);
177
178    cout.setf(std::ios::fixed); cout << setprecision(8);
179    cout << "Iterations: " << res.iters << "\n";
180    cout << "x*: [" << res.x[0] << ", " << res.x[1] << "]\n";
181    cout << "y*: [" << res.y[0] << "]\n";
182    cout << "s*: [" << res.s[0] << ", " << res.s[1] << "]\n";
183    cout << "Objective c^T x*: " << (c[0]*res.x[0] + c[1]*res.x[1]) << "\n";
184    return 0;
185}
186

This is a compact primal-dual interior point solver for small dense LPs in equality form. It implements Mehrotra’s predictor-corrector: compute an affine-scaling step (σ = 0), estimate μ_aff and σ = (μ_aff/μ)^3, then compute a corrector step using a modified complementarity residual. The KKT system is reduced to a symmetric positive definite Schur complement system (A D A^T) Δy = rhs and solved by a simple Cholesky factorization. Fraction-to-the-boundary rules keep x and s positive. The example problem has a known solution x* = [1, 0].

Time: Each iteration forms M = A D A^T in O(m^2 n) and solves it in O(m^3). Recovering Δx, Δs is O(n). Iteration counts are typically tens for good accuracy.Space: O(m n) to store A plus O(m^2) for the dense Cholesky factor and O(n) vectors for x, s.
#interior point method#logarithmic barrier#central path#primal-dual#kkt conditions#schur complement#newton method#line search#duality gap#self-concordant#mehrotra#fraction to the boundary#cholesky factorization#linear programming#quadratic programming