A dual Simplex-type algorithm for the smallest enclosing ball of balls and related problems

TL;DR

This paper introduces a dual Simplex-type algorithm for finding the smallest enclosing ball of balls, extending the classical simplex method to second order cone problems and related geometric optimization tasks.

Contribution

It proposes a novel dual algorithm that generalizes the simplex method for second order cone problems, addressing the smallest enclosing ball of balls and related issues.

Findings

Algorithm effectively solves the smallest enclosing ball of balls problem

Extends simplex method to second order cone optimization

Applicable to related geometric problems

Abstract

We define the notion of infimum of a set of points with respect to the second order cone. This problem can be showed to be equivalent to the minimum ball containing a set of balls problem and to the maximum intersecting ball problem, as well as others. We present a dual algorithm which can be viewed as an extension of the simplex method to solve this problem.

Tables

Table 1. Table 1: Computational results with the averages corresponding to 25 datasets with points randomly sampled from a standard normal distribution. † Matlab stopped responding.

Problems		Our dual algorithm			Gurobi (dual)		Gurobi (primal)
			S-pair	Time		Time		Time
$n$	$m$	Iters	updates	(secs)	Iters	(secs)	Iters	(secs)
$10$	${10}^2$	6.44	6.68	0.008	9.72	0.008	8.56	0.228
$10$	${10}^3$	8.12	8.76	0.010	11.56	0.094	9.16	1.627
$10$	${10}^4$	8.12	9.24	0.027	12.44	1.073	out of memory
${10}^2$	${10}^2$	15.56	15.56	0.031	9.56	0.127	8.88	7.455
${10}^2$	${10}^3$	20.96	21.16	0.145	could not solve†		9.36	87.601
${10}^2$	${10}^4$	25.28	25.72	0.375	could not solve†		out of memory
${10}^2$	${10}^5$	29.24	29.88	3.341	could not solve†		out of memory
${10}^3$	${10}^4$	77.32	77.44	8.641	could not solve†		out of memory
${10}^3$	${10}^5$	87.92	88.16	154.613	could not solve†		out of memory

Equations

$$\begin{array}[]{rl}Inf_{\mathcal{Q}}(\mathcal{P}):=\displaystyle\max_{x}&\langle{e_{1}},{x}\rangle\\ \displaystyle\text{ s.t. }&x\preceq_{\mathcal{Q}}p_{i},\,i=1,...,m\end{array}$$

$$\begin{array}[]{cl}\displaystyle\min_{y}&\displaystyle\sum_{i=1}^{m}\langle{p_{i}},{y_{i}}\rangle\\ \text{ s.t. }&\displaystyle\sum_{i=1}^{m}y_{i0}=1\\ &\displaystyle\sum_{i=1}^{m}\overline{y}_{i}=0\\ &y_{i}\succeq_{\mathcal{Q}}0,\quad i=1,...,m.\end{array}$$

$${\overline{y}}_i^{\ast }=\frac{y_{i0}^{\ast }}{p_{i0}-x_0^{\ast }}\left({\overline{x}}^{\ast }-{\overline{p}}_i\right)=y_{i0}^{\ast }\frac{{\overline{x}}^{\ast }-{\overline{p}}_i}{\Vert {\overline{p}}_i-{\overline{x}}^{\ast }\Vert }.$$

$${\overline{x}}^{\ast }\in \mathrm{conv}\left(\{{\overline{p}}_i:\Vert {\overline{p}}_i-{\overline{x}}^{\ast }\Vert =p_{i0}-x_0^{\ast }\}\right).$$

$$0=\sum _{i\in \mathcal{I}}{\overline{y}}_i^{\ast }=\sum _{i\in \mathcal{I}}\frac{y_{i0}^{\ast }}{p_{i0}-x_0^{\ast }}\left({\overline{x}}^{\ast }-{\overline{p}}_i\right)$$

$${\overline{x}}^{\ast }=\sum _{i\in \mathcal{I}}{\alpha }_i{\overline{p}}_i\text{with}{\alpha }_i=\frac{y_{i0}^{\ast }/p_{i0}-x_0^{\ast }}{{\sum }_{i\in \mathcal{I}}{y}_{i0}^{\ast }/p_{i0}-x_0^{\ast }}.$$

$${\overline{x}}^{\ast }=\sum _{i\in \mathcal{I}}{\alpha }_i{\overline{p}}_i,\text{with}{\alpha }_i\ge 0,i\in \mathcal{I},\text{and}\sum _{i\in \mathcal{I}}{\alpha }_i=1.$$

$${\alpha }_i=\frac{y_{i0}/p_{i0}-x_0^{\ast }}{{\sum }_j{y}_{i0}/p_{i0}-x_0^{\ast }},$$

