🤖 machine learning

Spherical Harmonic Optimal Transport: Application to Climate Models Comparisons

This paper establishes the theoretical convergence of heat kernel-based optimal transport on manifolds and introduces a fast, GPU-friendly Spherical Harmonic Sinkhorn algorithm with $\mathcal{O}(n^{3/2})$ complexity to efficiently compare global climate models on the 2-sphere.

Original authors: Pierre Houédry, Iskander Legheraba, Léo Buecher, Nicolas Courty

Published 2026-05-19

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Pierre Houédry, Iskander Legheraba, Léo Buecher, Nicolas Courty

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to compare two globes of the Earth. One globe represents the "real" weather patterns (like actual rainfall), and the other represents a computer simulation of the weather. Your goal is to figure out how different they are.

The Problem: The "Flat Map" Trap
Most computer programs compare these globes by flattening them out into a square map first. This is like trying to compare a basketball to a square piece of paper by squishing the ball flat. You distort the distances: places that are actually far apart on the globe might look close together on the flat map, and the shapes get stretched.

Furthermore, traditional math tools for comparing these globes are incredibly slow and expensive, like trying to move every single grain of sand on a beach one by one to see how the piles match up. They also assume the two globes have the exact same total amount of "stuff" (rain), which isn't true in real life—some models predict too much rain globally, while others predict too little.

The Solution: A Heat-Based "Diffusion" Trick
The authors, Pierre Houédry and his team, created a new, faster way to compare these spherical weather maps. They call it Spherical Harmonic Optimal Transport (SHOT).

Here is how they did it, using a few creative analogies:

The Heat Analogy (The "Warm Blanket"):
Instead of trying to calculate the exact distance between every single point on the globe (which is slow), they imagine wrapping the globe in a warm blanket. They ask: "If I drop a hot spot on this globe, how does the heat spread out over a tiny moment in time?"
Mathematically, this "heat spreading" (called the heat kernel) is much easier to calculate than the exact distance. The paper proves that if you let this heat spread for a very, very short time, it becomes a perfect stand-in for the actual distance between points. It's like using the smell of a perfume to guess how far away the bottle is, rather than walking the whole distance.
The Musical Analogy (The "Spherical Symphony"):
To make this calculation super fast, they use a mathematical tool called Spherical Harmonics. Think of the globe as a drum. When you hit a drum, it vibrates in specific patterns (notes).
- Old way: Trying to calculate how the drum skin moves by checking every single point on the skin individually.
- Their way: Breaking the movement down into musical notes (frequencies). They can calculate how the "heat" spreads by just adjusting the volume of these musical notes. This allows them to use powerful computer chips (GPUs) to do the math in parallel, making it thousands of times faster.
The "Unbalanced" Scale:
Traditional methods act like a strict scale that only works if you put exactly 1kg of apples on both sides. If one side has 1.2kg, the math breaks or gives a weird result.
The authors' method is like a smart scale that can say, "Okay, you have 20% more apples on this side. Let's compare the arrangement of the apples, but also account for the extra weight." This is crucial for climate models, which often predict the wrong total amount of rain globally.

What They Found (The "Climate Detective" Work)
They tested this new tool on real-world data: comparing 8 different computer climate models against actual weather observations (ERA5).

Better Rankings: When they used the old "strict scale" (balanced) method, some models looked good just because they got the total rain amount right by accident, even if the rain was falling in the wrong places. The new "smart scale" (unbalanced) revealed the truth: it separated models that got the location of rain right from those that just got the amount right.
The Arctic Mystery: They discovered something interesting about the Arctic. In winter, the models struggled to predict rain near the edge of the sea ice. The new tool showed a clear "bias map" (a gradient) that highlighted exactly where the models were wrong. It turned out the models were pushing the rain too far out, likely because they didn't simulate the sea ice edge correctly. In summer, this problem disappeared.
Seasonal Secrets: The tool showed that while some models were consistently bad at predicting tropical rain, others were specifically bad at predicting winter storms in the Northern Hemisphere.

In Summary
The paper presents a new mathematical "lens" that lets scientists compare weather globes quickly and accurately without flattening them. It uses the physics of heat diffusion and musical frequencies to solve a problem that was previously too slow to do at high resolutions. This allows climate scientists to pinpoint exactly where and why their weather models are failing, distinguishing between models that are just "lucky" with total rain amounts and those that truly understand the geography of the storm.

Technical Summary: Spherical Harmonic Optimal Transport for Climate Model Comparisons

Problem Statement

Optimal Transport (OT) offers a principled framework for comparing probability measures by respecting the geometry of their support, making it ideal for data residing on curved manifolds like the 2-sphere ( $S^2$ ). However, standard OT computations are prohibitively expensive ( $O(n^3)$ ) for large-scale applications. While entropic regularization (Sinkhorn algorithm) improves scalability, its implementation on manifolds faces two major hurdles:

Theoretical Limitations: The standard Gibbs kernel (based on squared geodesic distance) often fails to be positive definite on curved manifolds, invalidating desirable properties of Sinkhorn divergences.
Computational Bottlenecks: Efficient implementations in Euclidean space (e.g., FLASHSINKHORN) rely on inner-product structures unavailable for geodesic distances on $S^2$ . Existing manifold approaches often rely on generic numerical solvers that discard the manifold's harmonic structure or use sparse linear algebra (e.g., Cholesky decomposition) that is inefficient on GPUs.

Furthermore, in climate science, evaluating Global Climate Models (GCMs) requires comparing precipitation fields defined on $S^2$ . Standard metrics (e.g., RMSE) ignore spherical geometry and distributional properties, while existing OT-based metrics often assume balanced mass, failing to account for systematic global mean biases inherent in GCMs relative to observations.

