The NP-hard problem of computing the maximal sample variance over interval data is solvable in almost linear time with high probability

TL;DR

This paper improves an algorithm for computing the maximal sample variance over interval data, showing it runs in near-linear average time with high probability despite the problem's NP-hardness, and provides probabilistic bounds on worst-case instances.

Contribution

It introduces a more efficient version of an existing algorithm and analyzes its average-case complexity under a probabilistic data model, demonstrating near-linear expected runtime.

Findings

Algorithm's time complexity improved to $O(n ext{log} n + n 2^omega)$

Average complexity is $O(n^{1+epsilon})$ under the probabilistic model

Rare occurrence of instances requiring exponential time, with probability decaying as $e^{-n ext{log} ext{log} n}$

Abstract

We consider the algorithm by Ferson et al. (Reliable computing 11(3), p. 207-233, 2005) designed for solving the NP-hard problem of computing the maximal sample variance over interval data, motivated by robust statistics (in fact, the formulation can be written as a nonconvex quadratic program with a specific structure). First, we propose a new version of the algorithm improving its original time bound $O(n^2 2^\omega)$ to $O(n \log n+n\cdot 2^\omega)$, where $n$ is number of input data and $\omega$ is the clique number in a certain intersection graph. Then we treat input data as random variables as it is usual in statistics) and introduce a natural probabilistic data generating model. We get $2^\omega = O(n^{1/\log\log n})$ and $\omega = O(\log n / \log\log n)$ on average. This results in average computing time $O(n^{1+\epsilon})$ for $\epsilon > 0$ arbitrarily small, which may be considered as "surprisingly good" average time complexity for solving an NP-hard problem. Moreover, we prove the following tail bound on the distribution of computation time: hard instances, forcing the algorithm to compute in time $2^{\Omega(n)}$, occur rarely, with probability tending to zero at the rate $e^{-n\log\log n}$.

Equations

$$\underset{x\in {\mathbb{R}}^n}{\min }V(x)\text{s.t.}\underline{x}\le x\le \overline{x},$$

$$\underset{x\in {\mathbb{R}}^n}{\max }V(x)\text{s.t.}\underline{x}\le x\le \overline{x},$$

$$V(x):=\frac{1}{n}\sum _{i=1}^n{\left(x_i-\frac{1}{n}\sum _{j=1}^n{x}_j\right)}^2.$$

$$O(n\log n+n\cdot 2^{\omega }),$$

$$O(n^{1+ϵ})\text{for all}ϵ>0,$$

$$\mathsf{E}{2}^{\omega }=O(n^{\frac{1}{\log \log n}}),$$

$$\mathrm{Pr}[\omega \ge \delta n]\le e^{-n\log \log n}\text{ for every}\delta >0\text{ and a sufficiently large n}.$$

$$x_i^{\ast }:=\frac{1}{2}({\underline{x}}_i+{\overline{x}}_i),x_i^{\Delta }:=\frac{1}{2}({\overline{x}}_i-{\underline{x}}_i),{\mathbf{x}}_i^{1/n}:=[x_i^{\ast }-\frac{1}{n}{x}_i^{\Delta },x_i^{\ast }+\frac{1}{n}{x}_i^{\Delta }].$$

$$\mathbf{x}=\{x\mid \underline{x}\le x\le \overline{x}\}.$$

$$V(x)-V(x^{\prime })$$

$$\displaystyle\quad=\frac{1}{n}\Bigg{(}\overline{x}_{i}^{2}+\sum_{j\in J}x_{j}^{2}-\underline{x}_{i}^{2}-\sum_{j\in J}x_{j}^{2}$$

$$\displaystyle\qquad\qquad-\frac{1}{n}\bigg{(}\Big{(}\overline{x}_{i}+\sum_{j\in J}x_{j}\Big{)}^{2}-\Big{(}\underline{x}_{i}+\sum_{j\in J}x_{j}\Big{)}^{2}\bigg{)}\Bigg{)}$$

$$\displaystyle\quad=\frac{1}{n}\bigg{(}\overline{x}_{i}^{2}-\underline{x}_{i}^{2}-\frac{1}{n}\Big{(}\overline{x}_{i}^{2}-\underline{x}_{i}^{2}+2(\overline{x}_{i}-\underline{x}_{i})\sum_{j\in J}x_{j}\Big{)}\bigg{)}$$

$$\displaystyle\quad=\frac{1}{n}\bigg{(}4x^{*}_{i}x_{i}^{\Delta}-\frac{1}{n}(4x^{*}_{i}x^{\Delta}_{i})-\Big{(}4x^{\Delta}_{i}\sum_{j\in J}x_{j}\Big{)}\bigg{)}$$

$$\displaystyle\quad=\frac{4}{n}x^{\Delta}_{i}\bigg{(}x^{*}_{i}-\frac{1}{n}x^{*}_{i}-\frac{1}{n}\Big{(}-\overline{x}_{i}+\sum_{j\in\{1,\ldots,n\}}x_{j}\Big{)}\bigg{)}$$

$$<\frac{4}{n}{x}_i^{\Delta }\left(x_i^{\ast }-\frac{1}{n}{x}_i^{\ast }+\frac{1}{n}{x}_i^{\ast }+\frac{1}{n}{x}_i^{\Delta }-x_i^{\ast }-\frac{1}{n}{x}_i^{\Delta }\right)=0.\text{qed}$$

$$[a_0:=\mu [\underline{x}],a_1),(a_1,a_2),\dots ,(a_{m-1},a_m),(a_m,a_{m+1}:=\mu [\overline{x}]].$$