$$\left\{\begin{array}[]{l}\alpha_{1}(z_{1}+\ldots+z_{m})=z_{1}\\ \vdots\\ \alpha_{m}(z_{1}+\ldots+z_{m})=z_{m}\\ \sigma_{1}z_{1}+\ldots+\sigma_{m}z_{m}=1\end{array}\right.$$

$${\sigma }_1{\alpha }_1(z_1+\dots +z_m)+\dots +{\sigma }_m{\alpha }_m(z_1+\dots +z_m)=1$$

$$(z_1+\dots +z_m)=\frac{1}{{\sigma }_1{\alpha }_1+\dots +{\sigma }_m{\alpha }_m}.$$

$$z_i=\frac{{\alpha }_i}{{\sigma }_1{\alpha }_1+\dots +{\sigma }_m{\alpha }_m},$$

$$y_{i0}=\frac{{\alpha }_i(p_{i0}-x_0^{\ast })}{{\sum }_j{\alpha }_j(p_{j0}-x_0^{\ast })}.$$

$${\overline{y}}_i=\frac{y_{i0}}{p_{i0}-x_0^{\ast }}\left(\overline{x}-{\overline{p}}_i\right),$$

$$\overline{x}=\mathrm{ri}\mathrm{conv}(\{{\overline{p}}_i:y_{i0}^{\ast }>0\}).$$

$$\overline{x}\in \mathrm{aff}(\{{\overline{p}}_i:i\in \mathcal{I}\}),\text{for}\mathcal{I}=\{i:\Vert {\overline{p}}_i-\overline{x}\Vert =p_{i0}-x_0\}.$$

$$y_i=0,\text{for }i\notin \mathcal{I},$$

$$y_{i0}=\frac{{\alpha }_i(p_{i0}-x_0)}{{\sum }_j{\alpha }_j(p_{j0}-x_0)}\text{and}{\overline{y}}_i=\frac{y_{i0}}{p_{i0}-x_0}\left(\overline{x}-{\overline{p}}_i\right),\text{for }i\in \mathcal{I},$$

$$x_0^{\ast }=\min \left(p_{10},p_{20},\frac{p_{10}+p_{20}-\Vert {\overline{p}}_1-{\overline{p}}_2\Vert }{2}\right),{\overline{x}}^{\ast }=\frac{(p_{10}-x_0^{\ast }){\overline{p}}_2+(p_{20}-x_0^{\ast }){\overline{p}}_1}{(p_{10}-x_0^{\ast })+(p_{20}-x_0^{\ast })}.$$

$$\Vert {\overline{p}}_i-\overline{x}\Vert =p_{i0}-x_0,\forall p_i\in \mathcal{S};$$

$$\Vert {\overline{p}}_{i_j}-\overline{x}\Vert =p_{j_{i}0}-x_0,$$

$$\overline{x}\in \mathrm{conv}({\overline{\mathcal{S}}}^j\cup \{{\overline{p}}^{\ast }\}).$$

$$\overline{x}\in \mathrm{aff}({\overline{\mathcal{S}}}^j\cup \{{\overline{p}}^{\ast }\})$$

$$\displaystyle M=\left[\begin{array}[]{cccc}\overline{p}_{j_{2}}-\overline{p}_{j_{1}}&\overline{p}_{j_{3}}-\overline{p}_{j_{1}}&...&\overline{p}_{j_{s}}-\overline{p}_{j_{1}}\end{array}\right],$$

$$\displaystyle c=\left(\begin{array}[]{c}p_{{j_{2}}0}-p_{{j_{1}}0}\\ \vdots\\ p_{{j_{s}}0}-p_{{j_{1}}0}\end{array}\right),\,\text{and}\,\,b=\frac{1}{2}\left(\begin{array}[]{c}\left\|{\overline{p}_{j_{2}}}\right\|^{2}-p_{{j_{2}}0}^{2}-\left\|{\overline{p}_{j_{1}}}\right\|^{2}+p_{{j_{1}}0}^{2}\\ \vdots\\ \left\|{\overline{p}_{j_{s}}}\right\|^{2}-p_{{j_{s}}0}^{2}-\left\|{\overline{p}_{j_{1}}}\right\|^{2}+p_{{j_{1}}0}^{2}\end{array}\right).$$

$$M^T\overline{x}=b+x_{0}c,$$

$${\Vert {\overline{p}}_{j_1}-\overline{x}\Vert }^2={\left(p_{j_{1}0}-x_0\right)}^2,$$

$$x_0\le \underset{p_{j_i}\in {\mathcal{S}}^j}{\min }\{p_{j_{i}0}\}.$$

$$M^{+}=(M^{T}M{)}^{-1}{M}^T,$$

$$u=(M^{T}M{)}^{-1}(b-M^T{\overline{p}}_{j_1}),$$

$$v=(M^{T}M{)}^{-1}c,$$

$$w=-M^{+}({\overline{p}}^{\ast }-{\overline{p}}_{j_1}),$$

$$z=(I-M{M}^{+})({\overline{p}}^{\ast }-{\overline{p}}_{j_1}).$$

$$u=v=w=[]\text{and}z={\overline{p}}^{\ast }-{\overline{p}}_{j_1}$$

Full text

Keywords: Smallest enclosing ball, simplex-type methods, second order cone programming, computational geometry

1 Introduction

Let $\mathcal{Q}$ denote the second-order cone, $\mathcal{Q}:=\left\{x:=(x_{0};\overline{x})\in\mathbb{R}^{n}:\,\left\|{\overline{x}}\right\|_{2}\leq x_{0}\right\}$. Given a set of points $\mathcal{P}=\{p_{1},...,p_{m}\}\subset\mathbb{R}^{n}$, we consider the problem defined as follows

[TABLE]

with $e_{1}=(1,0,...,0)$. We define this problem as the infimum of $\mathcal{P}$ with respect to $\mathcal{Q}$ (given its resemblance with the problem of finding the minimum number of a set of numbers with a linear program). Throughout the article we shall refer to this problem as $Inf_{\mathcal{Q}}(\mathcal{P})$, or simply (P).

Denote by $x^{*}$ the optimal solution to (P). One possible geometric interpretation of the problem in question is to find $x^{*}$ as “high” as possible (when height is defined as the value of $x_{0}$) such that $x^{*}+\mathcal{Q}$ covers all points of set $\mathcal{P}$, in the sense that they either fall inside or on the boundary of $x^{*}+\mathcal{Q}$ (Figure 1). But more interestingly, it is possible to prove that (P) is in fact equivalent to several relevant problems in computational geometry, such as the smallest enclosing ball of balls or the largest intersecting ball.

Problem (P) is a Second-Order Cone Program (SOCP) and so it can be solved in polynomial time using interior point methods, such as the primal-dual path following algorithms [20, 21].These algorithms generate interior-points for the primal and dual problems that follow the so-called central path, which converges to a primal-dual optimal solution in the limit. Each iteration, in terms of computational work, basically consists on the solution of a linear system to compute the search direction resulting from the application of Newton’s method to the KKT conditions of a modified (P) with a logarithmic barrier function. The matrix in question is positive definite (see e.g. [1]), and usually its Cholesky factorization is used. In terms of our problem, the matrix has size $n\times n$ (note that the standard form cone LP is the dual), and would take $\mathcal{O}(mn^{2})$ basic arithmetic operations to be computed [25] resulting in a iteration with $\mathcal{O}({mn^{2}}+n^{3})$ computational complexity. Therefore, computing a solution with an error $\epsilon$ with an interior-point method would have an overall complexity of $\mathcal{O}((m+n)n^{2}\sqrt{m}\log{\frac{1}{\epsilon}}))$, though in practice it has been observed that the number of iterations is often very small and independent on $m$. Interior point methods have been used to solve (P) in particular in [14, 25]. Specifically, in [25] an algorithm that computes a $(1+\epsilon)$-approximation in $\mathcal{O}(\frac{mn}{\epsilon}+\frac{1}{\epsilon^{4.5}}\log\frac{1}{\epsilon})$ time is presented.

In this paper we introduce a dual algorithm for (P) that mechanically is analogous to the Simplex method for Linear Programming. Each (major) iteration of our algorithm starts with a primal solution corresponding to a dual feasible basic solution, but then, instead of a line search, our algorithm will perform a sequence of exact curve searches until it arrives to a new dual feasible basic solution with a better objective function value. We will define the notion of basis as well as dual feasible basic solution applied to our problem, which we will rename as support set and dual feasible S-pair, respectively. Mechanically speaking, our algorithm shares similarities with the dual active set algorithm for strictly convex QPs by Goldfarb and Idnani in [10], and with the work of Dearing and Zeck who proposed a dual simplex method for the minimum enclosing ball of points in [5].

Besides providing an exact solution, one advantage that a dual simplex algorithm has over interior point methods is that, after solving the problem, if it suffers small changes (e.g. adding an extra constraint), the dual simplex method will usually require a small number of iterations to calculate the new solution when it starts with the original solution. Another advantage of simplex-type algorithms in general is that they generate basic solutions. For instance in the case of the minimum enclosing ball of balls (and in particular of points too), that will be able to tell us which input balls (points) determine the smallest enclosing ball. On the other hand, we do not have any overall polynomial complexity guarantees for the algorithm.

The paper is organized as follows. In section 2 we show that problem (P) is equivalent to relevant problems in computational geometry involving hyperspheres. Section 3 presents important theoretical background to the algorithm, such as duality results and the definition of support set. The dual simplex algorithm for (P) is then introduced in section 4. We then briefly explain the implementation details in section 5 and show some computational results in section 6.

2 Equivalent geometric problems

We now show that problem (P) is equivalent to several relevant problems in computational geometry such as the smallest enclosing ball of balls. For what follows, denote by $B(c,r)$ a ball with the Euclidean norm with center at $c\in\mathbb{R}^{n-1}$ and radius $r\geq 0$.

The smallest enclosing ball problem. The classical and well studied problem of enclosing a set of Euclidean balls with an Euclidean ball of smallest radius can be reduced to (P) as a consequence of the fact that $B(c_{1},r_{1})\cap B(c_{2},r_{2})\neq\emptyset$ if and only if $\left\|{c_{2}-c_{1}}\right\|\leq r_{2}+r_{1}$. Thus, the problem of enclosing a set $\mathcal{B}$ of balls $B(c_{i},r_{i})$, $i=1,...,m$ with a ball $B(c,r)$ of minimum radius can be solved by (P) considering $\mathcal{P}=\{(-r_{i};c_{i}),\,i=1,...,m\}$. The optimal ball will then be given by the center $c=\overline{x}^{*}$ and radius $r=-x_{0}^{*}$. When $r_{i}=0$ for all $i$, the problem reduces to the smallest enclosing ball of points. Figure 2 illustrates this equivalence: the minimum enclosing ball of set $\mathcal{B}$ is the intersection of the cone $x^{*}+\mathcal{Q}$ with the plane $x_{0}=0$.

