% Revision: $Id$ ($Date$)
%
\documentclass[a4paper,twocolumn]{article}

%\usepackage{html}
\usepackage[dvips]{graphics,color,epsfig}
\usepackage{amssymb}
\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage{amsthm}
\usepackage{psfrag}

\newcommand{\N}{\ensuremath{\mathbb{N}}}
\newcommand{\F}{\ensuremath{\mathbb{F}}}
\newcommand{\Z}{\ensuremath{\mathbb{Z}}}
\newcommand{\D}{\ensuremath{T}}%\mathbb{D}}}
\newcommand{\R}{\ensuremath{\mathbb{R}}}
\newcommand{\Q}{\ensuremath{\mathbb{Q}}}
\newcommand{\C}{\ensuremath{\mathbb{C}}}
\newcommand{\M}{\ensuremath{\mathbb{M}}}
\newcommand{\X}{\ensuremath{\mathbb{X}}}
\newcommand{\E}{\ensuremath{\mathrm{E}}}

\newcommand{\zero}{\ensuremath{\mathbf{0}}}
\newcommand{\vol}{\ensuremath{\mathop{\rm vol}\nolimits}}
\newcommand{\lcm}{\ensuremath{\mathop{\rm lcm}\nolimits}}
\newcommand{\relint}{\ensuremath{\mathop{\rm relint}\nolimits}}
\newcommand{\diag}{\ensuremath{\mathop{\rm diag}\nolimits}}
\newcommand{\kernel}{\ensuremath{\mathop{\rm ker}\nolimits}}
\newcommand{\dom}{\ensuremath{\mathop{\rm dom}\nolimits}}
\newcommand{\rank}{\ensuremath{\mathop{\rm rank}\nolimits}}
\newcommand{\den}{\ensuremath{\mathop{\rm den}\nolimits}}
\newcommand{\im}{\ensuremath{\mathop{\rm im}\nolimits}}
\newcommand{\ulp}{\ensuremath{\mathop{\rm ulp}\nolimits}}
\newcommand{\ord}{\ensuremath{\mathop{\rm ord}\nolimits}}
\newcommand{\rnd}{\ensuremath{\mathop{\rm rnd}\nolimits}}
\newcommand{\eps}{\ensuremath{\mathop{\rm eps}\nolimits}}
\newcommand{\pr}{\ensuremath{\mathop{\rm Prob}\nolimits}}
\newcommand{\divi}{\ensuremath{\mathop{\rm div}\nolimits}}
\newcommand{\one}{\ensuremath{\mathop{\rm 1}\nolimits}}
\newcommand{\asso}{\ensuremath{\,\|\,}}
\newcommand{\lc}{\ensuremath{\mathop{\rm lc}\nolimits}}
\newcommand{\lm}{\ensuremath{\mathop{\rm lm}\nolimits}}
\newcommand{\rem}{\ensuremath{\mathop{\rm rem}\nolimits}}
\newcommand{\prem}{\ensuremath{\mathop{\rm prem}\nolimits}}
\newcommand{\quo}{\ensuremath{\mathop{\rm quo}\nolimits}}
\newcommand{\pquo}{\ensuremath{\mathop{\rm pquo}\nolimits}}
\newcommand{\spann}{\ensuremath{\mathop{\rm span}\nolimits}}
\newcommand{\QR}{\ensuremath{\mathop{\rm QR}\nolimits}}
\newcommand{\QNR}{\ensuremath{\mathop{\rm QNR}\nolimits}}
\newcommand{\aff}{\ensuremath{\mathop{\rm aff}\nolimits}}
\newcommand{\assoz}{\ensuremath{\,\|\,}}
\newcommand{\grad}{\ensuremath{\nabla}}%\mathop{\rm \underbar{grad}}\nolimits}}
\newcommand{\Hess}{\ensuremath{\mathop{\rm \underbar{Hess}}\nolimits}}
\newcommand{\trace}{\ensuremath{\mathop{\rm Tr}\nolimits}}
\newcommand{\conv}{\ensuremath{\mathop{\rm conv}\nolimits}}
\newcommand{\minp}{\ensuremath{\mathop{\rm mp}\nolimits}}
\newcommand{\argmax}{\ensuremath{\mathop{\rm argmax}\nolimits}}
\newcommand{\excess}{\ensuremath{\mathop{\rm ex}\nolimits}}