$$V(x^{\prime })=V_1+\frac{1}{n}((x^{\prime }{)}_i^2-x_i^2)-(V_2+\frac{1}{n}(x_i^{\prime }-x_i){)}^2.$$

$${\Phi }^{\ast }(\tilde{z})-{\Phi }^{\ast }(z)\le L(\tilde{z}-z)\text{ whenever }\tilde{z}>z.$$

$$\gamma :=\mathsf{E}[(x_i^{\Delta }{)}^{1+\epsilon }]<\infty .$$

$$O(\mathsf{E}[n\log n+n\cdot 2^{\omega }])$$

$$=O(n\log n+n\cdot n^{\frac{1}{\log \log n}})$$

$$=O(n^{1+ϵ})$$

$${\Phi }^{\ast }(q_j)=\frac{j}{n},j=1,\dots ,n-1.$$

$$p_j^n:=\mathrm{Pr}[{\mathbf{x}}_i^{1/n}\cap I_j\ne \varnothing ];$$

$${\phi }_{x_i^{\ast }}(z)\le L\forall z$$

$$\alpha :=1+2L\left(1+\frac{\gamma }{\epsilon }\right).$$

$$p_j$$

$$=\mathrm{Pr}[x_i^{\ast }\in I_j]+\underset{=:p_j^{-}}{\underbrace{\mathrm{Pr}[x_i^{\ast }+\frac{1}{n}{x}_i^{\Delta }\ge q_{j-1}\wedge x_i^{\ast }<q_{j-1}]}}$$

$$+\underset{=:p_j^{+}}{\underbrace{\mathrm{Pr}[x_i^{\ast }-\frac{1}{n}{x}_i^{\Delta }\le q_j\wedge x_i^{\ast }>q_j]}}.$$

$$\mathrm{Pr}[x_i^{\Delta }\ge z]\le \frac{\gamma }{z^{1+\epsilon }}.$$

$${\phi }_{x_j^{\prime }}(z)\le \frac{L}{n}\forall z.$$

$$p_j^{+}$$

$$=\mathrm{Pr}[x_i^{\Delta }\ge x_j^{\prime }\wedge x_j^{\prime }>0]$$

$$={\int }_0^{\infty }\mathrm{Pr}[x_i^{\Delta }\ge z|x_j^{\prime }=z]\cdot {\phi }_{x_j^{\prime }}(z)\mathrm{d}z$$

$$={\int }_0^{\infty }\mathrm{Pr}[x_i^{\Delta }\ge z]\cdot {\phi }_{x_j^{\prime }}(z)\mathrm{d}z$$

$$={\int }_0^1\mathrm{Pr}[x_i^{\Delta }\ge z]\cdot {\phi }_{x_j^{\prime }}(z)\mathrm{d}z+{\int }_1^{\infty }\mathrm{Pr}[x_i^{\Delta }\ge z]\cdot {\phi }_{x_j^{\prime }}(z)\mathrm{d}z$$

$$\le {\int }_0^1{\phi }_{x_j^{\prime }}(z)\mathrm{d}z+{\int }_1^{\infty }\frac{\gamma }{z^{1+\epsilon }}\cdot {\phi }_{x_j^{\prime }}(z)\mathrm{d}z$$

Full text

\colorlet

MAGENTAmagenta

\contactondrej.sokol@vse.cz

\submittedMay 10, 2022

\titleheadThe maximal variance over interval data

{subject} Quadratic programming, Geometric probability and stochastic geometry

1 Introduction and motivation

1.1 Problem formulation

[Ferson et al. (2005)] studied the pair of optimization problems

[TABLE]

where and $\underline{x}\leq\overline{x}\in\mathbb{Q}^{n}$ are given input data and

[TABLE]

It is obvious that \whole1 is a convex quadratic program (CQP) solvable in polynomial time, while \whole2 is easily proven to be NP-hard. It is worth noting that a general CQP solver yields a weakly polynomial algorithm for \whole1, but [Ferson et al. (2005)] introduced a strongly polynomial method.

They also introduced a method for solving \whole2 which works in exponential time in the worst case (not surprisingly). The method will be described in \whole2. Abbreviating the names of all authors (Ferson, Ginzburg, Kreinovich, Longpré and Aviles), we will refer to their method as FGKLA algorithm.

1.2 Summary of results

In this text we focus on the NP-hard case \whole2 and the FGKLA algorithm. Our contribution is twofold.

Improving the worst-case complexity of the FGKLA algorithm.

Firstly, we show that there exists an implementation of the FGKLA algorithm working in time

[TABLE]

where $\omega$ is the size of the largest clique in a certain intersection graph. The graph will be introduced in \whole2.2. This improves the bound $O(n^{2}\cdot 2^{\omega})$ from the original paper. For further discussion see \whole2.9.

Proving a “good” behavior in a probabilistic setting.

Secondly, we treat the input data $\underline{x},\overline{x}$ as random variables. We introduce a natural and fairly general probabilistic model (details are in \whole3), under which we show that

1. (i)

on average, the algorithm works in time

[TABLE]

which is surprisingly good considering the problem is NP-hard, 2. (ii)

the probability that the algorithm computes in time $2^{\Omega(n)}$ tends to zero faster than exponentially with $n\rightarrow\infty$. In other words, we show that “hard” instances occur indeed rarely.

More specifically: \whole(i) we prove that under the probabilistic model it holds

[TABLE]

where $\mathsf{E}[\cdot]$ stands for the expected value of $[\cdot]$. Combination of \whole5 with \whole3 yields \whole4 as $n\log n=O(n^{1+\epsilon})$ and $n^{\frac{1}{\log\log n}}=O(n^{\epsilon})$ for any $\epsilon>0$. In the entire text, “$\log$” stands for the natural logarithm.

