Source-linked AI summary

Efficient Inverse Maintenance and Faster Algorithms for Linear Programming

Yin Tat Lee, Aaron Sidford

arXiv:1503.01752v3cs.DSmath.NAmath.OC

TL;DR

The paper addresses the cost of solving sequences of slowly changing linear systems arising in optimization. It develops inverse-maintenance algorithms that combine numerical and sparsification techniques, achieving faster amortized maintenance and improved running times for linear programming and related problems. An open limitation remains for Dikin-walk sampling because determinant computation is not efficiently supported.

  • Problem

    Sequences of linear systems are a computational bottleneck in optimization, motivating efficient maintenance of solvers for slowly changing matrices A^T D^(k)A.

  • Method

    The paper combines low-rank updates, preconditioning, fast matrix multiplication, subspace embeddings, and spectral sparsification to maintain approximate solvers under stability assumptions.

  • Results

    ˜O(nnz(A)+d^2) amortized cost per iteration is achieved for inverse maintenance, with improved running times for linear programming and several related optimization problems.

  • Takeaways & Limitations

    The inverse-maintenance framework provides faster algorithms for linear programming, multicommodity flow, polytope rounding, and other applications.

  • Takeaways & Limitations

    For Dikin-walk sampling, the machinery improves linear-system solving but does not efficiently compute the determinant to sufficient accuracy, leaving the overall improvement open.

Abstract

from arXiv · show

In this paper, we consider the following inverse maintenance problem: given $A \in \mathbb{R}^{n\times d}$ and a number of rounds $r$, we receive a $n\times n$ diagonal matrix $D^{(k)}$ at round $k$ and we wish to maintain an efficient linear system solver for $A^{T}D^{(k)}A$ under the assumption $D^{(k)}$ does not change too rapidly. This inverse maintenance problem is the computational bottleneck in solving multiple optimization problems. We show how to solve this problem with $\tilde{O}(nnz(A)+d^ω)$ preprocessing time and amortized $\tilde{O}(nnz(A)+d^{2})$ time per round, improving upon previous running times for solving this problem. Consequently, we obtain the fastest known running times for solving multiple problems including, linear programming and computing a rounding of a polytope. In particular given a feasible point in a linear program with $d$ variables, $n$ constraints, and constraint matrix $A\in\mathbb{R}^{n\times d}$, we show how to solve the linear program in time $\tilde{O}(nnz(A)+d^{2})\sqrt{d}\log(ε^{-1}))$. We achieve our results through a novel combination of classic numerical techniques of low rank update, preconditioning, and fast matrix multiplication as well as recent work on subspace embeddings and spectral sparsification that we hope will be of independent interest.

1 Introduction

The paper formulates inverse maintenance for solving sequences of slowly changing systems A^T D^(k)A and develops faster algorithms under stability assumptions. These improvements yield faster linear-programming and related optimization algorithms.

  • Motivation: Solving sequences of linear systems is a computational bottleneck in interior point methods, polytope sampling, and grouped least-squares algorithms.Because the positive definite matrices change gradually, approximate inverses may reduce average iteration cost.
  • Approach: The approach combines low-rank updates, preconditioning, fast matrix multiplication, subspace embeddings, and spectral sparsification.It maintains a spectral sparsifier formed from a weighted subset of ˜O(d) rows and reduces the number and cost of updates.
  • Problem: The inverse maintenance problem maintains a solver for A^T D^(k)A across r rounds while keeping construction and application costs small.The input matrix A is fixed, while each round supplies a non-negative diagonal weighting vector d^(k).
  • Stability assumptions: The paper studies ℓ2 and weaker σ stability assumptions that bound how rapidly the diagonal weights change between rounds.Under σ stability, changes to important rows are controlled using leverage scores.
  • Contribution: The central goal is efficient inverse maintenance under the weaker σ stability assumption, which covers the applications considered by the paper.The σ assumption is strictly weaker than the ℓ2 assumption.

2 Preliminaries

The paper’s preliminaries establish notation for sparsity, positive definiteness, matrix approximation, leverage scores, and the sampling, solver, Woodbury, and subspace-embedding tools used later.

  • Notation: nnz counts nonzero entries, diagonal vectors represent diagonal matrices, and d typically denotes the smaller problem dimension.The paper uses these conventions throughout its algorithms and applications.
  • Matrix relations: Positive-definite matrix comparisons use Loewner order, while M ≈_ε N means e^-εN ⪯ M ⪯ e^εN.These relations formalize spectral approximation between symmetric matrices.
  • Leverage scores: Leverage scores measure row importance; weighted scores σ_A(d) are the scores of D^1/2A and lie in [0,1].Their sum is bounded by d, supporting score-based sampling and maintenance.
  • Core tools: Leverage-score sampling selects approximately O~(d) rows to obtain a spectral approximation with high probability.The sampling matrix is diagonal and independently retains rows according to prescribed probabilities.
  • Core tools: A T-time solver for A^TA yields approximate leverage scores in O~((nnz(A)+T)ε^-2 log(ε^-1)) time.This connects fast linear-system solving to efficient score estimation.
  • Core tools: The method combines Woodbury updates with subspace embeddings that preserve A^TA spectrally using O~(dε^-2) rows and O~(nnz(A)) sampling time.These results supply the algebraic and dimensionality-reduction machinery for inverse maintenance.

