: Selective \(G\)-Bispectra Made Practical

The usual way to make a network rotation-invariant is to pool. Each filter is applied in every orientation, and only the largest response is kept. The result is invariant, but it throws away much more than the orientation. The power spectrum, the textbook invariant, is no better, since it keeps the magnitudes of the Fourier coefficients and drops their relative phases.

The problem extends beyond rotations. For translations, reflections, or more general symmetry groups, we'd like a representation that preserves as much information as possible while ignoring the group action.

The \(G\)-bispectrum keeps the relationships between Fourier coefficients that the group action preserves and the power spectrum loses. The construction works across different groups, although computing the full bispectrum is often prohibitively expensive for neural networks. We show how to retain just a small, carefully chosen subset of its coefficients, with completeness guarantees for several symmetry groups.

Consider the Earth's land mask, represented using spherical harmonics up to degree 20, with 441 coefficients. Rotate the globe and the coefficients generally change, but its power spectrum, just 21 numbers, remains the same. Unfortunately, plenty of other signals have exactly the same power spectrum. The globe at the top right is one of them.

Now consider the 638 features in the lower panel, sums of products of three or four coefficients, chosen so that a rotation leaves every one unchanged. Unlike the power spectrum, they contain enough information for a numerical optimizer to reconstruct the original land mask, up to an arbitrary rotation.

Those features form a selective \(G\)-bispectrum, which retains a small subset of the higher-order invariants rather than computing all of them. Building on completeness results for several groups, we provide strong numerical evidence on the sphere and introduce a PyTorch library that makes these representations practical to use in neural networks.

Orbits and completeness

To say what an invariant should forget, let's write \(R \cdot f\) for the signal \(f\) rotated by \(R\). We call the rotation the signal's pose and the set \(\{R \cdot f\}\) of all its rotated copies its orbit. An invariant should forget the pose and nothing else. Formally, an invariant \(\Phi\) is constant on orbits, and it is complete if it also tells orbits apart, that is if \(\Phi(f) = \Phi(g)\) only when \(g = R \cdot f\) for some \(R\). Following the paper, we call the vector \(\Phi(f)\) the signal's invariant descriptor. Max pooling over the 8 rotations by multiples of 45° is far from complete, since it keeps one response out of 8.

Counting dimensions gives a lower bound on the size of a complete invariant that varies smoothly with the signal. A generic signal has no rotational symmetry, so its orbit is three-dimensional, and the orbits of signals with 441 coefficients form a space of dimension 441 − 3 = 438. The power spectrum, with 21 numbers, cannot be complete. The 638 features meet the bound, which is necessary for completeness but doesn't prove it.

The construction is easiest to follow on the ring first, where a rotation multiplies each Fourier coefficient by a phase.

Fourier coefficients on the ring

Let's take a periodic function \(r(\theta)\), drawn as a closed curve with radius \(r(\theta)\) at angle \(\theta\), and write it as a mean radius \(c_0\) plus 8 waves,

\[ r(\theta) = c_0 + 2 \sum_{k=1}^{8} |c_k| \cos(k\theta - \varphi_k). \]

The complex numbers \(c_k = |c_k|\, e^{i\varphi_k}\) are the Fourier coefficients of \(r\), conjugated so that a turn advances the phases. Rotating the curve by \(\alpha\) replaces \(\theta\) by \(\theta - \alpha\), which multiplies each coefficient by a phase,

\[ c_k \mapsto e^{ik\alpha}\, c_k, \]

so the phases advance at different rates, \(\varphi_k\) by \(k\alpha\), and the magnitudes don't change. The group acts diagonally on the coefficients, and the pose is a single phase shared by all of them in the ratio 1 : 2 : 3 : …

The power spectrum \(|c_k|^2\) drops the phases, so it is invariant, and it is also very far from complete. In step 2 of the figure below, a scramble keeps every \(|c_k|\), draws new phases at random, and gives a different curve with the same power spectrum. A rotation moves the phases \(\varphi_1, \dots, \varphi_8\) along a line, \(\varphi_k \mapsto \varphi_k + k\alpha\), whereas the scramble can move them anywhere. What we want are functions of the phases that are constant along such lines and tell different lines apart.

