Defining and implementing spheres from sampled points

Let $p_1,\dots,p_{n+1}\in \R^n$ be points with coordinates $p_i = (p_{i,1},\dots,p_{i,n})$, and $\R^n$ with coordinates $x_1,\dots,x_n$. It is clear that if theses $n+1$ points are in general position, then they define a unique $(n-1)$-sphere in $\R^n$, on which they all lie.

Guess: Every point $(x_1,\dots,x_n)$ on the unique $(n-1)$-sphere in $\R^n$ defined by $p_1,\dots,p_{n+1}$ satisfies
\begin{equation}
\det\begin{bmatrix}
\sum_{i=1}^n x_i^2 & x_1 & x_2 & \cdots & x_n & 1 \\
\sum_{i=1}^n p_{1,i}^2 & p_{1,1} & p_{1,2} & \cdots & p_{1,n} & 1 \\
\sum_{i=1}^n p_{2,i}^2 & p_{2,1} & p_{2,2} & \cdots & p_{2,n} & 1 \\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
\sum_{i=1}^n p_{n+1,i}^2 & p_{n+1,1} & p_{n+1,2} & \cdots & p_{n+1,n} & 1 \\
\end{bmatrix} = 0.
\end{equation}
Call this equation (1).

This guess is made based on the 2-sphere version presented in Zwillinger. It is immediate that every point $p_i$ satisfies this equation, as then the matrix has two rows with identical entries. From this guess, we may conclude the following.

Proposition: The radius of the $(n-1)$-sphere defined by $p_1,\dots,p_{n+1}$ in $\R^n$ is
\[
\sqrt{\sum_{j=2}^{n+1}\frac{A_{1,j}^2}{4A_{1,1}^2} +(-1)^n\frac{A_{1,n+2}}{A_{1,1}}},
\]
for $A_{i,j}$ the $(i,j)$-minor of the matrix in equation (1). Call this expression (2).

Proof: This follows by comparing two equations. Assume that these points define a sphere of radius $r$ centered at $(a_1,\dots,a_n)$. Then points on it satisfy
\[
r^2 = (x_1-a_1)^2+\cdots+(x_n-a_n)^2 = x_1^2-2a_1x_1+a_1^2+\cdots + x_n^2-2a_nx_n + a_n^2.
\]
Equation (1) may be expanded out along the first row as
\[
\sum_{i=1}^n x_i^2 A_{1,1} - x_1A_{1,2}+\cdots + (-1)^nx_nA_{1,n+1} + (-1)^{n+1}A_{1,n+2} =0,
\]
where none of the $A_{i,j}$ are in any of the $x_i$. Dividing by the leading factor and comparing coefficients of these two equations, we find
\begin{align*}
A_{1,1} & = 1, \\
-A_{1,2}/A_{1,1} & = -2a_1, \\
& \ \ \vdots \\
(-1)^nA_{1,n+1}/A_{1,1} & = -2a_n, \\
(-1)^{n+1}A_{1,n+2}/A_{1,1} &  = a_1^2+\cdots + a_{n^2} - r^2.
\end{align*}
The given expression follows by solving for $r$. $\square$

Since any collection of $k+1$ points in general position in $\R^n$ define a $k$-plane, it is natural to ask what would be the radius of the $(k-1)$-sphere in that $k$-plane defined by those $k+1$ points. One approach to answer this is to define a new coordinate system on $\R^n$ with the first $k$ vectors spanning the given $k$-plane, restrict to the first $k$-coordinates, and apply the proposition above. More precisely, subtract the first vector from the other $k$ vectors to define a new "origin," preform the Gram--Schmidt orthogonalization process on these shifted vectors, then take the QR-decomposition of this matrix of vectors whose inverse is the map from the standard basis to the new basis. In Sage code, this may be implemented as below.

# Returns the (i,j)-minor (determinant when ith row, jth col removed) of input matrix mat
def minor(mat,i,j):
    return mat.delete_rows([i]).delete_columns([j]).det()
    
# Returns the radius of an (n-1)-sphere defined by n+1 points in R^n
def sphere_radius(L,field=CDF):
    n = len(L)-1
    M1 = [[0]*(n+2)]
    for pt in L:
        tempL = [pt*pt]
        for pos in range(n):
            tempL.append(pt[pos])
        tempL.append(1)
        M1.append(tempL)
    M2 = matrix(field,n+2,n+2,M1)
    return sqrt(reduce(lambda x,y: x+y, map(lambda z: minor(M2,0,z-1)**2/(4*minor(M2,0,0)**2),
        range(1,n+2)))+(-1)**n*minor(M2,0,n+1)/minor(M2,0,0))
    
# Returns the radius of a (k-1)-sphere defined by k+1 points in R^n
def sphere_radius_general(L,field=CDF):
    k = len(L)-1
    n = len(L[0])
    L1 = []
    for vec in L[1:]:
        L1.append(vec-L[0])    
    M = matrix(field,k,n,L1)
    Q,R = M.transpose().QR()
    L2 = [vector(field,[0]*k)]
    Qinv = Q.inverse()
    for vec in L1:
        L2.append((Qinv*vec)[:k-n])
    return sphere_radius(L2)

Now in Mathematica.

(*Returns the (i,j)-minor of a an input matrix mat*)
minor[mat_,i_,j_] := Map[Reverse,Minors[mat],{0,1}][[i]][[j]]

(*Returns the radius of an (n-1)-sphere defined by n+1 points in R^n*)
SphereRadius[L_] := Module[{n, M},
    n = Length[L]-1;
    M = Join[{Array[0#&,n+2]},Table[Join[{Sum[L[[j]][[l]]^2,{l,1,n}]},L[[j]],{1}],{j,1,n+1}]];
    Sqrt[Sum[minor[M,1,j]^2/(4*minor[M,1,1]^2),{j,2,n+1}]+(-1)^n*minor[M,1,n+2]/minor[M,1,1]]]
    
(*Returns the radius of a (k-1)-sphere defined by k+1 points in R^n*)
SphereRadiusGeneral[L_] := Module[{n,k,Lv,L1,q,qq,qinv},
   n = Length[L[[1]]];
   k = Length[L]-1;
   Lv = Table[Unique["q"],{n}];
   L1 = L[[2;;]]-Table[L[[1]],{l,1,k}];
   q = QRDecomposition[Transpose[L1]][[1]];
   qq = Join[q,{Lv}]/.Solve[{q.Lv==0,Total[#^2&/@Lv]==1},Lv][[1]];
   qinv = Inverse[Transpose[qq]];
   SphereRadius[Join[{Array[0#&,k]},#[[;;-(n-k)-1]]&/@(qinv.#&/@L1)]]]

The variable L is a list of $(n+1)$-dimensional vectors of appropriate length. Both methods skip creating the first line of the matrix in (1), since it does not appear in the expression (2). The method in Sage is probably faster in practice, but less accurate. For example:

Note that the exact square root result is approximately $10.9682$, off by around $0.01$ from the Sage result.

References: Zwillinger (CRC Standard Mathematical Tables and Formulae, Section 4.8.1)