Methodology

The paper proposes Spherical Harmonic Optimal Transport (SHOT), a framework that leverages the heat kernel and spherical harmonic transforms (SHT) to compute Sinkhorn divergences efficiently on $S^2$ .

Theoretical Framework

Heat Kernel Regularization: Instead of the Gibbs kernel $e^{-d^2/\varepsilon}$ , the authors use the heat kernel $H_t$ of the manifold. They establish (Theorem 3.1) that the entropic unbalanced OT cost using the heat kernel converges to the classical OT cost as the time parameter $t \to 0$ (specifically $t = \varepsilon/4$ ).
Validity on $S^2$ : The authors prove (Proposition 3.2) that for the 2-sphere, the resulting Sinkhorn divergence retains key properties: symmetry, positive definiteness, smoothness, convexity, and the ability to metrize weak convergence.
Convolution Equivalence: A key theoretical insight (Proposition 3.4) is that for radial kernels (like the heat kernel) on $S^2$ , the integral operator reduces to a spherical convolution. This allows the use of the Spherical Harmonic Transform (SHT) to diagonalize the convolution operation.

Algorithm (SHOT)

The proposed algorithm replaces the dense matrix-vector multiplications of standard Sinkhorn iterations with fast spherical convolutions:

Discretization: The sphere is discretized using either the Clenshaw-Curtis (equiangular) or HEALPix (equal-area) grids.
Iteration: Each Sinkhorn iteration involves:
- Forward SHT of the dual variables.
- Pointwise multiplication by precomputed heat kernel weights $e^{-t\ell(\ell+1)}$ in the frequency domain.
- Inverse SHT to return to the spatial domain.
Complexity: This approach reduces memory requirements to $O(n)$ (no kernel matrix storage) and computational cost per iteration to $O(n^{3/2})$ (dominated by the SHT). Crucially, all operations are dense and parallelizable, making the algorithm fully GPU-friendly.

Key Contributions

Theoretical Justification: The paper rigorously proves that heat kernel-based Sinkhorn divergences converge to the true OT cost in both balanced and unbalanced regimes on $S^2$ , and that they satisfy the necessary conditions to define a valid metric-like loss function.
Efficient Algorithm (SHOT): The derivation of a fast Sinkhorn algorithm for $S^2$ that exploits the spectral properties of the Laplace-Beltrami operator. It achieves $O(n)$ memory and $O(n^{3/2})$ time complexity per iteration, enabling high-resolution computations on GPUs.
Climate Model Evaluation Framework: The application of this unbalanced OT framework to compare CMIP6 GCMs against ERA5 reanalysis data. The method decomposes model error into spatial (transport) and intensity (mass) components, addressing the global mean biases often present in climate models.

Experimental Results

Computational Efficiency

Scalability: On synthetic data, SHOT successfully handles resolutions up to $L=512$ (approx. 523,000 points) on both Clenshaw-Curtis and HEALPix grids without running out of memory, whereas methods relying on explicit kernel matrices (standard Sinkhorn, Cholesky-based approaches) fail at much lower resolutions.
Runtime: SHOT is competitive with sliced Wasserstein methods (LSSOT) and significantly faster than exact OT solvers. On Clenshaw-Curtis grids, it is the second-fastest method for high resolutions. On HEALPix, it scales well, though absolute speed is currently limited by the lack of optimized CUDA kernels in the s2fft library compared to torch-harmonics.

Climate Model Evaluation

Unbalanced vs. Balanced: The unbalanced formulation ( $\tau < \infty$ ) reveals significant differences in model rankings compared to the balanced case. For instance, the GFDL-ESM4 model improves from rank 5 to rank 1 when intensity biases are decoupled from spatial errors, while NorESM2-MM degrades, indicating its "good" balanced score was an artifact of compensating errors.
Regional Attribution: By analyzing the gradient of the Sinkhorn divergence, the authors identify specific geographic sources of bias.
- Seasonality: Tropical and ITCZ biases are found to be structurally season-independent.
- Arctic Bias: A strong seasonal signal is observed in the Arctic, where bias is significantly higher in winter (DJF) than summer (JJA).
- Mechanism: The gradient maps reveal that the winter bias is spatially organized around the sea-ice edge, suggesting a coupling between precipitation errors and sea-ice representation (e.g., ice-edge displacement). This mechanistic insight is not accessible via scalar metrics or standard difference maps.

Significance and Claims

The paper claims to provide a theoretically grounded and computationally efficient tool for comparing measures on the sphere. Its primary significance lies in:

Bridging Theory and Practice: It validates the use of heat kernels for OT on manifolds, overcoming the theoretical issues of Gibbs kernels on $S^2$ .
Enabling High-Resolution Analysis: By leveraging spherical harmonics and GPU acceleration, it makes optimal transport feasible for the high-resolution data typical of modern climate science.
Physical Interpretability: In the context of climate modeling, the method offers a "principled, geometrically faithful, and physically interpretable" framework. It moves beyond scalar rankings to provide gradient-based maps that attribute model errors to specific regions and physical mechanisms (e.g., sea-ice coupling), offering actionable insights for model development.

The authors remain modest, noting that numerical instabilities can occur for very small regularization parameters $\varepsilon$ due to the spiking behavior of the truncated heat kernel, and that the current climate study serves as an illustration of potential rather than a comprehensive assessment of all CMIP6 models. Future work is suggested to extend these methods to other manifolds and to use the differentiable divergence as a training objective for neural networks on spherical data.