Triple products

Phases add under multiplication, so a product of coefficients whose frequencies sum to zero is invariant. The smallest interesting ones are cubic:

\[ \beta_{k_1, k_2} = c_{k_1}\, c_{k_2}\, \overline{c_{k_1 + k_2}} \;\mapsto\; e^{i(k_1 + k_2 - (k_1 + k_2))\alpha}\, \beta_{k_1, k_2} = \beta_{k_1, k_2}. \]

The phase of \(\beta_{k_1, k_2}\) is \(\varphi_{k_1} + \varphi_{k_2} - \varphi_{k_1 + k_2}\), a relation between three phases, which a rotation preserves and the scramble doesn't. Three is the lowest order at which this happens, since a single coefficient is invariant only for \(k = 0\) and the only invariant products of two are \(c_k \overline{c_k} = |c_k|^2\). The table of all \(\beta_{k_1, k_2}\) is the bispectrum, so called because it has two frequency indices where the power spectrum has one. Its row \(k_1 = 0\) is \(c_0 |c_k|^2\), the power spectrum again.

The figure below draws the entries with \(k_1 \le k_2\) and \(k_1 + k_2 \le 8\) by their phases, each under the column of \(k_2\). Its second panel is the plane of \((\varphi_1, \varphi_2)\), in which a rotation moves the signal's point along a line of slope 2, a level set of \(\arg \beta_{1,1} = 2\varphi_1 - \varphi_2\).

The library's logo is this triangle of three coefficients whose phases cancel under rotation, set in the round counter of a b.

Fourier transforms on groups

Nothing above is specific to the ring. What made it work is that the Fourier basis diagonalises the group action: each coefficient transforms on its own, by a phase. For a general compact group the action can't always be diagonalised, but it can be block-diagonalised. The blocks are the irreducible representations of \(G\), or irreps: unitary matrix representations of \(G\) that cannot be split into smaller ones. The Fourier transform on \(G\) has one coefficient per irrep, and irreps play the role of frequencies. For a finite group, the coefficient at an irrep \(\rho\) of dimension \(d\) is a \(d \times d\) matrix,

\[ \hat f(\rho) = \sum_{g \in G} f(g)\, \rho(g)^\dagger. \]

Let's shift the signal by a group element \(h\), \((h \cdot f)(g) = f(h^{-1} g)\). Each coefficient is then multiplied on the right by a unitary matrix,

\[ \hat f(\rho) \mapsto \hat f(\rho)\, \rho(h)^\dagger, \]

and the power spectrum becomes the matrix \(\hat f(\rho)\, \hat f(\rho)^\dagger\), invariant because \(\rho(h)^\dagger \rho(h) = I\). For SO(2) the irreps are the \(1 \times 1\) matrices \(e^{ik\alpha}\), the transform is the usual one, and a rotation by \(\alpha\) multiplies \(\hat f(\rho_k)\) by \(e^{-ik\alpha}\). The ring's \(c_k\), which turn by \(e^{ik\alpha}\), are their complex conjugates. Conjugating every formula changes nothing else, and from here on we keep the paper's conventions.

Counting shows that the transform loses nothing. The irreps of a finite group satisfy \(\sum_\rho d_\rho^2 = |G|\), so the coefficients of the transform hold as many numbers as the signal has values, and the transform is a change of basis. The table below shows the transform for a few finite and infinite groups.

Commutative groups are the rows with only \(1 \times 1\) irreps. The sphere is the odd row, because it is a space that SO(3) acts on rather than the group itself. Each irrep therefore contributes a vector of length \(2\ell + 1\) instead of a square matrix, and the count is \(\sum_\ell (2\ell + 1)\), not \(\sum_\ell (2\ell + 1)^2\).

Non-commutative groups