The smallest enclosing ball of balls, and in particular of points, is a classical problem in computational geometry. This problem has been extensively studied, specially from a combinatorial point of view, in particular in the LP-type framework, see e.g. [6, 8, 16, 17, 23]. In particular, the minimum enclosing ball of points can easily be converted in a Quadratic Program (QP) and solved using off-the-shelf QP solvers. Gärtner and Schönherr [9] developed a generalization of the simplex method for QP with the goal of targeting geometric QPs, with one of the main applications being the MB problem, while later, Fischer and Gärtner [8] proposed an algorithm with a pivoting scheme resembling the simplex method for LP based on previous ideas from [12]. Using related ideas, Dearing and Zeck [5] developed a dual algorithm for the MB problem. This algorithm was further improved in [4]. Several approximation algorithms have also been developed focusing on finding an $\epsilon$-core set, [3], that is a subset of $\mathcal{S}\subset\mathcal{P}$ that has the property that the smallest ball containing $\mathcal{S}$ once expanded by $1+\epsilon$ covers $\mathcal{P}$. A surprising fact is the existence of an $\epsilon$-core set of size at most $\lceil\frac{1}{\epsilon}\rceil$, independent of the dimension $n$, for any point set $\mathcal{P}\subset\mathbb{R}^{n}$, [14, 2]. Several algorithms focused on finding $\epsilon$-core sets have been proposed [2, 3, 14, 15, 22, 24].

The smallest intersecting ball and the largest enclosed ball problems. Consider now the smallest intersecting ball problem of finding the Euclidean ball with smallest radius that intersects all balls $B(c_{i},r_{i})$, $i=1,...,m$. This problem also reduces to (P) given that $B(c_{1},r_{1})\cap B(c_{2},r_{2})\neq\emptyset$ if and only if $\left\|{c_{2}-c_{1}}\right\|\leq r_{2}+r_{1}$. Thus, considering $\mathcal{P}=\{(r_{i};c_{i}),\,i=1,...,m\}$ in (P), the smallest ball that intersects all balls has center $c=\overline{x}^{*}$ and radius $r=-x_{0}^{*}$. Figure 3 illustrates the equivalence. Balls $B(c_{i},r_{i})$, $i=1,...,m$, may all intersect, and so the smallest intersecting ball could be considered any point in the intersection. In such case, the smallest intersecting ball problem then looses its interest. In fact, when all balls intersect, the solution is such that $x_{0}^{*}>0$. This solution however is not meaningless, in fact, it is not difficult to see that, when $\bigcap B(c_{i},r_{i})\neq\emptyset$, ball $B(\overline{x},x_{0}^{*})$ is the largest radius ball that is enclosed in the intersection of the balls. We shall refer to this problem as the largest enclosed ball problem (see Figure 4). For previous work on this problem see [18, 19] and the references therein.

The smallest intersecting and enclosed ball problem. From the previous two sections, we conclude that the problem of finding a ball with smallest radius that simultaneously encloses balls $B(c_{i},r_{i})$, $i=1,...,m_{1}$, and intersects balls $B(c_{j},r_{j})$, $j=1,...,m_{2}$, can be solved by considering $\mathcal{P}=\{(-r_{i};c_{i}),\,i=1,...,m_{1}\}\cup\{(r_{j};c_{j}),\,j=1,...,m_{2}\}$. The optimal ball will then have center $\overline{x}^{*}$ and radius $-x_{0}^{*}$.

3 Preliminaries

Before we proceed, let us introduce/review some notation.

[TABLE]

3.1 Duality and optimality conditions

It is easy to see that the solution to (P) always exists and that it is unique. Moreover, it is possible to prove the following:

Lemma 1**.**

The solution to $\operatorname{Inf}_{\mathcal{Q}}(\mathcal{P})$ is $p_{k}\in\mathcal{P}$, for some $k=1,...,m$, if and only if $p_{k}\preceq_{\mathcal{Q}}p_{i}$ for all $i=1,...,m$ and $p_{k0}=\min\{p_{i0},\,i=1,...,m\}$.

The dual problem of $\operatorname{Inf}_{\mathcal{Q}}(\mathcal{P})$ is

[TABLE]

The solution to (D) may not be unique.

Lemma 2**.**

Both primal problem (P) and dual problem (D) are strictly feasible; i.e. there exists a primal-feasible vector $x$ such that $x\prec_{\mathcal{Q}}p_{i}$ for all $i=1,...,m$, and there exist dual-feasible $y_{1},...,y_{m}$ such that $y_{i}\succ_{\mathcal{Q}}0$ for all $i=1,...,m$.

Since both primal and dual problems are strictly feasible, the duality gap is zero, that is, strong duality holds and the Karush-Kuhn-Tucker conditions are also sufficient [1, p. 25], as Thoerem 3 states.

Theorem 3** (Optimality conditions).**

Let $x^{*}$ and $y_{i}^{*}$, $i=1,...,m$ be any points in $\mathbb{R}^{n}$. The pair $(x^{*},\{y_{i}^{*}\}_{i=1,...,m})$ is primal-dual optimal if and only if

• •

primal feasibility: $x^{*}\preceq_{\mathcal{Q}}p_{i},\,i=1,...,m$;

• •

dual feasibility: $\sum_{i=1}^{m}y_{i}^{*}=e_{1}$ and $y_{i}^{*}\succeq_{\mathcal{Q}}0,\,i=1,...,m$;

• •

complementary slackness: $\langle{p_{i}-x^{*}},{y_{i}^{*}}\rangle=0,\,i=1,...,m$.

The complementary slackness conditions imply

• •

if $y_{i}^{*}\succ_{\mathcal{Q}}0$ then $x^{*}=p_{i}$ (which can happen for a single $i$, and in that case $y_{i}^{*}=e_{1}$ and $y_{j}^{*}=0$ for all $j\neq i$);

• •

if $x^{*}\succ_{\mathcal{Q}}p_{i}$ then $y_{i}^{*}=0$;

• •

if $p_{i}-x^{*}\in\partial\mathcal{Q}$ and $y_{i}^{*}\in\partial\mathcal{Q}$ then

[TABLE]

A consequence of Theorem 3 is the following characterization of optimality:

Theorem 4**.**

$x^{*}$* is the optimal solution to problem (P) if and only if $x^{*}\preceq_{\mathcal{Q}}p_{i}$ for all $i=1,...,m$, and*

[TABLE]

Proof.

The theorem follows trivially from the optimality conditions in the case when $x^{*}=p_{k}$ for some $p_{k}\in\mathcal{P}$. Consider then, that is not the case.

First, consider $x^{*}$ is optimal. Let $\mathcal{I}=\{i:\,\left\|{\overline{p}_{i}-\overline{x}^{*}}\right\|=p_{i0}-x_{0}^{*}\}$, and $y_{i}^{*}$, $i=1,...,m$ be an optimal dual solution. From complementary slackness and dual feasibility we conclude

[TABLE]

since $y_{i}^{*}=0$ for $i\not\in\mathcal{I}$. From the previous equation we have

[TABLE]

Since $\sum_{i}y_{i0}^{*}=1$ there must be at least one $\alpha_{i}\neq 0$, and $y_{i}^{*}\succeq_{\mathcal{Q}}0$ implies $\alpha_{i}\geq 0$ for all $i$. Because $\sum_{i=1}^{m}\alpha_{i}=1$ we conclude that (2) holds.

Conversely, suppose we have $x^{*}$ that is primal feasible and such that $\overline{x}^{*}$ satisfies (2). Then

[TABLE]

We will now build a dual solution that is feasible and together with $x^{*}$ satisfy complementary slackness. Consider $y_{i0}$, $i=1,...,m$, such that

[TABLE]

Equations (3) together with $\sum_{i}y_{i0}=1$ give the following linear systems of equations

[TABLE]

with $z_{i}=y_{i0}/(p_{i0}-x_{0}^{*})$ and $\sigma_{j}=p_{i0}-x_{0}^{*}$. The last equation implies

[TABLE]

[TABLE]

Therefore

[TABLE]

that is,

[TABLE]

