Spherical Harmonics — the math foundation, visualized¶
Companion experiments for the research note 06-spherical-harmonics.md (tier 12, math foundations).
This notebook makes the mathematical facts tangible — the application to 3DGS color lives in tier 3B.
What we will see:
- The real SH basis $Y_{lm}$ as 3D lobes, band by band.
- Orthonormality: the Gram matrix $\int_{S^2} Y_i Y_j\,d\Omega$ is the identity.
- Truncating a sharp function on the sphere produces Gibbs-like ringing.
- A numerical check of the addition theorem $\sum_m Y_{lm}(x)Y_{lm}(y)=\tfrac{2l+1}{4\pi}P_l(x\cdot y)$.
- Per-band energy is rotation-invariant — the spectral fingerprint that makes SH useful for $SO(3)$.
Experiments 2 and 4 also serve as correctness checks of our hand-written SH implementation.
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from scipy.special import lpmv, eval_legendre, factorial
from mpl_toolkits.mplot3d import Axes3D # noqa: F401 (registers 3d projection)
import warnings
warnings.filterwarnings('ignore')
np.random.seed(42)
plt.rcParams.update({
'figure.dpi': 110,
'font.size': 11,
'axes.titlesize': 12,
'axes.labelsize': 11,
'figure.facecolor': 'white',
'axes.facecolor': '#f8f9fa',
'axes.grid': True,
'grid.alpha': 0.35,
})
print('Setup complete.')
Setup complete.
Helpers: a real spherical harmonic and a sphere-quadrature rule¶
We implement the real orthonormal SH directly from the associated Legendre function lpmv so the
notebook does not depend on a particular SciPy version of sph_harm. Convention: theta is the polar
angle in $[0,\pi]$ and phi is the azimuth in $[0,2\pi]$. The surface integral uses a uniform
$(\theta,\phi)$ grid with the $\sin\theta$ area weight.
def real_sph_harm(l, m, theta, phi):
"""Real, orthonormal spherical harmonic Y_lm(theta, phi).
theta: polar angle in [0, pi]; phi: azimuth in [0, 2*pi]."""
am = abs(m)
norm = np.sqrt((2 * l + 1) / (4 * np.pi) * factorial(l - am) / factorial(l + am))
P = lpmv(am, l, np.cos(theta)) # associated Legendre (with Condon-Shortley phase)
if m > 0:
return np.sqrt(2.0) * norm * P * np.cos(am * phi)
if m < 0:
return np.sqrt(2.0) * norm * P * np.sin(am * phi)
return norm * P
# Quadrature grid on S^2 (theta polar, phi azimuth) with sin-theta area weights.
N_THETA, N_PHI = 160, 320
_theta = np.linspace(0, np.pi, N_THETA)
_phi = np.linspace(0, 2 * np.pi, N_PHI, endpoint=False)
TH, PH = np.meshgrid(_theta, _phi, indexing='ij')
_dth = _theta[1] - _theta[0]
_dph = _phi[1] - _phi[0]
W = np.sin(TH) * _dth * _dph # area element dОmega
def sphere_integral(values):
"""Approximate the surface integral of `values` over S^2."""
return np.sum(values * W)
def dir_to_angles(d):
"""Unit vector(s) -> (theta, phi)."""
d = np.asarray(d, dtype=float)
theta = np.arccos(np.clip(d[..., 2], -1.0, 1.0))
phi = np.mod(np.arctan2(d[..., 1], d[..., 0]), 2 * np.pi)
return theta, phi
# Cartesian components of every grid direction (used later).
DX = np.sin(TH) * np.cos(PH)
DY = np.sin(TH) * np.sin(PH)
DZ = np.cos(TH)
print('Helpers ready. Grid shape:', TH.shape, '| total weight (should be 4pi):', round(sphere_integral(np.ones_like(TH)), 4))
Helpers ready. Grid shape: (160, 320) | total weight (should be 4pi): 12.566
Experiment 1: The SH basis as 3D lobes (bands $l=0,1,2$)¶
Each $Y_{lm}$ is a function on the sphere. The classic way to see it is to plot the surface $r(\mathbf d)=|Y_{lm}(\mathbf d)|$ and color it by the sign of $Y_{lm}$. Higher bands $l$ oscillate faster across directions — exactly the "angular frequency" intuition the eigenvalue $-l(l+1)$ promises.
bands = [(l, m) for l in range(3) for m in range(-l, l + 1)]
gt = np.linspace(0, np.pi, 80)
gp = np.linspace(0, 2 * np.pi, 160)
GT, GP = np.meshgrid(gt, gp, indexing='ij')
fig = plt.figure(figsize=(12, 8))
for (l, m) in bands:
Y = real_sph_harm(l, m, GT, GP)
r = np.abs(Y)
x = r * np.sin(GT) * np.cos(GP)
y = r * np.sin(GT) * np.sin(GP)
z = r * np.cos(GT)
idx = l * 5 + (3 + m) # center each band's row of 5 columns on column 3
ax = fig.add_subplot(3, 5, idx, projection='3d')
A = np.abs(Y).max() + 1e-12
facecolors = cm.coolwarm(0.5 + 0.5 * Y / A)
ax.plot_surface(x, y, z, facecolors=facecolors, rstride=2, cstride=2,
linewidth=0, antialiased=False, shade=False)
lim = r.max() + 1e-3
ax.set_xlim(-lim, lim); ax.set_ylim(-lim, lim); ax.set_zlim(-lim, lim)
ax.set_title(f'l={l}, m={m}', fontsize=10)
ax.set_axis_off()
fig.suptitle('Real spherical harmonics: radius = |Y_lm|, color = sign(Y_lm)', fontsize=13)
plt.tight_layout()
plt.show()
What you see¶
Band $l=0$ is a single round ball: the constant function, the "DC term" with no directional variation. Band $l=1$ gives three dumbbells aligned with the $x$, $y$, $z$ axes — the lowest-order directional variation (these are literally proportional to the Cartesian coordinates). Band $l=2$ adds five more intricate lobes with extra sign changes. The pattern: a degree-$l$ harmonic has $2l+1$ members, and the number of sign flips grows with $l$. That growing oscillation rate is the geometric meaning of the Laplace–Beltrami eigenvalue $-l(l+1)$: higher bands carry higher angular frequency, just as higher Fourier modes carry higher temporal frequency.
Experiment 2: Orthonormality — the Gram matrix is the identity¶
The single property that makes SH a basis (and lets us read off coefficients by projection) is orthonormality: $\int_{S^2} Y_i\,Y_j\,d\Omega = \delta_{ij}$. We compute that integral numerically for all 16 harmonics up to $l=3$ and display the Gram matrix. If our implementation is correct, it must be the identity.
basis = [(l, m) for l in range(4) for m in range(-l, l + 1)] # 16 harmonics, l = 0..3
Yvals = [real_sph_harm(l, m, TH, PH) for (l, m) in basis]
n = len(basis)
G = np.zeros((n, n))
for i in range(n):
for j in range(i, n):
G[i, j] = G[j, i] = sphere_integral(Yvals[i] * Yvals[j])
labels = [f'{l},{m}' for (l, m) in basis]
fig, ax = plt.subplots(figsize=(7.2, 6))
im = ax.imshow(G, cmap='RdBu_r', vmin=-1, vmax=1)
ax.set_xticks(range(n)); ax.set_xticklabels(labels, rotation=90, fontsize=7)
ax.set_yticks(range(n)); ax.set_yticklabels(labels, fontsize=7)
ax.set_title('Gram matrix $\\int_{S^2} Y_i Y_j\\, d\\Omega$ (l = 0..3)')
ax.set_xlabel('(l, m)'); ax.set_ylabel('(l, m)')
ax.grid(False)
fig.colorbar(im, ax=ax, label='inner product')
plt.tight_layout()
plt.show()
print('max |G - I| (off from identity):', float(np.max(np.abs(G - np.eye(n)))))
max |G - I| (off from identity): 0.0002277859206373245
What you see¶
A clean red diagonal of ones on a white (zero) background — the Gram matrix is the identity to within a tiny quadrature error (the printed deviation is on the order of $10^{-3}$ or smaller, shrinking if we refine the grid). This is the numerical proof that our $Y_{lm}$ are mutually orthogonal and unit-norm. The practical payoff: because the basis is orthonormal, the expansion coefficient of any function is just its inner product with the basis function, $a_{lm}=\int f\,Y_{lm}\,d\Omega$ — no linear system to solve, exactly like reading off a Fourier coefficient.
Experiment 3: Truncation and ringing (Gibbs on the sphere)¶
We take a spherical cap centered at the north pole — $f(\mathbf d)=1$ if $\theta<\theta_0$, else $0$ — a function with a sharp circular edge. Because it is azimuthally symmetric, only the $m=0$ harmonics contribute. We expand it in SH and reconstruct with increasing maximum degree $L$, then look at a cross-section along a meridian. A global basis cannot represent a jump cleanly, so we expect overshoot and oscillation near the edge.
theta0 = np.deg2rad(50.0)
f_grid = (TH < theta0).astype(float)
L_MAX = 25
coeffs0 = np.array([sphere_integral(f_grid * real_sph_harm(l, 0, TH, PH)) for l in range(L_MAX + 1)])
th_line = np.linspace(0, np.pi, 400)
ph_line = np.zeros_like(th_line)
def reconstruct(L):
s = np.zeros_like(th_line)
for l in range(L + 1):
s += coeffs0[l] * real_sph_harm(l, 0, th_line, ph_line)
return s
fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(np.rad2deg(th_line), (th_line < theta0).astype(float), 'k-', lw=2.5, alpha=0.4, label='target step')
ax.axvline(np.rad2deg(theta0), color='k', ls='--', lw=1, label='true edge ($\\theta_0=50°$)')
for L in [2, 5, 10, 25]:
ax.plot(np.rad2deg(th_line), reconstruct(L), lw=1.6, label=f'L = {L}')
ax.set_xlabel('polar angle $\\theta$ (degrees)')
ax.set_ylabel('reconstructed $f(\\theta)$')
ax.set_title('SH truncation of a spherical cap: ringing near the sharp edge')
ax.legend(loc='upper right', fontsize=9)
plt.tight_layout()
plt.show()
What you see¶
Low $L$ (e.g. 2) gives a smooth, badly blurred version of the step — the cap edge is smeared over tens of degrees. As $L$ grows the reconstruction tightens around the true edge, but overshoot and wiggles appear right next to the jump and do not disappear; they only get narrower. This is the spherical analogue of the Gibbs phenomenon in Fourier series, and it is a direct consequence of every $Y_{lm}$ having global support over the whole sphere. The lesson for applications: SH are excellent for smooth angular signals and poor for sharp ones (hard shadows, tight specular highlights) — pushing $L$ up fights the symptom at quadratic cost, never curing it.
Experiment 4: The addition theorem, checked numerically¶
The addition theorem collapses a whole band into a coordinate-free quantity:
$$\sum_{m=-l}^{l} Y_{lm}(\mathbf x)\,Y_{lm}(\mathbf y) = \frac{2l+1}{4\pi}\,P_l(\mathbf x\cdot\mathbf y).$$
We fix $\mathbf y$, sweep $\mathbf x$ along a great circle so the angle $\gamma$ between them runs from $0$ to $\pi$, and overlay the left-hand sum against the right-hand Legendre expression.
l = 3
y = np.array([0.3, -0.5, 0.8]); y = y / np.linalg.norm(y)
ty, py = dir_to_angles(y)
# Build x sweeping the great circle through y: x(s) = cos(s) y + sin(s) perp.
a = np.array([1.0, 0.0, 0.0])
perp = a - np.dot(a, y) * y
perp = perp / np.linalg.norm(perp)
s = np.linspace(0, np.pi, 200)
X = np.cos(s)[:, None] * y[None, :] + np.sin(s)[:, None] * perp[None, :]
tx, px = dir_to_angles(X)
lhs = np.zeros_like(s)
for m in range(-l, l + 1):
lhs += real_sph_harm(l, m, tx, px) * real_sph_harm(l, m, ty, py)
cos_g = X @ y
rhs = (2 * l + 1) / (4 * np.pi) * eval_legendre(l, cos_g)
fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(np.rad2deg(s), lhs, lw=3.2, label=r'$\sum_m Y_{lm}(x)\,Y_{lm}(y)$')
ax.plot(np.rad2deg(s), rhs, '--', lw=2, color='crimson', label=r'$\frac{2l+1}{4\pi}\,P_l(x\cdot y)$')
ax.set_xlabel(r'angle $\gamma$ between $x$ and $y$ (degrees)')
ax.set_ylabel('value')
ax.set_title(f'Addition theorem numerical check (l = {l})')
ax.legend(fontsize=11)
plt.tight_layout()
plt.show()
print('max |LHS - RHS|:', float(np.max(np.abs(lhs - rhs))))
max |LHS - RHS|: 7.216449660063518e-16
What you see¶
The two curves lie exactly on top of each other (printed difference near machine precision). The left-hand side was built from the absolute orientations of $\mathbf x$ and $\mathbf y$ in our chosen coordinate frame, yet it depends only on the relative angle $\gamma$ — that is the whole content of the theorem. It is the reproducing kernel of the degree-$l$ eigenspace, and it is why a band behaves as one indivisible, rotation-closed unit: rotating the frame reshuffles the individual $Y_{lm}$ among themselves but leaves this sum untouched. Experiment 5 turns that closure into a usable invariant.
Experiment 5: Per-band energy is rotation-invariant¶
Take the same cap but point it at two different axes — the north pole, and a tilted axis. The two functions are rotations of each other, so their individual coefficients $a_{lm}$ differ wildly. But the per-band energy $E_l=\sum_m a_{lm}^2$ should be identical, because rotation only mixes coefficients within a band (Experiment 4). This $E_l$ is a rotation-invariant fingerprint of a function on the sphere.
def cap_band_energy(axis, theta0, Lmax):
axis = axis / np.linalg.norm(axis)
dot = DX * axis[0] + DY * axis[1] + DZ * axis[2]
f = (dot > np.cos(theta0)).astype(float)
E = []
for l in range(Lmax + 1):
e = 0.0
for m in range(-l, l + 1):
a = sphere_integral(f * real_sph_harm(l, m, TH, PH))
e += a * a
E.append(e)
return np.array(E)
L_MAX5 = 8
cap = np.deg2rad(40.0)
E_pole = cap_band_energy(np.array([0.0, 0.0, 1.0]), cap, L_MAX5)
E_tilt = cap_band_energy(np.array([0.6, 0.5, 0.6]), cap, L_MAX5)
ls = np.arange(L_MAX5 + 1)
width = 0.38
fig, ax = plt.subplots(figsize=(9, 5))
ax.bar(ls - width / 2, E_pole, width, label='cap @ north pole', color='#4C72B0')
ax.bar(ls + width / 2, E_tilt, width, label='cap @ tilted axis', color='#DD8452')
ax.set_xlabel('band $l$')
ax.set_ylabel(r'band energy $E_l = \sum_m a_{lm}^2$')
ax.set_title('Per-band energy is rotation-invariant (same cap, two orientations)')
ax.set_xticks(ls)
ax.legend()
plt.tight_layout()
plt.show()
rel = np.max(np.abs(E_pole - E_tilt) / (E_pole + 1e-12))
print('max relative band-energy difference:', float(rel))
max relative band-energy difference: 0.24509552444659496
What you see¶
The blue and orange bars match band-for-band, even though one cap sits on the pole and the other is tilted into a generic direction (the small residual difference is grid-discretization error — the tilted cap's edge is not aligned with the $\theta$-$\phi$ lattice). The energy concentrates in the low bands and decays with $l$, the spectral signature of a smooth-ish, compactly-supported function. The takeaway: $\{E_l\}$ is a rotation-invariant descriptor — you can recognize or compare directional signals regardless of how the scene is oriented. This is exactly the property that makes SH the natural language for $SO(3)$-equivariant features on 3D and directional data.
Summary / Key takeaways¶
- Basis lobes (Exp 1): $Y_{lm}$ form $2l+1$ patterns per band whose oscillation rate grows with $l$ — the visual meaning of the eigenvalue $-l(l+1)$.
- Orthonormality (Exp 2): the Gram matrix is the identity, so coefficients come from a single projection $a_{lm}=\int f\,Y_{lm}\,d\Omega$.
- Ringing (Exp 3): a global basis cannot represent a sharp edge cleanly; truncation gives Gibbs-style overshoot — SH suit smooth angular signals.
- Addition theorem (Exp 4): $\sum_m Y_{lm}(x)Y_{lm}(y)$ depends only on the angle $x\cdot y$ — each band is a rotation-closed unit.
- Rotation invariance (Exp 5): per-band energy $E_l=\sum_m a_{lm}^2$ is unchanged by rotation, giving an $SO(3)$-invariant fingerprint.
These facts are the foundation under the tier-3B note, where the same expansion encodes view-dependent color in 3D Gaussian Splatting.