Ellipsoid method

Ellipsoid method is a weakly polynomial algorithm for solving linear programming. The preliminary version was introduced by the Russian mathematician Shor in 1977 for general convex optimization problems. Then it was applied to linear programming by Khachyan in 1979 [1]. Before interior point method, it was the first polynomial time algorithm for LP.

Feasibility implies Optimality

Now, consider the primal LP:

\mathcal{P}: \qquad \begin{aligned} &\max &c^T x & \\ &s.t. &Ax \le b & \\ & & x \ge 0 & \end{aligned}

whose dual is given by

\mathcal{D}: \qquad \begin{aligned} &\min &y^Tb \ \ & \\ &s.t. &y^TA \ge c^T & \\ & &y \ge 0 \ \ & \end{aligned}

The ellipsoid does not optimize the primal directly. Instead, it try to return a feasible point $(x, y)$ that satisfies the following constraints:

\begin{aligned} c^T x &= y^Tb & \\ Ax &\le b & \\ x &\ge 0 & \\ y^TA &\ge c^T & \\ y &\ge 0 \ \ & \end{aligned}

By weak duality, if such point exists, then $x$ is the optimal solution for the primal program. (Note that, if there is no feasible point for the third program, the primal could be either infeasible or unbounded.)

Ellipsoid

Ellipsoid is a unit ball that has been squashed or stretched along orthogonal directions. We illustrate this by an example in two dimension. The following picture shows an ellipse obtained by enlarging the unit ball by a facto of $2$ along the $v_1$ direction and squashing it ball by half along the $v_2$ direction, where $v_1$ and $v_2$ constitute an orthogonal base.

It is equivalent to first squash or stretch the ball along the directions of $e_1, e_2$ ,

then rotate the resulting ellipsoid to corresponding orthogonal directions.

The rotation can be done by left multiplying an orthogonal matrix

V = [v_1, v_2]

where $v_1, v_2$ are the orthogonal base for which the squash and stretch are performed.

To recover the unit ball, we can first perform the reverse rotation:

x \leftarrow V^T x

followed by a reverse scaling:

x \leftarrow \begin{bmatrix} 1/2 & \\ & 2 \end{bmatrix} x

Note that the reverse rotation is represented by matrix $V^T$ as $V V^T = I$ . Therefore, the ellipsoid is denoted as

\left( \begin{bmatrix} 1/2 & \\ & 2 \end{bmatrix} V^T x \right)^T \left( \begin{bmatrix} 1/2 & \\ & 2 \end{bmatrix} V^T x \right) \le 1

That is

x^T \left( V \begin{bmatrix} 1/2 & \\ & 2 \end{bmatrix}^2 V^T \right) x \le 1

In general, an ellipsoid is represented by

E(s, P) = \{ x : (x - s)^T P^{-1} (x - s) \le 1\}

where $s \in \mathbb{R}^n$ is the center of the ellipsoid and $P$ is a positive semi-definite matrix.

Recall that a positive semi-definite matrix $P$ can be rewritten as

P = V \Sigma V^T

and

P^{-1} = V \Sigma^{-1} V^T

where $V$ is an orthogonal matrix and $\Sigma$ is a diagonal matrix. Geometrically, $V$ are the directions and $\sqrt \Sigma$ are the ratios of the squashes or stretches. If we define $B = V \sqrt \Sigma$ , then $B^T = \sqrt \Sigma V^T$ , $B^{-1} = \sqrt \Sigma^{-1} V^T$ and $P = B B^T$ , $P^{-1} = (B^{-1} )^T B^{-1}$ .

The representation can be viewed as merely a recipe for converting the ellipsoid back to the unit ball:

First perform a translation $x \leftarrow (x - s)$ so that the center of the ellipsoid is now the center.
Rotate the ellipsoid, such that its axes are aligned with $e_1, e_2, ..., e_n$ . This is done by multiplying the matrix $V^T$ : $x \leftarrow V^T x$ .
Squash or stretch along the axes $e_1, e_2, ..., e_n$ : $x \leftarrow \sqrt \Sigma^{-1} x$ .
Therefore, the linear transformation for the recovery is

