# Simplicial Circumcenters

03 Nov 2017

Following up on the barycenter article I wrote I now want to go into computing the circumcenter of a simplex using barycentric coordinates.

### Definition of the circumcenter

The circumcenter is the center of the smallest hypersphere that contains a simplex. So for a sphere $S_{\mathbf{c}}^r$ paramterized by circumcenter $\mathbf{c}$ and circumradius $r$ we can define this closed sphere by

$\arg\min_{\mathbf{c},r} r \text{ s.t } \Delta \subset S_\mathbf{c}^r$

for a simplex $\Delta$. So long as we are computing the max minimizing $r$ here, because the radius is strictly positive we could actually minimize any monotonic increasing function of $r$. In particular $r^k$ for any positive $k$ is striclty increasing as the radius must be non-negative.

This document covers two methods for solving for this, one that discusses how to use the fact that we don’t care about evaluating the radius to compute the circumcenter itself with a change of variable that’s SPSD, and a second solution that translates the simplex to the origin to make $r = |c|$ which turns out to make the system SPD.

Given a circumcenter $\mathbf{c}$ The smallest radius that satisfies $\Delta \subset S_\mathbf{c}^r$ is given by

$r = \max_{p \in \Delta} \|p - \mathbf{c}\|.$

by convexity of distance, for a polygon this is maximized at the vertices.

$r = \max_i \|p_i - \mathbf{c}\|.$

#### Why simplices?

The reason why simplices are nice is that they have the following property: When the simplex in question is not degenerate, the circumradius is minimized when the boundary of the sphere touches the vertices of the simplex so this is equivalent to the circumcenter being equidistant to all of its vertices $\{v_i\}_{i \in \mathbb{Z}^n}$. Symbolically this means that

$r = \|v_i - \mathbf{c}\|, \forall i \in \mathbb{Z}_n..$

#### Barycentric circumcenter

By a simple triangle inequality argument it’s clear that $\mathbf{c}$ lies in the span of the simplex (otherwise we could find a sphere with smalleer radius slightly closer to the simplex). Therefore, this can be written in barycentric coordinates as

$\{c_i\} = \arg\min_c \max_i \left(\|v_i - \sum_{j=1}^n c_j v_j \|\right)^k.$

with the final circumcenter $\mathbf{c}$ computed by

$\mathbf{c} = \sum_i c_i v_i.$

### Change of Variable Solution

Begin by squaring things to make things quadratic we obtain

$r^2 = \|v_i - \mathbf{c}\|^2.$

By expanding this norm as inner products and then separating out terms we obtain

$2v_i^T\mathbf{c} + r^2 - \|\mathbf{c}\|^2 = \|v_i\|^2$

By applying a change of variable $\eta = r^2 -|\mathbf{c}|^2$ one obtains

$2v_i^T\mathbf{c} + \eta = \|v_i\|^2$

by substituting in barycentric coordinates for $c$ we obtain

$2\sum_j v_i^Tv_j c_j + \eta = \|v_i\|^2$

with an extra condition that $\sum_j c_j = 1$.

This linear system is SPSD and could be solved using Cholesky (LDLT) but always has a one dimensional kernel.

This is what is seen in PyDEC and as far as they can tell, is first derived here.

### Positive Definite Solution

We begin by moving the first vertex to the origin by $\tilde v_i = v_i - v_n$:

$r^2 = \|\tilde v_i - \sum_{j=1}^{n-1} c_j \tilde v_j \|^2, i\in \mathbb{Z}^n$

this guarantees the following relation:

$r^2 = \|\mathbf{\tilde c}\|^2.$

This also happens to “release” us of the constraint $\sum_{i=1}^n c_i = 1$ as we $c_n\tilde v_n=0$ so $c_n$ is free to take whatever value is necessary to satisfy this constraint.

The above equation can be expanded to

$r^2 = \|\tilde v_i\|^2 - \sum_{j=1}^{n-1} c_j 2\langle \tilde v_i,\tilde v_j\rangle + \|\tilde c\|^2, i\in \mathbb{Z}^{n}$

which by cancelling $r^2$ and $|\tilde {\mathbf{c}}|^2$ gives us

$\|\tilde v_i\|^2 = \sum_{j=1}^{n-1} c_j 2\langle \tilde v_i,\tilde v_j\rangle, i\in \mathbb{Z}^{n}.$

After noting that $\langle \tilde v_i, \tilde v_j \rangle$ is $0$ if $i$ or $j$ are the vertex removed we can reduce the dimension of the system. The barycentric coordinates, up to the vertex at the origin, is therefore given by

$\begin{bmatrix} \|\tilde v_1\|^2\\ \vdots\\ \|\tilde v_i\|^2\\ \vdots\\ \|\tilde v_{n-1}\|^2\\ \end{bmatrix} = 2V_{:,1:n-1}^TV_{:,1:n-1} c$

This system is SPD because the inner product matrix is clearly symmetric and the $v_i$ are assumed to not be degenerate. This lack of degenearcy implies that the inner product matrix has rank $n-1$, which finally implies that the resulting matrix has full rank $n-1$.

so this can all be implemented by

template <typename Derived>
auto circumcenter( Eigen::MatrixBase<Derived> & V) {
auto m = (V.rightCols(V.cols()-1).colwise() - V.col(0)).eval();
auto A = (2 * m.transpose() * m).eval();
auto b = m.colwise().squaredNorm().transpose().eval();
A.llt().solveInPlace(b);
return (V.col(0) + m*b).eval();
}


### A Broken Derivation that needs help

#### Augmented Lagrangian System

If we choose to use the minimization definition without assuming the existence of $r$, we can use the method of augmented Lagrangians. This formulation leads to a fairly straightforward augmented Lagrangian system:

$\mathcal{L}(\{c_i\}_{i=1}^n,\lambda) = \frac{1}{2}\sum_{i=1}^n \|v_i - \sum_{j=1}^n c_j v_j \|^2 + \lambda \left( \left[\sum_i c_i\right] - 1\right)$

where the $\lambda$ term is the constraint on the $c_i$ given by the definition of barycentric coordinates.

The partials of the above augmented Lagrangian are given by

$\partial_{c_i}\mathcal{L} = \left(v_i - \sum_{j=1}^n c_j v_j\right)^Tv_i + \lambda$

and

$\partial_{\lambda}\mathcal{L} = \sum_i c_i - 1$

The KKT conditions for this system are therefore given by

$\|v_i\|^2 = \sum_{j=1}^n v_i^Tv_jc_j + \lambda$ $\sum_i c_i = 1$

This has the same issue as the previous solution in that it has a one dimensional kernel. I derived this as my first attempt to re-derive the PyDEC code, and the minor difference between the solutions obtained implies to me that my math might just be wrong. )