---
title: 'Healpy fit_dipole: Cosmic Dipole Estimator'
url: https://www.emergentmind.com/topics/healpy-fit_dipole-algorithm
type: topic
---

# Healpy fit_dipole: Cosmic Dipole Estimator

Searching arXiv for the specified paper and closely related context.
arxiv_search.query({"search_query":"ti:\"Cosmic Multipoles in Galaxy Surveys II: Comparing Different Methods in Assessing the Cosmic Dipole\" OR id:2509.20651","start":0,"max_results":5,"sort_by":"submittedDate","sort_order":"descending"})
`healpy.fit_dipole` denotes, in the comparative framework of "Cosmic Multipoles in Galaxy Surveys II: Comparing Different Methods in Assessing the Cosmic Dipole" [2509.20651], a monopole+dipole least-squares estimator applied to a pixelized sky map of source number densities. In that treatment, the algorithm is used to determine the number count dipole from cosmological surveys and is analyzed as a representative estimator-based alternative to a Bayesian dipole analysis. The central result is that the estimator is unbiased in its native Cartesian parameterization of the dipole vector, but that its variance becomes large for low sky coverage, so practical reliability degrades in that regime. The same study further shows that apparent biases in dipole amplitude and direction can arise when outputs are transformed to spherical parameters and averaged there, even when the underlying Cartesian estimator remains unbiased [2509.20651].

## 1. Problem setting and fitted sky model

The algorithm is applied to a HEALPix map of pixelized number-count density,
\[
\mathcal N_i \equiv \text{number density / counts in pixel } i,
\]
with pixel \(i\) centered on the unit-vector direction
\[
\hat{\mathbf p}_i = (x_i,y_i,z_i).
\]
In the simulations used for comparison, the sky is pixelized with `nside = 64`, giving
\[
n_{\rm pix} = 49152
\]
equal-area pixels. The quantity of interest is the dipole modulation of the number-count map, represented internally by four coefficients,
\[
(A_0, A_1, A_2, A_3),
\]
where \(A_0\) is the monopole or mean number density and \((A_1,A_2,A_3)\) are the Cartesian dipole components of the modulated count field [2509.20651].

The assumed sky model is a linear monopole+dipole field on the sphere,
\[
\mathcal N_i \approx A_0 + A_1 x_i + A_2 y_i + A_3 z_i,
\]
or equivalently
\[
\mathcal N(\hat{\mathbf p}) = A_0 \left[1 + \mathbf{\mathcal D}\cdot \hat{\mathbf p}\right].
\]
The inferred dimensionless dipole vector is
\[
\mathbf{\mathcal D} = \left(\frac{A_1}{A_0}, \frac{A_2}{A_0}, \frac{A_3}{A_0}\right)
\equiv (\mathcal D_x,\mathcal D_y,\mathcal D_z).
\]
This formulation places the estimator squarely within the standard first-order dipolar modulation model for number-count dipoles.

The paper also situates the fitted field in the broader cosmological dipole problem. The expected kinematic number-count modulation is quoted as
\[
\frac{\Delta N(\mathbf n)}{N} = [2 + x(1+\alpha)]\,\frac{\mathbf v}{c}\cdot \mathbf n,
\]
where \(\mathbf n\) is the line-of-sight unit vector, \(\mathbf v\) is the observer velocity, \(c\) is the speed of light, \(x\) is the slope of cumulative counts \(N(>S)\propto S^{-x}\), and \(\alpha\) is the spectral index defined by \(S\propto \nu^{-\alpha}\). In this setting, the increase in discordance between the number count dipole and the CMB's kinematic dipole motivates a careful comparison of methods for dipole determination.

## 2. Least-squares formulation and regression structure

In the study, `fit_dipole` is explicitly interpreted as a least-squares estimator and, under Gaussian residual assumptions with uninformative priors, as the maximum of a Gaussian likelihood and the posterior mode. Its objective function is
\[
U = \sum_i \left[\mathcal N_i - (A_0 + A_1x_i + A_2y_i + A_3z_i)\right]^2,
\]
with the sum taken over unmasked pixels only. The stationarity equations are written as
\[
\sum_i \left[\mathcal N_i \mu_{j} - (A_0\mu_j + A_1\mu_j x_i + A_2\mu_j y_i + A_3\mu_j z_i)\right]=0,
\]
where \(\mu_j=\{1,x,y,z\}\). This is the normal-equations form of an ordinary linear least-squares problem [2509.20651].