3 Solving the Inverse Maintenance Problem Using ℓ2 Stability

Under ℓ2 stability, the paper maintains approximate inverses for slowly changing weighted systems by combining low-rank updates, fast multiplication, and subspace embeddings. The resulting solver has near-quadratic per-round cost, with improved total maintenance bounds.

  • Low-rank inverse maintenance: Theorem 9 maintains an implicit approximate inverse for B^(k)=A^TD^(k)A under bounded coordinate changes and spectral stability.It applies the inverse approximation to arbitrary vectors in O~(d^2 log β) time.
  • ℓ2-stable solver: Under ℓ2 stability, Algorithm 2 maintains an O~(nnz(A)+d^2)-time solver with high probability in total time stated by Theorem 10.The theorem’s total bound uses s=max{max_k nnz(d^(k)),d} and r rounds.
  • ℓ2-stable solver: The algorithm periodically restarts maintenance, using the dominance of the preprocessing term to obtain the stated total bound.The restart analysis compares the fast-multiplication terms across parameter ranges.
  • Proof strategy: The proof represents changes as low-rank modifications and computes the resulting products efficiently with fast matrix multiplication.Lemma 11 bounds weighted matrix-product maintenance, while Woodbury expresses the inverse update through a smaller matrix.
  • Low-rank inverse maintenance: Theorem 9 improves prior bounds by avoiding an additive O~(nd^(ω−1)) term that was prohibitively expensive for the target applications.The stated previous-best comparison concerns expected running time.
  • Proof strategy: Subspace embeddings replace an expensive explicit spectral computation, and the Woodbury identity converts the maintained low-dimensional correction into an inverse update.The construction maintains matrices associated with changed coordinates and uses fast multiplication for their products.

4 An algorithm for the σ Stable Case

For the stronger σ-stability setting, the paper stabilizes leverage scores, sparsifies the weighted system, and maintains its inverse efficiently. This yields a solver with O~(nnz(A)+d^2) per-round time and near-linear dependence on the number of rounds.

  • Main result: Theorem 13 maintains an O~(nnz(A)+d^2)-time solver with high probability in total time O~(d^ω+r(nnz(A)+d^2)).This is the section’s central σ-stability result.
  • Leverage-score stability: σ stability implies leverage scores remain stable in leverage-score norm, limiting how often the sparsifier requires low-rank updates.This stability lemma is identified as crucial to the inverse-maintenance analysis.
  • Sparsified maintenance: The algorithm samples rows using approximate leverage scores so that A^TH^(k)A spectrally approximates A^TD^(k)A.The maintained inverse is built for the sampled matrix and then used as a solver for the original weighted system.
  • Update complexity: Assuming independence between input changes, leverage-score errors, and the sampled matrix, the expected number of changed H-coordinates is O(r^2 log d).This bound controls the number of updates passed to the low-rank maintenance theorem.
  • Runtime analysis: The implementation restarts periodically and computes leverage scores from the maintained solver, producing the theorem’s total maintenance cost.The analysis combines restart accounting with approximate score computation and fast multiplication.
  • Randomness handling: The independence assumption is removed by transforming the solver into a statistically indistinguishable noisy solver with comparable running time.The paper defers this transformation to Section 5.

5 Hiding Randomness in Linear System Solvers

The paper addresses adversarial dependence on randomized solver outputs by adding carefully calibrated Gaussian noise. NoisySolver remains efficient and is statistically indistinguishable from an exact solver with suitable noise.

  • Motivation: Adaptive inputs or solver outputs correlated with internal randomness could invalidate the inverse-maintenance analysis.This motivates hiding the solver’s randomness from subsequent computation.
  • NoisySolver: NoisySolver combines approximate solves with a Gaussian perturbation to hide randomness while preserving the solver’s intended behavior.It transforms a linear solver into a nonlinear solver that resembles an exact solver plus Gaussian noise.
  • Construction: The construction uses two solver calls with separate accuracy parameters and a Gaussian vector sampled with identity covariance.The displayed algorithm specifies y1, y2, and the random perturbation procedure.
  • Correctness and cost: NoisySolver is an (nnz(A)+T)-time solver when its input solver for A^TA takes T time.Lemma 16 establishes the efficiency guarantee for every b and ε in the stated range.

6 Row Insertion and Removal

