COMPUTING THE SIGNED DISTANCE BETWEENOVERLAPPING ELLIPSOIDS∗SATORU IWATA , YUJI NAKATSUKASA , AND AKIKO TAKEDA †This version corrects an error in Sections 4.2, and revises Section 4.3 accordingly;
the arguments elsewhere are unaffected. The corrections are all colored.
Abstract. Computing the signed distance between two ellipsoids is a convex optimization
problem when the two ellipsoids have no intersection, but becomes non-convex when the ellipsoids
overlap. Efficient algorithms for convex optimization problems are thus not guaranteed to find thecorrect signed distance between overlapping ellipsoids. In this paper, we first show that computing
the signed distance is equivalent to minimizing the norm along the boundary of the Minkowski
difference. We then derive an algorithm with running time O(n6 ), where n is the dimension of theellipsoids, that obtains a global minimizer on the boundary of the Minkowski difference and hence
provides the exact signed distance. The algorithm first finds all the points that satisfy the KarushKuhn-Tucker (KKT) conditions, and then identifies a relevant KKT point with the smallest signed
distance. The primary difficulty in computing the KKT points is that they are the solutions of two
bivariate rational equations, whose poles are not known explicitly. Our key step is to convert therational equations into polynomial equations, which we do by constructing certain bivariate matrix
pencils whose zeros of determinants are the zeros of the rational functions. This reduces the problem
to a two-parameter quadratic eigenvalue problem, which can be solved via a single-parameter linear
eigenvalue problem of larger (squared) size, for which reliable algorithms are available. Thus weprovide the first algorithm for computing the signed distance between overlapping ellipsoids withpolynomial complexity.
Key words. ellipsoids, signed distance, non-convex optimization, KKT conditions, Minkowski
difference, two-parameter eigenvalue problem, Bézoutian
AMS subject classifications. 49M37, 65K05, 90C25, 90C30

1. Introduction. An ellipsoid E lying in an n-dimensional space centered atb ∈ Rn can be represented as the set of points x satisfying(x − b)> A−1 (x − b) ≤ 1,
where A ∈ Rn×n is a symmetric positive definite matrix. We denote such an ellipsoidby E(b, A). The signed distance between two ellipsoids E1 = E(b1 , A1 ) and E2 =E(b2 , A2 ) is defined by (see e.g. [19])(1.1)

dist(E1 , E2 ) = max

kwk=1



hw, b1 i −pphw, A1 wi − hw, b2 i − hw, A2 wi .

Here h·, ·i denotes the inner product, which for simplicity
we take to be hx, yi = x> y,pand k · k designates the norm defined by kxk = hx, xi.
Computing the signed distance between ellipsoids arises in various engineering
problems such as collision detection in graphics, motion planning in robotics, andbinary classification in machine learning [3, 25, 30].
∗ Version of June 6, 2017. This research is supported by Grants-in-Aid for Scientific Research(Challenging Exploring Rearch, No. 26540007) from Japan Society for Promotion of Science.
† Department of Mathematical Informatics, Graduate School of Information Science andTechnology, The University of Tokyo, Tokyo 113-8656, Japan (iwata@mist.i.u-tokyo.ac.jp,
nakatsukasa@mist.i.u-tokyo.ac.jp, takeda@mist.i.u-tokyo.ac.jp).

1

2

SATORU IWATA, YUJI NAKATSUKASA, AND AKIKO TAKEDA

When the two ellipsoids do not overlap, the signed distance dist(E1 , E2 ) coincides
with the perhaps more intuitive formulation of the distance (see Section 2)(1.2)

minimizekx1 − x2 k

subject to x1 ∈ E1 ,x2 ∈ E2 .

The equivalence can be verified by noting that for a fixed vector w, the quantity inside
the parenthesis of (1.1) represents the distance between the two parallel hyperplanes
that are tangent to the ellipsoids, see Figure 1.1 for an illustration. More specifically,
hw, b1 i − hw, b2 i is the distance between two parallelwith normal vector
p hyperplanes p
w that pass through the centers b1 and b2 , and hw, A1 wi and hw, A2 wi are thedistance between the centers b1 , b2 and the hyperplanes tangent to the ellipsoids andhave gradient parallel to w. To see this, note that the gradient at x1 on the boundary of
A−1>
1 (x1 −b1 )E1 is 2A−1
1 (x1 − b1 ), so its normalized vector w = kA−1 (x1 −b1 )k satisfies w (x1 − b1 ) =1phw, A1 wi, which ispthe distance between the two hyperplanes that pass through b1
and x1 . The term hw, A2 wi is analogous. Therefore the maximization problem(1.1) over w corresponds to finding parallel tangent hyperplanes of E1 , E2 with thesmallest distance.
When the ellipsoids overlap, much of the above interpretation remains valid; one
difference is that now dist(E1 , E2 ) is negative, whose absolute value quantifies theamount of overlap. In either case, dist(E1 , E2 ) is the smallest signed distance between
two parallel hyperplanes that are tangent to E1 , E2 such that the two ellipsoids lie on
opposite sides of each hyperplane. The optimal vector w in (1.1) is in fact parallel tothe gradients of both E1 and E2 at the points x1 and x2 .

Fig. 1.1. For non-overlapping ellipsoids E1 and E2 (left), the signed distance dist(E1 , E2 ) is
the distance between the two tangent planes indicated by the dashed lines. For overlapping ellipsoids(right) dist(E1 , E2 ) is negative: the distance between the tangent planes times −1.

The two cases of overlapping and non-overlapping ellipsoids as illustrated in Figure 1.1 lead to nontrivial mathematical differences. First, when the two ellipsoids
overlap the equivalence between (1.1) and (1.2) is lost: the signed distance (1.1) takesa negative value, whereas the minimum value of (1.2) is zero. In this case, the signed
distance provides information about how deeply the two ellipsoids cross with eachother.
The second, and computationally more important, difference is that the problem
of computing the signed distance becomes a non-convex optimization problem withmany local minima. By contrast, in the non-overlapping case, finding the distance canSIGNED DISTANCE BETWEEN ELLIPSOIDS3

be recast as a convex optimization problem and can be computed efficiently by solving
second order cone programming (SOCP). A more direct geometric method has been
proposed by Lin and Han [21], which computes the (standard) distance between two
non-overlapping ellipsoids. This algorithm, which crucially relies on the convexity ofthe problem, is used by the Ellipsoidal Toolbox [19]. Indeed the signed distance is notsupported by the toolbox: it returns the distance zero when the ellipsoids overlap.
Solving the non-convex optimization problem to compute the signed distance inthe overlapping case is nontrivial, and to our knowledge there is no method thatis guaranteed to provide a globally optimal solution. The goal of this paper is todevelop such an algorithm. Our approach is to introduce a dual problem of (1.1) withzero duality gap, and find all the points (x1 , x2 ) that satisfy the Karush-Kuhn-Tucker(KKT) conditions for its relaxation. We then identify a relevant KKT point withthe smallest dual objective value. The Lagrange multipliers corresponding to KKTpoints are characterized as the common zeros of two bivariate rational functions,whose poles are not known explicitly. A key step is to convert the rational equations
to polynomial matrix equations by considering the determinants of certain matrix
pencils. This results in solving a two-parameter quadratic eigenvalue problem, which
can be done via solving single-parameter linear eigenvalue problems of squared size,
for which reliable algorithms exist [16, 18, 23]. The main cost of our algorithm is insolving a 4n2 × 4n2 (generalized) eigenvalue problem, and the overall complexity isO(n6 ) time1 .
Since computing the signed distance is fundamentally more difficult in the overlapping case, it is recommended in practice to check in advance whether the two
ellipsoids overlap or not. This can be done via a quadratically constrained quadratic
programming problem. In the non-overlapping case, a convex optimization solver or amore specified geometric method [21] are more efficient than the method we describe
and hence recommended. The focus of this paper is the non-convex case, where such
efficient algorithms may fail to provide the global solution.
This paper is organized as follows. In Section 2, we introduce an equivalent dual
problem of (1.1). Section 3 derives the KKT conditions for the primal (1.1) and its
relaxed dual problem to give an outline of our algorithm. Section 4 discusses thesolution of bivariate determinantal equations resulting from the KKT conditions andshows that our algorithm works in generic cases with O(n6 ) time. In Section 5, we analyze the degenerate case and show that our algorithm still works by employing certain
preprocessing techniques. In Section 6, we describe how to check the local optimalityof KKT points in the primal problem. Finally, we provide numerical experiments in
Section 7 to illustrate the properties and performance of our algorithm.

2. Dual problem. Computing the signed distance by solving the optimization
problem (1.1) directly via finding all the KKT points appears to be difficult due tothe presence of the square roots in the objective functions. In this section, we show

1 Strictly speaking, the cost of an eigenvalue solver involves the precision  and our algorithm hascomplexity O(n6 log log( 1 )); however, since in most practical situations log log( 1 ) is a small constant(≤ 4 in double precision  ≈ 10−16 ) it is widely accepted that the eigenvalues of an N × N matrix
or matrix pencil can be computed in O(N 3 ) time [14, §7].

4

SATORU IWATA, YUJI NAKATSUKASA, AND AKIKO TAKEDA

by taking the dual that computing the signed distance is equivalent to(2.1)

kx1 − x2 k

minimize
x1 ,x2

subject to x1 ∈ E1 ,x2 ∈ E2 ,x1 − x2 ∈ bd(E1

E2 ).

Here
designates the Minkowski difference and bd(U ) indicates the boundary of aset U . For simplicity, defining U = E1 E2 , we write (2.1) as
minimize kxk

(2.2)

x

subject to x ∈ bd(U ).
Note that U is a convex set because E1 and E2 are convex.
The two ellipsoids are said to overlap if they share a common interior point.
Testing whether E1 and E2 overlap or not can be done by solving the optimization
problem:(2.3)

minimizeηsubject to(x − b1 )> A−1
1 (x − b1 ) ≤ η,
(x − b2 )> A−1
2 (x − b2 ) ≤ η.

The optimal value satisfying η∗ ≥ 1 indicates the non-overlapping case and η∗ < 1indicates the overlapping case. This is a convex quadratic optimization problem, andcan be solved efficiently by an interior-point method [7].
When the ellipsoids E1 , E2 have no overlap, U does not contain the origin in its
interior. In this case, even if bd(U ) is replaced by U in (2.2), the resulting problem(1.2) has an optimal solution on bd(U ). Therefore solving the convex problem (1.2)
is equivalent to computing the signed distance (1.1) when the two ellipsoids do not
overlap. The non-overlapping case of (2.2) includes the special case where the ellipsoids E1 , E2 are touching externally but do not overlap. In such a case, (1.2) has an
unique optimal solution satisfying x1 = x2 .
We mainly focus on the case where the ellipsoids E1 and E2 overlap. In this case,the origin lies inside the convex set U and (1.2) does not have an optimal solution
on bd(U ). Therefore, we need to deal with a non-convex optimization problem; theproblem (2.2) is equivalent to finding the minimum distance from a point in a convex
set U to the complement of U , for which a dual problem with strong duality is
formulated by Briec [8], building upon the work of Nirenberg [26]. Following Briec[8], we introduce the dual problem of (2.2) as(2.4)