Appendix A rewrites the fit in matrix form as
\[
M A = B,\qquad A=M^{-1}B,
\]
with
\[
A= \begin{bmatrix} A_0\\A_1\\A_2\\A_3 \end{bmatrix},
\]
\[
M= \begin{bmatrix}
\sum_i 1 & \sum_i x_i & \sum_i y_i & \sum_i z_i\\
\sum_i x_i & \sum_i x_i^2 & \sum_i x_i y_i & \sum_i x_i z_i\\
\sum_i y_i & \sum_i x_i y_i & \sum_i y_i^2 & \sum_i y_i z_i\\
\sum_i z_i & \sum_i x_i z_i & \sum_i y_i z_i & \sum_i z_i^2
\end{bmatrix},
\]
and
\[
B= \begin{bmatrix}
\sum_i \mathcal N_i\\
\sum_i \mathcal N_i x_i\\
\sum_i \mathcal N_i y_i\\
\sum_i \mathcal N_i z_i
\end{bmatrix}.
\]
This makes the connection to linear regression explicit: the fitted basis is \(\{1,x,y,z\}\), with design-matrix rows \(X_i=(1,x_i,y_i,z_i)\).

For a full, sufficiently symmetric sky, the paper notes the approximations
\[
\sum_i x_i \approx \sum_i y_i \approx \sum_i z_i \approx 0,
\qquad
\sum_i x_i y_i \approx \sum_i x_i z_i \approx \sum_i y_i z_i \approx 0,
\]
so that \(M\) becomes approximately diagonal. In that case,
\[
A_0 = \frac{\sum_i \mathcal N_i}{N},\quad
A_1 = \frac{\sum_i \mathcal N_i x_i}{\sum_i x_i^2},\quad
A_2 = \frac{\sum_i \mathcal N_i y_i}{\sum_i y_i^2},\quad
A_3 = \frac{\sum_i \mathcal N_i z_i}{\sum_i z_i^2}.
\]
Substituting the exact dipole template recovers the input dipole vector in the noiseless full-sky case, which clarifies why the method is unbiased in that limit.

## 3. Parameterization, coordinate transformation, and simulation setup

The template map used in the simulations is
\[
\mathcal N_i = \bar{\mathcal N}\left(1+\mathcal D \cos\theta_i\right),
\]
where \(\bar{\mathcal N}\) is the mean counts per pixel, \(\mathcal D\) is the dipole amplitude, and \(\theta_i\) is the angular separation between pixel \(i\) and the dipole direction. Writing the dipole direction as
\[
\hat{\mathbf d}=(x_0,y_0,z_0),
\]
and the pixel direction as
\[
\hat{\mathbf p}_i=(x_i,y_i,z_i),
\]
one has
\[
\cos\theta_i=\hat{\mathbf d}\cdot \hat{\mathbf p}_i = x_0x_i+y_0y_i+z_0z_i,
\]
so
\[
\mathcal N_i = \bar{\mathcal N}\left(1+\mathcal D x_0x_i+\mathcal D y_0y_i+\mathcal D z_0z_i\right).
\]
This identifies
\[
A_0=\bar{\mathcal N},\qquad
A_1=\bar{\mathcal N}\mathcal D x_0,\qquad
A_2=\bar{\mathcal N}\mathcal D y_0,\qquad
A_3=\bar{\mathcal N}\mathcal D z_0.
\]
The natural estimator outputs are therefore the Cartesian dipole components \((\mathcal D_x,\mathcal D_y,\mathcal D_z)\) rather than spherical amplitude-direction coordinates [2509.20651].

