In “More Notes on Calculating the Spherical Harmonics”, I used a simple iterative approach to converge on relatively accurate analysis of coefficients from a pixelized image on the sphere. In this article, we will look at the other method to improve the accuracy of approximating integration with finite sums by including so-called quadrature weight factors.
Introduction¶
When we discussed analysis of the sphere in “More Notes on Calculating the Spherical Harmonics”, we noted that discretization of the spherical harmonic transform (Eqn. $\ref{eqn:analysis_integral_sum}$) to a finite pixel grid necessarily approximates the transform.
This was most evident in the ability to recover the harmonic coefficients for the uniform
The first “fix” for this loss in precision was to reanalyze the residual and add to the previous harmonic coefficients, iterating until a desired level of convergence was achieved. While trivial to implement, iterating greatly increases the computational cost of performing the analysis transform because it actually requires multiple rounds of both synthesis and analysis to converge.
The alternative to iteration which was briefly discussed was to instead optimize the pixel grid so that the summation was compatible with a Gaussian quadrature1,2 rule, and we used the Gauss-Legendre nodes (and weights) for the colatitude coordinates instead of the uniformly spaced ECP grid. Using Gauss-Legendre nodes and weights greatly improved the accuracy of analysis to nearly machine precision, but it came at the cost of prescribing the location of the pixel rings.
To be compatible with arbitrary pixel grids (such as ECP and HEALPix4 pixelizations), the remainder of this article will demonstrate calculating quadrature weights that do not require controlling the node locations.
Quadrature Weight Formalism¶
Condition for quadrature weights on the sphere¶
Generically, we can approximate a definite integral of $f(\mat{x})$ with respect to the domain $\mat{x}$ as a finite sum of $N$ discrete points ${\mat{x}_i}$
where we are free to choose the quadrature weights $w_i$. For a Gaussian quadrature scheme, there are specific node and weight pairs $(\mat{x}_i, w_i)$ which guarantee the summation is exact when $f(\mat{x})$ is a polynomial of degree $2N-1$ or less1.
Deriving the nodes and weights for Guassian quadrature is a well-described procedure in the literature, so we won’t consider that case any futher here. Rather, our interest is in dropping the correspondance to exact integration and solving only for the [generic] quadrature weights given our choice of pixel locations (i.e. the integration nodes). We will also forego trying to derive any analytic derivation of weights, instead focusing only on a numerical procedure for solving for weight factors which minimize the error between the desired integration and concrete finite summation.
More specifically for spherical harmonic analysis, we want weights such that the error is minimized in the analysis transform
where $\mO = [\theta, \phi]^\top$ is vector of coordinates on the sphere and $i$ is an index over all discrete pixels. Since the left-hand side is an integral, a prudent step is to choose $f(\mO)$ to be a convenient function5 such that we know the analytic solution to the integral. The simplest such case is $f(\mO) = Y_{00}(\mO) = 1 / \sqrt{4\pi}$, since by orthonormality:
where $\delta_{ij}$ is the Kronecker delta function.
Substituting in our choice of function and simplifying, we arrive at the approximation
which formally provides the necessary conditions to solve for the desired weights $w_i$. The remaining challenge is to translate this expression into a practical implementation where solving for the weights can be performed efficiently.
Practical considerations for solving for quadrature weights¶
In order to keep the following math relatively compact, let us switch to the notation of linear operators. Let the synthesis of a map $\mat{m}$ on the sphere given spherical harmonic coefficients $\mat{a}$ be written as
where $\bfY$ is the spherical harmonic synthesis operator. This is a relatively direct translation of Eqn. 8 from Part I — each column6 of the matrix $\bfY$ corresponds to a particular mode $(\ell, m)$ (given some canonical ordering of the tuple of indices) and contains the pixel values for the function $Y_{\ell m}(\theta, \phi)$ over the sphere. Conversely, the reverse transform (i.e. analysis) is written using the the adjoint (conjugate-transpose) operator $\bfYt$:
where each row of the matrix is the complex conjugate $\overline{Y_{\ell m}(\theta, \phi)}$ for a particular $(\ell, m)$. Using this new notation, Eqn. $\ref{eqn:minimize_sum}$ is rewritten as
where $\mat{w}$ is the vector of quadrature weights and $\mat{a_0}$ is the complex vector $[1, 0, 0, \ldots]$ corresponding to the Kronecker delta product $\delta_{\ell0}\delta_{m0}$.
To actually solve the system of equations, a simple choice7 is to apply the conjugate gradient method, with the important caveat that the operator must be positive-definite for the CG method to apply. The transform $\bfYt$ is manifestly not positive-definite because it is not square. (The number of spherical harmonic modes need not equal the number of map pixels in general.) Thankfully, there is a simple solution — for an arbitrary operation $\mat A$, the operators $\mat{A}\mat{A}^\dagger$ and $\mat{A}^\dagger \mat{A}$ are at least positive-semidefinite.
In our case, it is convenient to create the operator $(\bfYt \bfY)$ by rewriting the weight map $\mat{w}$ as the synthesis of its corresponding harmonic coefficients:
where $\mat{w} \equiv \bfY \mat{\hat w}$ and $(\bfYt\bfY)$ is a positive-definite operator.
The final practical concern is the matter of turning the matrix of spherical harmonic coefficients into a vector. The natural choice is to simply unravel the matrix in its native ordering. For example, an $\ell_{\mathrm{max}} = 2$ set of coefficients unravels in Julia’s column-major order as:
(Remember that the spherical harmonics are defined for $|m| \le \ell$.) While this will work, unravelling the matrix results in almost twice as many entries as is strictly necessary since everything above the main diagonal is unused. The obvious improvement is to unravel only the lower-triangle of used elements into a vector, thereby reducing the effective dimensionality of the system of equations to solve.
Taken literally, this ordering of the harmonic coefficients also defines the order of the columns/rows in the synthesis/analysis operators, respectively. As we will see in the next section, though, we only ever invoke this ordering on the harmonic vectors.
Calculating Quadrature Weights¶
With the formalism all in place, we are now prepared to implement solving for appropriate quadrature weights for
our pixelization of choice on the sphere.
First, the spherical harmonic transforms as presented in Part III are assumed (or see the attached
sht.jl script which collects the implementation).
Furthermore, the unraveling of the spherical harmonic coefficients matrix to a vector and back requires new helper
functions, but because these are relatively straight forward, we will skip a detailed description of them.
In what follows, the unravelling operation is performed by packtril! which packs the lower triangle of a matrix
into a vector, and its inverse is unpacktril!.
The button to the right expands a hidden section with the implementation of these functions, plus a pair of
functions trillen and trildim which calculates the vector length and matrix dimensions given the other
sizing information, respectively.
Included in this aside are a number of utility helper functions which we will make use of in the rest of the article. Their implementation should be relatively straightforward and simple to rederive from their defined behavior alone, but we include them here in full for completeness and clarity.
The number of unique elements in the lower-triangle of an $n \times m$ matrix (where $m \le n$) is $L = \frac{n(n+1)}{2} - \frac{(n-m)(n-m+1)}{2}$.
"""
L = trillen(n::Int, m::Int = n)
Calculates the length `L` of a vector required to store the elements of an `n` × `m`
lower-triangular matrix.
"""
function trillen(n::Int, m::Int = n)
m <= n || throw(ArgumentError("require n ≤ m; got n = $n, m = $m"))
k = n - m
return n * (n + 1) ÷ 2 - k * (k + 1) ÷ 2
end
Conversely, the size of the matrix corresponding to $L$ packed elements depends on knowing something about the width of the matrix. If the matrix is square and we can assume $n = m$, then $n = \frac{\sqrt{8L + 1} - 1}{2}$. Otherwise if the width and height may be independent, we must be explicitly told the width $m$, and then the height is $n = \frac{2L + m(m - 1)}{2m}$.
"""
n, m = trildim(L::Int[, m::Int])
Calculates the dimensions `n` × `m` of a matrix capable of storing the `L` elements of a
lower-triangular matrix. If `m < n`, then `m` must be provided.
"""
function trildim end
function trildim(L::Int)
n′ = (sqrt(8L + 1) - 1) / 2
n = trunc(Int, n′)
n == n′ || throw(ArgumentError("length $L incompatible with inferred triangle dims ($n, $n)"))
return (n, n)
end
function trildim(L::Int, m::Int)
n′ = (2L + m^2 - m) / 2m
n = trunc(Int, n′)
n == n′ || throw(ArgumentError("length $L incompatible with inferred triangle dims ($n, $m)"))
return (n, m)
end
Transforming back-and-forth from packed vectors and unpacked matrices is then performed by relatively simple loops:
"""
vec = packtril!(vec::AbstractVector, mat::AbstractMatrix)
Packs the lower triangle of the matrix `mat` into the vector `vec`.
"""
function packtril!(vec::AbstractVector, mat::AbstractMatrix)
L, (n, m) = length(vec), size(mat)
L == trillen(n, m) || throw(DimensionMismatch(
"cannot pack size ($n, $m) triangular matrix into length $L vector"))
i₀, j₀, k₀ = firstindex(mat, 1), firstindex(mat, 2), firstindex(vec, 1)
kk = 0
@inbounds for jj in 0:m-1, ii in jj:n-1
vec[k₀+kk] = mat[i₀+ii, j₀+jj]
kk += 1
end
return vec
end
"""
mat = unpacktril!(mat::AbstractMatrix, vec::AbstractVector)
Unpacks the vector `vec` into the lower triangle of the matrix `mat`.
"""
function unpacktril!(mat::AbstractMatrix, vec::AbstractVector)
L, (n, m) = length(vec), size(mat)
L == trillen(n, m) || throw(DimensionMismatch(
"cannot unpack length $L vector into size ($n, $m) triangular matrix"))
i₀, j₀, k₀ = firstindex(mat, 1), firstindex(mat, 2), firstindex(vec, 1)
kk = 0
@inbounds for jj in 0:m-1, ii in jj:n-1
mat[i₀+ii, j₀+jj] = vec[k₀+kk]
kk += 1
end
return mat
end
For convenience, we can also provide versions which allocate the necessary output vector/matrix, inferring the necessary dimensions where possible.
"""
vec = packtril(mat::AbstractMatrix)
Packs the lower triangle of matrix `mat` into the vector `vec`.
"""
function packtril(mat::AbstractMatrix)
vec = fill!(similar(mat, trillen(size(mat)...)), zero(eltype(mat)))
return packtril!(vec, mat)
end
"""
mat = unpacktril(vec::AbstractVector[, m::Int])
Unpacks the vector `vec` into the lower triangle of a matrix `mat`, inferring the dimensions
based on the length `L = length(vec)` (and `m = size(mat, 2)` if necessary) with
[`trildim`](@ref).
"""
function unpacktril end
function unpacktril(vec::AbstractVector)
mat = fill!(similar(vec, trildim(length(vec))...), zero(eltype(vec)))
return unpacktril!(mat, vec)
end
function unpacktril(vec::AbstractVector, m::Int)
mat = fill!(similar(vec, trildim(length(vec), m)...), zero(eltype(vec)))
return unpacktril!(mat, vec)
end
To perform the conjugate gradient solve, we will make use of the cg method from the
IterativeSolvers.jl package8.
This requires expressing the transforms as a linear operator (applied via “matrix”-vector multiplication) — we
could do this by implementing the appropriate AbstractMatrix methods, but a more expedient option is to wrap
a simple function with LinearMaps.jl.
The following YᵀY function constructs a LinearMap object which corresponds to the
$\left(\bfYt \bfY\right)$ mathematical operator for a given pixelization and range of harmonic coefficients.
using LinearMaps
"""
YᵀY(maprings::MapRings, lmax::Int, mmax::Int = lmax)
Produces the successive synthesis-analysis linear operator \$(Y^\\dagger Y)\$ operating
upon a packed vector of harmonic coefficients up to \$(ℓ_{max}, m_{max})\$ for the
map described by `maprings`.
"""
function YᵀY(maprings::MapRings, lmax::Int, mmax::Int = lmax)
npack = trillen(lmax + 1, mmax + 1)
alms = zeros(ComplexF64, lmax + 1, mmax + 1)
function f!(alm_out::AbstractVector{<:Complex}, alm_in::AbstractVector{<:Complex})
# unpack the vector into matrix form
unpacktril!(alms, alm_in)
# do the round trip synthesis -> analysis on the provided alms
alms′ = analyze(maprings, synthesize(maprings, alms), lmax, mmax)
# then repack the new alms into the output vector
return packtril!(alm_out, alms′)
end
# wrap f! to make it a linear operator
return LinearMap{ComplexF64}(f!, f!, npack, npack,
issymmetric = true, isposdef = true, ismutating = true)
end
Solving for weights then proceeds by setting up the appropriate $\mat{a_0}$ vector, running the conjugate gradient method, and synthesizing the quadrature weight map $\mat w$ from the solution of harmonic coefficients $\mat{\hat w}$:
using IterativeSolvers: cg
"""
w = solve_weights(maprings::MapRings, lmax::Int, mmax::Int = lmax)
Solve for the quadrature weight map `w` (in the pixelization described by
`maprings`) which minimizes the the error in the equation
\$\$
Yᵀ w ≈ √(4π) δ_ℓ0 δ_m0
\$\$
for spherical harmonic modes up to `(lmax, mmax)`.
"""
function solve_weights(maprings::MapRings, lmax::Int, mmax::Int = lmax)
# a_ℓm = √(4π) δ_ℓ0 δ_m0
alm = zeros(ComplexF64, lmax + 1, mmax + 1)
alm[1, 1] = sqrt(4π)
a₀ = packtril(alm)
# run conjugte gradient descent to high relative convergence
op = YᵀY(maprings, lmax, mmax)
ŵ = cg(op, a₀, reltol = eps(sqrt(4π))^0.9)
# return the synthesize weight map
return synthesize(maprings, unpacktril!(alm, ŵ))
end
We can now finally solve for the analysis quadrature weights for our pixelization scheme of choice. Like the previous articles in this series, we’ll use a $50 \times 100$ pixel ECP grid as our test case.
nθ = 50
ringinfo = ecp_maprings(nθ, 2nθ)
test_alm = [1.0 0
0 1.0im]
test_map = synthesize(ringinfo, test_alm)
A new(ish) consideration is what $\ell_{\mathrm{max}}$ to choose when solving for the quadrature weights. The obvious condition is that the quadrature limits need to be at least as high as the highest multipoles used in any future analysis, but that only delays the question to what limits to use in analysis. A specific answer is outside the scope of this article9, but we can easily argue that an upper bounds is $\ell_\mathrm{max} + 1 \le n_\theta$ on the grounds of matching the $n=\ell_\mathrm{max}+1$ degrees of freedom in Legendre polynomial ($\lambda^\ell_0(x)$) amplitudes to the $n=n_\theta$ degrees of freedom (pixels) along a meridian of the map.
Proceeding with the maximum $\ell_\mathrm{max}$, we solve for the quadrature weights for our test map’s pixelization:
w = solve_weights(ringinfo, nθ - 1)
and the weight map for the ECP grid is easy to visualize as a matrix, so reshaping and plotting gives the image shown in Fig. 3.1.
wmap = reshape(w, nθ, 2nθ)
Notice that the weights appear to have two important symmetries — the weight map is (a) symmetric over the equator and (b) invariant to rotations in azimuth. This should not be a surprise since the ECP grid can both be flipped vertically and rotated in azimuth by any number of pixels without any observable change to the grid. This means that while in practice we have solved for $50\times100 = 5,000$ free parameters to construct the entire weight map, there are actually only $25$ independent values (corresponding to one half of a meridian), and the rest of the grid can be constructed from this subset. On larger grids (which corresponds to higher multipole ranges as well), it is prudent to take advantage of such symmetries to greatly reduce the dimensionality of the problem. For the sake of simplicity in this article, though, we will forego further optimization since the small example test grids already solve quickly.
Finally, let us compare the analyzed harmonic coefficients given the test map without quadrature weights to the case where quadrature weights are applied. First, the original case without quadrature weights…
# for highlighting non-zero elements in a matrix, exploit Julia's
# printing of sparse arrays
using SparseArrays
dropsmall(x) = droptol!(sparse(x), sqrt(eps(maximum(abs, x))))
analyze(ringinfo, test_map, 12, 3) |> dropsmall
13×4 SparseMatrixCSC{ComplexF64, Int64} with 13 stored entries:
1.00016+0.0im ⋅ ⋅ ⋅
⋅ -3.50502e-17+1.0im ⋅ ⋅
0.000368242+0.0im ⋅ ⋅ ⋅
⋅ -1.28227e-18-6.40155e-7im ⋅ ⋅
0.000495247+0.0im ⋅ ⋅ ⋅
⋅ 2.55144e-18-1.2748e-6im ⋅ ⋅
0.000597485+0.0im ⋅ ⋅ ⋅
⋅ 2.88304e-19-2.04774e-6im ⋅ ⋅
0.000686818+0.0im ⋅ ⋅ ⋅
⋅ -7.33496e-18-2.94784e-6im ⋅ ⋅
0.000768423+0.0im ⋅ ⋅ ⋅
⋅ 2.25251e-18-3.9715e-6im ⋅ ⋅
0.000845186+0.0im ⋅ ⋅ ⋅
…versus applying the quadrature weights to the map…
analyze(ringinfo, w .* test_map, 12, 3) |> dropsmall
13×4 SparseMatrixCSC{ComplexF64, Int64} with 2 stored entries:
1.0+0.0im ⋅ ⋅ ⋅
⋅ -3.59897e-17+1.0im ⋅ ⋅
⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅
In the first case without the quadrature weights, we again observe the response in higher degree multipoles than the
input harmonics, but application of the new quadrature weights [largely] eliminates the discrepancy from the known
input.
(More precisely without the truncation for presentation, the deviation between input and output is
$\operatorname{max}|\hat{a}_{\ell m} - a_{\ell m}| \approx 3\cdot 10^{-16}$ where $\hat{a}_{\ell m}$ and
$a_{\ell m}$ are the analyzed and known-input harmonic coefficients, respectively.)
As the final demonstration of how adding quadrature weights improves, we can repeat a similar comparison made in the two previous installments in this article series — we will draw a random realization of a pattern on the sphere which is consistent with a $C_\ell \propto \ell^{-2}$ power spectrum, and then we can compare the recovered power spectrum from analyzing the synthesized map both with and without appropriate quadrature weights.
function draw_realization(lmax)
alms0 = zeros(ComplexF64, lmax+1, lmax+1) # C_ℓ ∝ ℓ¯² ~ ⟨|alm|^2⟩
for ℓ in 1:lmax
alms0[ℓ+1,1] = randn(rng, Float64) / ℓ
alms0[ℓ+1,2:ℓ+1] .= randn.(rng, ComplexF64) ./ ℓ
end
return alms0
end
function reanalyze(pixelization, alms, weights = nothing)
lmax, mmax = size(alms) .- 1
map = synthesize(pixelization, alms)
if weights != nothing
map .*= weights
end
return analyze(pixelization, map, lmax, mmax)
end
alms0, alms_ecp, alms_ecpw = let nθ = 150
ringinfo = ecp_maprings(nθ, 2nθ)
alms0 = draw_realization(100) # lmax = 100 to match previous articles
quadw = solve_weights(ringinfo, nθ - 1)
alms_ecp = reanalyze(ringinfo, alms0)
alms_ecpw = reanalyze(ringinfo, alms0, quadw)
alms0, alms_ecp, alms_ecpw
end
# power spectra
Cl0 = alm2cl(alms0)
Cl_ecp = alm2cl(alms_ecp)
Cl_ecpw = alm2cl(alms_ecpw)
# fractional differences
δCl_ecp = (Cl_ecp .- Cl0) ./ Cl0
δCl_ecpw = (Cl_ecpw .- Cl0) ./ Cl0
Using quadrature weights during analysis clearly recovers the input spherical harmonics to much higher precision than the case without, and in contrast to the procedure discussed in Part II, no iteration is necessary. Near the right-hand side of the lower panel (around $\ell = 90$), small fluctuations in the dotted-green (with quadrature weights) line can be seen, though, so alternate pixelizations like the Gauss-Legendre grid should still be considered if near-perfect recovery is required. Nonetheless, generating quadrature weights is a relatively simple method to improve the analysis precision for arbitrary pixelization schemes when the particular pixelization may be chosen for reasons other than optimizing spherical harmonic transforms.
Conclusions¶
In this article, we took a closer look at the mathematics of approximating an integral with finite sums and how that impacts the ability to analyze maps on the sphere into spherical harmonic coefficients with relatively high accuracy. From previous articles, we knew that the precision of analysis could be greatly improved by modifying the pixel grid to align with Gaussian quadrature nodes, but here we considered the case where the pixel grid is unchanged (which often occurs in practice when the pixelization is chosen/restrained by other requirements).0
We then showed that the appropriate quadrature weights can be derived by solving a system of equations. While the system of equations in principle could be solved using standard linear algebra techniques, the potentially very large dimensions of any such operator motivated use of the conjugate gradient method to solve for the unknown quadrature weight map without explicitly forming the corresponding spherical harmonic transform operator as a matrix.
Finally, we implemented a quadrature weight solver in Julia and showed that on a couple of test cases of the ECP grid, we do in fact observe improved precision of analyzed maps when quadrature weighting is included. Importantly, we achieved higher precision results without needing to resort to iterative methods (using instead only a comparitively cheap multiplication of two vectors), which saves significant computational cost, especially on larger maps with higher multipole ranges.
“Gaussian Quadrature”. On: Wikipedia. url: https://en.m.wikipedia.org/wiki/Gaussian_quadrature ↩︎ ↩︎
“Legendre-Gauss Quadrature”. On: Wolfram MathWorld. url: https://mathworld.wolfram.com/Legendre-GaussQuadrature.html ↩︎
J.R. Shewchuk “An Introduction to the Conjugate Gradient Method Without the Agonizing Pain” (1994) url: https://www.cs.cmu.edu/~jrs/jrspapers.html ↩︎
K. M. Górski et al. “HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere”. In: Astrophys J. 662 (Apr 2005) ads: 2005ApJ…622..759G ↩︎ ↩︎
M. Gräf, S. Kunis, D. Potts. “On the computation of nonnegative quadrature weights on the sphere” In: Applied and Computational Harmonic Analysis, 27. (Jul 2009) doi: 10.1016/j.acha.2008.12.003 ↩︎
Note that while we use the terms “column” and “row” for matrices of finite dimensions for convenience, the operators and notation are well defined for the continuous case. See Hilbert spaces for a better discussion than will be presented here. ↩︎
We avoid traditional matrix-based linear solvers because methods such as LU or QR factorization make use of explicit matrices. In the case of high-resolution maps and high-$\ell$ spherical harmonic transforms, explicitly representing the operator $\bfYt$ would require prohibitively large amounts of RAM. Instead, the conjugate gradient method is one of the simplest examples of so-called “matrix-free” methods whereby the operator is implicitly defined in terms of its matrix-vector product, thereby avoiding the need to fill a very large, dense matrix. ↩︎
A good reference for a from-scratch implementation of a conjugate gradient solver is provided by Shewchuk3. ↩︎
A specific example of where the true bandwidth limit is lower than our generous upper bound is the HEALPix4 pixelization, for which the
anafastdocumentation states that grid supports signals up to $\ell_\mathrm{max} \sim \frac{3}{4} n_\theta$. ↩︎