x \leftarrow B^{-1} (x - s)

Conversely, $E(s, P) = E(s, V \sqrt \Sigma \sqrt \Sigma V^T)$ tells us how to get the ellipsoid from the unit ball:

Scale unit ball along the axes $e_1, e_2, .., e_n$ by factors that $\sqrt \Sigma_{1, 1}, \sqrt \Sigma_{2, 2}, ..., \sqrt \Sigma_{n, n}$ .
Rotate the ellipsoid by matrix $V$ .
Move the ellipsoid by $s$ .

That is, $x \leftarrow V \sqrt \Sigma y + s$ takes the unit ball $\{ y : y^T y \le 1 \}$ to the corresponding ellipsoid $\{ x : (x - s)^T V \Sigma^{-1} V^T (x - s) \le 1 \}$ .

Goal

Now return to our discussion of ellipsoid method. Without lose of generality, we consider only the following problem:

Given $F = \{x : Ax \le b\}$ , output any $x \in F$ or report $F = \emptyset$ , where $A \in \mathbb{R}^{m \times n}, x \in \mathbb{R}^n, b \in \mathbb{R}^m$ .

We assume further that

$\exists R \in \mathbb{R}$ , $F \subseteq B(0, R)$ . Besides, if $F \neq \emptyset$ , $\exists r \in \mathbb{R}$ , and $t \in \mathbb{R}^n$ , s.t., $B(t, r) \subseteq F$ . Here $B(0, R) = \{ x : \Vert x\Vert_2 \le R \}$ and $B(t, r) = \{ x : \Vert x - t\Vert_2 \le r \}$ .
Separation Oracle: Given a vertex $s \in \mathbb{R}^n$ , the separation oracle is a mechanism that asserts
- either $s \in F$
- or $s \notin F$ and returns a direction (which is a unit vector) $d \in \mathbb{R}^n$ , such that $d^T s \ge d^T x$ for $\forall x \in F$ .

The first assumption says that $F$ is neither too "large" nor too "small". Note that $B(t, r) \subseteq F$ also implies that $F$ has full dimension.

For LP, the separation is straightforward to implement: we just check whether $As \le b$ . If so, then $s \in F$ . Otherwise, $\exists i \in [m]$ , such that the $i^{th}$ constraint is violated: $(A s)_i > b_i$ . Just return $d^T = A_{i, \cdot}$ .

Question to ponder: the feasible region of an LP may not satisfy the first assumption? It is possible to reduce it to one satisfying the first assumption?

The Algorithm

The algorithm begins with a special ellipsoid $E_0 = B(0, R)$ . At each step, it queries the separation oracle and generates a smaller ellipsoid (with reduced volume). If $F$ is not empty, the ellipsoid at each step is guaranteed to contain it. The algorithm stops before the volume of the ellipsoid decreases below $ \mathtt{vol} (B(t, r))$.

The algorithmic flow is shown below:

$i = 0$ .
$E_0 = B(0, R)$ .
While $\mathtt{vol} (E_i(s_i, P_i)) \ge \mathtt{vol} (B(t, r))$ :
$\qquad$ if $s_i \in F$ , output $s_i$ .
$\qquad$ else,
$\qquad$ $\qquad$ Obtain a direction $d_i$ , s.t., $d_i^T s_i \ge d_i^T x$ for $x \in F$ *.
$\qquad$ $\qquad$ Calculate a ellipsoid $E_{i + 1} \supset E_{i} \cap \{x : d_i^T s_i \ge d_i^T x \}$ .
$\qquad$ $\qquad$ $i \leftarrow i + 1$ .
Output $F = \emptyset$ .

For the one dimension case, it is nothing more than a variant of binary search.

Bounding the Number of Iterations

The main result that characterizes the running time is

Theorem. It is guaranteed that for each $i$ ,
$\frac{ \mathtt{vol} ( E_{i + 1} ) } { \mathtt{vol} (E_i) } \le \exp \left( -\frac{1}{2 (n + 1) } \right).$