On the ring every irrep is \(1 \times 1\) because the group is commutative. Matrices that commute with each other can be diagonalised together, so a representation of a commutative group splits into lines. A non-commutative group has irreps larger than \(1 \times 1\). The smallest examples are the dihedral groups \(D_n\), the \(n\) turns and \(n\) flips of a regular \(n\)-gon, generated by a turn \(a\) and a flip \(x\) with \(x a x = a^{-1}\), so turning and then flipping is not the same as flipping and then turning. Two irreps of \(D_n\) are \(1 \times 1\) (four when \(n\) is even) and send each element to \(\pm 1\). The rest are \(2 \times 2\),

\[ \rho_k(a) = \begin{pmatrix} \cos \frac{2\pi k}{n} & -\sin \frac{2\pi k}{n} \\ \sin \frac{2\pi k}{n} & \cos \frac{2\pi k}{n} \end{pmatrix}, \qquad \rho_k(x) = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}, \qquad 0 < k < n/2. \]

The two matrices don't commute, so no change of basis diagonalises both. The coefficients \(\hat f(\rho_k)\) are \(2 \times 2\) matrices, and a shift mixes their columns.

The triple product needs a change of basis. On the ring, \(c_{k_1} c_{k_2}\) turns like the single frequency \(k_1 + k_2\), and one conjugate coefficient cancels it. Here the product of two coefficients is their tensor product \(\hat f(\rho_1) \otimes \hat f(\rho_2)\), a \(4 \times 4\) matrix that transforms by \(\rho_1 \otimes \rho_2\). That representation is reducible: \(\rho_i \otimes \rho_j\) splits into \(\rho_{i+j} \oplus \rho_{i-j}\). Indices are read modulo \(n\) and up to sign, and a block whose index is a multiple of \(n/2\) splits further into two \(1 \times 1\) irreps. For \(D_5\):

On the ring the same table has a single frequency in each cell, \(k_1 + k_2\). The split is a unitary change of basis \(C_{\rho_1, \rho_2}\), the matrix of Clebsch–Gordan coefficients:

\[ C_{\rho_1, \rho_2}^\dagger\, (\rho_1 \otimes \rho_2)(g)\, C_{\rho_1, \rho_2} = \bigoplus_{\rho \in \rho_1 \otimes \rho_2} \rho(g). \]

The \(G\)-bispectrum multiplies the product, in that basis, by the conjugate coefficients at the irreps it contains:

\[ \beta_{\rho_1, \rho_2} = \big[\hat f(\rho_1) \otimes \hat f(\rho_2)\big]\, C_{\rho_1, \rho_2} \Big[\bigoplus_{\rho \in \rho_1 \otimes \rho_2} \hat f(\rho)^\dagger\Big]\, C_{\rho_1, \rho_2}^\dagger, \]

where the middle factor is block-diagonal, with one block \(\hat f(\rho)^\dagger\) for each irrep in the split. Under a shift by \(h\) the first factor picks up \((\rho_1 \otimes \rho_2)(h)^\dagger = C \big[\bigoplus \rho(h)^\dagger\big] C^\dagger\) on the right, the block-diagonal factor picks up \(\bigoplus \rho(h)\) on the left, and the two cancel. On the ring every block is \(1 \times 1\), \(C = 1\), the split is the single frequency \(k_1 + k_2\), and the formula is \(\beta_{k_1, k_2}\) up to conjugation.

Kakarala (1992) proved that the \(G\)-bispectrum determines a signal on a finite group up to a shift when every Fourier coefficient is nonsingular. We call such signals generic. On the ring the condition is that every \(c_k\) is nonzero. In both cases the exceptions form a set of measure zero. On the sphere no such condition is known, so there generic just means outside a set of measure zero.

The recipe is the same for any compact group. Take the Fourier coefficients at two irreps, decompose their tensor product with Clebsch–Gordan coefficients, and multiply the result by the conjugate coefficients at the irreps it contains. Each module of the library is this construction for a different group and domain, with its own irreps and Clebsch–Gordan tables.

Sphere

