A spherical code, specified by the parameters $(n,N,t)$, is a set of $N$ coordinates on the $n$-dimensional unit hypersphere such that the set of dot products between any two unit vectors from the origin is larger than or equal to $t$.

What algorithms can be used to generate candidate spherical codes? So far my google skills have found suggestions to use some version of gradient descent, which I have managed to implement for the related Thomson problem, but I have not been able to find an explicit, detailed comparison of the suitability of different algorithms to generate spherical codes. I would like to implement the algorithm myself (in Python + numpy), so pointing me at e.g. a library of solvers is not sufficient.

A list of (putatively optimal) spherical codes is maintained here. Ideally the suggested algorithm(s) should be able to produce similar codes (with some suitable definition of 'similar').

I do not have access to a good library at the moment, but if such an exposition occurs in a book a reference would also be appreciated.


3 Answers 3


From the definition it seems that the optimization problem can be written as

\begin{align} &\min_{\mathbf{x}_i}& &1\\ &\text{subject to } & &\mathbf{x}_i \cdot \mathbf{x}_j \geq t\, , i\neq j\\ & & &\Vert \mathbf{x}\Vert = 1 \end{align}

I am not sure, though.

Thompson problem

In the case of the Thompson problem, you have

\begin{align*} &\text{minimize}\quad \sum_{i=1}^{N} \sum_{j>i}^{N} \frac{1}{\Vert \mathbf{x}_i - \mathbf{x}_j \Vert^p} \\ &\text{subject to}\quad \mathbf{x}_i \in \mathbb{R}^k,\, \Vert \mathbf{x}_i\Vert = 1\quad i=1,2\cdots,N \end{align*}

That have the following gradient

\begin{equation} \nabla_i U = p\sum_{j=1, j\neq i}^{N} \frac{\mathbf{x}_j - \mathbf{x}_i}{\Vert \mathbf{x}_i - \mathbf{x}_j \Vert^{p+2}} \enspace . \end{equation}

I have used SLSQP in Python for this problem for $k=3$ up to 100 points.

Another method that worked for me was to use a penalty method, when we add a function that quantifies the violation of the constraints. The most common option is to add a quadratic function of the constraints, i.e.

$$Q(\mathbf{x};\mu) = \sum_{i=1}^{N} \sum_{j>i}^{N} \frac{1}{\Vert \mathbf{x}_i - \mathbf{x}_j \Vert^p} + \frac{\mu}{2}\sum_{i=1}^N (\Vert \mathbf{x}_i\Vert^2 - 1)^2 \enspace ,$$

with the following gradient

$$\nabla_i Q(\mathbf{x};\mu) = p \sum_{j=1, j\neq i}^{N} \frac{\mathbf{x}_j - \mathbf{x}_i}{\Vert \mathbf{x}_i - \mathbf{x}_j \Vert^{p+2}} + 2\mu (\Vert \mathbf{x}_i\Vert^2 - 1)\frac{\mathbf{x}_i}{\Vert \mathbf{x}_i\Vert} \enspace ,$$

where $\mu$ is reduced in each iteration and an unconstrained optimization algorithm is used, such as BFGS.

Edit: 2020-04-06

According to @FedericoPolini comment, I think that the optimization problem could be rewritten as

\begin{align} &\min_{\mathbf{x}_i, t}& &t\\ &\text{subject to } & &\mathbf{x}_i \cdot \mathbf{x}_j \geq t\, , i\neq j\\ & & &\Vert \mathbf{x}\Vert = 1 \end{align}

  • 2
    $\begingroup$ I guess OP wants to obtain configurations with the larger possible value for $t$, so the first equation should probably read $\min t$. This makes it more appealing to use projected gradient methods for the first problem, without introducing those rational functions. $\endgroup$ Apr 6, 2021 at 9:37
  • $\begingroup$ @FedericoPoloni, I edited my answer. Please let me know if that's what you meant. $\endgroup$
    – nicoguaro
    Apr 6, 2021 at 19:21
  • $\begingroup$ @nicoguaro Yes, that is what I meant, thanks! $\endgroup$ Apr 6, 2021 at 19:28
  • $\begingroup$ @nicoguaro I also was thinking about approaches comparable to what you have discussed, but I was just worried about the rotational symmetries of the potential and the rotational ambiguity of the optimal point configuration. From your experience, is there a problem to implement the optimization "as is" without caring about the rotational invariance? It seems to me that the gradient is transverse to the orbits of the configurations, so maybe it is ok to go directly... $\endgroup$ Apr 7, 2021 at 13:57
  • $\begingroup$ @Futurologist, the approach I suggested does not consider rotational symmetry, and depends on the initial condition. For example, for $n=3, N=60, p=1$ I found a local minima, but if you place the points in a snub dodecahedron you get a lower value for the objective function. Regarding the gradient, yes, it should be a force pointing outside, that's why I start my points inside the sphere and not over it. $\endgroup$
    – nicoguaro
    Apr 7, 2021 at 17:21

