GEM Summer '22 tutorial: Spherical harmonic fitting¶

Vishal Upendran, IUCAA, Pune, India¶

Catch me at: uvishal@iucaa.in

In this notebook, we will look at how Spherical harmonics may be used to fit data present on a ( $\theta$ , $\phi$ ) grid. To motivate our general scheme of fitting this data, we shall first look at fitting a Fourier series to a 1-D data, and then generalizing it to Spherical harmonics.

The method presented here is generic to fitting any basis to a given dataset.

A sparse method for performing fourier transform¶

We consider the Fourier basis, and perform a sparse fitting of the Fourier basis to a given dataset. This is just to see how well the frequencies are retrieved by such a sparse computer.

FFT baseline¶

Let us perform FFT on this signal, and generate a "baseline" model.

In [4]:

fft = np.fft.fft(signal)/len(signal)
freq = np.fft.fftfreq(len(signal))

That's how a Discrete Fourier transform works. It decomposes the signal into discrete "modes", each of which has an associated coeffficient.

Typicall we can rewrite the decomposition of any signal into a "series" solution as: $s(n) = \sum_{i=0}^{N} a_i \cos(𝜔_i n/2N) + b_i \sin(𝜔_i n/2N),$

where we have accounted for Nyquist criterion to select the number of modes. The number $n$ stands for the time steps, and there are $\mathbf{2N}$ time steps in the signal. Generally, we have built in packages to give the max number of modes given a signal length.

Generally, $s$ is a time series. Let us denote time series variables as $\vec{s}$ . We can see clearly that each of the cosine and sine series is actually a time series vector. Hence, we can write the equation as:

$\vec{s} = \sum_{i=0}^{N} a_i \vec{c}(𝜔_i/2N) + b_i \vec{s}(𝜔_i/2N),$

where I have used a shorthand notation for cosine and sine time series. What we are basically doing is adding together a bunch of different time series with some weights. Hence, we can convert this into a matrix equation like:

$\vec{s} = a_0 + a_1\times[c_{1,0} c_{1,1} ... c_{1,2N}] + ...... +a_N\times[c_{N,0} c_{N,1} ... c_{N,2N}] + b_1\times[s_{1,0} s_{1,1} ... s_{1,2N}]+ .. b_N\times[s_{N,0} s_{N,1} ... s_{N,2N}],$

where $c_{a,b}$ represents the cosine value at frequency $a$ and time point $b$ .

This expression can be succintly written as a matrix equation as:

$\vec{s} = \begin{bmatrix}c_{0,1} & ... & c_{0,N} & s_{0,1} & ... & s_{0,N} \\ & ... & \\ c_{0,2N} & ... & c_{N,2N} & s_{0,2N} & ... & s_{N,2N} \end{bmatrix}\begin{bmatrix} a_1 \\ ... \\ a_N \\ b_1 \\ ... \\ b_N \end{bmatrix}$

The $\vec{a}$ is a set of coefficients, while the matrix is what we can call the "Basis matrix" $\mathbf{B}$ . Hence, the problem boils down to finding the vector $\vec{a}$ such that: $\vec{s} = \mathbf{B}\vec{a}.$

This is a typical "fitting" problem in linear algebra. We need to define a basis, invert the matrix $\mathbf{B}$ , and figure out $\vec{a}$ . Let us look at a few ways to perform this inversion.

Lasso regression.¶

Lasso, or Least absolute shrinkage regression, attempts to solve for $\vec{a}$ such that most of the elements of $\vec{a}$ become 0. In our particular case, it will try to make as many Fourier modes = 0 as possible. Let us see what solution we get from such an exercise.

Let us define the Lasso model, and fit it.¶

Lasso model has a regularization parameter $\alpha$ which can be tuned for the amount of shrinkage to be done.

In [10]:

model = Lasso(alpha=0.1,fit_intercept=False)
model.fit(basis,signal)

Out[10]:

Lasso(alpha=0.1, fit_intercept=False)

In [11]:

model.coef_

Out[11]:

array([ 0.        ,  9.85247823, 19.888655  , ...,  0.        ,
        0.        ,  0.        ])

Spherical harmonics¶

Now that we have seen how we may write a Fourier decomposition as a linear equation, we can look at a broader picture. In the equation for decomposition, $\vec{s} = \mathbf{B}\vec{a},$ the basis matrix $\mathbf{B}$ can be composed with any "good" basis functions. For periodic signals, the Fourier basis is good. If our signal $\vec{s}$ is a field on a sphere, the ideal set of basis would be the Spherical Harmonic functions.