To achieve \whole(ii): from \whole3 it follows that the computing time is exponential when $\omega\geq\delta n$ with $\delta>0$. We prove that

[TABLE]

1.3 Motivation from statistics

Problems \whole1 and \whole2 are studied in statistics; see e.g. [Antoch et al. (2010)] and references therein, and a pseudopolynomial method in [Černý & Hladík (2014)]. The statistical motivation is as follows: we are interested in sample variance $V(x)$ of a dataset $x=(x_{1},\dots,x_{n})^{\mathrm{T}}$. However, the data $x$ is not observable. What is available instead is a collection of intervals $\boldsymbol{x}_{i}\coloneqq[\underline{x}_{i},\overline{x}_{i}]$, $i=1,\dots,n$, such that $\underline{x}_{i}\leq x_{i}\leq\overline{x}_{i}$ (for example, instead of the exact values $x$ we have rounded versions only). Then, $V(x)$ cannot be computed exactly, but we can get tight bounds for $V(x)$ in the form \whole1 and \whole2. In econometrics, this phenomenon is sometimes called partial identification, see [Manski (2003)].

The problem is more general and is studied for various statistics in place of $V(x)$ in \whole1 and \whole2, see the reference books by [Kreinovich et al. (1998)] and [Nguyen et al. (2012)].

1.4 Related work

In general, this paper contributes to the analysis of complexity of optimization problems and algorithms when input data can be assumed to be random, drawn from a particular distribution or a class of distributions. As a prominent example recall the famous average-time analysis of the Simplex Algorithm by [Borgwardt (1982)] and [Spielman & Teng (2004)], where the phenomenon “exponential in the worst case but fast on average” has been studied since 1980’s.

The phenomenon is particularly interesting for NP-hard problems since the worst-case exponential time seems to be unavoidable. In the area of quadratic optimization, the simplex-constrained case has been studied by several authors (see e.g. [Bomze et al. (2017)] and references therein). It turns out that a quadratic form with entries randomly generated from “natural” distributions attains, with a high probability, its global maximum in a face of a small dimension. This property implies that the problem can be solved efficiently by enumerating faces.

Another nice example is the analysis of the average-case complexity of an NP-hard variant of the open shop scheduling problem by [Lu & Posner (1993)]. Their setup is similar to ours: they assume that input data (the job processing times) are generated from a certain class of probabilistic distributions and prove that the average complexity is polynomial in the number of jobs.

Finally, we mention the average-case complexity analysis of the NP-hard $k$-clique problem. [Rossman (2014)] derived bounds on average-case complexity on monotone circuits. [Fountoulakis et al. (2015)] then extended the analysis into a probabilistic setting, showing whether the “hard” instances occur frequently or rarely. Their results are, in a sense, analogous to ours: if the edges are sampled from a “natural” distribution, then $k$-clique can be solved in polynomial time with probability tending to one faster than polynomially.

Aside from the above mentioned general average-case results, the weak average-case complexity paradigm by [Amelunxen & Lotz (2017)] is worth mentioning; here, the word “weak” refers to omitting a “small” subset of hard instances from the set of all possible instances and performing the average-case complexity analysis on such a reduced class of instances. It turns out that for many interesting problems the removal of a subset with exponentially small measure is sufficient to making hard problems tractable.

2 FGKLA Algorithm

Recall that the input instance is given by the pair $\underline{x}=(\underline{x}_{1},\dots,\underline{x}_{n})^{\mathrm{T}}$ and $\overline{x}=(\overline{x}_{1},\dots,\overline{x}_{n})^{\mathrm{T}}$. Compact intervals will be denoted in boldface, e.g. $\boldsymbol{x}_{i}=[\underline{x}_{i},\overline{x}_{i}]$. For $i=1,\dots,n$ define

[TABLE]

The numbers $x_{i}^{*},x_{i}^{\Delta}$ are referred to as center and radius of $\boldsymbol{x}_{i}$, respectively, and $\boldsymbol{x}_{i}^{1/n}$ is called a narrowed interval (i.e., $\boldsymbol{x}_{i}$ shrunk by factor $n$ around its center). For $x\in\mathbb{R}^{n}$ we define $\mu[x]\coloneqq\frac{1}{n}\sum_{i=1}^{n}x_{i}$ (the mean of $x$).

Our version of the FGKLA algorithm is summarized as Algorithm \whole1. The main result of this section is \whole(b). In particular, it improves the worst-case complexity bound $O(n^{2}\cdot 2^{\omega})$ from [Ferson et al. (2005)] (see also \whole2.9). The proof of \whole(b) will be given in \whole2.1.

Theorem 2.1** (properties of the FGKLA algorithm (Algorithm \whole1)).**

• (a)

The FGKLA algorithm correctly solves \whole2*.*

• (b)

Let $G=(V=\{1,\dots,n\},E)$ be an undirected graph where $\{i,j\}\in E$ if and only if $\boldsymbol{x}^{1/n}_{i}\cap\boldsymbol{x}^{1/n}_{j}\neq\emptyset$ (here, $i\neq j$). Let $\omega$ be the size of the largest clique in $G$. Then, FGKLA algorithm works in time $O(n\log n+n\cdot 2^{\omega})$.

Definition 2.2**.**

The graph $G$ from \whole(b) is referred to as FGKLA intersection graph with data $\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{n}$.

2.1 Idea of the FGKLA algorithm

Since the quadratic form $V(x)$ is positive semidefinite, the maximum of \whole2 is attained in a vertex (an extremal point) of the feasible set

[TABLE]