We have that $y_{i0}\geq 0$ for $i\in\mathcal{I}$. Setting $y_{i0}=0$ for $i\not\in\mathcal{I}$, and considering

[TABLE]

for all $i=1,...,m$, we have that $(y_{i0};\,\overline{y}_{i})$ is feasible for the dual problem (D). Since $x^{*}$ is primal feasible and $\left\|{\overline{p}_{i}-\overline{x}^{*}}\right\|=p_{i0}-x_{0}^{*}$ for $i\in\mathcal{I}$, we have that, together with $y_{i}$, complementary slackness is satisfied, thus $x^{*}$ is optimal for (P). ∎

We now present two insightful conclusions from the proof of Theorem 4.

Observation 5*.*

The optimal solution $x^{*}$ to (P) is such that

[TABLE]

Note that $y^{*}_{i0}>0$ cannot be replaced by $y^{*}_{i}\in\partial\mathcal{Q}$ (when $\overline{x}^{*}=\overline{p}_{k}$ we have $y_{k}^{*}=e_{1}\not\in\partial\mathcal{Q}$).

Another conclusion from Theorem 4, which will be at the core of our algorithm, is the following corollary.

Corollary 6**.**

Consider $x$, not necessarily primal feasible, that satisfies

[TABLE]

Let $\alpha_{1},...,\alpha_{m}$ be the coefficients of the affine combination. There exists a dual solution given by

[TABLE]

[TABLE]

that satisfies the dual constraint $\sum_{i=1}^{m}y_{i}=e_{1}$, and, together with $x$, the complementary slackness conditions.

Note that, unless $\{\overline{p}_{i}:\,i\in\mathcal{I}\}$ is affinely independent, the coefficients of the affine combination are not unique and therefore $x$ may correspond to more than one such dual solution. If, additionally, $\overline{x}\in\operatorname{conv}(\{\overline{p}_{i}:\,i\in\mathcal{I}\})$, then $y_{i0}\geq 0$ for all $i$, and so $x$ corresponds to a dual feasible solution.

Finally, we can use the result from Theorem 4 to easily calculate algebraically the solution to $\operatorname{Inf}_{\mathcal{Q}}(\mathcal{P})$ when $\mathcal{P}$ only has two points.

Theorem 7**.**

The solution to $Inf_{\mathcal{Q}}(\{p_{1},p_{2}\})$ is given by

[TABLE]

3.2 The definition of support set and S-pair

We now define support set for the $Inf_{\mathcal{Q}}(\mathcal{P})$ problem the same way a basis is defined for an LP-type problem [7].

Definition 8** (Support set).**

A subset $\mathcal{S}$ is called a support set if no proper subset $\mathcal{S}^{\prime}$ of $\mathcal{S}$ is such that $Inf_{\mathcal{Q}}(\mathcal{S}^{\prime})=Inf_{\mathcal{Q}}(\mathcal{S})$. A subset $\mathcal{S}\subseteq\mathcal{P}$ is said to be an optimal support set if $\mathcal{S}$ is a support set and $Inf_{\mathcal{Q}}(\mathcal{S})=Inf_{\mathcal{Q}}(\mathcal{P})$.

Note that an optimal support set may not be unique.

Definition 9** (Dual feasible S-pair).**

Let $\mathcal{S}\subseteq\mathcal{P}$ and $x$ a vector. We say $(\mathcal{S},x)$ is a dual feasible S-pair if

1. (a)

the points of $\mathcal{S}$ lie on the boundary of $x+\mathcal{Q}$, that is,

[TABLE] 2. (b)

$\overline{x}\in\operatorname{ri}\operatorname{conv}(\overline{\mathcal{S}})$; 3. (c)

$\overline{\mathcal{S}}$ is affinely independent.

Comparing with LP, support set is analogous to the definition of basis, and dual feasible S-pair to the definition of dual basic feasible solution.

Using Theorem 4, it is now possible to prove the following equivalence.

Theorem 10**.**

$(\mathcal{S},x)$* is a dual feasible S-pair iff $\mathcal{S}$ is a support set and $x$ is the solution to $\operatorname{Inf}_{\mathcal{Q}}(\mathcal{S})$.*

As a consequence of the previous theorem, a support set has at least $1$ point and at most $n$ points. This is, in fact, a well known result: the minimum enclosing ball of balls in $\mathbb{R}^{n-1}$ is defined by at most $n$ balls. Another consequence, is the fact that if $(\mathcal{S},x)$ is a dual feasible S-pair and $x$ is primal feasible, then $x$ solves $\operatorname{Inf}_{\mathcal{Q}}(\mathcal{P})$.

Finally, Theorem 11 is a direct consequence of the definition of support set.

Theorem 11**.**

Let $\mathcal{S}$ be a support set. If a point $p^{*}$ is infeasible to $Inf_{\mathcal{Q}}(\mathcal{S})$, then $p^{*}$ belongs to an optimal support set for problem $Inf_{\mathcal{Q}}(\mathcal{S}\cup\{p^{*}\})$.

4 Algorithm description

We now outline the basic framework of the algorithm:

Initialization:

Suppose a dual feasible S-pair $(\mathcal{S}^{0},x^{0})$, $\mathcal{S}^{0}\subseteq\mathcal{P}$, is given.

Loop:

For $j=0,1,2,...$, do:

1. (a)

If $x^{j}$ is primal feasible, stop - $x^{j}$ is the optimal solution of (P). 2. (b)

Else, get a $p^{*}\in\mathcal{P}$ corresponding to a chosen violated constraint. 3. (c)

Obtain a new dual feasible S-pair $(\hat{\mathcal{S}}\cup\{p^{*}\},\hat{x})$, for $\hat{\mathcal{S}}\subseteq\mathcal{S}^{j}$ and $\hat{x}_{0}<x_{0}^{j}$;

Set $(\mathcal{S}^{j+1},x^{j+1})\leftarrow(\hat{\mathcal{S}}\cup\{p^{*}\},\hat{x})$.

Note that, as an initial solution, one can simply pick any point $p\in\mathcal{P}$, and consider the dual feasible S-pair $(\{p\},p)$.

The core of the algorithm is step (c) which will consist on a sequence of curve searches that keep dual feasibility. Assume, for now, that $\overline{\mathcal{S}}^{j}\cup\{\overline{p}^{*}\}$ is affinely independent. Since, at each iteration $j$, we start with a dual feasible S-pair $(\mathcal{S}^{j},x^{j})$, the set of constraints associated with $\mathcal{S}^{j}$ are active at $x^{j}$. As we shall see, the set of points $x$ that correspond to dual feasible solutions and that keep those constraints active define a curve (see section 4.1). A curve search basically consists on the following problem: starting at $x^{j}$, we “move” on the curve in the direction of decrease of $x_{0}$ until either dual feasibility is lost or $p^{*}$ becomes feasible, whichever happens first. As we will see, we will be able to calculate exactly the points where dual feasibility is lost and where $p^{*}$ becomes feasible, so the curve search will be exact. This then boils down to calculating what we will define as the partial step and the full step:

Partial step: the maximum step on the curve without violating dual feasibility. The point corresponding to that step, which we will denote by $x^{partial}$, is the first point on the curve where one of the dual variables, which vary non-linearly as we move on the curve, becomes zero.

Full step: the minimum step on the curve such that the constraint corresponding to $p^{*}$ is feasible, that is, when it becomes active. When it exists, we denote that point by $x^{full}$.

