Algorithms to generate spherical codes

Question

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.

Answer 1

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}

Answer 2

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.

Answer 3

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(q0.T)
print('')

print('deformed initial tetrahedron q: ')
print(q.T)
print('')
print('Gram matrix of the initial tetrahedron q: ')
print((q.T).dot(q))
print('')
print('Gram matrix of the regular tetrahedron q0: ')
print((q0.T).dot(q0))
print('')
print('')
print('tetrahedron after optimization: ')
print(q1.T)
print('')
print('check if its vertices are on the unit sphere')
print(np.linalg.norm(q1, axis=0))
print('')
print('Gram matrix of the regular tetrahedron q0: ')
print((q0.T).dot(q0))
print('')
print('Gram matrix of the optimizing tetrahedron q1: ')
print((q1.T).dot(q1))
print('')
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('')
print('number of iteration needed: ')
print(n_it)

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.