There are $2^{n}$ vertices in total. In a vertex $x$ we have $x_{i}\in\{\underline{x}_{i}=x^{*}_{i}-x^{\Delta}_{i},\ \overline{x}_{i}=x_{i}^{*}+x_{i}^{\Delta}\}$ for every $i$. FGKLA algorithm reduces the number of vertices to be examined from $2^{n}$ to $O(n\cdot 2^{\omega})$. The reduction is based on \whole2.3. A similar lemma was used in the original paper [Ferson et al. (2005)].

Lemma 2.3**.**

Let $x,x^{\prime}\in\boldsymbol{x}$ and let there exist $i\in\{1,\ldots,n\}$ such that

$x_{j}=x^{\prime}_{j}$* for all $j\neq i$ and* 2. 2.

one of the following is satisfied:

1. (a)

$x_{i}=\overline{x}_{i}$, $\mu[x]>x_{i}^{*}+\frac{1}{n}x_{i}^{\Delta}$ and $x^{\prime}_{i}=\underline{x}_{i}$, 2. (b)

$x_{i}=\underline{x}_{i}$, $\mu[x]<x_{i}^{*}-\frac{1}{n}x_{i}^{\Delta}$ and $x^{\prime}_{i}=\overline{x}_{i}$.

Then $V(x)<V(x^{\prime})$.

Proof 2.4**.**

We prove the claim for Assumption 2a, i.e. $x_{i}=\overline{x}_{i},\mu[x]>x_{i}^{*}+\frac{1}{n}x_{i}^{\Delta},x^{\prime}_{i}=\underline{x}_{i}$. The proof is analogous for 2b.

Let $J=\{1,\ldots,n\}\setminus\{i\}$. We want to prove $V(x)-V(x^{\prime})<0$. We have

[TABLE]

Corollary 2.5**.**

Let $x^{\prime}\in\boldsymbol{x}$ be a maximizer and $\mu^{\prime}=\mu[x^{\prime}]$. Let $X$ be the set of all vectors $x\in\boldsymbol{x}$ satisfying:

• (a)

$x_{i}=\underline{x}_{i}$* if $\mu^{\prime}>x_{i}^{*}+\frac{1}{n}x_{i}^{\Delta}$,*

• (b)

$x_{i}=\overline{x}_{i}$* if $\mu^{\prime}<x_{i}^{*}-\frac{1}{n}x_{i}^{\Delta}$, and*

• (c)

$x_{i}\in\{\underline{x}_{i},\overline{x}_{i}\}$* if $\mu^{\prime}\in\boldsymbol{x}_{i}^{1/n}$.*

Then $X$ contains a maximizer.

In cases (a) and (b) we say that variable $x_{i}$ (or index $i$) is determinable with respect to $\mu^{\prime}$; in case (c), variable $x_{i}$ (or index $i$) is free with respect to $\mu^{\prime}$.

Algorithm \whole1 works as follows. It builds the set $A$ (Line \whole1) containing all endpoints $x_{i}^{*}\pm\frac{1}{n}x_{i}^{\Delta}$ of the narrowed intervals $\boldsymbol{x}_{i}^{1/n}$, $i=1,\dots,n$ (here, $A$ acts as a set rather than a list, meaning that possible duplicities are removed), sorts them and denotes them by $a_{1}<\cdots<a_{m}$ (Line \whole2).

Consider the set $\{\mu[x]\mid x\in\boldsymbol{x}\}=[\mu[\underline{x}],\mu[\overline{x}]]$ of all possible means. The endpoints from $A$ divide the interval $[\mu[\underline{x}],\mu[\overline{x}]]$ into at most $2n+1$ regions

[TABLE]

Thanks to \whole2.3 and \whole2.5, every region $(a_{k},a_{k+1})$ contains means $\mu$ with the same set of free indices. For a region $(a_{k},a_{k+1})$, we denote this set by $I(a_{k},a_{k+1})$, i.e. $I(a_{k},a_{k+1})\coloneqq\{i\in\{1,\dots,n\}\mid\boldsymbol{x}_{i}^{1/n}\cap(a_{k},a_{k+1})\neq\emptyset\}$. The set $A$ of endpoints contains the worst possible mean values with respect to the number of free indices. More precisely: for every $a_{k}$, $k=1,\dots,m$, all indices from $I(a_{k-1},a_{k})\cup I(a_{k},a_{k+1})$ are free.

On Line \whole5, we examine the vertex $\overline{x}$. The value of $V(\overline{x})$ is computed and stored as $M$, the maximal value of $V$ found so far. Variables $V_{1}$ and $V_{2}$ will be useful in Line \whole2.

Then, Algorithm \whole1 takes means $a_{1},\ldots,a_{m}$ one by one. For every mean, say $a_{k}$, it takes the set $B_{a_{k}}=\{i\mid x_{i}^{*}-\frac{1}{n}x_{i}^{\Delta}=a_{k}\}$ of indices of narrowed intervals beginning in $a_{k}$ and inserts it to the set $L$ of free indices with respect to $a_{k}$ (Line \whole7). Indices $\{1,\dots,n\}\setminus L$ are determinable with respect to $a_{k}$. This yields $2^{\lvert L\rvert}$ candidate vertices that are examined by Algorithm \whole2, called on Line \whole8.

Then, indices from the set $E_{a_{k}}=\{i\mid x_{i}^{*}+\frac{1}{n}x_{i}^{\Delta}=a_{k}\}$ of narrowed intervals ending in $a_{k}$ are removed from $L$. Intervals with these indices will be fixed to the lower endpoint for every upcoming $k^{\prime}>k$ (Line \whole9 of Algorithm \whole1). The update of $V_{1}$ and $V_{2}$ will be explained later.

