Geometric Median in Nearly Linear Time

Michael B. Cohen; Yin Tat Lee; Gary Miller; Jakub Pachocki; Aaron Sidford

Source-linked AI summary

Geometric Median in Nearly Linear Time

Michael B. Cohen, Yin Tat Lee, Gary Miller, Jakub Pachocki, Aaron Sidford

arXiv:1606.05225v1Jun 16, 2016cs.DSmath.OC

TL;DR

The paper addresses how to compute (1+ε)-approximate geometric medians faster despite limited prior theoretical guarantees. It combines a long-step interior point method with stochastic subgradient descent, achieving nearly linear and sublinear running times.

Problem
The geometric median problem asks for a point minimizing the sum of Euclidean distances to n points in R^d, with faster approximate algorithms needed beyond prior guarantees.
Method
The paper develops a tailored long-step interior point method and uses sampling with stochastic subgradient descent for a separate algorithm.
Results
O(nd log^3(1/ε)) and O(dε^-2) compute a (1+ε)-approximate geometric median, with the latter requiring an oracle for sampling a random input point.
Takeaways & Limitations
The work demonstrates a nearly linear-time interior point algorithm for a canonical optimization problem traditionally requiring superlinear time.
Takeaways & Limitations
The paper leaves open the minimal self-concordance of a barrier for the relevant convex set and whether nearly linear-time interior point algorithms should generally be expected here.

Abstract

from arXiv · show