The corresponding spherical parameters are
\[
\mathcal D = \sqrt{\mathcal D_x^2+\mathcal D_y^2+\mathcal D_z^2},
\]
\[
\mathcal D_x = \mathcal D \cos b \cos l,\qquad
\mathcal D_y = \mathcal D \cos b \sin l,\qquad
\mathcal D_z = \mathcal D \sin b,
\]
with
\[
l = \operatorname{atan2}(\mathcal D_y,\mathcal D_x),\qquad
b = \arcsin\!\left(\frac{\mathcal D_z}{\mathcal D}\right).
\]
The paper stresses that Cartesian space is the correct parameterization in which to study bias and covariance, because the mapping to \((\mathcal D,l,b)\) is nonlinear.

The main benchmark simulations adopt the following settings: `nside = 64`, \(49152\) pixels, mean count per pixel
\[
\bar{\mathcal N}=50,
\]
fiducial dipole amplitude
\[
\mathcal D=0.007,
\]
and fiducial direction
\[
(l,b)=(264.021^\circ,48.253^\circ),
\]
corresponding to the CMB dipole direction. Noisy realizations are generated by Poisson sampling,
\[
P(\mathcal N_i\mid \lambda_i)=\frac{\lambda_i^{\mathcal N_i}e^{-\lambda_i}}{\mathcal N_i!},
\qquad
\lambda_i=\bar{\mathcal N}(1+\mathcal D\cos\theta_i).
\]
The fitted maps are therefore pixelized number-count maps contaminated primarily by Poisson shot noise, even though the least-squares objective corresponds to Gaussian residual assumptions.

## 4. Masking, coverage dependence, and frequentist behavior

The estimator is run on masked HEALPix maps by setting masked pixels to `UNSEEN`. The study tests Galactic plane masks with
\[
g_{\rm mask}=0^\circ,10^\circ,\dots,80^\circ,
\]
equatorial polar-cap masks with radii
\[
\Delta_s = 15^\circ,30^\circ,\dots,150^\circ,
\]
and a discontinuous survey mask comprising many small patches below declination \(-10^\circ\). For the discontinuous case, the mean count per pixel is increased to
\[
\bar{\mathcal N}\approx 2300,
\]
specifically to study sparse or discontinuous sky coverage at high number density. The estimator analysis uses
\[
100{,}000
\]
mock catalogs for each mask scenario [2509.20651].

The central frequentist result is that `healpy.fit_dipole` is unbiased in its native Cartesian dipole-vector parameterization. The paper states that
\[
\left<\mathbf{\mathcal D}\right>
=
\left<(\mathcal D_x,\mathcal D_y,\mathcal D_z)\right>
\]
is consistent with the true dipole vector across different noise levels, different sky coverages, and different mask shapes tested in the analysis. In the noiseless case, the appendix argues that the estimator recovers the true dipole exactly, with negligible variance, even with masks.

That conclusion is sharply qualified by the estimator variance. For sufficient sky coverage, the variance is small and the recovered amplitudes and directions are tightly centered on the truth. As sky coverage decreases, the variance grows strongly. The abstract emphasizes that, although the estimator gives unbiased results regardless of noise levels and sky coverage, low sky coverage leads to large variance and therefore practical unreliability. For Galactic masks, the discussion highlights a critical regime around
\[
g_{\rm mask}\sim 70^\circ,
\]
where the estimator variance is large. For equatorial caps, large \(\Delta_s\) values produce amplitude distributions shifted high and a bimodal direction distribution clustering near the two equatorial poles.

The discontinuous high-density benchmark is an important counterexample to a simplistic geometry-based interpretation. In that case, estimator outputs remain low-variance, both averaging methods agree with the truth, and the estimator is effectively unbiased and reliable. This suggests that sky geometry alone is not fatal; the decisive quantity is the effective information content, especially sky coverage together with count density.

## 5. Apparent bias in spherical coordinates

A recurrent misconception addressed in the paper is the claim that `fit_dipole` is intrinsically biased because average recovered amplitudes and directions can drift away from the injected truth. The study argues that this effect is largely a parameterization artifact. If one converts each realization from Cartesian dipole components to spherical parameters \((\mathcal D,l,b)\) and then averages amplitude and direction directly, apparent biases emerge: amplitudes shift high, and directions can shift toward structures associated with the mask boundaries, such as the Galactic equator or the equatorial poles [2509.20651].