Here is a heuristic approach, suitable if you want an approximate solution, or a good starting point for the optimization: take $M \gg N$ points uniformly at random on the $n$-sphere ($M=10N$ should be enough), apply the $k$-means algorithm with $k=N$, and use the resulting $N$ cluster centers.


I am posting some python toy code that I wrote hastily (so nothing sophisticated or guaranteed), after nicoguaro put my mind at ease, which encouraged me to give this a very simplistic try. I used the classical Coulomb's potential to construct the objective function, followed a straight forward gradient descend by first following the gradient of the objective function for very short period of time and then normalizing the result to force it back on the sphere, and finally I wrote a vectorized implementation):

import numpy as np
import math

the configuration matrix q should be of size [ dim ] x [ n_points ], e.g.
dim = 2:  q.shape = 2, n_points
dim = 3:  q.shape = 3, n_points
Dq is the matrix of size [dim] x [n_points] x [n_points] such that
Dq[:, i, j] = q[:,i] - q[:,j] which is the relative position vector,
Dq is a symmetric matrix with respect to axis=0 and axis=1
    The line:
    Dq = (D**3)*Dq 
    represents element-wise multiplication in the last two dimensions 
    Dq[0,:,:] = (D**3) * Dq[0,:,:] 
    Dq[1,:,:] = (D**3) * Dq[1,:,:]
Objective function is defined as the sum of Coulomb's repelling potentials 
between the points in the configuration matrix q [dim] x [n_points]
full_Grad() calculates the gradient of the objective function at the point 
configuration q in the ambient space 
def full_Grad(q):
    Dq = - q[:, :, np.newaxis] + q[:, np.newaxis, :]    
    D = np.linalg.norm(Dq, axis = 0) + np.identity(Dq.shape[1])
    D = 1/D
    D = D - np.identity(D.shape[1])     # column times row times elemntwise 
    Dq = (D**3)*Dq 
    return Dq.sum(axis=-1)
    sum the matrix along the last dimension, i.e. a tensor contraction 
    along the third index

Grad() projects the gradient of the objective function at the configuration 
q in the ambient space to the tangent planes of the dim-1 sphere, i.e. it 
extracts the gradient component that is tangent to the unit sphere 
def Grad(q):
    Gr = full_Grad(q)
    return  Gr - np.sum(Gr*q, axis=0)*q

Grad_flow() is the gradient flow, implemented by a simplistic Euler's method. 
Can be improved if necessary.
def Grad_flow(q, step):
    Q = q
    n_iter = 0
    while True:
        n_iter = n_iter + 1
        Q_prev = Q
        G = Grad(Q)
        G_max = np.max(np.linalg.norm(G, axis=0))
        h = min(1/(7*G_max), step)
        Q = Q - h * G
        Q = Q / np.linalg.norm(Q, axis = 0)
        if np.sum(np.abs((Q.T).dot(Q) - (Q_prev.T).dot(Q_prev))) < 1e-15:            
            return Q, n_iter

Toy test in dimension two with regular tetrahedral configuration

# vertices of the regular polyhedron        
A = np.array([math.sqrt(1 - (1/3)**2), 0, -1/3])
B = np.array([A[0]*math.cos(2*math.pi/3), A[0]*math.sin(2*math.pi/3), -1/3])
C = np.array([B[0], -B[1], -1/3])
D = np.array([0,0,1])

# configuration matrix of the regular polyhedron
q0 = np.array([A, B, C, D])

# random deformation of the regular polyhedron
dq0 = np.random.random((4,3))
q = q0 + dq0/2

q = q / np.linalg.norm(q, axis=1)[:, np.newaxis]

# transposition of the configuration matrices, bringing them in the format
# [dim] x [n_points]
q0 = q0.T
q = q.T

# Optimization, gradient flow
q1, n_it = Grad_flow(q, step=0.007)

print('templet regular tetrahedron q0: ')

print('deformed initial tetrahedron q: ')
print('Gram matrix of the initial tetrahedron q: ')
print('Gram matrix of the regular tetrahedron q0: ')
print('tetrahedron after optimization: ')
print('check if its vertices are on the unit sphere')
print(np.linalg.norm(q1, axis=0))
print('Gram matrix of the regular tetrahedron q0: ')
print('Gram matrix of the optimizing tetrahedron q1: ')
print('checking if the Gram matrices are almost the same, meaning that q1 and q0 are congruent: ')
print((q1.T).dot(q1) - (q0.T).dot(q0))
print('number of iteration needed: ')

Observe I generate the Gram matrices (which are basically the cosine of the spherical distances between the points from the configuration) of the regular tetrahedron q0 and the result of the optimization q1 (gradient descend along the surface of the sphere) and then measure the difference. Since the result of the optimization q1 is more likely to be an almost regular tetrahedron, but rotated with respect to the templet one q0, one can check that they are congruent by measuring the spherical distances between every pair of vertices, which is reflected in the Gram matrix, which contains the cosines of the said spherical distances.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.