If the full step happens before the partial step, then $x^{full}$ corresponds to a dual feasible variable and so it is the optimal solution to the subproblem $\operatorname{Inf}_{\mathcal{Q}}({\mathcal{S}}^{j}\cup\{p^{*}\})$. In this case, a major iteration is complete, and we go back to Step (a) after setting $(\mathcal{S}^{j+1},x^{j+1})\leftarrow(\mathcal{S}^{j}\cup\{p^{*}\},x^{full})$. Otherwise, if the partial step happens first, then the point (say $p_{k}\in\mathcal{S}^{j}$) corresponding to the dual variable that was about to become infeasible is dropped from $\mathcal{S}^{j}$. A new curve search starting at $x^{partial}$ is then performed. Eventually, after at most $|\mathcal{S}^{j}|$ curve searches where the value of $x_{0}$ either decreases or stays the same, the optimal solution of $\operatorname{Inf}_{\mathcal{Q}}(\mathcal{S}^{\prime}\cup\{p^{*}\})$, for some $\hat{\mathcal{S}}\subseteq\mathcal{S}^{j}$, with a strictly smaller objective function is found.

A curve search is only performed whenever $\overline{\mathcal{S}}^{j}\cup\{\overline{p}^{*}\}$ affinely independent. When that is not the case, a point is immediately dropped from $\mathcal{S}^{j}$ before a curve search is done (this step can actually be seen as a partial step with zero length). Section 4.1 defines the curve and section 4.2 explains the details of the curve search.

4.1 The curve

Consider iteration $j$. Consider

• •

$(\mathcal{S}^{j},x^{j})$, a dual feasible $\mathcal{S}$-pair for $\operatorname{Inf}_{\mathcal{Q}}({\mathcal{S}}^{j})$, with $\mathcal{S}^{j}\subseteq\mathcal{P}$;

• •

$p^{*}\in\mathcal{P}$ corresponding to an infeasible constraint at $x^{j}$;

• •

$\overline{\mathcal{S}}^{j}:=\{\overline{p}:\,p\in\mathcal{S}^{j}\}$ and $s=|\mathcal{S}^{j}|$.

The algorithm restricts the search for the next iterate to the set of primal solutions $x$ where dual feasibility and complementary slackness are maintained, that is

[TABLE]

while it decreases the objective function value $x_{0}$.

In this section we shall see that in general the set of points that satisfy (7) and

[TABLE]

constitute a manifold of dimension $1$, that is, a curve.

By squaring equations (7) it is easy to prove the following Lemma:

Lemma 12**.**

If $|\mathcal{S}^{j}|>1$, define

[TABLE]

Otherwise, when $|\mathcal{S}^{j}|=1$, $M=[\;]$, $b=[\;]$, and $c=[\;]$.

Conditions (7) are equivalent to the following conditions

[TABLE]

For the general case of $|\mathcal{S}^{j}|>1$, let us now define the following matrix and vectors

[TABLE]

Matrix $M^{+}$ is the Moore-Penrose inverse, or pseudo-inverse, of $M$. Since $M$ is full column rank, $M^{+}$ is a left-inverse ($M^{+}M=I$). When $|\mathcal{S}^{j}|=1$, simply consider

[TABLE]

Theorem 13**.**

Conditions (7) and (9) define points $(x_{0};\overline{x})$ such that

[TABLE]

for $x_{0}\leq\min_{p_{j_{i}}\in\mathcal{S}^{j}}\{p_{{j_{i}}0}\}$, $\alpha^{*}\in\mathbb{R}$, and $x_{0}$ and $\alpha^{*}$ such that

[TABLE]

Proof.

From (9) we know that there exist $\alpha_{1},\dots,\alpha_{s},\alpha^{*}$ such that

[TABLE]

that is,

[TABLE]

with $\alpha_{2:s}$ the vector with $\alpha_{2},\dots,\alpha_{s}$. Substituting in (10a), we have

[TABLE]

which yields

[TABLE]

with $u$, $v$, and $w$ as defined above. Combining (13) and (14) we obtain (11). Finally, by plugging (11) in (10b) we obtain (12). ∎

The curve: the affinely independent case

When $\overline{\mathcal{S}}^{j}\cup\{\overline{p}^{*}\}$ is affinely independent then we always have $z\neq 0$, thus the quadratic equation (12) can be solved for $\alpha^{*}$, the coefficient of the affine combination associated with $p^{*}$, obtaining

[TABLE]

We then get the points that satisfy (7,9) as a function of a single variable $x_{0}\in[-\infty,\min_{p_{j_{i}}\in\mathcal{S}^{j}}\{p_{{j_{i}}0}\}]$, thus defining a curve $\Gamma$. This is shown by Corollary 14.

Corollary 14**.**

Consider $\overline{\mathcal{S}}^{j}\cup\{\overline{p}^{*}\}$ affinely independent. Then, conditions (7,9) define a curve parameterized by

[TABLE]

such that

[TABLE]

Geometrically, the curve is either a hyperbola (Figure 5a) or a degenerate hyperbola, see Figure 5b). The curve is also symmetric with respect to a reflection through the hyperplane $\{(x_{0};\overline{x}):\,\overline{x}\in\operatorname{aff}(\overline{\mathcal{S}}^{j}\cup\{\overline{p}^{*}\}),\,x_{0}\in\mathbb{R}\}$ and intersects the boundary of $\operatorname{conv}(\overline{\mathcal{S}}^{j}\cup\{\overline{p}^{*}\})$ at two points (one of them in $\operatorname{conv}(\overline{\mathcal{S}}^{j})$).

When we are interested in the part of the curve that corresponds to a dual feasible solution, we need all coefficients of the affine combination to be non-negative. Therefore the portion of the curve we are interested on corresponds to $\alpha^{*}(x_{0})\geq 0$, that is

[TABLE]

Corollary 15 is a direct consequence of Theorem 13.

Corollary 15**.**

Consider $\overline{\mathcal{S}}^{j}\cup\{\overline{p}^{*}\}$ affinely independent. The points $x$ that satisfy (7,9) and for which $\alpha^{*}(x_{0})\geq 0$ can be written as an affine combination

[TABLE]

for $x_{0}\in[-\infty,\min_{p_{j_{i}}\in\mathcal{S}^{j}}\{p_{{j_{i}}0}\}]$, such that

[TABLE]

with $1_{s-1}$ a $(s-1)-$vector with entries all $1$, and $\alpha_{2:s}$ the vector with $\alpha_{2},\dots,\alpha_{s}$.

Using the formulas from Corollary 6, it is now possible to write explicitly the dual variables as functions of $x_{0}$ corresponding to each point $(x_{0};\Gamma^{+}(x_{0}))$ on the curve.

Case when $\overline{\mathcal{S}}^{j}\cup\{\overline{p}^{*}\}$ is affinely dependent

When $\overline{\mathcal{S}}^{j}\cup\{\overline{p}^{*}\}$ is affinely dependent, the linear system (10a) together with condition (9) define an affine space of dimension $1$, that is, a line, which intersects the manifold defined by (10b-10c), on a single point, $x^{j}$. That is because $x^{j}$ satisfies both (7) and $\overline{x}^{j}\in\operatorname{conv}(\overline{\mathcal{S}}^{j})=\operatorname{conv}(\overline{\mathcal{S}}^{j}\cup\{\overline{p}^{*}\})$ if and only if $x^{j}$ is the solution to $\operatorname{Inf}_{\mathcal{Q}}({\mathcal{S}}^{j})$, which is unique.

Theorem 16**.**

If $\overline{{\mathcal{S}}}^{j}\cup\{\overline{p}^{*}\}$ is affinely dependent then conditions (7,9) have as solution the single point $x^{j}$ that corresponds to a dual feasible solution.

When $\overline{\mathcal{S}}^{j}\cup\{\overline{p}^{*}\}$ is affinely dependent, Theorem 13 still holds with $z=0$. In particular (13) also still holds. But in this case $\alpha_{i}$, $i=1,...,s$, and $\alpha^{*}$ cannot be written as functions of $x_{0}$ as in Corollary 15, but they can be written instead in terms of $\alpha^{*}$. Since, in this case, $x_{0}(\alpha^{*})=x_{0}^{j}$ for any value of $\alpha^{*}$, from (14) we have

[TABLE]

4.2 The curve search