The explanation is the nonlinearity of the transformation
\[
(\mathcal D_x,\mathcal D_y,\mathcal D_z)\mapsto (\mathcal D,l,b).
\]
A roughly Gaussian cloud in Cartesian space becomes distorted in spherical coordinates, so in general
\[
\left< f(\mathbf{\mathcal D}) \right> \neq f\!\left(\left<\mathbf{\mathcal D}\right>\right)
\]
for nonlinear \(f\). The paper explicitly uses
\[
\mathcal D_z = \mathcal D \sin b
\]
to explain why high-amplitude solutions at low latitudes can still average to the correct \(z\)-component. The underlying estimator can therefore remain unbiased in Cartesian space even when averages in \((\mathcal D,l,b)\) appear displaced.

This distinction governs uncertainty propagation as well. The stated procedure is to estimate covariance in Cartesian components first and then transform to spherical coordinates via a Jacobian if needed. Directly reading uncertainty from transformed histograms in \((\mathcal D,l,b)\) does not correctly represent estimator uncertainty. In that sense, the paper treats Cartesian dipole-vector space as the native inferential space and spherical coordinates as a derived, potentially misleading summary representation.

## 6. Relation to Bayesian dipole analysis and methodological limits

The Bayesian comparison is based on the dipole model
\[
f(\hat{\mathbf p}_i)=1+\mathcal D\cos\theta_i,
\]
normalized over the observed sky as
\[
\hat f(\hat{\mathbf p}_i)= \frac{f(\hat{\mathbf p}_i)}{\sum_{j=1}^{n_{\rm pix}} f(\hat{\mathbf p}_j)}.
\]
The likelihood is
\[
\ln \mathcal L = \sum_{i=1}^{n_{\rm pix}} \mathcal N_i \ln \hat f(\hat{\mathbf p}_i)
= \sum_i \mathcal N_i \ln\left( \frac{f(\hat{\mathbf p}_i)}{\sum_j f(\hat{\mathbf p}_j)} \right),
\]
with \(M_0\) denoting a monopole-only model \(f=1\) and \(M_1\) a monopole+dipole model \(f=1+\mathcal D\cos\theta_i\). The Bayesian approach returns a posterior distribution over \((\mathcal D,l,b)\) and permits model comparison through
\[
\ln B_{10} = \ln \mathcal Z_1 - \ln \mathcal Z_0,
\]
with
\[
\ln B_{10}\gtrsim 2.3
\]
quoted as “strong support” for a dipole [2509.20651].

The comparison clarifies the role of `fit_dipole`. The estimator solves for a single best-fit vector, is naturally Cartesian, and does not provide built-in model comparison; uncertainty must be assessed by simulations or covariance estimation. By contrast, the Bayesian method uses an explicit likelihood for the count distribution on pixels, returns full posterior distributions, and can become inconclusive when the data are weak. This difference is most consequential at low sky coverage. There, `fit_dipole` still has unbiased mean in Cartesian space but can yield many high-amplitude realizations because of large variance, whereas Bayesian posteriors broaden and Bayes factors can fall to near zero or negative values. The study emphasizes that such inconclusive Bayesian outcomes safeguard against incorrect conclusions.

The paper also identifies several limitations and failure modes of the estimator. Low sky coverage is the dominant one, because it broadens amplitude distributions, distorts direction distributions, and undermines the reliability of single-realization estimates. Incomplete sky and higher multipoles can also induce power leakage from higher multipoles into the dipole amplitude, with worsening contamination as the masked fraction increases. Finally, `fit_dipole` is built for a single dipole fit to one map. A joint analysis across multiple catalogs would require a generalized likelihood or different estimator, written abstractly as
\[
\phi_{\rm MLE} = \arg\max_\phi \mathcal L_{\rm joint}(\phi).
\]
Accordingly, the paper’s practical recommendation is conditional rather than absolute: `healpy.fit_dipole` is acceptable and useful when sky coverage is sufficient and the analysis is conducted in Cartesian space, but Bayesian analysis is preferable when sky coverage is poor, when model comparison is required, or when a true multi-sample inference problem is being posed.

Source: https://www.emergentmind.com/topics/healpy-fit_dipole-algorithm