maximizew

min x> w
x∈U

subject to kwk = 1.
1/2

1/2

Since E1 = {b1 + A1 u1 : ku1 k ≤ 1} and E2 = {b2 + A2 u2 : ku2 k ≤ 1}, the optimal
objective value of (2.4) is equivalent to


1/2
1/2
(2.5)
max
min
w> (b1 − b2 ) + w> (A1 u1 − A2 u2 ) .
w:kwk=1 ku1 k,ku2 k≤1

Regarding the minimization with respect to u1 and u2 , the first term w> (b1 − b2 )is independent of u1 and u2 , so the minimum is clearly attained when u1 and u2 are

SIGNED DISTANCE BETWEEN ELLIPSOIDS

5

given by1/2

u1 = −

A1 w
1/2

kA1 wk

1/2

,

u2 =

A2 w
1/2

kA2 wk

.

We substitute this into (2.5) to obtain (1.1). Thus (2.2) is the dual problem ofcomputing the signed distance (1.1), and thus so is (2.1). The optimal objective valueof (2.4) is the negative of that of (2.1) and (2.2).
3. Finding the KKT points. Solving the dual problem (2.1) directly is highly
nontrivial. In fact, it is not an easy task to write down the KKT conditions, because
no explicit formula for bd(U ) is available.
We overcome this difficulty by introducing the following relaxation problem:(3.1)

minimize
x1 ,x2kx1 − x2 k

subject to x1 ∈ bd(E1 ),x2 ∈ bd(E2 ).

We find all the KKT points of (3.1), as described in Section 3.1. From among theobtained KKT points (x1 , x2 ), we select those having x1 − x2 on the boundary of U
as we explain in Section 3.2. It will be shown in Section 3.2 that these selected KKTpoints of (3.1) correspond to the KKT points of the primal problem (1.1). We canthen compute the signed distance by the objective values.
3.1. KKT conditions for relaxation. In this section, we present how to findthe KKT points of the relaxation problem (3.1), which can be recast as
minimize
x1 ,x2kx1 − x2 k2

subject to(3.2)(x1 − b1 )> A−1
1 (x1 − b1 ) = 1,

(3.3)

(x2 − b2 )> A−1
2 (x2 − b2 ) = 1.

The KKT conditions for the problem consist of these feasibility constraints and(3.4)

x1 − x2 = λA−1
1 (x1 − b1 ),(3.5)

x2 − x1 = µA−1
2 (x2 − b2 ),where λ and µ are Lagrange multipliers. These conditions imply(3.6)µ(x2 − b2 ) = −λA2 A1 −1 (x1 − b1 ),(3.7)

λ(x1 − b1 ) = −µA1 A2 −1 (x2 − b2 ).

Combining (3.4) and (3.6), we obtainµx1 − µb2 = λµA1 −1 (x1 − b1 ) − λA2 A1 −1 (x1 − b1 ).
Similarly, from (3.5) and (3.7), we obtain
λx2 − λb1 = λµA2 −1 (x2 − b2 ) − µA1 A2 −1 (x2 − b2 ).
Defining
N (λ, µ) = λµI − µA1 − λA2 ,6

SATORU IWATA, YUJI NAKATSUKASA, AND AKIKO TAKEDA

we can write the two equalities as(3.8)µ(b1 − b2 ) = N (λ, µ)A1 −1 (x1 − b1 ),(3.9)

λ(b2 − b1 ) = N (λ, µ)A2 −1 (x2 − b2 ).

The two equations (3.8) and (3.9) can be solved for x1 , x2 to express them in terms of
λ, µ, which can be substituted into (3.2), (3.3) to yield two bivariate rational equations
with respect to λ, µ. However, solving them for λ, µ is difficult as the poles of therational equations are unknown.
Instead, we formulate a pair of matrix equations that provide the appropriate
Lagrange multipliers: we introduce a pair of matrices M1 (λ, µ) and M2 (λ, µ) defined
by




µA1
N (λ, µ)λA2
N (λ, µ)(3.10)M1 (λ, µ) =, M2 (λ, µ) =,N (λ, µ)µD
N (λ, µ)λD
where D = dd> , d := b1 − b2 . We show in Lemma 3.2 that the Lagrange multipliersmust satisfy the bivariate determinantal equations(3.11)

det(M1 (λ, µ)) = det(M2 (λ, µ)) = 0.

To prove Lemma 3.2, we utilize the following lemma, which will be used repeatedly
throughout this paper.
Lemma 3.1. Let X be a real matrix in the formΛ
Ψ
(3.12)X=,Ψ dd>
where d is a vector, Λ, Ψ are n × n diagonal, Λ is nonsingular and Ψ is singular. Then
X is nonsingular if and only if there is no null vector of Ψ that is orthogonal to d.
Proof. Defining


Λ Ψ
O
d ,
X̂ =  Ψ O
O d> −1
we have det X̂ = − det X, so it suffices to examine when X̂ is singular. Recallingthat Λ and Ψ are both diagonal with Λ nonsingular, by the determinant expansion
we easily see that X̂ is singular if Ψii = di = 0 for some i. Similarly, X̂ is singular if Ψ
has two or more diagonal elements equal to zero. It remains to treat the case where Ψ
has exactly one zero diagonal element Ψii = 0 with the corresponding element di 6= 0.
In this case there is only one nonzero term in the Leibniz determinant formula, whichis indeed nonzero, hence det X̂ 6= 0.
Lemma 3.2. For every pair of x1 ∈ bd(E1 ) and x2 ∈ bd(E2 ) that satisfies the
KKT conditions with multipliers λ and µ, we have det M1 (λ, µ) = det M2 (λ, µ) = 0.
Proof. By (3.2)–(3.5), if one of λ or µ is zero, then so is the other. In this case,M1 (λ, µ) = M2 (λ, µ) = O, the zero matrix. Hence we may assume that neither λ norµ is zero.
By the conditions (3.8) and (3.9), we see that d = b1 − b2 must belong to
Im N (λ, µ). Since A1 , A2 are positive definite, the pair (A1 , N (λ, µ)) can be simultaneously diagonalized by congruence, so Lemma 3.1 implies that M1 (λ, µ) is singular.
The same argument shows that M2 (λ, µ) is also singular.

SIGNED DISTANCE BETWEEN ELLIPSOIDS7

Now suppose that N (λ, µ) is nonsingular. Note that




µA1
N (λ, µ) O
λA2N (λ, µ) OOµd  , M20 (λ, µ) = N (λ, µ)O
λd M10 (λ, µ) = N (λ, µ)
>
>
O
µd
−µ
O
λd
−λ
satisfy det M10 (λ, µ) = −µ det M1 (λ, µ) and det M20 (λ, µ) = −λ det M2 (λ, µ). Since

−1 

µA1
N (λ, µ)O
N (λ, µ)−1=,
N (λ, µ)O
N (λ, µ)−1 −µN (λ, µ)−1 A1 N (λ, µ)−1
we also havedet M10 (λ, µ) = (−1)n det N (λ, µ)2 {µ3 d> N (λ, µ)−1 A1 N (λ, µ)−1 d − µ}.
Thus we obtaindet M1 (λ, µ) = (−1)n−1 det N (λ, µ)2 {µ2 d> N (λ, µ)−1 A1 N (λ, µ)−1 d − 1}
= (−1)n−1 det N (λ, µ)2 {(x1 − b1 )> A1 −1 (x1 − b1 ) − 1},
where the last equality follows from (3.8). Similarly, using (3.9), we obtaindet M2 (λ, µ) = (−1)n−1 det N (λ, µ)2 {λ2 d> N (λ, µ)−1 A2 N (λ, µ)−1 − 1}
= (−1)n−1 det N (λ, µ)2 {(x2 − b2 )> A2 −1 (x2 − b2 ) − 1}.
It then follows from (3.2) and (3.3) that det M1 (λ, µ) = 0 and det M2 (λ, µ) = 0.
Lemma 3.2 suggests computing all possible pairs of Lagrange multipliers λ and µ
for the KKT points by solving the bivariate determinantal equations (3.11). We will
discuss how to solve (3.11) in Section 4. In practice, the computed pairs (λ, µ) may
contain nonreal pairs and those at infinity. These can immediately be removed fromconsideration because (λ, µ) must be real at the KKT points.
For each pair of λ and µ thus obtained, one can compute x1 and x2 by solving
the system of linear equations (3.8) and (3.9). If N (λ, µ) is nonsingular, then x1 andx2 are uniquely determined, and they naturally satisfy the feasibility conditions (3.2)
and (3.3). If N (λ, µ) is singular, among all the solutions of (3.8) and (3.9), we selectthose that satisfy (3.2) and (3.3), and verify that (3.4) and (3.5) hold. This uniquelydetermines x1 , x2 if the null space of N (λ, µ) is 1-dimensional. If dim Ker N (λ, µ) > 1,
then we have infinitely many solutions of x1 , x2 satisfying (3.2) and (3.3). Specifically,
denoting by N0 the null space of N (λ, µ), they can be written as x1 = x1,∗ +N0 v1 , x2 =x2,∗ + N0 v2 , where x1,∗ , x2,∗ are any vectors satisfying (3.8) and (3.9); for example
the least-squares solution. We then solve (3.2) and (3.3) for v1 , v2 . More specifically,
substituting the form of x1 into (3.2) we obtain(3.13)(x1,∗ + N0 v1 − b1 )> A−1
1 (x1,∗ + N0 v1 − b1 ) = 1

with unknown v1 . We can solve this for example as follows: first solve the convex
unconstrained optimization problem minimizev1 (x1,∗ + N0 v1 − b1 )> A−1
1 (x1,∗ + N0 v1 −b1 ), and check if the optimal value is below or above 1. If it is above 1 then there is no
solution for (3.13). If it is below 1, then we can find v1 by v1 = v1,∗ + αw where α isan arbitrary vector. We then solve (3.13) for α, which must have a real solution. We
treat v2 analogously. Note that such solutions x1 , x2 may not be unique. Fortunately,
it can be shown that the corresponding objective function values are the same forevery solution, so it suffices to compute just one of them.

8

SATORU IWATA, YUJI NAKATSUKASA, AND AKIKO TAKEDA

3.2. KKT points for the primal problem. The KKT conditions for the
primal problem (1.1) consist of(3.14)−(b1 − b2 ) + p

A2 w
+p
= ϕw,
w > A1 w
w > A2 w
A1 w

and the feasibility constraint w> w = 1, where ϕ is a Lagrange multiplier. Then theobjective value of (1.1) at a KKT point w coincides with
(3.15)