\newcommand{\MEL}{\ensuremath{\mathop{\textnormal{\textsc{mel}}}\nolimits}}
\newcommand{\ov}[1]{\overline{#1}}

\newtheorem{theorem}{Theorem}
\newtheorem{alg}[theorem]{Algorithm}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{dilemma}[theorem]{Dilemma}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{proposition}[theorem]{Proposition}

\bibliographystyle{plain}

\begin{document}
\title{Appoximate smallest enclosing ellipse---an \\
implementation for {\sc Cgal}}
\author{Kaspar Fischer \\ ETH Z\"urich}
\maketitle

\begin{center}
  \begin{sc}
    Warning: this document does not fully apply to the implementation
    in Cgal; for instance, the Cgal implementation does not run the
    algorithm with exact arithmetic (as is ``documented'' here).
  \end{sc}
\end{center}

\begin{abstract}
  This report discusses an implementation of Khachiyan's algorithm
  \cite{k-rprnmc-96} for computing an $(1+\epsilon)$-approximation to
  the smallest enclosing ellipsoid $\MEL(S)$ of a point set $S$ in
  $d$-dimensional space.

  We <<have>> implemented the algorithm in C++ in a generic way such that
  it can either be run using floating-point (i.e., double) or exact
  arithmetic.  (Exact arithmetic means here that all intermediate
  operations are free of rounding errors.)  While both approaches seem
  to work in practice for sufficiently large $\epsilon$ (i.e.,
  $\epsilon>10^{-4}$), floating-point arithmetic cannot safely be used
  for very small $\epsilon$ because numerical errors might detain the
  algorithm from terminating.  (However, the algorithm has the very
  nice property that given a floating-point solution $E$ we can easily
  determine the \emph{exact} $\varepsilon$ for which the ellipsoid $E$ is a
  $(1+\epsilon)$-approximation.)

  Since our experminents show that Khachiyan's algorithm is very slow
  when used with exact arithmetic on rational input, we formulate the
  underlying operations in such a way that the growth of the numbers
  in intermediate results is lowered somewhat (without resorting to
  gcd-computations).  Using the resulting algorithm as a subroutine of
  a simple heuristic to find a $(1+\epsilon)$-approximation to
  $\MEL(S)$, we can solve instances of $10,\!000$ points in 3-space in
  less than ... seconds on a modern PC.
\end{abstract}

\section{Introduction}

The rest of this report is organized as follows.  In
Section~\ref{sec:khachiyan} we briefly describe Khachiyan's algorithm
for computing an approximation to the smallest enclosing ellipsoid of
a \emph{centrally symmetric} point set and mention how this helps in
finding an approximation to $\MEL(S)$ for arbitrary input points.  In
Section~\ref{sec:exact} we revise the algorithm's update step:
assuming rational input points, we will be able to store all
intermediate results $\alpha/\beta$ as pairs $(\alpha,\beta)$, where
$\alpha$ and $\beta$ only have ``intrinsic'' common factors, the
``avoidable'' common factors being removed \emph{without}
gcd-computations.  Based on this revised algorithm we implement a
simple heuristic (similar to the one in~\cite{kmy-ccsasehhd-03}) to
find a $(1+\epsilon)$-approximation to $\MEL(S)$.  Finally, we
discusses our experimental findings in Section~\ref{sec:exp}.

\section{Khachiyan's algorithm for centrally symmetric input}
\label{sec:khachiyan}

We start with the case where the input points $P\subset \R^d$ are
\emph{centrally symmetric}, i.e., where $P=-P$.  Khachiyan's algorithm
is based on the following program in the variables
$x=(x_1,\ldots,x_n)$:
\begin{equation}
  \label{eq:dual}
  \begin{array}{lll}
    \mbox{(D)} & \mbox{max}  & \ln \det(\sum_{i=1}^n x_i p_i p_i^T) \\
    & \mbox{s.t.} & \sum_{i=1}^n x_i = 1 \\
    &             & x_1,\ldots,x_n \ge 0
  \end{array}
\end{equation}
The program helps in finding (an approximation to) the minellipse
because for any solution $x$ of (D), the matrix $M(x):= \sum_{i=1}^n
x_i p_i p_i^T$ defines an ellipsoid
\begin{equation}
  \label{eq:sol_ell}
  E_x:= \{ x \in \R^d \mid x^T M(x)^{-1} x \le 1 \}
\end{equation}
with the following nice properties.

\begin{lemma}[Khachiyan \cite{k-rprnmc-96}]
  \label{thm:rounding}
  Let $P\subset\R^d$ be a set of $n$ points with $\aff(P)=\R^d$.  If
  for some $\epsilon\ge 0$ and all $j\in\{1,\ldots,n\}$
  \begin{equation}
    \label{eq:relaxed_conds}
    \excess_j(x):= p_j^T M(x)^{-1} p_j \le (1+\epsilon)\,d
  \end{equation}
  then the ellipsoid $E_x$ satisfies 
  \begin{equation}
    \label{eq:rounding}
    E_x \subseteq \conv(P) \subseteq \sqrt{(1+\epsilon)\,d}\,E_x
  \end{equation}
  and $\vol(\sqrt{(1+\epsilon)\,d}\,E_x) \le (1+\epsilon)^{d/2}
  \vol(\MEL(P))$.
\end{lemma}
%
In other words, if you have any solution $x$ to program (D), and
$\epsilon$ is a large enough number such that the equations
\eqref{eq:relaxed_conds} hold, then you know that (i) $E_x$ is an
ellipsoid inscribed in $\conv(P)$, (ii) that $E'_x:=
\sqrt{(1+\epsilon)\,d}\,E_x$ is circumscribing and (iii) that $E'_x$
is a $(1+e)^{d/2}$-approximation of the smallest enclosing ellipsoid
of~$P$.

Khachiyan's algorithm exploits the above lemma by maintaining the
conditions \eqref{eq:relaxed_conds} for some progressively smaller
$\epsilon_0>\epsilon_1>\ldots$ until the desired approximation ratio
is achieved.

To start with, we choose the feasible solution
\begin{equation}
  \label{eq:feasible}
  x^{(0)} := (1/n,\ldots,1/n)
\end{equation}
and iteratively update $x^{(k)}$ to a feasible solution $x^{(k+1)}$
which is better in the sense that it has a higher objective value %
$w(x^{(k+1)})$, for $w(x):= \ln \det M(x)$.

\begin{lemma}
  \label{lemma:start}
  The start value $x^{(0)}$ represents a proper ellipsoid, i.e.,
  $w(x^{(0)})>-\infty$.
\end{lemma}

%% \begin{proof}
%%   The matrix $M(x^{(0)}) = 1/n\, \sum_j p_j p_j^T$ is positive
%%   semidefinite; it thus suffices to show that its determinant doesn't
%%   vanish.  So assume for a contradiction that there are coefficients
%%   $\lambda_i$, not all zero, such that $\sum_i \lambda_i m_i = 0$, the
%%   $m_i$ being the columns of $M(x^{(0)})$.  Then
%%   \begin{eqnarray*}
%%     0 &=& \sum_{i=1}^d \lambda_i m_i \\
%%       &=& \sum_{i=1}^d \lambda_i \frac{1}{n} \sum_{j=1}^n p_j p_{ij} \\
%%       &=& \sum_{j=1}^n p_j
%%           \left( \frac{1}{n}\sum_{i=1}^d \lambda_i p_{ij} \right)
%%   \end{eqnarray*}
%% \end{proof}

\begin{lemma}
  Let $x$ be a solution to (D) satisfying the $j$th optimality
  condition \eqref{eq:relaxed_conds} with equality (i.e.,
  $\excess_j(x)=(1+\epsilon)\,d$ for some $\epsilon$).  Then
  the objective function in direction $e_j$
  \begin{equation}
    \label{eq:lin_search}
    w_j(\lambda):= w((1-\tau)\,x + \tau e_j)
  \end{equation}
  attains its maximum for $\tau = \epsilon / (\excess_j(x) - 1)$.
\end{lemma}

The following lemma proves lemma~\ref{lemma:start} since our resulting
ellipsoid always contains the origin (so that the input points
together with the origin affinely span the whole space).
%
\begin{lemma}
   Let $p_1,\ldots,p_n\in\R^d$ be a set of $n\ge d+1$ points.  Then
   the $p_i$ linearly span $\R^d$ if and only if $\det(\sum_{i=1}^n
   p_i p_i^T)>0$.
\end{lemma}
%
\begin{proof}
  For convenience, we set $P = [p_1,\ldots,p_n]$ and $Q=\sum_{i=1}^n
  p_i p_i^T$ for the duration of this proof.  Observe first that $Q$
  is positive semidefinite and hence $\det Q \ge 0$.  So it suffices
  to prove that the points $P$ affinely span $\R^d$ if and only if
  $\det Q \not=0$.
  
  The $p_i$ spanning the whole space is equivalent to $Px=q$ having a
  solution $x$ for every given $q\in\R^d$.  This, however, is
  equivalent to $P$ having full rank.  It follows that $PP^T$ is
  regular\footnote{To see this, write $P^T = QR$ for an orthogonal
  $(n\times n)$-matrix $Q$ and a upper triangular $(d\times d)$-matrix
  $R$ of the form
  \[ R = \left[
  \begin{array}{c}
  R_0 \\ 0
  \end{array}
  \right].
  \]
  We must have $\det R_0\not=0$, because $P$ has full rank.  It
  follows that $PP^T = R^TQ^T Q R = R^T R$ and clearly,
  $\det(PP^T)\not=0$.} and from $PP^T = \sum_{i=1}^n p_i p_i^T$ the
  claim follows.
\end{proof}

\section{The update step revised}
In each iteration of the algorithm we need to determine the
``direction'' $e_j$, $j\in\{1,\ldots,n\}$, which violates the
optimality conditions \eqref{eq:relaxed_conds} most.  Since for this
we need the inverse of the matrix $M_k:= M(x^{(k)})$, we maintain
$M_k^{-1}$ explicitly in the code.  As Khachiyan suggests in his
paper, $M_{k+1}^{-1}$ can easily be obtained from $M_k^{-1}$.

\begin{lemma}
  \label{lm:rank-1}
  Suppose we know $M(x)^{-1}$ for some feasible solution $x$ to (D)
  (and the matrix exists).  If
  \[
  x' = (1-\tau)\, x + \tau e_k,
  \]
  where $\tau = \epsilon/((1+\epsilon)d - 1)$ and $p_k^T M(x)^{-1} p_k
  = (1+\epsilon)d$ for some $\epsilon\in\R$ then the inverse of
  $M(x')$ reads
  \[
    \Big(
      1+\frac{\epsilon}{(d-1)(1+\epsilon)}
    \Big) M(x)^{-1} 
    -
    \frac{\epsilon}{(d-1)(1+\epsilon)^2} l l^T,
  \]
  where $l:= M(x)^{-1} p_k$.
\end{lemma}

\section{Using exact arithmetic}

In this section, we want to formulate Khachiyan's algorithm in a way
which is better suited for computation with exact numbers (i.e., not
with double arithmetic).  Our experiments have shown that running
Khachiyan's algorithm with an exact number type (as provided by the
{\sc gnu mp} package) results in enormously high running times; the
reason for this is not really clear to us, either it is the numbers
growing very fast or the many gcd-computations done by the abritrary
precision number type.  We try to overcome (?) these problems in the
following.

For this, we assume that our input points have coordinates from an
integral domain $T$ (i.e., a ring in which $xy=0$ implies that $x$ or
$y$ is zero).  The matrix $M:= M(x)$ has then entries in $T$ and its
inverse has, by Cramer's rule, rational entries over $T$ with common
denominator $\det M$.
%
\begin{lemma}
  Let $A'=\alpha A + \beta uu^T$ for some regular symmetric matrix
  $A\in T^{d\times d}$, a vector $u\in T^d$, and $\alpha,\beta\in
  T$.  Then $\det A' = \alpha^{d-1}(\alpha+\beta u^Tv) \det A$, and
  \begin{equation}
    \label{eq:exrank1inverse}
    \frac{1}{\alpha} A^{-1} -
    \frac{\beta}{\alpha\,(\alpha+\beta u^T v)}vv^T
  \end{equation}
  is, for $v = A^{-1} u$, the inverse of $A'$ (if it exists).
\end{lemma}
%
\begin{proof}
  For the determinant formula, you walk to room B-forty-someting and
  ask Bernd G\"artner.

  The given update formula for the inverse of $A'$ is easily verified
  by multiplying out the product $A' J$, $J$ being the
  matrix~\eqref{eq:exrank1inverse} from the claim.
\end{proof}
%
For the following lemma we recall from Cramer's Rule that the inverse
of $A$ can be written as $A^{-1} = \hat{A} / \det{A}$ for some matrix
$\hat{A}$ with entries from $T$ (actually, $\hat{A} = A^{-1}
\det{A}$).
%
\begin{lemma}
  \label{lemma:rank1update}
  In the context of the previous lemma, the inverse of $A'$ reads
  \begin{equation}
    \label{eq:exrank1inv_sep}
    \frac{(\delta \hat{A} - \beta\hat{v}\hat{v}^T)/\det{A}}%
	 {\delta\alpha}
  \end{equation}
  where $\hat{v} = \hat{A} u$ and $\delta = \alpha \det{A} + \beta
  u^T\hat{v}\in T$.  Further, the division in the nominator of
  \eqref{eq:exrank1inv_sep} is proper in the integral domain~$T$ and
  $\det{A'} = \alpha^{d-1} \delta$.
\end{lemma}
%
\begin{proof}
   From the previous lemma we know that
   \[ \det{A'} = \alpha^{d-1} (\alpha\det{A}+\beta u^T\hat{v})
               = \alpha^{d-1} \delta.
   \]
   so the formular for the determinant holds.  Also, the lemma tells
   us that with $v = A^{-1} u$,
   \begin{eqnarray*}
   A'^{-1} &=& \frac{\hat{A}}{\alpha \det{A}} -
               \frac{\beta vv^T}{\alpha\,(\alpha+\beta u^T v)} \\
           &=& \frac{(\alpha + \beta u^T v) \hat{A} - \beta vv^T\det{A}}%
               {\alpha\,(\alpha+\beta u^T v) \det{A}} \\
           &=& \frac{(\alpha + \beta u^T v) \hat{A} -
               \beta vv^T\det{A}}{\delta\alpha} \\
	   &=& \frac{(\delta \hat{A} - \beta
	       \hat{v}\hat{v}^T)/\det{A}}{\delta\alpha}.
   \end{eqnarray*}
   On the other hand, we know from \eqref{eq:exrank1inverse} that the
   entries of $A'^{-1}$ can be written as rationals over $T$ with
   common denominator $\delta\alpha$; it follows that the division by
   $\det{A}$ in the above nominator is without remainder.
\end{proof}
%
%
\paragraph{The revised update step.}
The situation is like in Lemma~\ref{lm:rank-1}; we have already at
hand the matrix $M(x)^{-1}$ for some feasible solution $x\in Q(\D)^d$
to (D) and we want to find the matrix $M(x')^{-1}$ where
\[
x' = (1-\tau) x+ \tau e_k \quad \mbox{for} \quad
\tau = \frac{\epsilon}{(1+\epsilon)d-1}.
\]
We know that $k=\argmax_{1\le j\le n} p_j^T M(x)^{-1} p_j$, and $\epsilon$
satisfies the relaxed optimality conditions with equality, i.e.,
$p_k^T M(x)^{-1} p_k = (1+\epsilon)d$.  For convenience, we set $m:=
p_k M(x)^{-1} p_k$.

Observe first that $\tau$ is a rational over $\D$ and therefore the
new $x'$ need not lie in $\D$.  So although the algorithm starts with
$x^{(0)}\in\D^n$, the subsequent $x^{(i)}$ are rationals over $T$, and
consequently, the entries of $M(x)$ are fractions, too.  Therefore, we
represent $M(x)$ as
\[
M(x) = \frac{N}{\nu}
\quad \mbox{for $N\in \D^{d\times d}$, $\nu\in T$}.
\]
By Cramer's rule, its inverse is then of the form $M(x)^{-1} = \nu
\hat{N}/\det{N}$ for some matrix $\hat{N}\in \D^{d\times d}$
(actually, $\hat{N} = N^{-1}\det{N}$).  This means that in general, we
simply cannot avoid a denominator of magnitude $|\det{N}|$, even if we
actually preferred a smaller value in order to avoid growth of
numbers.  Luckily, we will be able to use a smaller denominator than
$\det{N}$, and for the moment we just call it $\den(N)$; we will have
$\det{N} = \alpha \den{N}$ for some $\alpha\in T$.  So our inverse is
of the form $M(x)^{-1} = \nu\hat{N}/\den{N}$ for a number $\den{N}\in
T$ such that $\hat{N}:= N^{-1}\den{N}$ has entries in $T$ only.  And
since Khachiyan's algorithm only needs this latter inverse, we thus
store the tuple $(\hat{N},\nu,\den{N},\det{N},\alpha)$ in order to
represent $M(x)^{-1}$.

By definition of the function $M(x)$ we have
\[ 
M(x') = \sum_{i=1}^n x'_i p_ip_i^T = (1-\tau) M(x) + \tau p_kp_k^T.
\]
So $M(x')$ is obtained from $M(x)$ by means of a rank-one update.
Therefore, we would like to apply Lemma~\ref{lemma:rank1update} to get
a formula for $M(x')^{-1}$.  However, the numbers $1-\tau$ and $\tau$
are not from the integral domain~$T$ and thus we cannot immediately
invoke the lemma.  We will therefore express these two numbers as
rationals over $T$ and deal with their common denominator at a later
stage.  To this end we set $\ov{m}:= p_k \hat{N} p_k$, so that
$(1+\epsilon)d = \nu \ov{m}/\den{N}$.  From this, we get
\[
\epsilon = \frac{v\ov{m}-d\den{N}}{d \den{N}}.
\]
Plugging this into the formular for $\tau$ yields
\begin{eqnarray*}
  \tau   &=& \frac{\epsilon\den{N}}{\nu\ov{m} - \den{N}}
          =  \frac{\nu\ov{m}-d\den{N}}{d\,(\nu\ov{m} - \den{N})}, \\
  1-\tau &=& \frac{\nu\ov{m}\,(d-1)}{d\,(\nu\ov{m} - \den{N})}.
\end{eqnarray*}
It follows that
\begin{eqnarray*}
M(x') &=& \frac{\nu \ov{m}\,(d-1) N + \nu\,(\nu\ov{m}-d\den{N})\, p_kp_k^T}%
               {\nu d\,(\nu \ov{m}-\den{N})} \\
      &=& \frac{\ov{m}\,(d-1) N + (\nu\ov{m}-d\den{N})\, p_kp_k^T}%
               {d\,(\nu \ov{m}-\den{N})} \\
      &=:& \frac{\alpha N + \beta p_kp_k^T}{\nu'}.
\end{eqnarray*}
In order to invoke Lemma~\ref{lemma:rank1update}, we need to have
$N^{-1}\det{N} = \hat{N} \alpha$.  This gives
\[
M(x')^{-1} = \nu'\, \frac{(\alpha\delta\hat{N}-\alpha^2\beta
\hat{v}\hat{v}^T)/\det{N}}{\delta\alpha}
\]
with $\hat{v} = \hat{N}p_k$, where the division in the nominator is
proper and $\alpha^{d-1}\delta$ is the determinant of $N':= \alpha N + \beta
p_kp_k^T$.  Further, the lemma tells us that
\begin{eqnarray*}
\delta &=& \alpha \det{N} + \alpha\beta p_k^T\hat{N}p_k %\\
%       &=& \ov{m}(d-1) \det{N} + (\nu\ov{m}-d\den{N})\,\ov{m} \\
%       &=& \ov{m}\left((d-1) \det{N} + \nu\ov{m}-d\den{N}\right) \\
%       &=& \ov{m}\left( d\,(\det{N}-\den{N}) - \det{N} + \nu \ov{m}\right)
%       %&=& (\nu\ov{m}-\det{N})\,\ov{m}
\end{eqnarray*}
[{\sc todo:} If $\den{N}=\det{N}$ then
$\delta=(\nu\ov{m}-\det{N})\,\ov{m}$.  Compare with $\den{N}:=
\det{N}$?]

Consequently, our implementation maintains $M(x)^{-1}$ as
$\nu\hat{N}/\den{N}$, i.e., as a tuple
\[
(\hat{N},\nu,\den{N},\det{N})
\]
with the invariant that $M(x)=N/\nu$ for $N,\nu$ over $\D$ and
$\den{N}\in T$ such that $\hat{N} = N^{-1}\den{N}$ has entries in $T$
only.  In the update step from $M(x)^{-1}$ to $M(x')^{-1}$, we know
that $M(x')=N'/v'$ and thus compute (see formulas above)
\begin{eqnarray}
  \ov{m}  &=& p_k^T \hat{N} p_k, \nonumber\\
  \alpha  &=& \ov{m}\,(d-1), \label{eq:alpha}\\
  \beta   &=& \nu\ov{m}-d\den{N}, \nonumber\\
  \hat{v} &=& \hat{N}p_k, \nonumber\\
  \nu'    &=& d\,(\nu\ov{m}-\den{N}), \nonumber\\
  \delta  &=& \alpha\det{N} + \beta\ov{m}\label{eq:delta}.
\end{eqnarray}
Using these quantities we first calculate the integral matrix $\hat{N}':=
(\delta \hat{N}-\beta \hat{v}\hat{v}^T)/\det{N}$ which by the lemma
satisfies $\hat{N}' = N'^{-1}\den{N'}$ for $\den{N'}:= \delta\alpha$.
Since $\alpha^{d-1}\delta = \det{N'}$, we finish with the tuple
$(\hat{N'},\nu',\delta\alpha,\alpha^{d-1}\delta)$.

Our implementation initially chooses $\nu$, $\hat{N}$ and $\den{N}$ in
such a way that $\den{N}>0$ and $\nu>0$.  From the positive
definiteness of $\hat{N}$ and equation~\eqref{eq:alpha} we already
have $\alpha>0$.  Furthermore,
\begin{eqnarray*}
  \beta &=& \nu\ov{m} - d\den{N} \\
        &=& \den{N}\left(\frac{\nu\ov{m}}{\den{N}} - d\right) \\
        &=& \den{N} \left( (1+\epsilon) d - d \right) > 0
\end{eqnarray*}
(because $\epsilon\in(0,1)$ in Khachiyan's algorithm). Similarly,
\begin{eqnarray*}
  \nu' &=& d\den{N}\left(\frac{\nu\ov{m}}{\den{N}} - 1\right) \\
       &=& d\den{N}\left((1+\epsilon) d - 1\right) \\
       &>& d\den{N} \cdot \epsilon d > 0.
\end{eqnarray*}
Consequently, $\delta>0$ by equation~\eqref{eq:delta} hence both
$\nu'$ and $\den{N'}$ are strictly positive.

\paragraph{Updating the excesses.}
In order to speed things up, Khachiyan's algorithm maintains the
excesses $\excess_x(p_j) = p_j^T M(x)^{-1} p_j$ for all points $p_j$
explicitly (with $x$ being the current solution to program~(D)).  In
an iteration we therefore need to update the excesses of the input
points to reflect the change from $x$ to $x'$.

Our implementation stores the excess $\excess_x(p_j)$ as
$\hat{\excess}_x(p_j):= p_j^T \hat{N} p_j$, that is, we have
\begin{equation}
  \excess_x(p_j) = p_j^T \frac{\nu \hat{N}}{\den{N}} p_j
               = \frac{\nu}{\den{N}}\,\hat{\excess}_x(p_j).
\end{equation}
(In particular, since $\nu/\den{N}>0$, we can search for the maximum
among the $\hat{\excess}_x(p_j)$ when we need to find the index $j$
for which $\excess_x(p_j)$ is maximal.)

Since $\hat{N}$ updates to $\hat{N}'$ by means of $\hat{N}' =
\alpha\,(\delta\hat{N}-\alpha\beta\hat{v}\hat{v}^T)/\det{N}$ (where
the division is proper), we get
\begin{eqnarray*}
  \hat{\excess}_{x'}(p_j) &=& p_j^T \hat{N}' p_j \\
       &=& \alpha\,\frac{\delta \hat{\excess}_x(p_j) -
                         \alpha\beta (p_j^T \hat{v})^2}%
                {\det{N}},
\end{eqnarray*}
where again the divison is without remainder because we know that
$p_j^T \hat{N'} p_j=\hat{\excess}_{x'}(p_j)$ has entries in $T$ only.

\paragraph{Finding the initial inverse.}
The algorithm starts with $x:=(1/n,\ldots,1/n)$ and we have to find
$M(x)^{-1}$.  We have $M(x) = N/\nu$ with $\nu=n$ and $N =
\sum_{i=1}^n p_ip_i^T$.  What is $\hat{N}=N^{-1} \det{N}$?

We use a technique called
\emph{$Q$-pivoting}~\cite{em-nsme-97,g-ealcc-98}.  We start with the
identity matrix (the inverse of which we happen to know) and replace
one column after the other by the respective column of our desired
matrix $N$, in each step updating the inverse.

\begin{lemma}
  Let $A$ be some invertible matrix and $u$ some vector over $T$.
  Denote by $A_{(i;u)}$ the matrix $A$ with the $i$th column replaced
  by $u$.  Set $v=A^{-1} u$ and
  \begin{equation}
    \label{eq:q-pivot-w}
    w = [-v_1, \ldots, -v_{i-1}, 1, -v_{i+1}, \ldots, -v_d ]^T / v_i.
  \end{equation}
  Then $\det A_{(i;u)} = v_i \det{A} = \hat{v}_i$ and
  \[
  A_{(i;u)}^{-1} = \frac{I_{(i;w')} \hat{A}^{-1} /\det{A}}{\hat{v}_i},
  \]
  with the division in the nominator being proper.
\end{lemma}
%
\begin{proof}
  We have $A_{(i;u)} = A I_{(i;v)}$ where $I_{(i;v)}$ denotes the
  identity matrix with the $i$th row replaced by $v$.  Clearly, $\det
  I_{(i;v)} = v_i$, which shows the formula for the determinant of
  $A_{(i;v)}$.

  The inverse of $I_{(i;v)}$ is the matrix $I_{(i;w)}$ where $w$ is as
  in \eqref{eq:q-pivot-w}; this is easiliy verified by multiplying out
  $I_{(i;v)} I_{(i;w)}$.  Since $v_k = \hat{v}_k/\det{A}$ for all $k$,
  we can write the entries of $w$ as rational numbers over $T$ with
  common denominator $\hat{v}_i$; that is, we have $w = w'/\hat{v}_i$
  for
  \[
  w' = [-\hat{v}_1,\ldots,\hat{v}_{i-1},\det{A},-\hat{v}_{i+1},-\hat{v}_d]^T.
  \]
  So 
  \[
  A_{(i;u)}^{-1} = \frac{I_{(i;w')} A^{-1}}{\hat{v}_i} = 
  \frac{I_{(i;w')} \hat{A}^{-1} /\det{A}}{\hat{v}_i},
  \]
  and since the entries of $A_{(i;u)}^{-1}$ can by Cramer's rule be
  written as rationals over $T$ with common denominator
  $\det{A_{(i;u)}}$, the division in the nominator must be proper.
\end{proof}
%
\bibliography{mel}
\end{document}