In this paper we provide faster algorithms for solving the geometric median problem: given $n$ points in $\mathbb{R}^{d}$ compute a point that minimizes the sum of Euclidean distances to the points. This is one of the oldest non-trivial problems in computational geometry yet despite an abundance of research the previous fastest algorithms for computing a $(1+ε)$-approximate geometric median were $O(d\cdot n^{4/3}ε^{-8/3})$ by Chin et. al, $\tilde{O}(d\exp{ε^{-4}\logε^{-1}})$ by Badoiu et. al, $O(nd+\mathrm{poly}(d,ε^{-1})$ by Feldman and Langberg, and $O((nd)^{O(1)}\log\frac{1}ε)$ by Parrilo and Sturmfels and Xue and Ye. In this paper we show how to compute a $(1+ε)$-approximate geometric median in time $O(nd\log^{3}\frac{1}ε)$ and $O(dε^{-2})$. While our $O(dε^{-2})$ is a fairly straightforward application of stochastic subgradient descent, our $O(nd\log^{3}\frac{1}ε)$ time algorithm is a novel long step interior point method. To achieve this running time we start with a simple $O((nd)^{O(1)}\log\frac{1}ε)$ time interior point method and show how to improve it, ultimately building an algorithm that is quite non-standard from the perspective of interior point literature. Our result is one of very few cases we are aware of outperforming traditional interior point theory and the only we are aware of using interior point methods to obtain a nearly linear time algorithm for a canonical optimization problem that traditionally requires superlinear time. We hope our work leads to further improvements in this line of research.

1 Introduction

The paper studies the geometric median problem and develops faster algorithms for computing (1+ε)-approximate solutions, including a nearly linear-time interior point method and a sublinear stochastic method.

Problem and motivation: Given n points in d dimensions, the geometric median minimizes their sum of Euclidean distances and supports facility location, clustering, and robust estimation.Its robustness includes an optimal breakdown point of 0.5, tolerating arbitrary corruption of up to half the input data.
Main results: O(nd log^3(1/ε)) computes a (1+ε)-approximate geometric median, improving the prior running times.The paper also extends the results to the more general Weber’s problem.
Main results: O(dε^-2) provides another (1+ε)-approximation when a random-point sampling oracle is available.This algorithm combines sampling to locate a containing Euclidean ball with stochastic subgradient descent inside that ball.
Algorithmic approach: The nearly linear-time algorithm is a long-step interior point method tailored so that every iteration can be implemented in nearly linear time.Its analysis improves a starting polynomial-time interior point method and uses a non-standard approach relative to interior point literature.
Open questions: The approach leaves open whether a nearly linear-time interior point algorithm should generally be expected for geometric median, and asks for the minimal self-concordance of a relevant barrier.These boundaries reflect limitations of standard self-concordance analysis and the specialized nature of known improvements.

2 Notation

The paper defines the geometric median objective, its approximation criterion, and notation for the penalized objective and central path. It also introduces weights, residual directions, and the minimum-eigenvalue direction used later.

The geometric median minimizes f(x)=∑_i∈[n] ∥a(i)−x∥2 over x∈R^d.
A point is (1+ε)-approximate when f(x)≤(1+ε)f(x∗).
The paper replaces f with penalized objectives f_t and denotes their minimizers by x_t.The central path is {x_t:t>0}.
The weight w_t(x) is defined from the quantities g_t^(i)(x), while u^(i)(x) is the unit direction of x−a(i).
The smallest Hessian eigenvalue is λ_d(∇²f_t(x)), with corresponding eigenvector v_t(x).

3 Properties of the Central Path

The central-path analysis establishes stability properties that support long multiplicative increases in the path parameter. It then connects sufficiently large t to a (1+ε)-approximate geometric median.

The penalized objective is stable inside an O(1/t) Euclidean ball, enabling rapid local minimization near a central-path point.The Hessian changes by at most a constant in this region, so Newton’s method can find the nearby center.
The Hessian and its associated weights change controllably when the path parameter changes multiplicatively.
Except in the smallest-eigenvector direction, the central path moves little over a constant multiplicative increase in t.
At t=2n/(εf(x∗)), solving the penalized problem yields a (1+ε)-approximate geometric median.
The algorithm’s multiplicative schedule therefore requires O(log(1/ε)) iterations.

4 Nearly Linear Time Geometric Median

The nearly linear-time algorithm follows the central path by computing an approximate minimum-eigenvector direction, locally centering, and line searching across multiplicative path-parameter steps. The resulting algorithm outputs an (1+ε)-approximate geometric median in O(nd log^3(1/ε)) time with constant probability.

The algorithm first computes a 2-approximate median, then repeatedly advances to a larger path parameter while tracking the central path.
Each iteration computes an ε_v-approximate minimum eigenvector of the penalized objective’s Hessian.The ApproxMinEig procedure provides both an eigenvalue and eigenvector approximation.
The computed direction is sufficiently aligned with the central path’s bad direction, while orthogonal directions have bounded alignment.
LineSearch uses local centering and one-dimensional convex minimization to find a point close to the next central-path point.
LocalCenter and LineSearch achieve their stated proximity guarantees with high probability under the listed path and initialization conditions.
O(nd log^3(1/ε)) time yields an (1+ε)-approximate geometric median with constant probability.

A Properties of the Central Path (Proofs)

The proofs establish derivative, Hessian, stability, and central-path movement properties for the penalized objective. Together, these results explain why local Newton-style centering and long-step continuation are effective.

The proof develops derivative identities and basic Hessian facts for the penalized objective along the central path.
The function components and Hessian remain stable when points move within a sufficiently small Euclidean neighborhood.
The Hessian admits structural bounds derived from its component matrices and rank-one form.
The central path is almost straight in directions orthogonal to the minimum-eigenvalue direction over a long step.
The proofs relate central-path optimality conditions to the geometric median objective through gradient and path-derivative calculations.

A.6 Simple Lemmas

This section establishes geometric and analytic lemmas used to control the smoothed objective and its Hessian. It also relates objective values to distances from the minimizer.

Technical role: The technical lemmas quantify the inequalities needed to control algorithmic operations and accuracy.These lemmas are introduced to simplify and apply bounds from Section 4.
Distance and objective bounds: For any x, the distance to the minimizer xt satisfies ∥x − xt∥2 ≤ f(x).The proof uses that xt is a convex combination of the input points.
Smoothed-objective bounds: Lemma 20 bounds the smoothed objective and its Hessian for all x ∈ R^d and t > 0.The supplied statement gives the lemma’s scope, while surrounding passages provide the proof steps and resulting inequalities.

B Nearly Linear Time Geometric Median (Proofs)

This section collects proofs, algorithms, and technical lemmas supporting the nearly linear-time geometric-median analysis.

Section purpose: The section provides proofs, algorithms, and technical lemmas from Section 4.It serves as a repository for the results used in the later analysis.

B.1 Eigenvector Computation and Hessian Approximation

This section develops power-method eigenvector computation and uses it to approximate the Hessian of the smoothed objective. The guarantees exploit eigenvalue gaps and stable rank.

Power method: The power method applies repeated matrix-vector multiplication from a Gaussian start to produce a normalized approximate top eigenvector.Algorithm 5 takes a symmetric PSD matrix A and k iterations, then returns u = y/∥y∥2.
Power-method guarantee: In time O(nnz(A) · log(ns/ε)), PowerMethod outputs u with ⟨v1(A), u⟩2 ≥ 1 − ε and u⊤Au ≥ (1 − ε)λ1(A) with high probability in n/ε.The runtime uses the stable rank s rather than the ambient dimension d in the logarithmic factor.
Hessian approximation: ApproxMinEig uses the approximate top eigenvector to obtain Hessian information for the smoothed objective.The Hessian relation is ∇2ft(x) = t2 · wt(x) − A.
Gap-dependent accuracy: When A has a constant multiplicative eigenvalue gap and constant stable rank, the computed vector approximates the top eigenvector with squared correlation at least 1 − ε.The paper identifies g = Ω(1) and s = O(1) as the relevant conditions.
Central-path stability: For t′ ∈ [t, (1 + 1/600)t], the central-path points satisfy ∥xt′ − xt∥2 ≤ 1/(100t).This controls the movement between nearby smoothing parameters.

B.2 Line Searching

This section develops line searching along approximate Hessian directions. It establishes convexity and Lipschitzness of the one-dimensional objective, then derives approximate central-path updates and their costs.

Section overview: The line-search analysis uses centering and one-dimensional minimization as its main tools for approximating the central path.The section explicitly identifies these as the main tools for computing central-path approximations.
Local centering: LocalCenter computes x(k) with high probability in n/ε, using ApproxMinEig and iterative local updates.The runtime analysis assigns O(nd log(n/ε)) to ApproxMinEig and O(nd) per iteration.
Line-search objective: For t > 0, y ∈ R^d, and unit v, gt,y,v is convex and nt-Lipschitz.The function is defined by minimizing ft over a ball translated along the line y + αv.
Search interval: Under the stated proximity and parameter conditions, a minimizer of gt′,y,v reaches ft′(xt′) within α ∈ [−6f(x∗), 6f(x∗)].The conditions include t′ ≤ (1 + 1/600)t and y sufficiently close to xt.
Line-search guarantee: LineSearch outputs x′ satisfying ∥x′ − xt′∥2 ≤ ε/t′ with high probability in n/ε.Its stated cost is O(nd log2(n/(ε˜ε∗))) with O(log(n/(˜ε∗ε))) calls to LocalCenter.
Same-parameter search: A simpler same-parameter line search requires only ∥x − xt∥2 ≤ 1/(100t) initially and returns y with ∥y − xt∥2 ≤ ε/t.The guarantee holds for any vector u.

B.3 Putting It All Together

The paper combines a constant-factor initialization based on geometric-median robustness with algorithms that refine the solution to a (1+ϵ)-approximation. The resulting procedures provide both nearly linear-time and dimension-dependent guarantees.

Main result: O(nd log^3(1/ϵ)) time yields a (1+ϵ)-approximate geometric median with constant probability.Theorem 1 analyzes the algorithm through progressively refined iterates and bounds the final objective using the triangle inequality.
Main result: O(dϵ^-2) time computes a (1+ϵ)-approximate geometric median using sampling and stochastic subgradient descent.The method relies on geometric-median stability to obtain a constant approximation before improving accuracy.
Initialization: The robustness lemma allows up to half the points to be treated as modified while relating the candidate’s distance scale to the geometric median.Its proof compares distances to unmodified points and uses optimality of the geometric median.
Initialization: A constant approximation is obtained by sampling two independent subsets and selecting a candidate using distance percentiles.The procedure CrudeApproximateK computes percentile-based scores for sampled candidates and outputs a selected point together with its scale estimate.
Initialization: The sampled percentile estimates succeed with probability 1−e^-Θ(K) and support the candidate-selection guarantee.Independence of the two samples and a union bound provide simultaneous control over sampled candidates.

C.2 A 1 + ϵ Approximation of Geometric Median

This section improves a constant approximation to a (1+ϵ)-approximation by restricting stochastic subgradient descent to a bounded search domain. The algorithm achieves high-probability accuracy in O(d/ϵ^2) time.

Method: The algorithm samples an index from [n] at each iteration and uses the corresponding point to form the update.The implementation begins by sampling i_k and defining the sampled update quantities.
Guarantee: With probability 1−e^-Θ(1/ϵ), ApproximateMedian returns the claimed approximate geometric median.The theorem states the success probability for the algorithm’s output.
Guarantee: O(d/ϵ^2) time is sufficient for the algorithm to compute a (1+ϵ)-approximate geometric median.The runtime guarantee accompanies the approximation theorem.
Method: The method applies stochastic subgradient descent to the geometric-median objective after computing an initial point and scale parameter.The sampled gradient is a subgradient in expectation, while the initial scale bounds the search domain.

D Derivation of Penalty Function

The paper derives a penalized objective whose minimizers approach the geometric-median solution as the penalty parameter increases. Supporting lemmas control curvature, eigenvectors, line searches, and noisy one-dimensional minimization.

Penalty construction: Increasing the penalty parameter makes minimizers of the penalized objective converge to a solution of the geometric median problem.The construction minimizes auxiliary variables and removes terms that do not affect the minimizing x.
Penalty construction: The Hessian of the penalized objective can be spectrally characterized, enabling high-precision application of its inverse to a vector.A matrix approximation lemma is used to analyze the Hessian structure.
Spectral tools: Spectral similarity bounds how accurately top eigenvectors must be computed when approximating the penalized objective’s curvature.The analysis relates eigenvectors of spectrally similar positive-semidefinite matrices and bounds their alignment error.
Spectral tools: Approximate transitivity of large inner products transfers eigenvector alignment across two successive approximations.If each adjacent pair has squared inner product at least 1−ϵ, the endpoints have squared inner product at least 1−4ϵ.
Line search: A noisy one-dimensional convex minimization routine achieves an additive-error solution using logarithmically many oracle calls.The routine assumes a Lipschitz convex function on an interval and an additive evaluation oracle.

F Weighted Geometric Median

The weighted geometric median extends the objective by assigning positive weights to input points. The paper reduces weighted instances to bounded-sum integer weights and simulates the unweighted algorithms without explicitly expanding repeated points.

Problem: The weighted problem minimizes the weighted sum of Euclidean distances from a point to the input points.This is the Weber-problem formulation used in the section.
Weighted guarantee: O(nd log^3(1/ϵ)) time computes a (1+ϵ)-approximate weighted geometric median.Lemma 35 gives the weighted runtime guarantee in d dimensions for n points.
Reduction: Integer weights with total at most nϵ^-1 represent an equivalent unweighted instance consisting of repeated points.The algorithms simulate this implicit multiset rather than constructing it explicitly.
Algorithm selection: For ϵ>n^-1/2, ApproximateMedian runs in O(dϵ^-2)=O(nd), while smaller ϵ uses AccurateMedian with O(nd log^3(ϵ^-1)) total time.Weighted sampling uses the alias method in the first regime, and implicit repeated-point computations support the second.