Algorithm \whole2 consecutively traverses all $2^{\lvert L\rvert}$ vertices of $\boldsymbol{x}$ resulting from fixing either $x_{i}=\underline{x}_{i}$ or $x_{i}=\overline{x}_{i}$ for the free indices $i\in L$. For every such vertex, say $x$, the variance $V(x)$ is computed. To make these computations cheap, the traversal of $L$ is performed in a way that two successive vertices $x,x^{\prime}$ differ in just one component. Then \whole2.6 shows how to get $V(x^{\prime})$ from $V(x)$ with $O(1)$ arithmetic operations. The variance is stored indirectly as variables $V_{1}$ and $V_{2}$; they can be easily updated when $\underline{x}_{i}$ is switched to $\overline{x}_{i}$, or vice versa.

Lemma 2.6**.**

For $x\in\mathbb{R}^{n}$, we have $V(x)=V_{1}-V_{2}^{2}$, where $V_{1}=\mu[(x_{1}^{2},\ldots,x_{n}^{2})]$ and $V_{2}=\mu[x]$. Furthermore, if $x^{\prime}$ differs from $x$ in just one component, say $i$th, then

[TABLE]

Algorithm \whole2 is an adaptation of the algorithm from Rohn (2006, p. 37) for enumeration of elements of the set $\{\pm 1\}^{\ell}$ for a given $\ell$. The enumeration can in general start from an arbitrary element. The proof of correctness can be found therein. In our variant, the variable $s\in\{\pm 1\}^{\lvert L\rvert}$ indicates the current vertex. In every iteration of while cycle, some $s_{i}$ is set to $-s_{i}$. The $i$th index $L_{i}$ is taken from $L$ (here we consider $L$ as a list rather than a set) and $x_{L_{i}}$ is switched to the other endpoint. For this new vertex, $V_{1}$ and $V_{2}$ are updated (Line \whole8) and the resulting variance $V$ is compared to the best value found so far (Line \whole9).

The following property of Algorithm \whole2 is crucial for the correctness and complexity of FGKLA algorithm: When Algorithm \whole2 ends, then $s=(1,\dots,1)$. The next proposition immediately follows.

Lemma 2.7**.**

Let $V_{1}^{*},V_{2}^{*}$ be the values of the global variables $V_{1},V_{2}$ when Algorithm \whole2 starts and let $V_{1}^{**},V_{2}^{**}$ be the values of $V_{1},V_{2}$ when Algorithm \whole2 ends. Then $V_{1}^{*}=V_{1}^{**}$ and $V_{2}^{*}=V_{2}^{**}$.

In particular, this means that when entering Line \whole8 of Algorithm \whole1, we always start examining the free indices with $x_{i}\in\overline{x}_{i}$ for each $i\in L$.

Finally, Line \whole9 of Algorithm \whole1 removes intervals ending in $a_{k}$ from $L$. These intervals are going to be fixed to their lower endpoints in the following iterations. Since they are at the upper endpoint now, Line \whole9 updates $V_{1}$ and $V_{2}$ accordingly.

Proof 2.8**.**

Proof of \whole(b)

1. a)

Correctness.* Let $x\in\mathbb{R}^{n}$ be a maximizer of \whole2**. Since the maximum is attained in a vertex of the feasible set $\boldsymbol{x}$, we can assume $x_{i}\in\{\underline{x}_{i},\overline{x}_{i}\}$ for all $i$. Moreover, thanks to \whole2.5**, we can assume $x_{i}=\overline{x}_{i}$ for every $i$ such that $\mu[x]<x_{i}^{*}-\frac{1}{n}x_{i}^{\Delta}$ and $x_{i}=\underline{x}_{i}$ for every $i$ such that $\mu[x]>x_{i}^{*}+\frac{1}{n}x_{i}^{\Delta}$. Put all other indices to set $L^{\prime}$, i.e. $L^{\prime}=\{i\in\{1,\ldots,n\}\mid\mu[x]\in\boldsymbol{x}_{i}^{1/n}\}$. Set $k=\arg\max_{k\in\{1,\ldots,m\}}\lvert\mu[x]-a_{k}\rvert$. Consider the set $L$ processed by Algorithm \whole2** in $k$th iteration of Algorithm \whole1. By construction, $L^{\prime}\subseteq L$. Hence, the maximizer $x$ is among the examined vertices.* 2. b)

Complexity.* On Line \whole2**, the algorithm sorts $2n$ numbers with complexity $O(n\log n)$. Algorithm \whole2** is called at most $m$ times, where $m\leq 2n=O(n)$. Recall that $\omega$ is the size of the maximal clique of the FGKLA intersection graph. In the $k$th iteration of the for cycle on Lines \whole6 to \whole10** of Algorithm \whole1** we have $\lvert L\rvert=\lvert\{i\mid a_{k}\in\boldsymbol{x}^{1/n}_{i}\}\rvert$. Thus $\lvert L\rvert\leq\omega$.*

Algorithm \whole2* performs exactly $2^{\lvert L\rvert}$ iterations of the while cycle on Lines \whole2** to \whole11**. Inside its iteration, there is the for cycle on Lines \whole3** to \whole6**. The amortized time complexity of this for cycle is $O(1)$, because in its iteration it either sets some nonzero $z_{i}$ to [math] or stops iterating. Since $z_{i}$ is set to a nonzero value only $2^{\lvert L\rvert}$ times, the overall time of all courses of the for cycle is $O(2^{\lvert L\rvert})$.*

Computing time in the remaining steps is negligible. In particular, note that since $B_{a_{1}},\ldots,B_{a_{m}}$ are pairwise disjoint sets (the same holds true for $E_{a_{1}},\ldots,E_{a_{k}}$), the total number of iterations of for cycles on Lines \whole7* to 9 is at most $n$ during the whole course of FGKLA algorithm.*