At the beginning of a curve search step, let $\mathcal{S}^{j}$ be the support set and $x^{j}$ the current iterate. If we are starting a major iteration (we come from step (b)) then $(\mathcal{S}^{j},x^{j})$ is a dual feasible S-pair, otherwise, if we come from a previous curve search, it is not. Regardless of the case, at the beginning of a curve search we always have the following

• •

$\overline{\mathcal{S}}^{j}\cup\{\overline{p}^{*}\}$ is affinely independent;

• •

$\overline{x}\in\operatorname{conv}(\overline{\mathcal{S}}^{j}\cup\{\overline{p}^{*}\})$.

This fact will be proved in section 4.4.

The curve search consists on, starting at $x^{j}$, moving on the curve defined by $\Gamma^{+}(x_{0})$ in the direction of decrease of $x_{0}$, until either dual feasibility is lost or the violated constraint becomes feasible. Using the conclusions from the previous section, that corresponds to solving the following problem

[TABLE]

Note that since $\overline{x}^{j}\in\operatorname{conv}(\overline{\mathcal{S}}^{j}\cup\{\overline{p}^{*}\})$ we know that $\alpha_{i}(x_{0}^{j})\geq 0$, for all $i=1,...,s$, and $\alpha^{*}(x_{0}^{j})\geq 0$.

Therefore, the solution to the problem above, is always smaller or equal to $x_{0}^{j}$ which always satisfies the last inequality.

Unless stated otherwise, for the remainder of this section, we shall consider $|\mathcal{S}^{j}|>1$ and $\overline{\mathcal{S}}^{j}\cup\{\overline{p}^{*}\}$ affinely independent.

The partial step

The partial step consists on finding the point corresponding to the maximum step on the curve that can be made without losing dual feasibility. When $\overline{{\mathcal{S}}}^{j}\cup\{\overline{p}^{*}\}$ is affinely independent, Corollary 15 describes the dual variables as we move on the curve. Therefore, the point where dual feasibility is lost corresponds to the value of $x_{0}$ that solves

[TABLE]

Note that, if we consider the definitions of $\alpha_{i}$, $i=1,...,m$ given by (17), the constraint $\alpha^{*}(x_{0})\geq 0$ is already implied. The solution to (19) is then found by simply solving the $s$ equations (17b-17c) subject to the inequality $x_{0}\leq x_{0}^{j}$ and picking the largest of the solutions. These are radical equations that can be easily solved by isolating the square root term in one side, squaring both sides, solving the squared equation, and finally discarding any extraneous solutions. Note that, for some $i$ the corresponding equation may not have a solution, since the curve may not intersect the supporting hyperplane of the corresponding facet of $\operatorname{conv}(\overline{\mathcal{S}}^{j}\cup\{\overline{p}^{*}\})$.

The partial step procedure is summarized in Algorithm 1. Note that this procedure will only be done when $|\mathcal{S}^{j}|>1$, as we shall see ahead.

The full step

The full step procedure finds the point corresponding to the minimum step on the curve for which $p^{*}$ is feasible, that is, the first point on the curve where the primal constraint corresponding to $p^{*}$ is feasible is active:

[TABLE]

The system of equations (10a,10b,20) is equivalent to the system of equations (10a,10b,21), with (21) being

[TABLE]

such that

[TABLE]

and for $x_{0}$ satisfying (10c) and $x_{0}\leq p^{*}_{0}$. Therefore, the full step consists on solving

[TABLE]

Note that the solution to (22) to may not exist. It is easy to see that equations (7,9) together with (20) seek the position of the vertex of a second-order cone that has points $\overline{\mathcal{S}}^{j}\cup\{p^{*}\}$ on its boundary, which is not guaranteed to exist. Additionally, such point may not be unique.

To solve (22), we could plug $\overline{x}=\Gamma^{+}(x_{0})$ in (20), and solve for $x_{0}$. However we realized this approach would result in rather intricate calculations involving square roots. So, instead, we go back to the results of Lemma 12 and Theorem 13 and proceed as we explain next. We start by plugging (11) in (21), and solve for $\alpha^{*}$. We then obtain $\alpha^{*}$ in terms of $x_{0}$:

[TABLE]

Now, we could solve (12) for $x_{0}$. However, the following is easier: we plug (23) back in (11), which now allows $\overline{x}$ to be written solely as a linear function of $x_{0}$ simply as

[TABLE]

with

[TABLE]

and we solve (10b). This results in a much simpler quadratic equation on $x_{0}$:

[TABLE]

Note that, in this process, we never imposed that $\alpha^{*}(x_{0})\geq 0$. If the solution(s) to (24) are real, we only keep the one(s) that satisfy both

[TABLE]

for $\alpha^{*}(x_{0})$ given by (23). If there are more than one such solution, we pick the one with the maximum value of $x_{0}$ (the one that “occurs first” as we move the curve), which will become $x_{0}^{full}$, the solution of (22). If we are left with no solutions or the solutions were not real, then, as explained previously, that means that such point does not exist. In this case, we shall consider $x_{0}^{full}=-\infty$. The details of this procedure are summarized in Algorithm 2.

Taking a step

If $x_{0}^{partial}<x_{0}^{full}$, then the $p^{*}$ becomes feasible while all dual variables associated with $\mathcal{S}^{j}$ are feasible too. That implies that $x^{full}$ is the solution to $\operatorname{Inf}_{\mathcal{Q}}({\mathcal{S}}^{j}\cup\{p^{*}\})$. The algorithm now starts a new major iteration with a dual feasible S-pair $(\mathcal{S}^{j}\cup\{p^{*}\},x^{full})$.

On the other hand, if $x_{0}^{partial}\leq x_{0}^{full}$, then a dual variable became [math] before $p^{*}$ was primal feasible. Geometrically, that means the curve hits a facet of $\operatorname{conv}(\overline{\mathcal{S}}^{j}\cup\{\overline{p}^{*}\})$ before hitting the translated cone $\overline{p}^{*}-\mathcal{Q}$. The point $p_{k}$ that results from applying (26) is an opposite point to that hit facet. When this happens, $p_{k}$ is dropped from $\mathcal{S}^{j}$: $\mathcal{S}^{j}=\mathcal{S}^{j}\setminus\{p_{k}\}$, and the algorithm then performs a new curve search with the new $\mathcal{S}^{j}$ starting at $x^{partial}$. It may happen that there may be more than one dual variable that became [math] simultaneously, that is, at line 6 of Algorithm 1, $x_{0}^{partial}$ belongs to more than one $\mathcal{A}_{i}$. When that happens, after dropping one of such points at the first curve search, the next curve search will consist on a partial step with zero length, that is, the objective function value will be maintained, and a point corresponding to one of those dual variables that are [math] will be picked to be removed from the support set. This is done as many times as necessary.

4.2.1 $\overline{{\mathcal{S}}}^{j}\cup\{\overline{p}^{*}\}$ is affinely dependent

When $\overline{{\mathcal{S}}}^{j}\cup\{\overline{p}^{*}\}$ is affinely dependent, we have seen that conditions (7,9) do not define a curve. That means, a movement in the primal space is not possible. However, it is possible to move in the dual space to different dual solutions corresponding to the same $x^{j}$. Once one the dual variables becomes [math], we drop that point from $\mathcal{S}^{j}$, fixing the issue of the affine dependence. In order to do that, we will solve the following problem

[TABLE]

for $\alpha_{i}(\alpha^{*})$ defined as in (18). This problem finds an affine combination of $x^{j}$ where $\alpha^{*}>0$ and $\alpha_{i}=0$ for some $i=1,...,s$. Problem (25) can easily be solved using a minimum ratio rule, see Algorithm 3.

Let $p_{k}\in\mathcal{S}^{j}$ be the point resulting from (26). In section 4.4 it will be proved that $\overline{\mathcal{S}}^{j}\setminus\{\overline{p}_{k}\}\cup\{\overline{p}^{*}\}$ is affinely independent, and so, a curve search can then be performed.

4.2.2 The special case when $p^{*}$ is the solution and the relation to the $|\mathcal{S}^{j}|=1$ case