The paper extends inverse maintenance to changing matrix rows under the K stability assumption, using reductions to low-rank updates and intermediate steps. It obtains a high-probability solver with amortized near-quadratic maintenance cost for sufficiently small row changes.

  • K Stability Assumption: K stability permits rows of A to be revealed or hidden according to nonzero entries of d^(k), with at most K coordinate changes between rounds.This setting addresses applications where rows are added or removed while diagonal weights also change.
  • Algorithmic Extension: The algorithm avoids requiring the entire matrix A explicitly, unlike previous inverse-maintenance methods.It modifies Algorithm 3 to support row addition and removal under the K stability assumption.
  • Running Time: For K ≤ d^(3−ω)/(ω−1), a variant of Algorithm 3 maintains a solver in total time ˜O(d^ω + r(nnz(A) + d^2)) with high probability.Restarting every d^(ω−2) steps yields amortized ˜O(d^2) maintenance cost under this condition.
  • Stability Reduction: Intermediate steps split insertion or removal changes so consecutive systems remain within a constant-factor spectral approximation.For insertion, O(log(n)) steps suffice; removal is handled analogously.
  • Update Complexity: The expected total number of row changes is O(r^2 + rK), accounting for both small diagonal-weight changes and row resampling.Interactions between the two update cases increase the count by only a constant factor.
  • Running Time: Using ω < 2.37287, the method supports d^0.4568 row additions and removals in amortized ˜O(d^2) time.This specializes the general K bound to the current matrix-multiplication constant.

7 Applications

The paper applies inverse maintenance to linear programming, regression, rounding ellipsoids, multicommodity flow, and related convex optimization problems. These applications reuse slowly changing linear systems to obtain improved running times under stated input and stability conditions.

  • 7.1 Linear Programming: The linear-programming result applies to bounded polytopes represented by equality constraints and coordinate bounds, given an interior initial point.Theorem 21 assumes A ∈ R^{d×n}, d ≤ n, and an interior point x ∈ R^n.
  • 7.1 Linear Programming: The inverse maintenance solver improves linear-programming algorithms by solving slowly changing systems A^T D(k)A to 1/poly(n) accuracy.The resulting algorithm uses d log(U/ϵ) rounds with additional amortized ˜O(nnz(A) + d^2) time per round.
  • 7.2 ℓ1 and ℓ∞ Regression: The method also yields ℓ1 and ℓ∞ regression algorithms, while leaving further improvements for p ∉ {1, 2, ∞} as an open question.The regression corollary specifically covers p = 1 or p = ∞.
  • 7.3 Rounding Ellipsoid for Polytopes: 65%?
  • 7.3 Rounding Ellipsoid for Polytopes: The ellipsoid application computes E satisfying E ⊂ {x ∈ R^d : Ax ≥ b} ⊂ 100d · E from an initial feasible point.The proof combines a lemma characterizing a suitable ellipse with an algorithm for computing it.
  • 7.4 Multicommodity Flow: The framework is also applied to multicommodity flow and yields an exact minimum-cost maximum-flow algorithm for weighted directed graphs.The stated corollary runs in ˜O(|V|^2.5 log^O(1)(U)) time for integer capacities and costs bounded by U.
  • 7.5 Convex Optimization: The inverse-maintenance machinery further supports a separation-oracle algorithm for finding a point in a bounded convex set or certifying the absence of a sufficiently large ball.Theorem 29 gives time O(dT log(dR/ϵ) + d^3 log^O(1)(dR/ϵ)).

8 Open Problem: Sampling from a Polytope

The paper examines sampling from explicitly represented polytopes and explains why faster linear-system solves do not by themselves improve Dikin walk. Determinant computation remains the unresolved bottleneck.

  • Sampling random points from polytopes has applications in optimization, counting, learning, and rounding, with Hit-and-Run and Dikin walk as key methods.
  • Dikin walk is expensive because each iteration solves a linear system and computes determinants for importance sampling.
  • The paper leaves open whether determinant computation can be circumvented to improve Dikin walk's running time.An earlier claimed improvement was withdrawn because it accounted only for linear-system solving.

A Relationships Between T -time Linear Solver and Inverse Matrix

This appendix establishes near-equivalence between maintaining approximate inverses and maintaining linear solvers. A spectral approximation can be converted into a solver using preconditioned iteration, and a solver yields an inverse approximation.

  • A spectral approximation to M^-1 that can be applied in O(T) time yields a linear solver for positive-definite M.The solver runs in linear time in nnz(M) + T, up to logarithmic accuracy factors.
  • The construction uses preconditioned gradient descent, with the number of iterations chosen to achieve the target accuracy.
  • Conversely, a linear solver S(b, ϵ) = Qϵb yields a spectral approximation to the inverse of M.

B Remarks for Figure 1.1

Figure 1.1 compares previous linear-programming running times across algorithmic approaches. The listed bounds differ in their dependence on n, d, and the linear-system solver cost.

  • Figure 1.1 summarizes the previous fastest algorithm for linear programming min c^T x subject to Ax ≥ b.
  • Vaidya's interior-point method achieves one of the displayed bounds, while Karmarkar acceleration is applied to a short-step path-following method.
  • Using the best current linear-system solver with short-step path following gives a bound of ˜O(√n(z + d^2.38)).
  • Vaidya's cutting-plane method gives a displayed bound of ˜O(d(z + d^2.38)).
Loading 1503.01752v3…