The overall complexity is $O(n\log n+n\cdot 2^{\omega})$.∎

Remark 2.9**.**

Aside of the implementation details (which are important for the reduced time complexity bound), our formulation of the algorithm differs from the original paper Ferson et al. (2005) also for another reason. The original formulation can lead to complexity $O(n^{2}\cdot 4^{\omega})$, for example if $\omega=\ell$ and if there are $\ell$ narrowed intervals ending in some $a_{k}$ and further $\ell$ narrowed intervals starting in $a_{k+1}$. However, a minor modification of the original formulation would be sufficient to achieve the time $O(n^{2}\cdot 2^{\omega})$.

3 A probabilistic model

This section is devoted to the main probabilistic result: on average, FGKLA algorithm works in “almost” linear time.

Here we use the statistical motivation of the problem as described in \whole1.3. Namely, in statistics, data are often assumed to form a random sample from a certain distribution. This is exactly our probabilistic model: we assume that both centers of the intervals and their radii form two independent random samples from fairly general classes of distributions.

{assumption}

[the probabilistic model]

1. (A)

The centers $x^{*}_{1},\dots,x^{*}_{n}$ are independent and identically distributed (“i.i.d.”) random variables with a Lipschitz continuous cumulative distribution function (“c.d.f.”) $\Phi^{*}(z)$. That is, there exists a constant $L>0$ such that

[TABLE] 2. (B)

The radii $x^{\Delta}_{1},\dots,x^{\Delta}_{n}$ are i.i.d. nonnegative random variables with a finite moment of order $1+\varepsilon$ for some $\varepsilon>0$. In other words, we assume

[TABLE] 3. (C)

The random variables $x_{i}^{*}$, $x_{i}^{\Delta}$ are independent.

Theorem 3.1**.**

Denote by $\omega$ the size of the largest clique of the FGKLA intersection graph with data $[\underline{x}_{i}\coloneqq x_{i}^{*}-x^{\Delta}_{i},\ \overline{x}_{i}\coloneqq x_{i}^{*}+x^{\Delta}_{i}]_{i=1,\dots,n}$. If $n$ is sufficiently large, then

1. (a)

$\mathsf{E}2^{\omega}\leq 1+n^{\frac{1}{\log\log n}}$* and $\mathsf{E}\omega\leq\frac{3}{2}\left(1+\frac{\log n}{\log\log n}\right)$,* 2. (b)

$\Pr[\omega\geq\delta n]\leq e^{-n\log\log n}$* for any $\delta>0$.*

Remark 3.2**.**

Proof of \whole3.1 will be given in \whole3.1. Statement (b) should be understood more precisely as follows: for every $\delta>0$ there exists $n_{\delta}$ such that $\Pr[\omega\geq\delta n]\leq e^{-n\log\log n}$ if $n\geq n_{\delta}$.

Corollary 3.3** (main result).**

The average computing time is

[TABLE]

for an arbitrarily small $\epsilon>0$. Moreover, the computing time is $2^{\Omega(n)}$ when $\omega$ is linear in $n$ and this event occurs with probability as small as $O(e^{-n\log\log n})$.

Remark 3.4**.**

\whole

B* on the distribution of radii is very mild; indeed, we need just something a little more than existence of the expectation (we even do not need finite variance). On the other hand, Lipschitz continuity of $\Phi^{*}$ (\wholeA) is unavoidable; we will show what can happen without Lipschitz continuity in \whole4. We will also discuss there what happens when we relax the independence assumption (\wholeC) and what is the cost for dependence paid by existence of higher-order moments.*

3.1 Proof of 3.1

Notation.

For a random variable $X$, its probability density function (“p.d.f.”) is denoted by $\varphi_{X}$. Denote by $q_{j}$ the $j$th $n$-quantile of the distribution of centers: i.e., let $q_{1},\dots,q_{n-1}$ satisfy

[TABLE]

Let $I_{1},I_{2},\ldots,I_{n-1},I_{n}$ stand for the intervals $(-\infty,q_{1}],[q_{1},q_{2}],\ldots,$ $[q_{n-2},q_{n-1}],[q_{n-1},\infty)$. Consider the probabilities

[TABLE]

observe that $p^{n}_{j}$ does not depend on $i$ by the i.i.d. assumptions. Probabilities $p^{n}_{1},\dots,p^{n}_{n}$ may differ due to the different shape of $\varphi_{x^{*}_{i}}$ around $q_{1},\ldots,q_{n-1}$. In the upcoming lines, we utilize the fact

[TABLE]

implied by Lipschitz continuity of $\Phi^{*}(z)$ to derive the following upper bound:

Lemma 3.5**.**

For every $j=1,\ldots,n$ we have $p^{n}_{j}\leq\frac{\alpha}{n}$, where

[TABLE]

Remark 3.6**.**

The value of $\alpha$ depends on the properties of distributions of centers and radii. The “hard” cases are those with $\alpha\gg 1$. Indeed, the difficult case is when $\varepsilon$ is close to zero (“radii can be large with a high probability”), $\gamma\gg 0$ (“radii are large on average”) and $L\gg 0$ (“the density of centers can have high peaks”, or “many centers can be close to one another”).

Proof 3.7** (of lemma 3.5).**

Let $n$ and $i$ be fixed; we omit the index $n$ for brevity. For $j=2,\ldots,n-1$, we can decompose $p^{n}_{j}\equiv p_{j}$ as

[TABLE]

Observe that probability $p^{-}_{j}$ vanishes for $j=1$, as well as probability $p^{+}_{j}$ does for $j=n$.

By definition of $I_{j}$, we have $\Pr[x^{*}_{i}\in I_{j}]=\frac{1}{n}$.

