GaussianRandomFieldS2 samples are not isotropic: orders m > 0 get twice the power of m = 0

Open
  • Is something wrong?
#349 0 comments 0 reactions 0 assignees View on GitHub

Maintainers usually reply within 1 day

Nobody has claimed this yet.

Assessment

Difficulty
3/5
Estimated time
1-2 days
Newbie friendliness
62/100
Issue type
Bug
Clarity
Clearly specified
Activity status
Active
Tech stack
python, pytorch
Domain
machine-learning

Research direction

Start by reading GaussianRandomFieldS2.forward and the inverse transform behavior at order m = 0; the issue identifies both as the source of the variance difference. Apply the proposed scaling to orders m > 0 for generated and supplied xi, then update test_spectral_amplitudes_are_computed_in_requested_dtype and add the stated variance check across the five grids. Done means the tested variance ratio is 1 within the reported numerical spread.

Written by the indexing model from the issue text.

Description

Problem

GaussianRandomFieldS2 defines its covariance through a power spectrum C_l. Such a covariance is rotation invariant, so the variance of the field should be the same at every point. Instead it depends on the latitude:

grid variance, ring nearest the pole variance, equator ratio
legendre-gauss, 24x48 0.004320 0.006850 0.631
equiangular, 33x66 0.003334 0.006434 0.518

These use alpha=1.5, tau=3.0, sigma=1.0 and float64. They are exact rather than sampled: the map from noise to field is linear, so pushing one unit per noise component through it and summing the squares gives the variance. Sampling agrees: 4000 fields on legendre-gauss 48x96 give a ratio of 0.589.

import torch
import torch_harmonics as th
from torch_harmonics.random_fields import GaussianRandomFieldS2

field = GaussianRandomFieldS2(th.as_grid("legendre-gauss", nlat=24, nlon=48), alpha=1.5, tau=3.0, sigma=1.0, dtype=torch.float64)
lmax, mmax = field.isht.lmax, field.isht.mmax
count = 2 * lmax * mmax
noise = torch.eye(count, dtype=torch.float64).reshape(count, lmax, mmax, 2)
variance = field(count, xi=torch.view_as_complex(noise)).pow(2).sum(dim=0).mean(dim=-1)
print(variance[0] / variance[12])  # 0.631, expected 1

Cause

forward draws xi with unit-variance real and imaginary parts at every order and calls self.isht(xi * self.sqrt_eig). The inverse transform drops the imaginary part at m = 0, while a coefficient of order m > 0 also stands for its conjugate at -m. So E|a_lm|^2 is C_l at m = 0 and 2 C_l at m > 0. For comparison, healpy.synalm with C_l = 1 gives E|a_16,0|^2 = 0.997 and E|a_16,8|^2 = 0.991 over 4000 draws.

Proposed fix

Scale orders m > 0 by 1 / sqrt(2) in forward, which also covers a user-supplied xi. With that change the ratio is 1.0000 on the two grids above, with a spread of about 2e-15, and I have a test for the spread on five grids (legendre-gauss with two longitude counts, equiangular, lobatto, fejer1).

This changes behavior: samples for a given seed or xi change, and the variance at the equator about halves. One existing assertion, in test_spectral_amplitudes_are_computed_in_requested_dtype, encodes the current formula and would need to change.

A question on the overall scale

This is not part of the proposed fix. After it, the variance equals sum_l C_l (2l + 1) / (4 pi)^2, not / (4 pi). The transform uses norm="backward", which legendre.py treats as four-pi, so a factor of 4 pi separates the two: 0.00354 against 0.0445 for the parameters above. The docstring says the default sigma makes the variance O(1), while the defaults alpha=2.0, tau=3.0 give 0.0058. Is that scale intended?

Would you take a PR for the fix above? This was found and reproduced with AI assistance (Claude Opus 5.5). I can follow up on review.

Dominant language
Jupyter Notebook
Stars
726
Forks
78
Avg merge
3d 13h
Merged PRs (30d)
24

Getting set up

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. Open a pull request that references the issue number.

More from NVIDIA/torch-harmonics

All issues in NVIDIA/torch-harmonics

Similar issues

More Machine Learning issues

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.