w> (b1 − b2 ) −p

w > A1 w −p

w> A2 w = −ϕ.

Given a KKT point of (1.1), one can obtain a KKT point of (3.1) as follows.
Lemma 3.3. Suppose that w is a KKT point of (1.1) with a Lagrange multiplierϕ. Then defining(3.16)(3.17)

A1 w,
x1 := b1 − p
w > A1 w
A2 w
x2 := b2 + p,
w > A2 wp
λ := ϕ w> A1 w,
pµ := ϕ w> A2 w,

we obtain a KKT point (x1 , x2 ) of (3.1) with Lagrange multipliers λ, µ.
Proof. It is easy to see from (3.16) and (3.17) that x1 and x2 satisfy the feasibility
constraints (3.2) and (3.3), respectively. By (3.16), we have
λA−1
1 (x1 − b1 ) = − pλww > A1 w= −ϕw.

Similarly, by (3.17), we also have
µw
µA−1
= ϕw.
2 (x2 − b2 ) = p >w A2 w
It follows from (3.14), (3.16) and (3.17) that x1 − x2 = −ϕw holds. Hence x1 and x2
satisfy (3.4) and (3.5).
Note that the Lagrange multipliers λ and µ thus obtained must both have thesame sign as ϕ. Intuitively, such a KKT point (x1 , x2 ) with Lagrange multipliers
having the same sign places x1 − x2 on the boundary of E1 E2 , because at such−1
points the gradients 2A−1
1 (x1 − b1 ) of (3.2) and 2A2 (x2 − b2 ) of (3.3) point towards
exactly the opposite directions.
We now intend to claim conversely that any KKT point of (3.1) with nonzero
Lagrange multipliers with the same signs leads to a KKT point of (1.1).
Lemma 3.4. Suppose (x1 , x2 ) is a KKT point of the relaxed dual problem (3.1)
with nonzero Lagrange multipliers λ, µ having the same sign, say σ ∈ {+, −}. Thendefining(3.18)

w := −

σ(x1 − x2 ),kx1 − x2 k

ϕ := σkx1 − x2 k,

we obtain a KKT point w of (1.1) with a Lagrange multiplier ϕ.

SIGNED DISTANCE BETWEEN ELLIPSOIDS

9

Proof. The definitions (3.18) of w and ϕ imply that x1 − x2 = −ϕw. By (3.4)
and (3.2), we haveσA1 (x1 − x2 )
A w
p 1
= −p
(x1 − x2 )> A1 (x1 − x2 )w > A1 w
σλ(x1 − b1 )=− q
= −(x1 − b1 ).
|λ| (x1 − b1 )> A−1
(x
−
b
)
1
1
1
Similarly, by (3.5) and (3.3), we also haveA2 w
pw > A2 w=p

σA2 (x2 − x1 )(x2 − x1 )> A2 (x2 − x1 )σµ(x2 − b2 )q=
= x 2 − b2 .
|µ| (x2 − b2 )> A−1
(x
−
b
)
2
2
2

Combining these equalities, one can easily verify that w and ϕ satisfy (3.14).
It will be demonstrated later in Section 7 that (3.1) admits KKT points with
Lagrange multipliers of different signs. Such KKT points correspond to interior pointsof E1 E2 , which are infeasible in the original dual problem (2.1).
For a KKT point (x1 , x2 ) with Lagrange multipliers λ, µ, if one of λ or µ is
zero, then so is the other by (3.6) and (3.7). In this case, x1 = x2 must hold.
Conversely, any point x ∈ bd(E1 ) ∩ bd(E2 ) provides a KKT point (x1 , x2 ) := (x, x)
with Lagrange multipliers λ, µ := 0. Thus, in the overlapping case, the relaxed dual
problem (3.1) has infinitely many KKT points with zero Lagrange multipliers, butthey are irrelevant to the signed distance problem. Fortunately, such KKT points are
removed automatically in the algorithm we describe in Section 4.1.
If (3.1) has a KKT point (x1 , x2 ) with Lagrange multipliers λ, µ both negative,then the objective value at the corresponding KKT point w of (3.1) must be positive,
which means the signed distance is positive. This can happen only if the two ellipsoids
are disjoint.
Conversely, since the signed distance between two non-overlapping ellipsoids must
be nonnegative, (3.15) and Lemma 3.3 imply that there is a KKT point (x1 , x2 ) for(3.1) with Lagrange multipliers λ, µ ≤ 0. Note that (x1 , x2 ) remains to be a KKT pointwith the same Lagrange multipliers if we replace the feasibility constraints by x1 ∈ E1
and x2 ∈ E2 (including the interior). The resulting convex optimization problem (1.2)
has a unique KKT point, which is optimal. This establishes the uniqueness of the
KKT point for (3.1) with nonpositive Lagrange multipliers with the same signs in the
non-overlapping case, and the signed distance is obtained by dist(E1 , E2 ) := kx1 −x2 k.
3.3. Outline of the algorithm. For computing the signed distance between
overlapping ellipsoids, Lemmas 3.3 and 3.4 together with (3.15) suggest finding all
the KKT points (x1 , x2 ) for (3.1) with positive Lagrange multipliers and choosing theone that attains the minimum value of kx1 − x2 k among them. In fact, Lemma 3.4
and (3.15) ensure that the selected (x1 , x2 ) provides a KKT point for (1.1) with theobjective value equal to −kx1 − x2 k. On the other hand, it follows from Lemma 3.3
and (3.15) that any KKT point for (1.1) with the objective value −ϕ yields a KKTpoint (x1 , x2 ) for (3.1) with positive Lagrange multipliers such that kx1 − x2 k = ϕ.
This means if there were a better KKT point for the primal problem, then there

10

SATORU IWATA, YUJI NAKATSUKASA, AND AKIKO TAKEDA

should have been a better one in the relaxed dual than the selected one. Thus thealgorithm computes the signed distance correctly in the overlapping case.
We now summarize the whole algorithm for computing the signed distance by apseudocode.
Algorithm 3.1 Compute the signed distance between two ellipsoids E1 and E2 .
1: Test whether E1 and E2 overlap via the problem (2.3). If they do not overlap,
compute the distance as a convex optimization problem [21].
2: Solve the bivariate determinantal equations (3.11).
3: Among the pairs with λ > 0, µ > 0, find the corresponding x1 , x2 .
4: For each pair obtained in Step 3, compute and compare kx1 − x2 k. The pair withthe smallest one gives the signed distance by dist(E1 , E2 ) := −kx1 − x2 k.
4. Solving the bivariate determinantal equations. We now discuss how
to solve the bivariate determinantal equations (3.11) for λ and µ. The algorithm isshown to run in O(n6 ) time, which is the dominant cost in Algorithm 3.1.
4.1. Reduction to univariate linear eigenvalue problems. Observe that
(3.11) is a quadratic two-parameter eigenvalue problem, for which recent studies describe numerical solutions by reducing it to a single-parameter eigenvalue problem[16, 18, 23]. In our case, (3.11) has a special form that can be expressed as(4.1)

det(F00 + λF10 + µF01 + λµF11 ) = 0,

(4.2)

det(G00 + λG10 + µG01 + λµG11 ) = 0.

In the general quadratic problem there are two additional terms F20 λ2 + F02 µ2 in(4.1) and G20 λ2 + G02 µ2 in (4.2), but in our case we do not have them. We cantake advantage of this structure to reduce the computational cost: following [18], the2n × 2n quadratic two-parameter eigenvalue problem (4.1), (4.2) can be solved via
two (2n)2 × (2n)2 quadratic single-parameter eigenvalue problems2(4.3) det Q1 (λ) := det λ2 (F11 ⊗ G10 − F10 ⊗ G11 ) + (F01 ⊗ G00 − F00 ⊗ G01 )+ λ(F11 ⊗ G00 − F00 ⊗ G11 − F10 ⊗ G01 + F01 ⊗ G10 ) = 0,(4.4) det Q2 (µ) := det µ2 (F11 ⊗ G01 − F01 ⊗ G11 ) + (F10 ⊗ G00 − F00 ⊗ G10 )+ µ(F11 ⊗ G00 − F00 ⊗ G11 + F10 ⊗ G01 − F01 ⊗ G10 ) = 0.
The solutions (λ, µ) for (3.11) match those of (4.3) and (4.4). Moreover, if the eigenvectors of (3.11) are v1 and v2 , i.e., M1 (λ, µ)v1 = 0 and M2 (λ, µ)v2 = 0, then those of(4.3) and (4.4) are v1 ⊗ v2 , with eigenvalues λ and µ, respectively [18]. That is, eachsolution of (4.3) is paired with one solution of (4.4). In our case, defining J = OI OI ,we have




O
−A2
A1 −A1(4.5)
F00 = O, F10 =, F01 =, F11 = J,
−A2
O
−A1
D
2 By solving det M1 (λ,µ)λµ

M2 (λ,µ)= 0 for λ, µ 6= 0 we can obtain a linear two-parameterλµ
1
1and µ , whose standard solution [2] gives (4.7) and (4.8). This
λ

= det

eigenproblem with respect tobypasses the quadratic eigenproblems (4.3), (4.4). Here we have chosen to work directly with theoriginal Lagrange multipliers λ, µ and exhibit the quadratic eigenproblems since the Bézoutian connection in Section 4.2 is most explicit with respect to Q1 (λ), Q2 (µ).

SIGNED DISTANCE BETWEEN ELLIPSOIDS


(4.6)

G00 = O, G10 =

A2
−A2



−A2
O
, G01 =
D
−A1

11


−A1
, G11 = J.
OGenerally, we can solve the single-parameter quadratic eigenvalue problem (4.3), (4.4)
by linearization, which results in a linear (generalized) eigenvalue problem of doubled
size. However, here we have F00 = G00 = O, which can be taken advantage of. All
the terms involving F00, G00 in (4.3) and (4.4) are zero, and in particular the constantterms (F01 ⊗ G00 − F00 ⊗ G01 ) and (F10 ⊗ G00 − F00 ⊗ G10 ) are both zero. It follows
that instead of solving the quadratic eigenvalue problems, we can simply solve the
4n2 × 4n2 linear generalized eigenvalue problems(4.7) det L1 (λ) := det λ(F11 ⊗ G10 − F10 ⊗ G11 ) + (F01 ⊗ G10 − F10 ⊗ G01 ) = 0,(4.8) det L2 (µ) := det µ(F11 ⊗ G01 − F01 ⊗ G11 ) + (F10 ⊗ G01 − F01 ⊗ G10 ) = 0.
This problem is twice smaller than the 8n2 × 8n2 eigenvalue problem one obtains by
reducing the quadratic eigenvalue problems (4.3), (4.4) by the standard approach of
linearization [31]. This size reduction has the effect of removing the “trivial” solutions
λ = µ = 0 for det M1 = det M2 = 0.
4.2. Connection with Bézoutians. Here we mention a connection of theabove process to Bézoutians, which is also used later. Forming Q1 , Q2 in (4.3),(4.4) from M1 , M2 is related to constructing the Bézoutians for matrix polynomials [28]. Specifically, this conversion from a two-parameter eigenvalue problem to asingle-parameter eigenvalue problem is equivalent to taking the Bézoutian for the two
matrix polynomials I2n ⊗ M1 (λ, µ) and M2 (λ, µ) ⊗ I2n , which are of size 4n2 × 4n2 .
Here the Kronecker products are taken to achieve commutativity, which makes it
convenient for forming the Bézoutian for matrix polynomials [20].
Matrix polynomials P1 and P2 are said to commute if P1 (ξ)P2 (ξ) = P2 (ξ)P1 (ξ)
holds for every value of ξ. The Bézoutian for commuting regular matrix polynomials
P1 , P2 of size n × n and degree k is defined by the bivariate matrix polynomial(4.9)