Corollary. The number of iterations of the ellipsoid method is bounded by
$2n(n + 1) \ln \frac{R}{r}.$

To simplify the notation further, we write $E_i = E(s, P)$ , $d_i = d$ and $E_{i + 1} = E(s', P')$ . The key lemma for the theorem is:

Lemma. Given an ellipsoid $E(s, P) = \{ x \in \mathbb{R}^n : (x - s)^T P^{-1} (x - s) \le 1\}$ and the half space $H(s, d) = \{ x \in \mathbb{R}^n : d^T (x - s)\le 0 \}$ , we can find an ellipsoid $E(s', P')$ , such that

$E(s', P')$ contains the intersection of $E(s, P)$ and $H(s, d)$ :

$E(c, P) \cap H(s, d).$

The volume of $E(s', P')$ is at most $\exp \left( -\frac{1}{2 (n + 1) } \right)$ times that of $E(s, P)$ .

In particular, $s'$ and $P'$ is given by

\begin{aligned} s' &= s - \frac{1}{n + 1} \frac{P d}{\sqrt{d^T P d} } \\ P' &= \frac{n^2}{n^2 - 1} \left( P - \frac{2}{n + 1} \frac{Pd}{\sqrt{d^T P d} } \left( \frac{Pd}{\sqrt{d^T P d} } \right)^T \right) \end{aligned}

Proof. The transformation $y = B^{-1} (x - s)$ maps $E(s, P)$ to $E(0，I)$ , and

\{x \in R^n: d^T B (B^{-1} (x - s)) \le 0\}

H(0, B^T d) = \{y \in R^n: d^T B y \le 0\}.

Let

e' = \frac{B^T d}{\sqrt{d^T B B^T d} } = \frac{B^T d}{\sqrt{d^T P d} }