Note that the curve search always assumes that the solution to $\operatorname{Inf}_{\mathcal{Q}}({\mathcal{S}}^{j}\cup\{p^{*}\})$ will satisfy (7) for a subset of $\mathcal{S}^{j}$. But that is not the case when $p^{*}$ happens to be the solution. In fact, as we show next, this case fails to be identified by the curve search procedure.

Consider that $p^{*}$ is the solution to $\operatorname{Inf}_{\mathcal{Q}}({\mathcal{S}}^{j}\cup\{p^{*}\})$. After a series of partial steps where a point from $\mathcal{S}^{j}$ is dropped each time, the algorithm finally performs a curve search with $|\mathcal{S}^{j}|=1$, in which, the full step essentially consists on solving the following system of equations

[TABLE]

When $p^{*}$ is the solution to $Inf_{\mathcal{Q}}(\{p_{j_{1}}\}\cup\{p^{*}\})$, the solution to equations (27) does not exist in general (unless it happens that $\left\|{\overline{p}_{1}-\overline{p}^{*}}\right\|=p_{10}-p^{*}_{0}$), and so from the full step we would have $x_{0}^{full}=-\infty$. On the other hand, the partial step consists of solving $\alpha_{1}(x_{0})=0$ given simply by

[TABLE]

(note that $M=[\,\,]$) for $x_{0}\leq p_{{j_{1}}0}$. This yields

[TABLE]

The partial step in this case does not get the correct value of $x_{0}$, which would be $p^{*}_{0}$. That happens because, in the curve search we made the assumption that $\left\|{\overline{p}_{1}-\overline{x}}\right\|=p_{10}-x_{0}$. Note that the result of the partial step is independent of whether $p^{*}$ is the solution or not.

So, when $p^{*}$ is the solution to $Inf_{\mathcal{Q}}(\{p_{j_{1}}\}\cup\{p^{*}\})$, following the partial and full step procedures as explained before, we would have $x_{0}^{partial}<x_{0}^{full}$. As a consequence, a partial step would be taken, point $p_{j_{1}0}$ would correctly be dropped from $\mathcal{S}^{j}$, but the next iterate would have the incorrect value of $x_{0}^{j}$!

This issue conducted us to add an extra step after the Optimality Check where we check whether $p^{*}$ is the solution, that is,

[TABLE]

If $p^{*}$ is not the solution, then a curve search is performed. This not only avoids issues with the partial step when $\mathcal{S}^{j}$ has a single point, but it also avoids potentially $n$ curve search procedures where a point would be removed from the support set in each one, only then to find out that $p^{*}$ is the solution.

This extra step does not fix the fact that, when $\mathcal{S}^{j}$ has a single point, the partial step will not work correctly. However that will never be an issue once it is known that $p^{*}$ is not the solution, because when that is the case, at the full step we will obtain

[TABLE]

as per Theorem 7. The value of $x_{0}^{partial}$ will be as in (28), and so $x_{0}^{full}>x_{0}^{partial}$ as a consequence of the fact that

[TABLE]

Alternatively, and for the sake of simplicity, whenever $|\mathcal{S}^{j}|=1$ one can simply use the formulas given by Theorem 7 instead of doing a curve search.

4.3 Pseudo-code

We now aggregate the results/discussion of the previous sections on Algorithm 4.

4.4 Finiteness and correctness of the algorithm

For the proof of the correctness of the algorithm presented in Algorithm 4, we will follow a similar approach as in [10]. We first introduce the following definition:

Definition 17**.**

The triple $(x,\mathcal{S},p)$ is said to be a (violated) V-triple if the following four conditions hold:

1. (a)

$\overline{\mathcal{S}}\cup\{\overline{p}\}$ is affinely independent, 2. (b)

$\left\|{\overline{p}-\overline{x}}\right\|\geq(p_{0}-x_{0})$, 3. (c)

$\left\|{\overline{p}_{i}-\overline{x}}\right\|=(p_{i0}-x_{0})$, for $p_{i}\in\mathcal{S}$, 4. (d)

$\overline{x}\in\operatorname{conv}(\overline{\mathcal{S}}\cup\{\overline{p}\})$.

We now prove a series of theorems that will culminate in the correctness of the algorithm.

Theorem 18**.**

Given a V-triple $(x^{j},\mathcal{S}^{j},p^{*})$, if the solution given by the full step, $x^{full}$, exists and is dual feasible, then it is the optimal solution to $\operatorname{Inf}_{\mathcal{Q}}(\mathcal{S}^{j}\cup\{p^{*}\})$. Moreover, $x_{0}^{full}\leq x_{0}^{j}$, and, if $x_{0}^{full}>x_{0}^{partial}$ then $(\mathcal{S}^{j}\cup\{p^{*}\},x^{full})$ is a dual feasible S-pair.

Proof.

The statement is a direct consequence of the fact that the full step finds a point $x^{full}$ such that

[TABLE]

whenever the above conditions are feasible. If, additionally, $\overline{x}^{full}\in\operatorname{conv}(\overline{\mathcal{S}}^{j}\cup\{\overline{p}^{*}\})$ then $x^{full}$ solves $\operatorname{Inf}_{\mathcal{Q}}(\mathcal{S}^{j}\cup\{p^{*}\})$ as per Corollary 4. The fact that $x_{0}^{full}\leq x_{0}^{j}$ is a consequence of the full step definition.

Now assume $x_{0}^{full}>x_{0}^{partial}$. It is easy to see that $(\mathcal{S}^{j}\cup\{p^{*}\},x^{full})$ is a dual feasible S-pair. Clearly $\overline{\mathcal{S}}^{j}\cup\{\overline{p}^{*}\}$ is affinely independent. Additionally, $\overline{x}^{full}\in\operatorname{ri}\operatorname{conv}(\overline{\mathcal{S}}^{j})$, because otherwise there would have been a coefficient of the convex combination that would be [math], implying that $x_{0}^{full}=x_{0}^{partial}$, which contradicts the assumption. ∎

From Theorem 18 we conclude that, before the “Optimality Check” procedure, we always have a dual feasible S-pair $(\mathcal{S}^{j},x^{j})$. After picking a $p^{*}$ corresponding to an infeasible constraint, if $\overline{\mathcal{S}}^{j}\cup\{\overline{p}^{*}\}$ is affinely independent, then we have a V-triple. Otherwise, Theorem 20, proves that Algorithm 3 returns a V-triple. Either way, the V-triple at this stage is always such that $\left\|{\overline{p}^{*}-\overline{x}^{j}}\right\|>(p_{0}^{*}-x_{0}^{j})$. Combined with Theorem 19, we conclude that before a curve search we always have a V-triple.

Theorem 19**.**

Suppose $x_{0}^{partial}\geq x_{0}^{full}$. Given a V-triple $(x^{j},\mathcal{S}^{j},p^{*})$, the partial step returns another V-triple $({x}^{partial},\mathcal{S}^{j}\setminus\{p_{k}\},p^{*})$ such that ${x}_{0}^{partial}\leq x_{0}^{j}$.

Proof.

We need to prove properties (a-d) from Definition 17 for $({x}^{partial},\mathcal{S}^{j}\setminus\{p_{k}\},p)$. Property (a) follows from the fact that $\overline{\mathcal{S}}^{j}\cup\{\overline{p}^{*}\}$ is affinely independent. Properties (c) and (d) and the fact that ${x}_{0}^{partial}\leq x_{0}^{j}$ are a direct consequence of the partial step definition. Property (b) follows from a continuity argument since $\left\|{\overline{p}^{*}-\overline{x}^{j}}\right\|\geq(p_{0}^{*}-x_{0}^{j})$ and the fact that $x_{0}^{partial}\geq x_{0}^{full}$, implying that $\left\|{\overline{p}^{*}-\overline{x}^{partial}}\right\|\geq(p_{0}^{*}-x_{0}^{partial})$. ∎

Theorem 20**.**