On the sphere, under SO(3), the irreps are indexed by a degree \(\ell = 0, 1, 2, \dots\), and the corresponding basis functions are the spherical harmonics \(Y_\ell^m\), \(-\ell \le m \le \ell\). The coefficients of degree \(\ell\) form a vector \(\mathbf F_\ell = (a_\ell^{-\ell}, \dots, a_\ell^{\ell}) \in \mathbb C^{2\ell + 1}\), and a rotation \(R\) acts on it by a unitary matrix of size \(2\ell + 1\), the Wigner D-matrix \(D^\ell(R)\):

\[ \mathbf F_\ell \mapsto D^\ell(R)\, \mathbf F_\ell. \]

The figure below shows a rotation on all coefficients at once: a block-diagonal matrix, in which degrees never mix but the coefficients within a degree do. Rotations about the pole are the exception, since they multiply each \(a_\ell^m\) by a phase, \(m\) times the angle, exactly as on the ring. In practice a signal is kept up to some degree \(L\), its band-limit, which leaves \((L + 1)^2\) coefficients. The power spectrum is \(\|\mathbf F_\ell\|^2\), one number per degree. It is invariant because \(D^\ell\) is unitary, and as incomplete as on the ring.

As on \(D_n\), the product of two coefficient vectors doesn't transform like a single degree. The tensor product \(\mathbf F_{\ell_1} \otimes \mathbf F_{\ell_2}\) transforms by \(D^{\ell_1} \otimes D^{\ell_2}\), which splits into one irrep of every degree from \(|\ell_1 - \ell_2|\) to \(\ell_1 + \ell_2\). Clebsch–Gordan coefficients again perform the split, and for SO(3) they are the numbers \(C^{\ell m}_{\ell_1 m_1 \ell_2 m_2}\). Projecting the product onto degree \(\ell\) and taking its inner product with \(\mathbf F_\ell\) gives the spherical bispectrum,

\[ \beta_{\ell_1, \ell_2, \ell} = \sum_{m_1 + m_2 = m} C^{\ell m}_{\ell_1 m_1 \ell_2 m_2} \; a_{\ell_1}^{m_1}\, a_{\ell_2}^{m_2}\, \overline{a_\ell^{m}}, \]

one number for each triple of degrees with \(|\ell_1 - \ell_2| \le \ell \le \ell_1 + \ell_2\). The projected product and \(\mathbf F_\ell\) transform by the same unitary \(D^\ell(R)\), so their inner product is invariant. On the ring a pair of frequencies couples to the single frequency \(k_1 + k_2\), on \(D_n\) a pair of \(2 \times 2\) irreps couples to two, three or four, and on the sphere a pair of degrees couples to a whole range of them.

Cost

The bispectrum has been known since 1953 and applies to any compact group, but it has rarely been used in networks, mostly because of its size. Take a 256 × 256 image under cyclic shifts. The FFT gives 65,536 coefficients, and the bispectrum has an entry for every pair of frequencies: 2,147,516,416 numbers, about 17 GB at 8 bytes each (our arithmetic), per image and per forward pass. On any finite group the full bispectrum has \(O(|G|^2)\) scalars, quadratic in the size of the signal.

Ring and finite groups

Most of those entries are redundant. On the ring, \(c_0\) is real for a real signal, so \(\beta_{0,0} = c_0^3\) gives \(c_0\) and \(\beta_{0,1} = c_0 |c_1|^2\) gives \(|c_1|\). Each entry \(\beta_{1,k} = c_1 c_k \overline{c_{k+1}}\) then gives the next coefficient,

\[ c_{k+1} = \overline{\beta_{1,k} / (c_1 c_k)}. \]

Step 4 of the ring figure runs this chain. No entry fixes the phase of \(c_1\), and step 5 shows that any value of it reproduces the signal, rotated. For a generic signal none of the divisions is by zero, so on a ring of \(n\) points this chain, \(n\) of the \(n(n+1)/2\) distinct entries, determines the signal up to its pose, and every other entry can be recomputed from it.

A selective bispectrum is a subset of this kind: enough entries to determine a generic signal up to the group action, and a number of them that grows linearly with the size of the signal rather than quadratically. For the image that is 65,536 numbers, one per pixel. Mataigne et al. (2024) select \(O(|G|)\) entries on any finite group by a search over its table of tensor products, like the one for \(D_5\) above, and prove that the selection is complete for commutative groups, the dihedral groups and the rotations of a cube, with or without reflections. For the 24 rotations of a cube it keeps 172 of the 576 scalars.