For $p^{+}_{j}$, we derive the upper bound 12. Note that the bound does not depend on $j$. Note that the bound holds also true for $p^{-}_{j}$ by symmetry.

By Markov’s inequality and B we get

[TABLE]

Setting $x_{j}^{\prime}:=n(x^{*}_{i}-q_{j})$, from 7 we have

[TABLE]

Now

[TABLE]

Finally,

[TABLE]

Let us introduce indicator variables for all $i,j=1,\dots,n$:

[TABLE]

Note that $W_{ij}$ is alternatively distributed with parameter $p^{n}_{j}$ and that the variables

[TABLE]

are independent.

Hence, the sum of these variables follows binomial distribution, i.e.

[TABLE]

By introducing

[TABLE]

the estimate $p^{n}_{j}\leq\tfrac{\alpha}{n}$ from 3.5 implies

[TABLE]

If $\omega\geq\kappa$ for some $\kappa$, at least $\kappa$ intervals $\boldsymbol{x}^{1/n}_{i}$ have to share a common intersection. This common intersection belongs to one of intervals $I_{1},\ldots,I_{n}$, hence there exists $j$ such that $E_{j}\geq\kappa$. More precisely,

[TABLE]

Fact 14 allows us to use an estimate based on Penrose’s tail bound for binomial distribution.

Lemma 3.8** (Tail bound for the binomial distribution ((Penrose, 2003, p. 16))).**

Let

[TABLE]

If $Z\sim\mathrm{Bi}(n,\pi)$ and $\kappa\geq n\pi$, then

[TABLE]

If $\kappa\geq\alpha,$ we can extend the estimate \whole15 to the form

[TABLE]

In 18 we used the union bound $\Pr[Q_{1}\vee\cdots\vee Q_{n}]\leq\sum_{i=1}^{n}\Pr[Q_{i}]$ for any events $Q_{1},\dots,Q_{n}$. In 19 we used 14 and the tail bound from 3.8.

Let

[TABLE]

If $n$ is sufficiently large, we can estimate

[TABLE]

where $K\coloneqq\log\frac{8\alpha}{c}$. Clearly, $1+n^{\frac{1}{\log\log n}}=O(n^{\epsilon})$ for every $\epsilon>0$. Inequalities \whole22 and \whole23 follow from \whole3.9 showing basic properties of the sequence $u_{\ell}$. Namely, it shows that it decreases exponentially fast. The estimate $\zeta_{n}\leq\frac{1}{4}$ if $n$ is sufficiently large $n$ follows from the observation that $1-c<0$; thus $\zeta_{n}\stackrel{{\scriptstyle n\rightarrow\infty}}{{\longrightarrow}}0$.

Lemma 3.9**.**

Let $k_{n}\geq 4e\alpha$ (here, $e=\exp(1)$).

1. (a)

$\sum_{\ell=\lfloor k_{n}\rfloor+1}^{n}u_{\ell}<2u_{\lfloor k_{n}\rfloor+1}$. 2. (b)

$u_{\lfloor k_{n}\rfloor+1}\leq 2\zeta_{n}$.

Proof 3.10**.**

To prove (a) we show that $u_{\ell}<\frac{1}{2}u_{\ell-1}$ for

[TABLE]

It follows that

[TABLE]

since $\sum_{\ell=\lfloor k_{n}\rfloor+1}^{n}u_{\ell}$ can be bounded by the sum of a geometric sequence with quotient $\frac{1}{2}$. We have

[TABLE]

and

[TABLE]

the last inequality follows from 25. For (b) we use the fact that $H(\xi)=1-\xi+\xi\log\xi$ is a nondecreasing function for $\xi\geq 1$. Thus

[TABLE]

Remark 3.11**.**

Note that a huge $n$ might be needed to achieve $\zeta_{n}\leq\frac{1}{4}$; the particular value depends on $\alpha$. However, if we admit a greater $c$, e.g. $c=8\alpha$, we get $K=0$ and $\zeta_{n}\leq n^{1-8\alpha+8\frac{\alpha}{e}}\approx n^{1-5\alpha}$, which tends to zero fast, so condition $\zeta_{n}\leq\frac{1}{4}$ is not at all restrictive even from the practical viewpoint. On other hand, the exponent in $n^{\frac{c\log 2}{\log\log n}}+1$ becomes a bit worse.

This shows that at the cost of a pair of worse constants, the method behaves well even for small $n$.

In order to complete the proof of \whole3.1(a), we need to estimate $\mathsf{E}\omega$. Using Jensen’s inequality we get

[TABLE]

The proof of \whole3.1(b) is a corollary of the above theory. Indeed, using the notation from \whole17 – \whole20, definition of $H(\xi)$ from \whole16 and \whole3.8, we have

[TABLE]

because the term $(\star)$ is of the order $n\log n$ and dominates all other terms in the limit. The proof of \whole3.1 is complete. ∎

Remark 3.12**.**

The same proof method can be easily generalized to estimate, for example, the probability that the clique is as large as $n^{\eta}$ for a fixed $0\leq\eta\leq 1$ (i.e., this is the event “the computing time exceeds $2^{n^{\eta}}$”). In this case we get $\Pr[\omega\geq n^{\eta}]\leq\exp(-n^{\eta}\log\log n)$.

4 Concluding remarks and comments

4.1 “Unfriendly” distributions for the FGKLA algorithm: Why Lipschitz continuity (A) is unavoidable

We show that if we drop the Lipschitz continuity assumption, we can get $\mathsf{E}\omega\geq\pi n$ for some $\pi>0$ and thus exponential computing time on average (using the fact that average computing time $\geq\mathsf{E}2^{\omega}\geq 2^{\mathsf{E}\omega}=2^{\Omega(n)}$).