so that $H(0, B^Td)$ is equivalent to $H(0, e')$ .

It suffices to find an ellipsoid that contains $E(0, I) \cap H(0, e')$ . Then applying the reverse transformation $x = By + s$ will give us an ellipsoid that contains $E(s, P) \cap H(s, d)$ .

We can further reduce the problem to finding an ellipsoid that covers $E(0, I) \cap H(0, e_1)$ as we can recover the ellipsoid by rotating $e_1$ to $e'$ .

Claim. The ellipsoid that covers $E(0,I) \cap H(0, e_1)$ is given by
$E \left( -\frac{1}{n + 1}e_1, \frac{n^2}{n^2 - 1} \left( I - \frac{2}{n + 1}e_1 e_1^T \right) \right)$
and has minimum volume $\exp ( -\frac{1}{2(n + 1)} )$ .

To under the ellipsoid, we check how it transforms the unit ball to an ellipsoid:

It squashes the ball along the $e_1$ axis, by a factor of $\sqrt{ \frac{n^2} {n^2 - 1} \frac{n - 1} {n + 1} } = \frac{n}{ n + 1 }$ .
Then it stretches all other directions $e_2, e_3, ..., e_n$ by a factor of $\sqrt \frac{n^2}{n^2 - 1}$ .
Finally, move the center of the ellipsoid to $-\frac{1}{n + 1}e_1$ .

Note that the ellipsoid has volume

\begin{aligned} \left(1 - \frac{1}{n + 1} \right) \left( 1 + \frac{1}{n^2 -1} \right)^{(n - 1) / 2} &\le \exp \left( - \frac{1}{n + 1} + \frac{1}{n^2 -1} (n - 1) / 2 \right) \\ &= \exp \left( - \frac{1}{2(n + 1)} \right). \end{aligned}

This is also the ratio of the ellipsoid over that of the unit ball, whose value is 1.

Before we prove the claim, we show how to transform the ellipsoid back to $E(s', P')$ that covers $E(s, P) \cap H(s, d)$ . First, by symmetry, the ellipsoid that covers $E(0, I) \cap H(0, e')$ is given by

E \left( -\frac{1}{n + 1}e', \frac{n^2}{n^2 - 1} \left( I - \frac{2}{n + 1}e' (e')^T \right) \right)

Applying the mapping $B y + s = x$ :

E \left( -\frac{1}{n + 1} e', \frac{n^2}{n^2 - 1} \left( I - \frac{2}{n + 1}e' (e')^T \right) \right) \\ \longrightarrow E \left( s - \frac{1}{n + 1} B^T e', \frac{n^2}{n^2 - 1} B \left( I - \frac{2}{n + 1}e' (e')^T \right) B^T \right).

As linear transformation preserves the ratio of volumes, the ratio of the volume of the final ellipsoid over that of $E(s, P)$ is also $\exp(- \frac{1}{2 (n + 1)} )$ .

Question to ponder:

Show that $\frac{n^2}{n^2 - 1} (I - \frac{2}{n + 1}e_1 e_1^T)$ is positive definite.

Show that $\left( \frac{n^2}{n^2 - 1} (I - \frac{2}{n + 1}e_1 e_1^T) \right)^{-1} = \frac{n^2 - 1}{n^2} I + \frac{2n + 2}{n^2} e_1 e_1^T$ .

Show that $\frac{n^2}{n^2 - 1} B (I - \frac{2}{n + 1}e' (e')^T) B^T$ is positive definite.

Proof of the Claim. As both $E(0,I)$ and $H(0, e_1)$ are symmetric to $e_1$ , so should the ellipsoid covering the intersection of them. Therefore it should be of the form

a_1 (x_1 - \lambda) + \sum_{i = 2}^n a_i x_i^2 \le 1

where $a_i$ 's and $\lambda$ ( $0 < \lambda < 1$ ) are the parameters to be determined. The goal is to minimize its volume:

\sqrt{\frac{1}{a_1} \prod_{i = 2}^n \frac{1}{a_i} } \mathtt{vol} [ E(0, I) ].

where $ \mathtt{vol} [ E(0, I) ] = 1$ is the volume of unit ball in $R^n$ . It is equivalent to maximize

a_1 \prod_{i = 2}^n a_i.

As the ellipsoid covers $E(0,I) \cap H(0, e_1)$ and should has as small volume as possible, the points $e_i$ ( $i = 1, 2, ..., n$ ) should be on its boundary. Therefore,

\begin{aligned} a_1 (1 - \lambda)^2 + \sum_{i = 2}^n a_i 0^2 &= 1 &\rightarrow a_1(1 - \lambda)^2 = 1 &\\ a_1 (0 - \lambda)^2 + a_i 1^2 &= 1 &\rightarrow a_i = 1 - a_1 \lambda^2 &\quad \forall i = 2, 3, ..., n \\ \end{aligned}

Replacing $a_i$ with $1 - a_1 \lambda^2$ and then $a_1$ with $\frac{1}{(1 - \lambda)^2}$ , the goal becomes

\max \quad \frac{1}{(1 - \lambda)^2} \left[ 1 - \frac{\lambda^2}{(1 - \lambda)^2} \right]^{n - 1}

Substitute $\frac{1}{1 - \lambda}$ by a variable $x > 1$ , and let

f \doteq x^2 [1 - (x - 1)^2]^{n - 1}

Its derivative is

f' = 2x [1 - (x - 1)^2]^{n - 1} - 2 (n - 1)(x - 1) x^2 [1 - (x - 1)^2]^{n - 2}

Setting the derivative to zero, we obtain

\begin{aligned} &[1 - (x - 1)^2] - (n - 1)(x - 1)x = 0 \\ &\Rightarrow 1 - x^2 + 2x - 1- (n -1 )x^2 + (n -1 )x = 0 \\ &\Rightarrow - n x^2 + (n + 1) x = 0 \\ &\Rightarrow x = \frac{n + 1}{n} \end{aligned}

It follows that $a_1 = (\frac{n + 1}{n})^2$ , $\lambda = \frac{1}{n + 1}$ and $a_i = \frac{n^2 - 1}{n^2}$ for $\forall i = 2, 3, ..., n$ .

Reference

[1]. Michel X. Goemans, 7. Lecture notes on the ellipsoid algorithm, 18.433: Combinatorial Optimization, MIT