Sphere, degree by degree

On the sphere the full bispectrum grows more slowly than on a finite group, as \(L^3\) for \((L + 1)^2\) coefficients, but each entry is a sum of Clebsch–Gordan terms. At \(L = 15\), the setting of the paper's experiments, a signal has 256 coefficients and the full bispectrum 1,124 entries. At \(L = 64\) its 70,257 entries add up 267 million such terms. Our selective version keeps 380 numbers at \(L = 15\) and 5,610 at \(L = 64\), for 25 times less work.

On a real signal the full list repeats itself, since entries that differ only in the order of their three degrees are proportional and some are always zero. Counted once each, there are 352 distinct entries at \(L = 15\) and 22,945 at \(L = 64\). Against those the selective version saves little at small band-limits, and at \(L = 64\) it does 8 times less work.

On the sphere there is no recursion as simple as the ring's, since each degree brings \(2\ell + 1\) unknowns at once. The paper's construction goes degree by degree. Once the lower degrees are known, an entry \(\beta_{a, b, \ell}\) with \(a, b < \ell\) is linear in \(\mathbf F_\ell\), so the entries kept at degree \(\ell\) form a system in its \(2\ell + 1\) real unknowns. Each degree keeps up to \(2\ell + 1\) such triple products, a few self-couplings \(\beta_{\ell, \ell, \ell'}\), and a few product energies \(\|(\mathbf F_{\ell_1} \otimes \mathbf F_{\ell_2})|_\ell\|^2\), the squared norms of Clebsch–Gordan projections. They are chosen so that the system has full rank, and the total is \(\Theta(L^2)\). The paper calls the triple products bispectral coefficients and the product energies Clebsch–Gordan power coefficients. Since the energies are not entries of the bispectrum, it calls the whole descriptor the augmented selective spherical bispectrum.

Testing completeness

The figure below tests completeness numerically. A solver gets only the invariant descriptor and tries to recover the signal. It runs at \(L\) = 12, as do the classifier figure and the playground further down, so that a phone can solve it in a moment. At that band-limit the descriptor has 254 numbers (196 triple products and 58 product energies) for a signal of 169 coefficients. The paper's, at \(L = 15\), has 380. (The globe at the top runs at \(L\) = 20, to draw the coastlines more sharply.)

Modules

Here is the ring's chain again, as a module:

import torch
from bispectrum import CnonCn, OctaonOcta, SO2onDisk

f = torch.randn(16, 8)        # 16 signals on a ring of 8 points
bsp = CnonCn(n=8)             # selective by default
beta = bsp(f)                 # (16, 8): the chain, n entries
f_rec = bsp.invert(beta)      # complex: f up to the chain's free phase

# the same interface for other groups
beta_o = OctaonOcta()(torch.randn(16, 24))         # the 24 rotations of a cube
beta_d = SO2onDisk(L=32)(torch.randn(16, 32, 32))  # images on a disk

bispectrum is a PyTorch library of seven such modules, one for each group and the domain it acts on, named {Group}on{Domain}:

• CnonCna ring of \(n\) points, turned in \(n\) steps
• SO2onS1a circle, sampled at \(n\) points
• TorusOnTorusimages and volumes under cyclic shifts
• DnonDnthe \(2n\) turns and flips of a polygon
• SO2onDiskan image on a disk, rotated in the plane
• SO3onS2a signal on the sphere, rotated in space
• OctaonOctathe 24 rotations of a cube

Each is a torch.nn.Module without parameters. Calling it returns the selective bispectrum, selective=False returns the full one where it is implemented, and invert recovers the signal where an inversion is known. The modules are differentiable, so they can sit anywhere in a network, and two placements are natural.

Example usage

In front of an MLP. When the input is itself a signal on the sphere, the selective bispectrum can be the whole feature extractor, and the classifier on top sees the pose only through the grid's sampling error. This is the model of the paper's Spherical MNIST experiment:

import torch
from torch import nn
from bispectrum import SO3onS2

bsp = SO3onS2(lmax=15, nlat=64, nlon=128)  # 380 numbers, real or imaginary
head = nn.Sequential(
    nn.Linear(760, 256), nn.BatchNorm1d(256), nn.ReLU(),
    nn.Linear(256, 128), nn.BatchNorm1d(128), nn.ReLU(),
    nn.Linear(128, 10),
)

beta = bsp(x)                              # x: (batch, 64, 128)
z = torch.cat([beta.real, beta.imag], dim=-1)
z = z.sign() * z.abs().log1p()             # signed log: they span many scales
logits = head(z)

In place of pooling. A group convolution applies each filter in every pose of a group, say the 8 rotations by 45°, so its feature maps carry an extra axis with one slot per group element. Rotating the input by 45° rotates the feature maps and shifts them cyclically along that axis. Equivariant networks usually reach an invariant by taking a max over the axis, which keeps one response out of 8, the max pooling we started from. The axis is a signal on a ring of 8 points, so CnonCn can take the place of the max:

import torch
import torch.nn.functional as F
from bispectrum import CnonCn

# x: output of a group convolution, (batch, channels, 8 poses, height, width)
before = x.max(dim=2).values               # keep the largest of the 8

bsp = CnonCn(n=8)                          # keep everything but the pose
B, C, G, H, W = x.shape
# one ring of 8 per channel and pixel
beta = bsp(x.permute(0, 1, 3, 4, 2).reshape(-1, G))
z = torch.cat([beta.real, beta.imag], dim=-1)
z = z.sign() * z.abs().log1p()
z = z.reshape(B, C, H, W, -1).permute(0, 1, 4, 2, 3).reshape(B, -1, H, W)
after = F.relu(proj(z))                    # proj = nn.Conv2d(16 * C, C, 1)

Speed

Two implementation details make the modules fast enough to use as layers. The first is sparsity. Most Clebsch–Gordan coefficients are zero, so when a module is built it lists every nonzero term of every kept entry once: the indices of its three coefficients, its weight and its output entry. The forward pass is then a gather, a product and a scatter-add, with no loop over degrees. From so3_on_s2.py, with names shortened:

prods = flat[:, i1] * flat[:, i2] * cg * torch.conj(flat[:, i3])
_scatter_add_complex(result, entry, prods)

At \(L = 64\) the list has 10.8 million terms, against 267 million for the full bispectrum.

The second is launch overhead. At batch 16 a forward pass is a few dozen small kernels, and the time goes into launching them from Python rather than running them. With autograd off, each module records its forward pass as a CUDA graph on the first input it sees and replays it for every later input of the same shape, with nothing to configure. An input of another shape runs without the graph.

Every module in the table takes under a millisecond with autograd on, from 0.07 to 0.67 ms, and 0.03 to 0.20 ms at inference.

Spherical MNIST

The figure below runs three small classifiers on the same spherical digit, all trained on upright digits only: an MLP on the selective bispectrum, a CNN on pixels, and an MLP on the power spectrum.

The paper runs this experiment at \(L = 15\), over five seeds and against stronger baselines, so its numbers differ from the figure's in the last digit. The bispectrum with its MLP reaches 95.8% on rotated digits, not significantly different from a spherical CNN of matched size (96.1%, \(p = 0.14\)) and a quarter of a point below a tensor-product network (\(p = 0.006\)), while running 4 times faster than the first and 364 times faster than the second. The pixel CNN falls from 99.1% to 13.4%. The power spectrum reaches 77.7%, and 80.1% with an MLP 48 times larger.

As a pooling layer

The paper puts the bispectrum into two equivariant networks in place of the pooling layer and compares it with up to five alternatives, 46 comparisons in all, and it is never significantly worse. On tissue images under 8 rotations (PatchCamelyon) everything ties, in tests that three seeds per model make sensitive only to gaps of about 0.025 AUC. On 3D organ scans under the cube's 24 rotations (OrganMNIST3D) it beats max pooling by 5 points (78.5% against 73.6%). Its smallest model, with 0.46 million parameters, beats max pooling's largest, with 6.0 million, by 2.6 points. The gain over max pooling disappears with 250 training examples or fewer and is not significant at the largest width. Our guess is that the more a max discards, the more there is to recover. A max over the cube's 24 rotations discards more than one over 8.

Try it

Draw something on the pad: a letter, a face, a continent of your own. It starts with a 0, so clear it first if it's in the way. When you lift the pen, the page computes the 254 numbers of your drawing's selective bispectrum and gives only those to the solver. The reconstruction appears on the right.

Very symmetric drawings, such as a single dot or a centred ring, are the hard case. They are unchanged by every rotation about some axis, so they are far from generic. If the solver gives up, press Reconstruct for other random starts.

References

The results quoted here come from the paper, bispectrum: Selective G-Bispectra Made Practical, which has the proofs, the protocols and the tables. The asides, the small classifiers behind the figures and the rerun of the benchmark are additions of this page. The library is on GitHub and installs with pip install bispectrum on Python 3.12. The numbers quoted here are those of the camera-ready version, and an earlier version on arXiv may differ in places. If you use the library, please cite:

@article{mathe2026bispectrum,
  title   = {bispectrum: Selective $G$-Bispectra Made Practical},
  author  = {Mathe, Johan and {Myers Lantow}, Adele and Mataigne, Simon
             and Miolane, Nina},
  journal = {arXiv preprint arXiv:2605.07270},
  year    = {2026}
}

1. J. W. Tukey (1953). The spectral representation and transformation properties of the higher moments of stationary time series. Reprinted in The Collected Works of John W. Tukey, vol. 1, 165–184, Wadsworth, 1984.
2. R. C. Jennison (1958). A phase sensitive interferometer technique for the measurement of the Fourier transforms of spatial brightness distributions of small angular extent. Monthly Notices of the Royal Astronomical Society 118(3), 276–284.
3. Z. Kam (1980). The reconstruction of structure from electron micrographs of randomly oriented particles. Journal of Theoretical Biology 82(1), 15–39.
4. A. W. Lohmann, G. Weigelt and B. Wirnitzer (1983). Speckle masking in astronomy: triple correlation theory and applications. Applied Optics 22(24), 4028–4037.
5. G. B. Giannakis (1989). Signal reconstruction from multiple correlations: frequency- and time-domain approaches. Journal of the Optical Society of America A 6(5), 682–697.
6. R. Kakarala (1992). Triple Correlation on Groups. PhD thesis, University of California, Irvine.
7. R. Kondor (2007). A novel set of rotationally and translationally invariant features for images based on the non-commutative bispectrum. arXiv:cs/0701127.
8. R. Kakarala (2012). The bispectrum as a source of phase-sensitive invariants for Fourier descriptors: a group-theoretic approach. Journal of Mathematical Imaging and Vision 44(3), 341–353.
9. T. S. Cohen and M. Welling (2016). Group equivariant convolutional networks. ICML.
10. T. S. Cohen, M. Geiger, J. Köhler and M. Welling (2018). Spherical CNNs. ICLR.
11. B. S. Veeling, J. Linmans, J. Winkens, T. Cohen and M. Welling (2018). Rotation equivariant CNNs for digital pathology. MICCAI. The PatchCamelyon benchmark.
12. J. Yang, R. Shi, D. Wei, Z. Liu, L. Zhao, B. Ke, H. Pfister and B. Ni (2023). MedMNIST v2: a large-scale lightweight benchmark for 2D and 3D biomedical image classification. Scientific Data 10, 41. The source of OrganMNIST3D.
13. S. Mataigne, J. Mathe, S. Sanborn, C. Hillar and N. Miolane (2024). The selective G-bispectrum and its inversion: applications to G-invariant networks. NeurIPS 37, 115682–115711.
14. A. Myers Lantow and N. Miolane (2026). Selective disk bispectrum: a complete and rotation invariant image descriptor. arXiv:2511.19706.