# Algorithms to generate spherical codes

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.

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*}

$$$$\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 .$$$$

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 ,$$

$$\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}

• Thanks you, that is helpful! Apr 6 at 8:23
• 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. Apr 6 at 9:37
• @FedericoPoloni, I edited my answer. Please let me know if that's what you meant. Apr 6 at 19:21
• @nicoguaro Yes, that is what I meant, thanks! Apr 6 at 19:28
• @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... Apr 7 at 13:57

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]
configuration q in the ambient space
'''
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
'''

'''
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
'''
return  Gr - np.sum(Gr*q, axis=0)*q

'''
Can be improved if necessary.
'''
Q = q
n_iter = 0
while True:
n_iter = n_iter + 1
Q_prev = 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

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.

• Thanks, this looks promising! Apr 9 at 6:00