B(s, t) =P1 (s)P2 (t) − P2 (s)P1 (t)=s−tk−1,k−1X

Bij si tj

i,j=0

in s and t. Here Bij is the n × n coefficient matrix corresponding to the term si tj inB(s, t). Then the block symmetric nk × nk matrix


B0,0
···
B0,k−1..
..
B =  ...
.
.
.
Bk−1,0···Bk−1,k−1

is called the Bézout matrix.
Lemma 4.1 ([20, Theorem 1.1]). The Bézout matrix B is singular if and only
if P1 and P2 share an eigenpair (ξ, v), i.e., a scalar ξ and a vector v 6= 0 such that
P1 (ξ)v = P2 (ξ)v = 0.
More generally, the null space of the Bézoutian is related to the so-called common
restriction [13, 20].
We now claim that det(Q1 (λ)) = 0 in (4.3) solves the equation det(B(λ)) = 0 for
λ, where B(λ) is the Bézout matrix between P1 (µ) := I2n ⊗ M2 (λ, µ) and P2 (µ) :=M1 (λ, µ) ⊗ I2n with respect to µ. More specifically, define the Bézoutian
1
B(s, t) =(I2n ⊗ M2 (λ, s))(M1 (λ, t) ⊗ I2n ) − (M1 (λ, s) ⊗ I2n )(I2n ⊗ M2 (λ, t)) .
s−t

12

SATORU IWATA, YUJI NAKATSUKASA, AND AKIKO TAKEDA

P1 P1Then (writing B(s, t) = i=0 j=0 si tj Bij we have Q1 (λ) = B10 = B01 (which is amatrix polynomial in λ), the coefficient matrix of B(s, t) with linear terms in s, or
equivalently t. Noting that B11 = O,)←B(s, t) is actually a scalar here, so the Bézout
matrix isB00
Q1 (λ)(4.10)
B(λ) = ()←−Q1 (λ),
Q1 (λ)
Oso solving det Q1 (λ) = 0 is equivalent to finding the values of λ = λ∗ for which theBézout matrix is singular. By Lemma 4.1, this is in turn equivalent to finding λ∗such that (I2n ⊗ M2 (λ∗ , µ∗ ))v = (M1 (λ∗ , µ∗ ) ⊗ I2n )v = 0 for some µ∗ and v 6= 0. The
discussion for Q2 (µ) is completely analogous, in which we start with the bivariate
matrix polynomial
B(s, t) =
1(I2n ⊗ M2 (s, µ))(M1 (t, µ) ⊗ I2n ) − (M1 (s, µ) ⊗ I2n )(I2n ⊗ M2 (t, µ)) .
s−tThis process is a generalization from scalars to matrices of solving a pair of bivariate polynomial equations via Bézoutians [5, 24].
4.3. Relation between Mj and Lj . The Bézoutian connection just describedreveals the relation between Mj and Lj , namely we can prove that the regularity of
L1 , L2 is precisely the condition for det M1 (λ, µ) = det M2 (λ, µ) = 0 to have a finite
number of solutions.
Lemma 4.2. Let M1 (λ, µ), M2 (λ, µ) and L1 (λ), L2 (µ) be as defined in (3.10),(4.7), (4.8). Then the pencils L1 (λ) and L2 (µ) are regular if and only if the numberof solutions (λ, µ) for which det M1 (λ, µ) = det M2 (λ, µ) = 0 is finite.
Moreover, if det M1 (λ, µ) = det M2 (λ, µ) = 0 has infinitely many solutions, thenfor every fixed λ (or µ) there exists a value of µ (or λ) for which det M1 (λ, µ) =det M2 (λ, µ) = 0.
Proof. The key is to use the connection to Bézoutians. As seen in Section 4.2,
values of λ for which det Q1 (λ) = 0 are those for which the Bézout matrix B(λ)
is singular between P1 (µ) = I2n ⊗ M2 (λ, µ) and P2 (µ) = M1 (λ, µ) ⊗ I2n viewed as
matrix polynomials in µ, regarding λ as a hidden variable. Hence by Lemma 4.1,
if P1 (µ), P2 (µ) are regular matrix polynomials at λ = λ0 , they share an eigenpair
2(µ0 , x), that is, P1 (µ0 )x = P2 (µ0 )x = 0 for a nonzero x ∈ R4n , and in this caseL1 (λ0 ) has the same null vector x. The discussion for L2 is entirely analogous.
To invoke this result it is crucial to show that L1 , L2 correspond to the Bézout
matrices of two regular matrix polynomials. Now note that for any fixed value ofλ∈/ λ(A1 ), both M1 (λ, µ) and M2 (λ, µ) are regular matrix polynomials in µ, in that
there exist values of µ for which they are nonsingular.
Note that care is needed when λ is an eigenvalue of A1 ; in this case M1 is still
regular but M2 can be singular. Nonetheless, for each such λ, det M1 (λ, µ) = 0 has nsolutions and hence finite, so the statement is unaffected by λ ∈ λ(A1 ).
It follows that Q1 (λ) is always (the off-diagonal (1,2)-block part of)←−1 times
the Bézout matrix (4.10) between two regular matrix polynomials P1 (µ), P2 (µ), andsince L1 (λ) = Q(λ)λ , L1 (λ0 ) is a singular matrix if and only if P1 (µ), P2 (µ) share an
eigenpair (µ0 , x) at λ0 .
This condition is equivalent to det M1 (λ0 , µ0 ) = det M2 (λ0 , µ0 ) = 0. To see this,note that P1 (µ0 ) is singular if and only if M2 (λ0 , µ0 ) is singular, and if M2 v2 = 0 then
P1 (µ) has null space w ⊗ v2 for any vector w ∈ R2n . Indeed, we have P1 (µ)(w ⊗ v2 ) =SIGNED DISTANCE BETWEEN ELLIPSOIDS

13(I2n ⊗ M2 )(w ⊗ v2 ) = w ⊗ M2 v2 = 0. Similarly, P2 (µ) is singular if and only if
M1 (λ0 , µ0 ) is singular, and if M1 v1 = 0 then P2 (µ) has null space v1 ⊗ w for any
w ∈ R2n . Hence taking x := v1 ⊗ v2 we have P1 (µ0 )x = P2 (µ0 )x = 0, or equivalently
L1 (µ0 )x = L2 (µ0 )x = 0.
To establish the last statement of the lemma we see that if L1 (λ) is singularthen for every fixed λ = λ0 the matrix L1 (λ0 ) has a null vector, and by the Bézoutian
connection this means the matrix pencils M1 (λ0 , µ) and M2 (λ0 , µ) share an eigenvalueµ. The same argument proves the counterpart for any fixed µ = µ0 .
Note that the proof indicates that λ ∈ λ(A1 ) needs some care: in our algorithm wetreat it separately by taking λ to be each eigenvalue of A1 and solving det(M1 (λ, µ)) =0, then check if any of the solutions satisfy det(M2 (λ, µ)) = 0; we keep those that doin the candadates for the KKT multipliers. This process costs O(n4 ), and negligible
in the overall algorithm.
4.4. The number of solutions is finite in the generic case. For the process
in Section 4.1 to correctly find the KKT points, we need to show that the solution
to (3.11) is zero-dimensional, that is, the number of pairs (λ, µ) satisfying (3.11) isfinite. In fact, if we can show it is finite, then from the linear eigenvalue problems(4.7) and (4.8), (or by Bézout’s theorem) we see that it is in fact bounded by 4n2 .
Lemma 3.2 shows that every pair (λ∗ , µ∗ ) corresponding to a KKT point satisfies
det M1 (λ∗ , µ∗ ) = det M2 (λ∗ , µ∗ ) = 0. In addition, if M1 (λ, µ)v1 = 0 and M2 (λ, µ)v2 =0, then L1 (λ)(v1 ⊗ v2 ) = L2 (µ)(v1 ⊗ v2 ) = 0 ([16], see also Section 4.3). Thus it
suffices to show that the matrix pencils L1 (λ), L2 (µ) are regular. We now claim thatgenerically this property holds. First we make precise the meaning of “generic” in our
problem.
Definition 4.3. We call the signed distance problem with (A1 , A2 , b1 , b2 ) genericif at least one of A1 and A2 has no eigenvector orthogonal to d = b1 − b2 .
Calling this generic is justified because the property holds for any positive definite
matrices A1 , A2 except for a set of measure zero, as long as d 6= 0.
We now show that L1 , L2 are regular pencils for generic (A1 , A2 , b1 , b2 ).
Theorem 4.4. The matrix pencils L1 (λ), L2 (µ) in (4.7), (4.8) are regular, i.e.,there exist (λ, µ) such that det L1 (λ) 6= 0 and det L2 (µ) 6= 0, if at least one of A1 andA2 has no eigenvector orthogonal to d.
Proof. We write L1 (λ) = λX1 + Y1 , L2 (µ) = µX2 + Y2 and proceed by examining
the linear terms X1 , X2 , whose nonsingularity is a sufficient condition for L1 (λ), L2 (µ)
to be regular. We have 
A2 −A2O−A2
X1 = (F11 ⊗ G10 − F10 ⊗ G11 ) = J ⊗
−
⊗ J.
−A2
D
−A2
O
To simplify the expression, we reformulate the problem in the coordinate system
e1 := Q> A1 Q, eb1 =determined by the eigendecomposition A2 = QΛ2 Q> , in which A>>>>
e := Q DQ = d˜d˜ . Then we define the matrix
Q b1 , eb2 = Q b2 , d˜ = eb1 − eb2 , and D
e1 unitarily similar to X1 by
X

 

Λ2 −Λ2
O
−Λ2
e
e
e
e
e
X1 = (F11 ⊗ G10 − F10 ⊗ G11 ) = J ⊗⊗ J.
e − −Λ2O−Λ2
D
We can write this as"#e
2 ×2n2
O
X
2n
e1 =
X
,
e
X
O2n2 ×2n2