Non-continuous distributions.

First consider $\Phi^{*}(z)$, the c.d.f. of $x_{i}^{*}$, with a discontinuity point $z_{0}$. Then $\pi\coloneqq\Pr[x_{i}^{*}=z_{0}]>0$. Setting

[TABLE]

we get $\mathsf{E}U_{i}=\Pr[U_{i}=1]\geq\pi$ and $\omega\geq\sum_{i=1}^{n}U_{i}$ a.s. Thus

[TABLE]

Continuous non-Lipschitz distributions.

We show that the misbehavior of the non-continuous distribution from the previous paragraph can be “simulated” by a non-Lipschitz continuous distribution. Let $z_{0}$ be a discontinuity point of $\Phi^{*}(z)$ from the last paragraph, let $\Phi_{0}\coloneqq\lim_{z\nearrow z_{0}}\Phi^{*}(z)$ and $\eta\coloneqq\Phi^{*}(z_{0}+1)-\Phi_{0}$. Clearly $\eta>0$. Consider another distribution of $x_{i}^{*}$ with c.d.f.

[TABLE]

with $\varepsilon>0$ arbitrarily small. Now $\widetilde{\Phi}^{*}(z)$ is continuous (if there are more discontinuity points of $\Phi^{*}(z)$ outside $[z_{0},z_{0}+1)$, a similar construction can be done in each of them). If $x_{i}^{\Delta}=1$ a.s. and $U_{i}$ has the same meaning as in \whole26, we get

[TABLE]

and thus

[TABLE]

Taking $\varepsilon$ close to zero, we get a clique with average size arbitrarily close to the order $n$.

4.2 The independence assumption (C) is also essential

If we relax the independence assumption, we can get only a weaker estimate on $p^{n}_{j}$ than the bound $p^{n}_{j}=O(n^{-1})$ from \whole3.5. Said informally, we needed $p^{n}_{j}=O(n^{-1})$ in \whole3.8 to satisfy $np^{n}_{j}=O(1)$. Then, since $k_{n}$ grows unboundedly (although slowly), we were able to apply the tail bound for $n$ sufficiently large.

But in the dependent case we can derive only the bound

[TABLE]

resulting in $np^{n}_{j}=O(n^{\frac{1}{2}})$. Then, $k_{n}$ would have to grow faster than $n^{\frac{1}{2}}$ to be able to apply the tail bound and we would get

[TABLE]

or even something worse. Then, the average computation time bound would be as poor as $2^{\sqrt{n}}$. This is a high price for dependence. For specific extremal distributions, the situation can indeed be so bad, as shown in \whole4.3; but for “usual” distributions with enough moments the situation is much better, as explained in \whole4.4.

Let us show \whole27 without the assumption of independence of $x^{*}_{i}$ and $x_{i}^{\Delta}$. By Markov’s inequality we have

[TABLE]

similarly as in \whole10; recall that we have only assumed the existence of a finite moment of order $(1+\varepsilon)$ with value $\gamma$. We have

[TABLE]

as the bounds for $p^{+}_{j}$ (and $p^{-}_{j}$, similarly) can be obtained by setting $x^{\prime}_{j}=n(x_{i}^{*}-q_{j})$ and using 11:

[TABLE]

4.3 An extremal distribution

Unfortunately, the bounds from the previous sections cannot be generally improved. We show an example where \wholeA and B are satisfied, \wholeC is violated and a slightly weaker form of \whole28 holds true — the clique is as large as $n^{\frac{1}{2}-\varepsilon}$ on average, for an arbitrarily small $\varepsilon>0$. Thus we can push the average computation time of FKGLA algorithm arbitrarily close to $2^{\sqrt{n}}$.

Let $x_{1}^{*},\dots,x_{n}^{*}\sim\text{Unif}(0,1)$ independent. Then, clearly, \wholeA is satisfied. Let $0<\varepsilon<1$ (a choice with $\varepsilon$ close to zero is interesting). Define

[TABLE]

\whole

B is satisfied: indeed, the moment of order $1+\varepsilon$ is finite, since

[TABLE]

and $x_{1}^{\Delta},\dots,x_{n}^{\Delta}$ are independent.

For $i=1,\dots,n$ define

[TABLE]

We have

[TABLE]

Obviously, $\omega\geq\sum_{i=1}^{n}U_{i}$ a.s. Thus

[TABLE]

which is close to $n^{\frac{1}{2}}$ if $\varepsilon$ is small.

4.4 FGKLA algorithm can benefit from high-order moments: A trade-off between dependence and existence of such moments

Note that the problem with dependence of the input random variables sketched in \whole4.2 is closely related to the value of $\varepsilon$. Here we show that

if we assume the existence of high-order moments, we can push the bound on $p^{n}_{j}$ close to the “desired” order $O(n^{-1})$ and get good computation time of the FGKLA algorithm even in the dependent case. Indeed, if $\widetilde{\gamma}\coloneqq\mathsf{E}[(x_{i}^{\Delta})^{d}]<\infty$ for some $d$, then Markov’s inequality gives us $\zeta_{n}\coloneqq\Pr[x_{i}^{\Delta}\geq n^{\frac{1}{1+d}}]\leq\widetilde{\gamma}n^{-\frac{d}{1+d}}$. Now we have

[TABLE]

which is close to $n^{-1}$ if $d$ is large.

{acknowledge}

All authors were supported by the Czech Science Foundation pro-ject 22-19353S. O. Sokol was also supported by the Internal Grant Agency (project F4/19/2019) of Faculty of Informatics and Statistics, Prague University of Economics and Business.