The spherical harmonics functions $Y_{nm} (\theta,\phi)$ are the angular parts of the solution to the Laplace equation on a spherically symmetric structure. These functions are defined as: $Y_{nm}(\theta,\phi) := \sqrt{\frac{2n+1}{4\pi} \frac{(n-m)!}{(n+m)!}} e^{i m \theta} P^m_n(\cos(\phi)).$

Here, $P^m_n(\cos(\phi))$ are the associated Legendre polyynomials.

Any scalar function $f(\theta,\phi)$ can then be decomposed into the Spherical harmonic basis (like a Fourier basis) as:

$f(\theta,\phi) = \sum_{n=0}^{\infty} \sum_{m=-n}^{n} a_{nm} Y_{nm}(\theta,\phi)$

Similar to how we approximated the Fourier series decomposition of our signal till the maximum number of modes (or frequency) N, we can limit the maximum number of modes we would like to use for fitting the 2-D function. Hence, the function now becomes: $f(\theta,\phi) \approx \sum_{n=0}^{N} \sum_{m=-n}^{n} a_{nm} Y_{nm}(\theta,\phi),$

where N is the largest mode in our expansion.

Now we will have to define the correct Basis matrix $\mathbf{B}$ . To do so, we will need to transform the 2-D indexing to a 1-D indexing.

Consider a mode with n = $n_0$ . If we consider all the modes from $n=0$ to $n=n_0$ , we will have a total of $(n_0+1)^2$ number of modes. Let the index $i_0$ correspond to the mode with $n=n_0$ and $m=0$ . Then, we may write $i_0$ as: $i_0 = (n+1)^2-n-1\implies i_0=n^2+n,$ where an extra 1 is subtracted since $m=0$ . With this, we may write the index for any $m$ for $n=n_0$ as: $i = n^2+n+m.$

Thus, we can "unroll" the 2D summation into 1D.

We will also need to put this back together in 2D. For this, observe that for $m=-n$ , we have: $i_{n,-n} = n^2 \implies n = \sqrt{i_{n,-n}}.$ If we now consider the $m=n$ mode, we would get $i_{n,n}=n^2+2n$ . This $i$ is not a perfect square, and we would get the next perfect square for $i_{n+1,-n-1} = (n+1)^2$ . This, for a given $i$ , we have: $n =floor(\sqrt{i}), m = i-n^2-n.$

With these two transformations, we can move back and forth between the 2D and 1D arrays.

With this theory, let us try to perform a Lasso fit to a toy example, and see if we recover the modes correctly.

Acknowledgement: Aspects of this notebook have been adapted from a notebook written by Mark Cheung, LMSAL.

Recipe for fixing N and $\alpha$ ¶

There are a couple of ways of fixing N and $\alpha$ . I will outline 3 methods below:

N constrained by your computational resources, and $\alpha$ from the reconstruction error.
Constrain N, $\alpha$ through cross validation.
Constrain N, $\alpha$ through validation.

We shall look at cross validation in this tutorial.

So let us sweep over $N_{max}$ and $\alpha$ . Let us run $N_{max}=$ {5,10,15,20}, and $\alpha=$ {10,1,0.1,0.01}.

In [29]:

# Let us consider
Nmax_list = [5,10,15,20]
alpha_list = [10,1,0.1,0.01]
alpha,Nmax = np.meshgrid(alpha_list,Nmax_list)
comb_param = list(zip(alpha.ravel(),Nmax.ravel()))

In [30]:

np.asarray(comb_param).shape

Out[30]:

(16, 2)

In [31]:

from sklearn.model_selection import cross_val_score

GEM Summer '22 tutorial: Spherical harmonic fitting¶

Vishal Upendran, IUCAA, Pune, India¶

A sparse method for performing fourier transform¶

FFT baseline¶

Lasso regression.¶

Let us define the Lasso model, and fit it.¶

Let's check the coefficients for different modes.¶

Spherical harmonics¶

Let us define a toy example function using Scipy.¶

Define a sparse basis to generate the spherical harmonic decomposition¶

Check the error due to decomposition¶

Recipe for fixing N and $\alpha$ ¶

Cross validation: The validation of validations ☕¶