When $\overline{\mathcal{S}}^{j}\cup\{\overline{p}^{*}\}$ is affinely dependent, $(x^{j},\mathcal{S}^{j}\setminus\{p_{k}\},p^{*})$ for $p_{k}$ satisfying (26) is a V-triple.

Proof.

In order to prove that $(x^{j},\mathcal{S}^{j}\setminus\{p_{k}\},p^{*})$ we need to prove the four properties from Definition 17. Properties (b) and (c) are trivial since $(\overline{\mathcal{S}}^{j},x^{j})$ is an S-pair and $p^{*}$ corresponds to a primal infeasible constraint at $x^{j}$. Now, in order to prove (a), suppose, by contradiction, that $\overline{\mathcal{S}}^{j}\setminus\{\overline{p}_{k}\}\cup\{\overline{p}^{*}\}$ is affinely dependent. Then, since $\overline{\mathcal{S}}^{j}\setminus\{\overline{p}_{k}\}$ is affinely independent, there exists $\beta_{i}$, $i=1,...,s$ and $i\neq k$ s.t.

[TABLE]

With loss of generality consider $k\neq 1$. Since $z=Mw+(\overline{p}-\overline{p}_{1})=0$ (a consequence of the affine dependence of $\overline{\mathcal{S}}^{j}\cup\{\overline{p}^{*}\}$), we have

[TABLE]

for $\delta_{1}=1-\sum_{i}w_{i}$ and $\delta_{i}=w_{i}$ for $i=2,...,s$. As a consequence

[TABLE]

We have that $\sum\frac{\beta_{i}-\delta_{i}}{\delta_{k}}=1$, and note that $\delta_{k}\neq 0$. That is, $\overline{p}_{k}$ is an affine combination of $\overline{\mathcal{S}}\setminus\{\overline{p}_{k}\}$, which contradicts the assumption.

Finally, consider $\beta=\rho-\frac{\rho_{k}}{\sigma_{k}}\sigma$ and $\beta^{*}=\frac{\rho_{k}}{\sigma_{k}}$. It is easy to see that $\beta\geq 0$ and, in particular, $\beta_{k}=0$. We now prove that $\overline{x}^{j}=M\beta_{2:s}+\beta^{*}(\overline{p}-\overline{p}_{1})+\overline{p}_{1}$, that is property (d), by observing that

[TABLE]

since $z=0$. ∎

With the previous theorems, we conclude that, starting from a V-triple $(x^{j},\mathcal{S}^{j},p^{*})$ for which $\left\|{\overline{p}^{*}-\overline{x}^{j}}\right\|>(p_{0}^{*}-x_{0}^{j})$, one can obtain a dual feasible S-pair $(\hat{\mathcal{S}}\cup\{p\},\hat{x})$ such that $\hat{\mathcal{S}}\subseteq\mathcal{S}^{j}$ in at most $|\mathcal{S}^{j}|-|\hat{\mathcal{S}}|\leq n$ partial steps and a full step. Moreover, $\hat{x}_{0}<x_{0}^{j}$, since, even though when taking a partial step or a full step the value of $x_{0}$ may be maintained, we know that the value of $x_{0}$ must decrease either in one of the partial steps taken or in the full step, because otherwise we would have that $\hat{x}=x^{j}$, that is, $Inf_{\mathcal{Q}}(\mathcal{S}^{j})=Inf_{\mathcal{Q}}(\hat{\mathcal{S}}\cup\{p^{*}\})$ contradicting the fact that $p^{*}$ is infeasible to $Inf_{\mathcal{Q}}(\mathcal{S}^{j})$. Therefore, since the value of $x_{0}$ strictly decreases at each major iteration the same S-pair can never reoccur. Since the number of possible S-pairs is finite we conclude:

Theorem 21**.**

The proposed dual algorithm solves problem $Inf_{\mathcal{Q}}(\mathcal{P})$ in a finite number of iterations.

An observation on degeneracy

It is easy to see that in a primal algorithm cycling would be a possibility. The current primal feasible solution may correspond to different support sets, that is, different dual solutions. Thus, it could happen that, after a sequence of adding/dropping points from the support set, the algorithm did not move in the primal space and ended up in a support set visited previously. That is a consequence of the lack of freedom of which point to enter the support set in the primal setting, it has to be one corresponding to a primal constraint that is active. In the dual algorithm, the equivalent situation, of when a movement is not possible because a dual variable is zero at the current iterate, is easily dealt by the algorithm by removing the point corresponding to that dual variable from $\mathcal{S}^{j}$ (either at a partial step or when $\overline{\mathcal{S}}^{j}\cup\{\overline{p}^{*}\}$ is affinely dependent). This is done as many times as the number of dual variables that are zero, after which a movement will then be possible in the next iteration. Therefore, as far as degeneracy is concerned no special procedure is required in our algorithm to prevent cycling.

5 Implementation details

The main computational work that is required in the algorithm happens before each curve search, when three linear systems need solved in order to get vectors $u,v$ and $w$. These linear systems all have the same matrix $M^{T}M$. Our implementation is based on the QR factorization of matrix $M$ of size $n\times(s-1)$ such that $s-1\leq n$, where $s=|\mathcal{S}^{j}|$. We have that $M=QR$, with $Q$ a $n\times n$ orthogonal matrix and $R$ a $n\times(s-1)$ upper triangular matrix. Let

[TABLE]

with $Q_{M}$ with size $n\times(s-1)$ and $Q_{\perp}$ with size $n\times(n-s+1)$, and let, and $R_{\triangle}$ a $(s-1)\times(s-1)$ matrix. As a consequence, $M^{T}M=R^{T}R$.

Using the QR factorization of $M$ allows the following linear systems

[TABLE]

to be reduced to two triangular linear systems, which can efficiently be computed using Back Substitution. Moreover, $M^{+}=R_{\triangle}^{-1}Q_{M}^{T}$, and so

[TABLE]

can be obtained by solving the triangular linear system

[TABLE]

Matrix $M$ is given by the points of the support set $\mathcal{S}^{j}$ which is updated throughout the algorithm either by adding a point or removing a point. Next we explain how to update $M$ in those circumstances.

Adding a point to $\mathcal{S}^{j}$. Whenever a point, $p_{i_{k}}$ is added to $\mathcal{S}^{j}$, we append column $\overline{p}_{i_{k}}-\overline{p}_{i_{1}}$ to matrix $M$.

Removing a point from $\mathcal{S}^{j}$. Whenever a point, $p_{i_{k}}$ is removed from $\mathcal{S}^{j}$, there are two cases. If $i_{k}\geq 1$ then we need to remove the $(k-1)$-th column from $M$. When $i_{k}=1$, we need to remove the first column from $M$, obtaining $\hat{M}$, and get the new $M$ by adding the rank-$1$ matrix $(p_{i_{1}}-p_{i_{k}})^{T}1_{n}$ to $\hat{M}$.

Given the above, every time $M$ needs to be updated, we can use the QR factorization of the old matrix $M$ to calculate efficiently the QR factorization of the new $M$ [11, §12.5]. This is accomplished by using Givens rotations, which, in the case when $M$ is $n\times n$, has $\mathcal{O}(n^{2})$ computational complexity.

6 Computational Results

We have implemented our algorithm in MATLAB as explained in Section 5 choosing the most infeasible point at each iteration. We also compared it with solving both the primal (P) and dual (D) problems with the interior point method solver of Gurobi [13] (version 8.0.0) using its MATLAB interface. All Gurobi parameters were kept at their default values. Our experiments were conducted using MATLAB R2016a (version 9.0) on a Mac with an Intel Core i5 1.6 GHz processor, with 8GB RAM, running Mac OS X version 10.11.6. The results are shown in Table 1.

One reason for the good performance of the dual algorithm is that it does not add many points to the support set that are not in the final support set. This can be seen by the fact that the number of iterations (number of points added to the support set) is very close to the number of dual S-pair updates (number of curve searches performed).