e = diag(Ge 10 + λ1 J, G
e 10 + λ2 J, . . . , G
e 10 + λn J),
X

14

SATORU IWATA, YUJI NAKATSUKASA, AND AKIKO TAKEDA

e1 is singular if and only ifwhere λi is the ith eigenvalue of A2 . It follows that X
e
G10 + λi J is singular for some i. For each i we haveΛ2
λi I − Λ2e 10 + λi J =G.
λi I − Λ2dede>Note that the off-diagonal blocks have a zero diagonal in the (i, i) entry.
e 10 + λi J) is singular if and only if dei = 0. Since the argument
By Lemma 3.1, (Gholds for every i, we conclude that X1 is singular if and only if dei = 0 for some i. In
the original coordinate system, this condition asserts that A2 has an eigenvector thatis orthogonal to d. An analogous argument shows that X2 is singular if and only if
A1 has an eigenvector that is orthogonal to d.
When just one of A1 and A2 has an eigenvector v orthogonal to d, one of L1 or
L2 has a nonsingular linear coefficient, so we can compute 4n2 values of λ (or µ) andsolve det M1 = det M2 = 0 for the other variable. Note that for λ fixed to any value,det M1 and det M2 are not identically zero as µ varies. Hence det M1 = det M2 = 0
yields a pair of regular eigenvalue problems. Hence for each λ there are finitely manysolutions. It follows that det M1 = det M2 = 0 has a finite number of solutions, andby Lemma 4.2 it follows that L1 , L2 are both regular.
We have thus proved that L1 and L2 are regular matrix pencils in the genericcase. It remains to consider the situation where both A1 and A2 have eigenvectorsthat are orthogonal to d.
Note that indeed L1 and L2 can be singular. Suppose that A1 and A2 share an O >  Q O eigenvector v orthogonal to d. For simplicity, we consider M̂1 := Q
M1 O Q ,O Q
 Q O >  Q O 
M̂2 := O Q M2 O Q , where Q is an orthogonal matrix whose last column is v.
Then we see that for every λ, there exists µ such that M̂1 , M̂2 both have a zero column
in the second block columns, which shows they are clearly singular. In such cases some
remedy is needed for the algorithm to compute the signed distance correctly, as we
describe in Section 5.
4.5. Complexity analysis. We now examine the computational cost. The dominant cost lies in solving the linear generalized eigenvalue problems det L1 (λ) = 0,
det L2 (µ) = 0 of size (2n)2 = 4n2 . Since the standard QZ algorithm for computing the eigenvalues of an N × N linear generalized eigenvalue problem requires
about 30N 3 floating point operations [14, §7.7.7], the computational cost is about30(4n2 )3 ≈ (1.9 × 103 )n6 flops.
Since a ballpark estimate of a feasible arithmetic cost for a current standard
desktop machine is about 1013 flops, a realistic limit for our algorithm is n . 50.
5. Degenerate case analysis. Algorithm 3.1 successfully computes the signed
distance for almost all pairs of E1 and E2 , including virtually any randomly generated
example. However, there remain cases where the matrix pencils L1 , L2 are singular,
and so the algorithm may fail to compute all the KKT points. In this section, we
analyze such situations and show that Algorithm 3.1 can still be used with appropriate
treatments to compute the signed distance reliably.
5.1. Conditions for L1 , L2 to be singular. We have seen in Section 4.4 thata necessary condition for L1 , L2 to be singular is that both A1 and A2 have an eigenvector that is orthogonal to d. We next present a sufficient condition for L1 , L2 to be
singular. Recall that a linear subspace W is called an invariant subspace of a squareSIGNED DISTANCE BETWEEN ELLIPSOIDS

15

matrix A if AW ⊆ W . By the same notation, we may express an orthonormal basisof the subspace.
Proposition 5.1. Suppose that A1 and A2 have a shared invariant subspace W
that is orthogonal to d. Then det M1 (λ, µ) = det M2 (λ, µ) = 0 has infinitely manysolutions with N (λ, µ) singular, and L1 (λ), L2 (µ) are singular matrix pencils.
Proof. Let k = dim W . By the assumption, W > d = 0 and there exist k × k
matrices X, Y such that A1 W = W X, A2 W = W Y . HenceN (λ, µ)W = λµW − µA1 W − λA2 W = W (λµI − µX − λY ),and for each fixed value of µ, there exists λ such that µXx = λ(µI − Y )x with x 6= 0,
i.e., λ is a generalized eigenvalue of the matrix pencil µX − λ(µI − Y ). Hence for the>nonzero vector v :=
 x ∈ W we have N (λ, µ)v = 0 and v d = 0, so M1 and M2
 W
have a null vector v0 , hence det M1 = det M2 = 0. It follows that det L1 (λ) = 0 for
every value of λ, and hence, by Lemma 4.2, L1 is singular. The discussion for L2 isanalogous.
Note that most of the solutions (λ, µ) for det M1 = det M2 = 0 with N singular as
in the above proposition are irrelevant for our problem, as they do not satisfy the KKT
conditions; indeed |λ|, |µ| are necessarily bounded at the KKT points, see Section 6.2.
Below we discuss remedies for this issue of having infinitely many solutions.
Although we have derived both a necessary condition and a sufficient condition
for L1 , L2 to be singular (and we derive another sufficient condition in Section 5.4.2),
identifying the precise necessary and sufficient condition appears to be highly nontrivial. We leave this as an open problem.
5.2. Regularity test for L1 , L2 . We describe a process to test, and guarantee
when successful, that the matrix pencils L1 , L2 are regular. One direct method is to
use the GUPTRI algorithm [10] for L1 or L2 , but this can cost O(n8 ) operations in
the worst case as the matrices have size O(n2 ). We look for cheaper alternatives.
By the necessary condition, the first test to guarantee regularity of L1 , L2 is tosee if both A1 , A2 have an eigenvector orthogonal to d. This can be done with O(n3 )cost. Generic cases pass this test, in which we proceed with solving (4.7), (4.8).
If this first test fails, we employ another test: fix λ to an arbitrary number andcheck whether det M1 (λ, µ) = det M2 (λ, µ) = 0 has a solution µ by solving two linear
eigenvalue problems of size 2n, which can be done with O(n3 ) cost. If there is no
solution µ, then we can safely conclude that L1 , L2 are regular. On the other hand,
if the test fails for more than 4n2 distinct values of λ, we may assert that L1 , L2
are singular. In this case, we take the additional preprocessing steps described below
before invoking Algorithm 3.1. Thus the test for regularity requires O(n5 ) cost in theworst case.
5.3. Preprocessing by dimension reduction. Our next goal is to rectify
Algorithm 3.1 so that it can deal with pathological cases where L1 , L2 are singular.
Proposition 5.1 implies that a sufficient condition for L1 , L2 to be singular is for
A1 and A2 to have a common eigenvector orthogonal to d. Such cases may arise
naturally in practice, an extreme example being A1 = A2 = I. In general, if theintersection of eigenspaces of A1 and A2 is of dimension two or larger, then there
must be a common eigenvector orthogonal to d.
Fortunately, there is a remedy in such cases. Let W1 and W2 be the eigenspacesof A1 and A2 corresponding to the eigenvalues ρ1 and ρ2 , respectively, and let W :=W1 ∩ W2 be nonzero. We now define subspaces Ŵ and W̄ as follows. If W > d 6= 0,

16

SATORU IWATA, YUJI NAKATSUKASA, AND AKIKO TAKEDA

let W̄ be the set of vectors in W orthogonal to d. If W > d = 0, we take W̄ to be
an arbitrary subspace of W with one dimension less than W . Then we define Ŵ as
the orthogonal complement of W̄ in Rn . Using the same notation, we may designate
orthonormal bases of these subspaces. We can simultaneously block-diagonalize A1and A2 by the orthogonal matrix [Ŵ W̄ ], so that




Â1 O
Â2 O
>
>
(5.1)
[Ŵ W̄ ] A1 [Ŵ W̄ ] =, [Ŵ W̄ ] A2 [Ŵ W̄ ] =
.
O Ā1
O Ā2
We then solve the problem of reduced size with d := Ŵ > d, A1 := Â1 , A2 := Â2 . We
can repeat the process when more than one eigenspaces are shared.
Algorithm 5.1 Preprocessing by dimension reduction.
1: Compute the eigenvalue decompositions of A1 , A2 , and look for the intersection
W of eigenspaces with dim W ≥ 2.
2: If such W exists, then take Ŵ as in (5.1) and call Algorithm 5.1 with b1 := Ŵ > b1 ,A1 := Ŵ > A1 Ŵ , b2 := Ŵ > b2 , A2 := Ŵ > A2 Ŵ .
3: If there exists no such W , then return (E1 , E2 ) with E1 := E(b1 , A1 ) and E2 :=E(b2 , A2 ).
Below we verify that this process yields a pair of ellipsoids having the same
 signed
distance as the original one. Recalling the definition (1.1) and writing w = ŵw̄ withŵ ∈ Ŵ and w̄ ∈ W̄ , we havedist(E1 , E2 ) =max

kŵk2 +kw̄k2 =1

f (ŵ, w̄)where using W̄ > d = 0 we can writeqq
f (ŵ, w̄) = hŵ, Ŵ > di − hŵ, Â1 ŵi + hw̄, Ā1 w̄i − hŵ, Â2 ŵi + hw̄, Ā2 w̄i.
The following result shows that the dimension reduction employed in Algorithm 5.1
results in an equivalent signed distance problem of reduced size.
Proposition 5.2. Under the above notation, for every pair (ŵ, w̄) with kŵk2 +2kw̄k = 1 and w̄ 6= 0, there exists a vector ŵ0 with kŵ0 k = 1 such that f (ŵ0 , 0) ≥f (ŵ, w̄).
Proof. Note that there exists a common eigenvector v̂ such that Â1 v̂ = ρ1 v̂ and>Â2 v̂ == 1. Put ĉ = v̂ > ŵ andp ρ2 v̂. We choose v̂ 0so that hv̂,0 W di ≥ 0 and kv̂k
0
0 2
2
2
ĉ = ĉ + kw̄k . Then ŵ = ŵ + (ĉ − ĉ)v̂ satisfies kŵ k = kŵk2 + kw̄k2 = 1. In
addition, we havehŵ0 , Âj ŵ0 i = hŵ, Âj ŵi + hw̄, Āj w̄ifor j = 1, 2, and hencef (ŵ0 , 0) − f (ŵ, w̄) = hŵ0 , W̄ > di − hŵ, W̄ > di = (ĉ0 − ĉ)hv̂, W̄ > di ≥ 0,
as required.
The preprocessing step in Algorithm 5.1 results in an equivalent signed distanceproblem with lower-dimensional Aj , bj . The reduced problem belongs to the genericcase (no eigenvector is orthogonal to d) for some instances including those with A1 =α1 I, A2 = α2 I, and d 6= 0 in the original setting.

17

SIGNED DISTANCE BETWEEN ELLIPSOIDS

5.3.1. Determining shared eigenspace. In practice, numerically the sharedeigenspace is not exactly shared, and instead of (5.1) we haveÂ1 E >
Â2 F >
(5.2)
[Ŵ W̄ ]> A1 [Ŵ W̄ ] =, [Ŵ W̄ ]> A2 [Ŵ W̄ ] =
,
E Ā1
F Ā2
where kEk, kF k are small but nonzero. The question arises therefore of when to
determine the eigenspaces are shared and set E, F to the zero matrix and invoke
Algorithm 5.1.
We introduce a criterion based on conditioning analysis. In Section 5.4.3 we showthat the condition number of the signed distance with respect to perturbation in Aj is√ 1+ √ kF k.
, so setting E, F ← O results in error bounded by √ kEk
2

σmin (Aj )2σmin (A1 )2σmin (A2 )On the other hand, when (5.2) holds, defining w̄ toin the nearly an eigenvector
 beshared eigenspace orthogonal to d and setting v = w̄0 ⊗ w̄0 , we see that L1 (λ)v andL2 (µ)v are both O(kEk + kF k) in norm for all λ, µ of O(1). Thus L1 , L2 are nearly
singular if kEk, kF k are small. Such near singularity of a pencil is known [29, § 15.1]
to inversely affect the conditioning of the eigenvalues, and we expect error of size
(kA1 k+kA2 k)
in the computed eigenvalues. Hence the computed x1 , x2 by (3.8),(3.9)
kEk+kF k
1 k+kA2 k)kdk
have error estimate of size (kAkEk+kF. These observations suggest that we choosekto set E, F = O when it gives the smaller error estimate, that is, when below holds:

(5.3)

kF k
(kA1 k + kA2 k)kdk
kEk
p
+ p
≤
.
kEk + kF k
2 σmin (A1 ) 2 σmin (A2 )√√Note that (5.3) is invariant under the scaling A1 , A2 , b1 , b2 ← cA1 , cA2 , cb1 , cb2
for a scalar c > 0, which (recalling (1.1)) is the scaling that results in an equivalent
problem. Roughly, when A1 , A2 are well-conditioned
and O(1) in norm, (5.3) means√
we set E, F to zero when they are of size .
5.4. Perturbing the ellipsoids to avoid singular L1 , L2 . There are situations in which even with the preprocessing step just described, Algorithm 3.1 does
not result in L1 , L2 regular. An example is when W > d = 0 in the above setting.
Another extreme example is when b1 = b2 , i.e., when the center of the ellipsoids
are the same. In this case d = 0, and it is easy to see that (3.11) holds for every
pair (λ, µ) such that N (λ, µ) is singular. Even worse, one can show that even the
Kronecker canonical form of the linear problems L1 , L2 do not provide the KKTpoints λ, µ. Handling this case analytically appears to be surprisingly difficult.
To illustrate the distinct cases, Figure 5.1 shows the representative configurationof the ellipsoids in the plane n = 2: (i) the generic case, (ii) the degenerate casetreated in Section 5.3, and (iii) the degenerate case treated in this subsection.
5.4.1. Perturbing a center. A simple remedy for the difficulty arising fromd = 0, and also the case where W := W1 ∩ W2 as above satisfies W > d = 0, is tointroduce a small perturbation in the center bj , thus perturbing d.
Specifically, the following procedure constructs a perturbed problem such thateigenvectors of A1 corresponding to simple eigenvalues all become non-orthogonal tod. The perturbed problem belongs to the generic case for example when the original
Aj had only simple eigenvalues.
5.4.2. Perturbing a matrix. It is sometimes necessary to also perturb thematrices A1 , A2 . This happens when for every fixed λ, (even after the dimension

18

SATORU IWATA, YUJI NAKATSUKASA, AND AKIKO TAKEDA

1
1

1

0

0

0

−1

−1
−1

−1

0

1

2

3

−1

0

1

2

−2

−1

0

1

2

Fig. 5.1. Two-dimensional illustration of generic and degenerate cases. Left: generic casewhere at least one of A1 or A2 has no eigenvector (shown by arrows) orthogonal to d, shown by
the dotted line. Center: degenerate case of shared invariant subspace caused by multiple eigenvaluesof A1 and A2 . This case can be resolved by the process of dimension reduction as described in
Section 5.3. Right: degenerate case of shared eigenspace caused by eigenvectors for distinct eigenvalues being orthogonal to d. In this case dimension reduction cannot be employed; indeed thetangent planes corresponding to the signed distance (dashed black lines) are not perpendicular to thehorizontal axis. In this case our solution is to slightly perturb the center b1 .

Algorithm 5.2 Preprocessing by perturbing bj .
1: Let Aj = QΛj Q> be the eigenvalue decomposition with an orthogonal matrix Q.
2: For each column vector qi of Q corresponding to a simple eigenvalue, if qi> d = 0,then set bj := bj + i qi with an arbitrary nonzero number i of magnitude O(),where  is the machine precision.

reduction) there exists µ such that N (λ, µ) = λµI − µA1 − λA2 has a null space ofdimension ≥ 2. An n = 4 example is00A1 = 2I + 
1
0

0
0
0
0

1
0
0
1


0
0
,
1
0



0
1

A2 = 2I + 
0
0

1
0
0
1

0
0
0
0


0
1
.
0
0

In such cases, L1 , L2 are singular regardless of d.
Proposition 5.3. Suppose that for every fixed λ ∈ R, there exists µ ∈ R suchthat N (λ, µ) = λµI − µA1 − λA2 has a null space of dimension ≥ 2. Then L1 , L2
in (4.7), (4.8) are singular matrix pencils.
Proof. Let λ ∈ R any fixed nonzero finite number. By the assumption, there
exists a finite µ ∈ R such that N (λ, µ) has a null space of dimension two or larger,
which means the nullcontains a nonzero vector x with x> d = 0. Then we easily
 0 space2nsee that the vector x ∈ R is a null vector of both M1 (λ, µ) and M2 (λ, µ), so (λ, µ)is a solution of det M1 (λ, µ) = det M2 (λ, µ) = 0. Hence by Lemma 4.2 we conclude
that L1 , L2 are singular.
This is another sufficient condition for L1 , L2 to be singular, which is neither asubset nor a superset of the condition in Proposition 5.1.
Since the statement of the proposition makes no reference to d, this result indicatesthat perturbing the ellipsoid centers is sometimes not enough. If we do not pass theregularity test in Section 5.2 after having perturbed the center bj as in Algorithm 5.2,
then as a last resort we perturb the matrix Aj as well.
The matrix Aj thus obtained has only simple eigenvalues, and no eigenvector is
orthogonal to d, resulting in a generic situation of Definition 4.3.

SIGNED DISTANCE BETWEEN ELLIPSOIDS

19

Algorithm 5.3 Preprocessing by perturbing Aj and bj .
1: Let Aj = QΛj Q> be the eigenvalue decomposition with an orthogonal matrix Q.
2: For each column vector qi of Q corresponding to a multiple eigenvalue, if qi> d = 0,then set bj := bj + i qi with an arbitrary nonzero number i of magnitude O().
3: Update Aj := Q(Λj + ∆)Q> , where ∆ is a diagonal matrix with nonzero distinct

diagonals of magnitude O().

5.4.3. Sensitivity analysis under perturbation. Perturbing the inputs clearly
perturbs the output dist(E1 , E2 ), which is undesirable if the change in the output
is significantly magnified. We now examine the sensitivity of the signed distancedist(E1 , E2 ) under perturbation in bj and Aj .
Recalling (1.1), we first examine the term hw, bj i, the only term that depends on
bj . Defining b̂j = bj + δbj , for any vector w with hw, wi = 1, we have|hw, b̂j i − hw, bj i| ≤ max |hw, δbj i| = kδbj k.
hw,wi=1

Therefore, perturbing bj by  can only perturb the signed distance by at most , thatis, the signed distance is a well-conditioned function of bj , so slightly perturbing theellipsoid centers is always harmless.
Perturbation in Aj is a bit more subtle. Defining Âj = Aj + δAj and assuming
the perturbation is small enough so that kδAj k ≤ σmin (Aj ), in which we take thespectral norm kAk := σmax (A), we haveqqqqhw, Âj wi − hw, Aj wi ≤hw, Aj wi − kδAj k − hw, Aj wikδAj kphw, Aj wi − kδAj k + hw, Aj wikδAj kkδAj k≤ p.
. p
2 hw, Aj wi
2 σmin (Aj )=p

Strict inequality holds here up to O(kδAj k2 ). Using these bounds we obtainkδA1 kkδA2 k(5.4) |dist(E1 , E2 ) − dist(Ê1 , Ê2 )| . kδb1 k + kδb2 k + p+ p.
2 σmin (A1 ) 2 σmin (A2 )Note that equality can be attained in (5.4) up to O(kδAj k2 ): denoting by v theeigenvector for σmin (Aj ), suppose w = v and take δAj = vv > . The above argumentshows that the conditioning with respect to perturbation in Aj is proportional to thesmallest eigenvalue of Aj , so care is needed when Aj is ill-conditioned. For a matrix
A, we denote by κ(A) its condition number, i.e., κ(A) := σmax (A)/σmin (A). Between
A1 and A2 , it is advised to perturb Aj with smaller condition number. If κ(Aj )  1,then the change in the solution may be unacceptably large, and one may be advised
to use higher precision arithmetic. If κ(Aj ) = O(1), which is easy to verify, thenperturbing Aj can be done safely.
We note that while the conditioning of the signed distance can be analyzed asabove, that of the corresponding points (x1 , x2 ) on the ellipsoids cannot be boundedin general. For an extreme example, let E1 , E2 both be unit spheres centered at theorigin. Then the signed distance is clearly −2, attained at (u, −u) for an arbitrary

20

SATORU IWATA, YUJI NAKATSUKASA, AND AKIKO TAKEDA

point u with kuk = 1. However, by a small  perturbation in b1 , the corresponding points become uniquely determined with signed distance  − 2, illustrating thesensitivity of the points (x1 , x2 ) to perturbation in bj .
This does not necessarily indicate the weakness of the algorithm, but rather thedifficulty of the problem of finding xj . Indeed, small perturbation in Aj , bj can be
regarded as a small backward error [15, §1.8], so the backward stability is not lost.
In practice, perturbing b1 to introduce a small nonzero d still results in det M1 (λ, µ),det M2 (λ, µ) ≈ 0 whenever det N = 0, indicating numerical instability may√be an issue. One remedy for this is to apply a diagonal congruence scaling diag( cI, √1c I)
for some c > 0 to balance the norms of the diagonal blocks: we redefinec2 λA2 N (λ, µ)
c1 µA1 N (λ, µ), M2 (λ, µ) =M1 (λ, µ) :=.
µλN (λ, µ)N (λ, µ)
c1 D
c2 D
q
q
kDk
kDk
An appropriate choice is c1 ≈ kA
,
c
≈
2
k
kA2 k . Clearly their determinants are1
the same as those in (3.10), but the numerical behavior is improved significantly. Thisscaling is recommended whenever the norms of kDk and kAj k differ, so we employ it
by default in our algorithm.
5.5. The linearization is regular. In summary, our overall algorithm proceeds
as follows:
Algorithm 5.4 Overall algorithm to compute the signed distance between ellipsoids.
1: If the regularity test is not passed, then reduce dimension by Algorithm 5.1.
2: If the regularity test is not passed, then select j ∈ {1, 2} as follows: set j if Aj has
only simple eigenvalues. Otherwise, take j with smaller condition number κ(Aj ),and perturb bj by Algorithm 5.2.
3: If the regularity test is not passed, then perturb Aj , bj by Algorithm 5.3.
4: Invoke Algorithm 3.1 to compute the signed distance.
The rationale behind the choice of j is that if Aj has only simple eigenvalues
then Algorithm 5.2 alone is sufficient, and otherwise we choose j so that perturbing
Aj by Algorithm 5.3 is more stable. After the preprocessing steps in Algorithms
5.1, 5.2 and 5.3, L1 , L2 are regular matrix pencils, and det M1 = det M2 = 0 has afinite number of solutions. As a result, one can reliably compute the eigenvalues of
L1 , L2 , which give a finite number 4n2 of candidate points for the KKT conditions,and hence Algorithm 3.1 computes the signed distance between overlapping ellipsoids
with O(n6 ) cost.
6. KKT points and local optimality. We now turn to analyzing the KKTpoints; in particular we discuss how to check the local optimality of the KKT points.
Although this is not directly necessary for Algorithm 3.1, which simply compares
the objective values at the KKT points, investigating the local optima reveals thenonconvex nature of the problem and deepens our understanding.
6.1. Checking local optimality of the primal KKT points. For a KKTpoint w of the primal problem (1.1) obtained through (3.18), we would like to identify
whether it is a local optimum or a saddle point. For simplicity, we introduce theequivalent problem


p
p
(6.1)
minimize −hw, b1 i + hw, A1 wi + hw, b2 i + hw, A2 wi
kwk=1

21

SIGNED DISTANCE BETWEEN ELLIPSOIDS

and its Lagrangian function(6.2)

L(w, ϕ) = −hw, b1 i +phw, A1 wi + hw, b2 i +p
ϕhw, A2 wi + (1 − kwk2 ).
2Classical results in constrained optimization [27, §12.5] implies that if a KKT point ŵ
with a Lagrange multiplier ϕ̂ is a local minimum then the projected Hessian of L at ŵ
onto the tangent plane T (ŵ) = {v : ŵ> v = 0} is positive semidefinite. Conversely, ŵ
is a strict local minimum if the projected Hessian of L onto T (ŵ) is positive definite.
Direct computation shows that the Hessian is




1
A1 ww> A1
1
A2 ww> A2
∂2L
A1 −
+p
A2 −
− ϕI.
=p
(6.3)
∂w2
w > A1 w
w > A2 w
w > A1 w
w > A2 w
Recall that our algorithm obtains the KKT point w for (1.1) through (3.18), and thecorresponding KKT point (x1 , x2 ) and Lagrange multipliers λ, µ for (3.1) are alreadyavailable. Using these we can rewrite the right-hand side of (6.3) as
 ϕ

∂2L
ϕ
A1 − (x1 − b1 )(x1 − b1 )> +A2 − (x2 − b2 )(x2 − b2 )> − ϕI.
=
∂w2
λ
µ
2

∂ L
We can then test the local optimality by checking if the matrix Q>w ∂w2 Qw is positive
definite, where Qw is the n × (n − 1) matrix such that [ŵ, Qw ] is orthogonal. This isa basic problem in linear algebra requiring O(n3 ) operations, which can be done byeither attempting the Cholesky factorization or computing the eigenvalues.
There is some subtlety in testing local optimality by this process: when the∂2L
projected Hessian Q>w ∂w2 Qw is positive semidefinite but not definite. In this case
we can neither guarantee that the KKT point is locally optimal, nor exclude it as asaddle point. Fortunately, this rarely happens. Under the nondegeneracy assumptionthat the projected Hessian is nonsingular at every KKT point, we can enumerate all
the locally optimal solutions.

6.2. Bounds for Lagrange multipliers. We show that at the KKT points the
Lagrange multipliers λ, µ are bounded in absolute value.
Proposition 6.1. At every KKT point (x1 , x2 ) of (2.1) with Lagrange multipliers
λ and µ satisfying (3.2)–(3.5), the absolute values of λ and µ are bounded byp
p
kA1 k(kb1 k + kb2 k + kA1 k + kA2 k)
p
|λ| ≤
(6.4)
,
σmin (A1 )p
p
kA2 k(kb1 k + kb2 k + kA1 k + kA2 k)
p
(6.5)
|µ| ≤
.
σmin (A2 )Proof. By (3.2) and (3.3) we have
p
p
(6.6)
σmin (Ai ) ≤ kxi − bi k ≤ σmax (Ai ),i ∈ {1, 2}.
pThese also imply bounds on kx1 k, kx2 k: we have kxi k ≤ kbi k + kAi k. Now by (3.4)
and (3.5) we have
p
p
(6.7)
|λ|kA−1
kA1 k + kA2 k.
1 (x1 − b1 )k = kx1 − x2 k ≤ kb1 k + kb2 k +22

SATORU IWATA, YUJI NAKATSUKASA, AND AKIKO TAKEDA−1
Since kA−1
1 (x1 − b1 )k ≥ σmin (A1 )kx1 − b1 k, from (6.7) we obtainpp
p
p
kb1 k + kb2 k + kA1 k + kA2 k
kA1 k(kb1 k + kb2 k + kA1 k + kA2 k)|λ| ≤≤kx1 − b1 k
kA−1
1 (x1 − b1 )kp
p
kA1 k(kb1 k + kb2 k + kA1 k + kA2 k)p,≤σmin (A1 )where we used the lower bound in (6.6) for the last inequality. This proves (6.4). The
bound (6.5) can be obtained entirely analogously.
The bounds (6.4), (6.5) can be used to exclude some of the computed solutions
for det L1 (λ) = 0, det L2 (µ) = 0 that violate them.
7. Numerical experiments. This section presents numerical experiments toillustrate Algorithm 3.1. All the experiments were conducted in Matlab 2013A.
Unless otherwise mentioned, we generate the ellipsoids by setting Ai = Xi> Xi with
Xi being n×n random matrices and by taking random vectors as bi , using the standard
normal distribution randn.
7.1. Illustration in two dimensions. For ease of visualization we first consider
the two-dimensional case n = 2.
Recall that in order to compute the signed distance our algorithm finds the KKTpoints for the relaxed dual problem (3.1). As shown in Section 3.2, the computed
KKT points contain those of the primal problem (1.1) (those with λµ > 0) and more(those with λµ < 0).
We illustrate the situation by visualizing the set x1 − x2 with the conditions
x1 ∈ bd(E1 ), x2 ∈ bd(E2 ), then plotting the KKT points of (3.1) in the left figures
below. Plotting the set {x1 − x2 | x1 ∈ bd(E1 ), x2 ∈ bd(E2 )} directly is nontrivial, so
we fix x1 ∈ bd(E1 ) on a boundary of E1 and plot x1 − x2 where x2 moves along x2 ∈bd(E2 ), then repeat this for many points on x1 ∈ bd(E1 ) (black curves). The convex
set U = E1 E2 in Section 2 is the interior of the outer boundary of the black region{x1 − x2 | x1 ∈ bd(E1 ), x2 ∈ bd(E2 )}. The KKT points computed by Algorithm 3.1
are indicated by the dots (green, blue and red); those on the boundary bd(E1 E2 )correspond to the KKT points for the primal problem (1.1). In particular, the blue
dots correspond to the local solutions as verified by the procedure in Section 6.1, andthe red dot to the global solution.
In the center figures, we draw the ellipsoids together with the KKT points, witheach pair of KKT points connected by black lines: the solid thick line corresponds tothe solution, the dashed line connects x1 , x2 at the KKT points with λµ > 0 (criticalpoints for the primal problem (1.1)), and the dot-dashed lines show the other KKTpoints with λµ < 0 (critical points for the relaxation problem (3.1), with x1 − x2 in
the interior of U = E1 E2 ).
Observe from Figures 7.1–7.4 that the number of local minima for (1.1) in theoverlapping case differs depending on the relative positioning of E1 and E2 . Our
experiments suggest that the number of local minima is 1, 2, 3 or 4. Moreover, weobserve that each local minimum for (1.1) is also a local minimum for the dual problem(2.1). The total number of KKT points for the relaxation problem (3.1) also depends
largely on the relative positioning of the ellipsoids.
Figure 7.5 is a typical illustration of a non-overlapping (convex) case, for whichmore efficient algorithms are available but Algorithm 3.1 is still applicable.

23

SIGNED DISTANCE BETWEEN ELLIPSOIDS
2.5

1.5

1.5

2

1

1
1.5

0.5

−0.5

µ

0.5

0

−1

0

−0.5

−1.5
−2

−1

0

1

2

kkt points
local minimum
solution

1
0.5

0

−0.5

0

0.5

1

1.5

−0.5
−0.5

2

0

0.5

1

1.5

2

2.5

λ

Fig. 7.1. Overlapping case with four local minima. For the figures in this subsection, the threeplots illustrate the following: Left: The Minkowski difference {x1 − x2 |x1 ∈ bd(E1 ), x2 ∈ bd(E2 )}and the KKT points. Center: ellipsoids and the KKT points. Right: (λ, µ)-values satisfying theKKT conditions.
6
2

2

5

1.5

4

1

3

kkt points
local minimum
solution

µ

1
0
−1

0.5

2

0

1
0

−0.5
−2
−2

−1

0

1

2

3

−2

−1

0

−1
−1

1

0

1

2

3

4

λ

Fig. 7.2. Overlapping case with three local minima.

2

8

2

1

6

1

kkt points
local minimum
solution

µ

0
−1
−2

0

4

−1

2

−3
0
−2

−2
−4

−2

0

2

4

−3

−2

−1

0

1

2

−1

0

1

2

3

4

λ

Fig. 7.3. Overlapping case with two local minima.
83
3

kkt points
local minimum
solution

6
2

2

µ

4

1

1

0

2
0

0
−1

−2
−1

−2
−4

−2

0

−1

0

1

2

3

−4
−2

0

2

4

6

8

10

λ

Fig. 7.4. Overlapping case with one local minimum equal to the global solution.

7.2. KKT points for larger dimension. We next run Algorithm 3.1 for larger
dimension n and illustrate by plotting the Lagrange multipliers (λ, µ) corresponding
to the KKT points for the relaxation problem (3.1). As before, the KKT points forthe primal problem are those for which λ and µ have the same sign; see Figure 7.6.
In the non-overlapping case, we generally observe that there is at least one point

24

SATORU IWATA, YUJI NAKATSUKASA, AND AKIKO TAKEDA

1

6

10

0

5

kkt points
solution

4

µ

−1
−2

5

3
−3
0

2

−4

1

−5
−4

−2

0

2

−4

−2

0

−2

2

0

2

4

λ

Fig. 7.5. Non-overlapping case. There is just one local minimum equal to the global solution.
No Overlap

Overlap

80

40

60
40

50

30
20

30

µ

µ
0

20

10

10

−20
0

−40
−60
−100

kkt points
solution40

kkt points
solutionµ20Overlap
60

kkt points
solution

−50

0

λ

50

100

−10
−10

0

0

10

20

λ

30

40

50

−10
−20

0

20

40

60

80

λ

Fig. 7.6. Values of (λ, µ) for the KKT points, non-overlapping case (left) and overlapping case(center and right).

in each quadrant of the (λ, µ)-plane and the unique point with λ, µ < 0 corresponds
to the global solution, verifying our results in Section 3.2.
In the overlapping case there are generally many more KKT points, reflecting thenon-convexity of the problem. We also observe that the solution in the overlappingcase generally corresponds to one of the KKT points closest to the origin with λ, µ > 0.
This is analogous to the point-ellipsoid distance problem, in which the solution corresponds to the KKT point with the smallest Lagrange multiplier [12, 17]. However,
the right plot of Figure 7.6 suggests that the solution may not have the smallest valueof λ or µ. Experiments suggest it may hold that the pair of Lagrange multipliers(λ∗ , µ∗ ) for the global solution admits no other positive Lagrange multipliers (λ, µ)for KKT points with λ > λ∗ and µ > µ∗ ; a precise characterization appears to be
nontrivial and we leave it as an open problem.
We also used the process in Section 6.1 to examine the number of local optima
among the KKT points. Fortunately, the projected Hessian was always either strictlyindefinite or positive definite, so we accurately counted the number of local optima.
Experiments suggest that in the overlapping case there are at most four local optima,
regardless of the dimension n. Note that for the special case of point-ellipsoid distance problem it is known that there can be at most two local minima [12, 22], andour experiments indicate that an analogous result may hold in our ellipsoid-ellipsoid
problem. Making this observation precise is also an open problem.
7.3. Runtime with increasing dimension. The dominant cost of our algorithm is in computing the eigenvalues of a 4n2 × 4n2 linear generalized eigenvalue
problem, for which the complexity is O(n6 ). To verify this numerically we generated
random examples of varying dimension n = 2, 3, . . . , 40 and examined the runtime.
Figure 7.7 confirms that the runtime scales as O(n6 ). In all cases with n ≥ 10, atleast 90% of the runtime was spent on the eigensolver.

SIGNED DISTANCE BETWEEN ELLIPSOIDS

25

4

10

runtime(s)

2

10

0

10

Algorithm
Trend O(n6)

−2

10

5

10

20

40

n

Fig. 7.7. Runtime as the dimension n varies.

8. Conclusion and discussion. We have proposed the first algorithm for computing the signed distance between two ellipsoids with polynomial complexity. The
algorithm computes the KKT points for a relaxation of the dual problem, for which weanalyze connections with the primal KKT points. A crucial component is to convert
the KKT conditions into a two-parameter eigenvalue problem, which is then reduced
to linear generalized eigenvalue problems of size O(n2 ). We analyze when the conversion results in a regular pencil, and propose remedies for the pathological cases.
We present numerical experiments, which, in addition to illustrating the ellipsoids
and the algorithm performance, pose several open problems. Besides the applications
mentioned in the introduction, the signed distance may be of use as a measure ofoverlap in the ellipsoid packing proglem [32].
Let us conclude with remarks on possible future work. First, the O(n6 ) complexity
can be a bottleneck when n is large, so a natural question is to design a more efficientalgorithm. In addition, our algorithm perturbs the ellipsoids in the pathological cases,
and it would be nice to have an approach that does not require such treatments.
Finally, an active area of research in non-convex optimization is the Celis-DennisTapia (CDT) problem, which arises as a subproblem of a nonlinear optimization
problem [9, 11]. In the CDT problem one is asked to minimize a non-convex quadratic
function over the intersection of two ellipsoids. A number of recent studies (e.g. [1, 6])
examine and analyze the CDT problem, but to our knowledge there is no deterministic
polynomial-time algorithm for solving it3 .
This paper has focused on the signed distance problem, but we suspect some of theideas that we have introduced may be useful for the CDT problem. Specifically, bothproblems have two ellipsoidal constraints, which lead to KKT conditions involving
two Lagrange multipliers. In this paper we have computed the KKT points for thesigned distance problem by solving an eigenvalue problem, and we conceive that thisframework may be applicable also to the CDT and other optimization problems.
Acknowledgements. We are very grateful to Vanni Noferini for showing us
the connection between the linearizations (4.7), (4.8) and Bézoutians, and we thankMartin Lotz for a fruitful discussion.
REFERENCES
3 A very recent preprint [4] addresses this question and describes a polynomial-time algorithm for
the CDT problem under some solvability assumptions on a system of quadratic equations. However,the polynomial degree appears to be unclear and the complexity depends on the accuracy  as log 1 .
Our approach, if successful, will have complexity O(n6 log log( 1 )).

26

SATORU IWATA, YUJI NAKATSUKASA, AND AKIKO TAKEDA[1] W. Ai and S. Zhang, Strong duality for the CDT subproblem: a necessary and sufficientcondition, SIAM J. Optim., 19 (2009), pp. 1735–1756.
[2] F. V. Atkinson, Multiparameter Eigenvalue Problems, Academic Press, NY, USA, 1972.
[3] C. Bhattacharyya, Second order cone programming formulations for feature selection, J.
Mach. Learn. Res., 5 (2004), pp. 1417–1433.
[4] D. Bienstock, A note on polynomial solvability of the CDT problem, ArXiv 1406.6429, (2014).
[5] D. A. Bini and A. Marco, Computing curve intersection by means of simultaneous iterations,
Numerical Algorithms, 43 (2006), pp. 151–175.
[6] I. M. Bomze and M. L. Overton, Narrowing the difficulty gap for the celis-dennis-tapia
problem, Preprint, Department of Statistics and Operations Research, Universitaet Wien,
Vienna, Austia, (2013).
[7] S. P. Boyd and L. Vandenberghe, Convex Optimization, Cambridge University Press, 2004.
[8] W. Briec, Minimum distance to the complement of a convex set: duality result, J. Optim.
Theory Appl., 93 (1997), pp. 301–319.
[9] M. R. Celis, J. E. Dennis, and R. A. Tapia, A trust region strategy for nonlinear equality
constrained optimization, in Numerical Optimization, P. T. Boggs, R. H. Byrd, and R. B.
Schnabel, eds., SIAM Philadelphia, 1985, pp. 71–82.
[10] J. Demmel and B. Kågström, The generalized Schur decomposition of an arbitrary pencil
A-zB: robust software with error bounds and applications. Part I: theory and algorithms,
ACM Trans. Math. Soft., 19 (1993), pp. 160–174.
[11] M. El-Alem, A global convergence theory for the celis-dennis-tapia trust-region algorithm forconstrained optimization, SIAM J. Numer. Anal., 28 (1991), pp. 266–290.
[12] G. E. Forsythe and G. H. Golub, On the stationary values of a second-degree polynomial
on the unit sphere, J. SIAM, 13 (1965), pp. 1050–1068.
[13] I. Gohberg, P. Lancaster, and L. Rodman, Matrix Polynomials, SIAM, Philadelphia, USA,2009. Unabridged republication of book first published by Academic Press in 1982.
[14] G. H. Golub and C. F. Van Loan, Matrix Computations, The Johns Hopkins University
Press, 4th ed., 2012.
[15] N. J. Higham, Accuracy and Stability of Numerical Algorithms, SIAM, Philadelphia, USA,
2nd ed., 2002.
[16] M. E. Hochstenbach, A. Muhič, and B. Plestenjak, On linearizations of the quadratic
two-parameter eigenvalue problem, Linear Algebra Appl., 436 (2012), pp. 2725–2743.
[17] S. Iwata, Y. Nakatsukasa, and A. Takeda, Global optimization methods for extended Fisher
discriminant analysis, in Proc. Seventh AISTATS, JMLR W&CP 33, 2014, pp. 411–419.
[18] E. Jarlebring and M. E. Hochstenbach, Polynomial two-parameter eigenvalue problems andmatrix pencil methods for stability of delay-differential equations, Linear Algebra Appl.,
431 (2009), pp. 369–380.
[19] A. A. Kurzhanskiy and P. Varaiya, Ellipsoidal toolbox, Tech. Report EECS-2006-46, EECS,UC Berkeley, 2006.
[20] L. Lerer and M. Tismenetsky, The Bézoutian and the eigenvalue-separation problem for
matrix polynomials, Integral Equations and Operator Theory, 5 (1982), pp. 386–445.
[21] A. Lin and S.-P. Han, On the distance between two ellipsoids, SIAM J. Optim., 13 (2002),
pp. 298–308.
[22] J. M. Martı́nez, Local minimizers of quadratic functions on Euclidean balls and spheres, SIAM
J. Optim., 4 (1994), pp. 159–176.
[23] A. Muhič and B. Plestenjak, On the quadratic two-parameter eigenvalue problem and its
linearization, Linear Algebra Appl., 432 (2010), pp. 2529–2542.
[24] Y. Nakatsukasa, V. Noferini, and A. Townsend, Computing the common zeros of two
bivariate functions via Bézout resultants, Numer. Math., (2014), pp. 1–29.
[25] J. S. Nath and C. Bhattacharyya, Maximum margin classifiers with specified false positiveand false negative error rates., in Proc. Seventh SDM, SIAM, 2007, pp. 35–46.
[26] L. Nirenberg, Functional Analysis, Academic Press, 1961.
[27] J. Nocedal and S. J. Wright, Numerical Optimization, Springer New York, 2nd ed., 1999.
[28] V. Noferini. private communication, 2013.
[29] B. N. Parlett, The Symmetric Eigenvalue Problem, SIAM, Philadelphia, 1998.
[30] A. Takeda, H. Mitsugi, and T. Kanamori, A unified classification model based on robust
optimization, Neural Computation, 25 (2013), pp. 759–804.
[31] F. Tisseur and K. Meerbergen, The quadratic eigenvalue problem, SIAM Rev., 43 (2001),
pp. 235–286.
[32] Caroline Uhler and Stephen J. Wright, Packing ellipsoids with overlap, SIAM Rev., 55(2013), pp. 671–706.