数理解析研究所講究録
第 1606 巻 2008 年 18-29

18

Fast Spherical Harmonic Transform Algorithm
based on Generalized Fast Multipole Method
Reiji Suda
the University of Tokyo / CREST, JST

1

Introduction

Spherical harmonic transform is the most important orthogonal function transform only except
Fourier transform, and is used not only for climate simulation and signal processing but also for abase of several numerical algorithms. Fast Fourier Transform (FFT), which runs in time $O(N\log N)$
is quite well known, but, for spherical harmonic transform, there is no fast algorithm which is assimple as FFT.
We have proposed a fast transform algorithm for spherical harmonic transform using Fast Multipole Method (FMM)[17, 25], and its program is publicly available as FLTSS[24, 20]. Our algorithm
computes the transform using polynomial interpolations, and to make it possible, we have introduced split Legendre functions[25], which nearly double the computational costs, and also affect thenumerical stability[21].
To reduce computational costs, we tried fast transform without using split Legendre functiors[22,23] with the help of generalized FMM[29, 28]. This algorithm (we call this new algorithm and theformer old algorethm) runs faster than the old algorithm with improved numerical stability.
This paper discusses some details of the new algorithm and its implementation.

2

Spherical Harmonic Transform

2.1

Spherical harmonic transformA spherical harmonic function
and a trigonometric function as$Y_{n}^{m}(\lambda, \mu)$can be decomposed into an associated Legendre function$Y_{n}^{m}(\lambda, \mu)=P_{n}^{m}(\mu)e^{im\lambda}$.

Thus evaluation of a spherical harmonic expansion$g( \lambda, \mu)=\sum_{m=0}^{T}\sum_{n=m}^{T}g_{n}^{m}Y_{n}^{m}(\lambda, \mu)$can also be decomposed into associated Legendre function transform and Fourier transform as,$g^{m}(\mu)$$=$$\sum_{n=m}^{T}g_{n}^{m}P_{n}^{m}(\mu)$$g(\lambda, \mu)$$=$$\sum_{m=0}^{T}g^{m}(\mu)e^{im\lambda}$

(1)

.

(2)

The inverse spherical harmonic transform is a computation of, while the sphericalfromharmonic transform is a computation of, where $1\leq j\leq J,$ $1\leq k\leq K$ , andfrom
and$K$ are usually determined by the so-called alias-free condition$g(\lambda_{j}, \mu k)$$g_{n}^{m}$$J\geq 3T+1$ ,$g(\lambda_{j}, \mu k)$$K\geq(3T+1)/2$ .

$g_{n}^{m}$

$J$

(3)

19
For asymptotic complexity analysis, we can simply assume that $J=O(T)$ and $K=O(T)$ . The
inverse transform consists of the associated Legendre function transform and Fourier transform,
and the Fourier transform (2) can be computed in time $O(T^{2}\log T)$ by using FFT for each.
Conventional method of the associated Legendre function transform is the direct summation of
the form of (1), thus requires time $O(T^{3})$ . Let us call this method direct $\omega mputation$. By using
direct computation for associated Legendre function transform, the spherical harmonic transform
requires time $O(T^{3})$ . Thus, by accelerating associated Legendre function transform, we can reduce thecomputational complexity of spherical harmonic transform up to $O(T^{2}\log T)$ , which is the complexityof FFT. It is well known that the forward transform algorithm is easily derived from the inverse
transform algorithm. So we focus on the inverse associated Legendre function transform (1) in thefollowing discussion.
$\mu k$2.2

Related works

Some related research works are reviewed.
Driscoll and Healy[6] proposed a fast spherical harmonic transform using FFT. Their algorithm runs in time $O(T^{2}\log^{2}T)$ , and is precise (under infinite precision arithmetic operations).
Their algorithm is quite well known, but unfortunately it is numerically unstable. Several modified
versions[7, 14] were proposed to remedy the instability, but they increase computational costs. The
research is still continuing[10, 9] but a tight complexity bound for the stabilized algorithm is notfound yet.
Mohlenkamp[ll] proposed two stable and fast algorithms for the spherical harmonic transformusing fast wavelet transform. His algorithms run in time $O(\tau^{5/2}\log T)$ and $O(T^{2}\log^{2}T)$ , but theformer runs faster in practical sizes.
Rokhlin and Tygert[15] proposed another fast spherical harmonic transform algorithm by using
FMM. The complexity of their algorithm is $O(T^{2}\log T)$ , but their implementation in that paper is
not faster than the direct computations for $T\leq 4096$ , so it requires a better implementation.
The following works are of the Legendre polynomial transform, which is a part of sphericalharmonic transform (only for $m=0$ ). To clarify the difference from the spherical harmonic transform,
we will use $N$ instead of .
Orszag[13] invented an algorithm of Legendre polynomial transform that uses the WKB approximation and FFT. Orszag’s algorithm runs in time ( $N\log^{2}$ N/log log $N$ ). We have improved hisalgorithm using the asymptotic expansion by Stieltjes[16], with the resulting computational complexity of $O(N\log N)[12]$ .
Alpert and Rokhlin[l] proposed a fast algorithm for Legendre polynomial transform and runs intime $O(N\log N)$ . They provide a pair of transforms between the coefficients of a Legendre polynomial
expansion and those of a Chebyshev expansion that represents the same function. Accelerating thosetransforms as the same manner as the FMM, they can show their complexity as $O(N)$ . Then byusing FFT, one can compute the Legendre polynomial transforms (forward and backward) in time
$O(N\log N)$ .
Beylkin, Coifman and Rokhlin showed that the transform between a Legendre polynomial expansion and the corresponding Chebyshev expansion can be computed in time $O(N)$ by using fast
wavelet transform. Thus their algorithm also runs in time $O(N\log N)$ .
The followings are not of spherical harmonic transforms, but our old algorithm owes them the
key idea.
Boyd[3] showed that a function represented by an associated Legendre function expansion can beinterpolated in linear time with a help of FMM.
Jakob-Chien and Alpert[8] proposed a fast algorithm for the spherical filtering, which removes thehigher frequency components than a given cut-off frequency, in time $O(T^{2}\log T)$ , also using FMM.
$T$

$O$

20
3

The new fast transform algorithm

3.1

Fast interpolation of associated Legendre function expansions

Legendre function can be represented with ultraspherical polynomial$P_{n}^{m}(\mu)=c_{n}^{m}P_{m}^{m}(\mu)q_{n-m}^{m}(\mu)$$q_{n-m}^{m}(\mu)$

as
(4)

,

is a constant.
is a polynomial of degree $n-m$ and thus $g^{m}(\mu)/P_{m}^{m}(\mu)$ (whereis defined in (1) is also a polynomial of degree $T-m$ . Now we have Lagrange interpolation
formula for polynomial

where$P_{n}^{m}(\mu)/P_{m}^{m}(\mu)$$c_{n}^{m}$$g^{m}(\mu)$$)$$p( \mu)=\omega(\mu)\sum_{i=1}^{D}\frac{1}{\mu-\mu_{i}}\frac{p(\mu i)}{\omega_{i}(\mu i)}$with $D=\deg(p)+1$ and$\omega(\mu)=\prod_{i=1}^{D}(\mu-\mu_{i})$,$\omega_{i}(\mu)=\frac{\omega(\mu)}{\mu-\mu_{i}}$.

Thus letting $M=T-m+1$ we can evaluate $g^{m}(\mu k)$ on a set of arbitrary points
sampling points
as$\{\mu k\}$from $M$$\{g^{m}(\mu_{i})\}_{i=1}^{M}$$g^{m}( \mu k)=P_{m}^{m}(\mu k)\omega(\mu k)\sum_{i=1}^{M}\frac{1g^{m}(\mu_{i})}{\mu k-\mu_{i}P_{m}^{m}(\mu_{i})\omega_{i}(\mu_{i})}$.

(5)

The computation (5) can be computed in three steps.
1. Point-wise multiplication

2. Summation$\overline{g}_{i}=g^{m}(\mu_{i})/P_{m}^{m}(\mu_{i})\omega_{i}(\mu_{i})$$\overline{g}k=\sum_{i=1}^{At}\overline{g}_{i}/(\mu k-\mu_{t})$3. Point-wise multiplication

for $1\leq i\leq M$ .

for $1\leq k\leq K$ .

$g^{m}(\mu k)=\overline{g}kP_{m}^{m}(\mu k)\omega(\mu k)$

for $1\leq k\leq K$ .

, then steps
If we have fixed the sets
and precomputeandand
1 and 3 are computed in time $O(M)$ and $O(K)$ . By applying FMM[5], we can compute step 2 intime $O(M+K)$ , as is pointed out by Boyd[3]. Thus the interpolation (5) can be computed in time
linear to the input size $(i.e., O(M+K))$ .
$\{\mu_{i}\}$3.2$\{\mu k\}$$P_{m}^{m}(\mu_{i})\omega_{i}(\mu_{i})$$P_{m}^{m}(\mu k)\omega(\mu k)$Base scheme of fast associated Legendre function transform

Using the fast interpolation of associated Legendre function expansion discussed in the previous
section, we can accelerate the associated Legendre function transform in the following way.

1. Compute values on the sampling points$g^{m}( \mu_{i})=\sum_{n=m}^{T}g_{n}^{m}P_{n}^{m}(\mu_{i})$for $1\leq i\leq M$ .

2. Compute values on the target points $g^{m}(\mu k)$ for $1\leq k\leq K$ by interpolation from$g^{m}(\mu_{i})$.

The computational complexity of step 1 is clearly $O(M^{2})$ . By choosing the sampling points
from the target points
and using fast interpolation, the computational complexity of step
2 is $O(K)$ . Remembering $M=T-m+1$ and the alias-free condition (3), one can see that theasymptotic computational costs of the above scheme is 4/9 of the direct computation.
$\{\mu i\}_{i=1}^{M}$$\{\mu k\}_{k=1}^{K}$

3.3

The divide-and-conquer scheme

To reduce the computational costs furthermore, we apply the above scheme to step 1 recursively. Let
us restate step 1:$g^{m}( \mu_{i})=\sum_{n=m}^{T}g_{n}^{m}P_{n}^{m}(\mu_{i})$for $1\leq i\leq M$ .

21
We divide the summation into two roughly equal parts.
,(6)$g^{m}(\mu_{i})$$=$$g_{0}^{m}(\mu_{i})+g_{1}^{m}(\mu_{i})$$g_{0}^{m}(\mu i)$$=$$\sum_{n=m}^{n_{O}-1}g_{n}^{m}P_{n}^{m}(\mu_{i})$,(7)$g_{1}^{m}(\mu_{i})$$=$$\sum_{n=no}^{T}g_{n}^{m}P_{n}^{m}(\mu_{i})$.

(8)is a nearest integer of $m+(T-m)/2$ .
where
Now the fast interpolation scheme is applicable. To accelerate computation of (7), we choose
an appropriate set of $n_{0}-m$ points
from the set of $M$ points, computevaluesfor. and interpolate those values onto the other points.
This is accomplished exactly as stated in the previous section. The computational complexity is$n_{0}$$\mathcal{M}_{0}=\{\mu_{i_{0}}\}$$g_{0}^{m}(\mu_{i_{(}})$$\mathcal{M}=\{\mu_{i}\}$$\mu_{i_{()}}\in \mathcal{M}_{0}$$\mu i\in \mathcal{M}_{0}-\mathcal{M}$$O((n_{0}-m)^{2}+M)$ .

We want to do the same thing for (8). Let us choose $T-n_{0}+1$ points.
from(Note thatandcan be different sets of points.) Because the degree of freedom ofis$T-n_{0}+1$ , the values
on the remaining points
can be determined from thoseon. However,is a polynomial of degree $M=T-m+1$ , not $T-n_{0}+1$ . For thisreason, we cannot use the polynomial-based interpolation (5).
In the old algorithm[17, 25], we solved this problem by introducing split Legendre functions:$\mathcal{M}_{1}=\{\mu_{i_{1}}\}$$\mathcal{M}_{0}$$\mathcal{M}_{1}$$\mathcal{M}$$g_{1}^{m}(\mu)$$g_{1}^{m}(\mu_{i})$$\mu_{i}\in \mathcal{M}-\mathcal{M}_{1}$$g_{1}^{m}(\mu)/P_{m}^{m}(\mu)$$\mathcal{M}_{1}$By splitting the expansion$P_{n}^{m}(\mu)$$=$$P_{n,\nu}^{m,0}(\mu)+P_{n,\nu}^{m,1}(\mu)$,$P_{n,\nu}^{m,l}(\mu)$$=$$q_{n,\nu}^{m,l}(\mu)P_{\nu+l}^{m}(\mu)$$l=0,1$ .

$g_{1}^{m}(\mu)$accordingly$g_{1}^{m}(\mu)$$=$$g_{1,\nu}^{m,0}(\mu)+g_{1,\nu}^{m,1}(\mu)$$g_{1,\nu}^{m,l}(\mu)$$=$$\sum_{n=n_{0}}^{T}g_{n}^{m}P_{n,\nu}^{m,l}(\mu)$$l=0,1$ ,becomes a polynomial of degree $|T-\nu+l-1|$ , and thus interpolation scheme (5)becomes applicable. This method requires two interpolations for $l=0$ and 1, thus the computationalcosts are doubled. Also we have additional computational costs to change the split point , and thenumerical stability is affected 21 .
In the new algorithm[22], we use generalized FMM[29, 28]. Let
be an $M\cross(T-n_{0}+1)$ matrix
whose $(i,j)$ element is, and. By lettingasbe a column vector
a column vector withas its elements, we have matrix-vector representation for (8) as$g_{1,\nu}^{m,l}(\mu)/P_{\nu+l}^{m}(\mu)$

$\nu$

$[$

$]$

$P_{1}$$P_{n_{0}+g}^{m}(\mu_{i})$$(g_{n_{t1}}^{m}, g_{n_{|)}+1}^{m}, \cdots, g_{T}^{m})^{T}$$g_{1}$$v_{1}$$g_{1}^{m}(\mu_{i})$$v_{1}=P_{1}g_{1}$.

Let us consider the set of sampling points
which has $T-n_{0}+1$ points. Collect the set of rows
fromandthat correspond toand represent them as
and . Then we have$\mathcal{M}_{1}$$P_{1}$$\overline{P}_{1}$$\mathcal{M}_{1}$$v_{1}$$\overline{v}_{1}$$\overline{v}_{1}=\overline{P}_{1}g_{1}$.

Also collect the remaining rows and represent them as$\hat{v}_{1}=\hat{P}_{1}g1$

Since each element of$v_{1}$

is either in

an appropriate permutation matrix$\overline{v}_{1}$$\Pi_{1}$

or in
as

$\hat{v}_{1}$$\hat{P}_{1}$and$\hat{v}_{1}$, then we have.

, we can obtain$v_{1}$

by aggregating$v_{1}=\Pi_{1}(\begin{array}{l}\overline{v}_{1}\hat{v}_{1}\end{array})=\Pi_{1}(\begin{array}{l}\overline{P}_{1}\hat{P}_{1}\end{array})g1$.

$\overline{v}_{1}$and$\hat{v}_{1}$by using

22
Note thatis a square matrix, and assume that it is non-singular. Then we have$\overline{P}_{1}$$v_{1}$$=$$\Pi_{1}(\begin{array}{l}IQ_{1}\end{array})\overline{P}_{1}g_{1}$$Q_{1}$$=$$\hat{P}_{1}\overline{P}_{1}^{-1}$

(9)

,

.

represents the interpolations of the elements of
The matrix
from . We apply generalized. Then the computationis accelerated.
FMM (details described in [28]) to
For the problems to which the analytical FMM is applicable, generalized FMM works no worse
represents a polynomial interpolation, then the complexity
than the analytical FMM. Thus if
is linear to the sum (instead of the product) of the numbers of the rowsof computing. But
is surely not a polynomial interpolation. So we lose the complexity
and the columns ofbound. However, we can avoid the split Legendre functions that doubled the computational costs inthe old algorithm. Thus we may have some gain. Let us evaluate it experimentally in section 4.
Before closing this subsection, we evaluate the computational complexity on the assumption thatgeneralized $FMM$ runs in linear time. We will discuss about the preprocessing phase somewhat later,and here let us consider only transform computations assuming that appropriate preprocessing has
been done. Let us compute the computational costs in a backward manner. The last step of thecomputation is the polynomial interpolation (5). As is stated in the previous section, its complexity
is $O(K)$ . The second last step is the summation (6). The complexity is clearly $O(M)$ . The third last
step is the interpolationand the corresponding interpolation of the other half. From theassumption, their computational costs are summed up as $O(M)$ . The step before that is summations
similar to (6). They are two summations of half sizes, thus their costs are summed up as $O(M)$ .
Then we have four interpolations of half sizes, which also amounts $O(M)$ . After $O(\log M)$ steps ofrecursion, we have $O(M)$ problems of a constant size. Summing up them, we have the computational
complexity $O(K+M\log M)$ . This complexity is for each $m$ . Note that $M=T-m+1$ and $K=O(T)$ .
Summing up for all $m$ , the computational complexity is $O(T^{2}\log T)$ .
$Q_{1}$$\hat{v}_{1}$$Q_{1}$$\overline{v}_{1}$$Q_{1}\overline{v}_{1}$$Q_{1}$$\hat{v}_{1}=Q_{1}\overline{v}_{1}$$Q_{1}$$Q_{1}$$\hat{v}_{1}=Q_{1}\overline{v}_{1}$3.4

Error control and sampling point selection

As stated in the previous subsection, the interpolation matrixin (9) is accelerated by using
generalized FMM. Because the FMM is an approximate algorithm, we have to consider the errors at
the output ( in (9)) incurred by the approximation. The following methodology of error analysis
and control is introduced in [26, 27] and used also for generalized FMM[28].
Before discussion, note that (9) does not reach the end of the computation. The resultingareinterpolated into the remaining points by the interpolation (5). Let us represent the interpolation(5) by a matrix
and consider the computation$Q_{1}$$v_{1}$$v_{1}$$Q$)$w_{1}=Q\Pi_{1}(\begin{array}{l}IQ_{1}\end{array})\overline{P}_{1}g_{1}$.

be the matrix representation of the approximate interpolation accelerated by generalized
FMM. Then the actual computation isLet$\tilde{Q}_{1}$$\overline{w}_{1}=Q\Pi_{1}(\begin{array}{l}I\tilde{Q}_{l}\end{array})\overline{P}_{1}g1\cdot$

Thus the error is$\tilde{w}_{1}-w_{1}=Q\Pi_{1}(\begin{array}{l}0\tilde{Q}_{1}-Q_{1}\end{array})\overline{P}_{1}g_{1}$Let us ignore

$0$

above

$\tilde{Q}_{1}-Q_{1}$, and we want to control$e_{1}=Q\hat{\Pi}_{1}(\tilde{Q}_{1}-Q_{1})\overline{P}_{1}g_{1}$where$\hat{\Pi}_{1}$represents the part of$\Pi_{1}$that corresponds to$Q_{1}$

.

,

.

23
Let

$\Vert\cdot\Vert$

represent some norm both for vectors and matrices that satisfies the consistency inequality:$\Vert AB\Vert\leq\Vert A\Vert\Vert B\Vert$

,

(10)

$\Vert Av\Vert\leq\Vert A\Vert\Vert v\Vert$

for any matrices
and $B$ and any vector (but of course the number of columns of
number of rows of $B$ and should be the same).
Then we have
$A$

$v$

$A$

and the$v$$\Vert e_{1}\Vert\leq\Vert Q\hat{\Pi}_{1}\Vert\Vert\tilde{Q}_{1}-Q_{1}\Vert\Vert\overline{P}_{1}\Vert\Vert g_{1}\Vert$

.

(11)

Note that if we use double precision, then we cannot expect approximation of higher precision than
the round-off error:$\epsilon\approx 10^{-16}$.

$\frac{\Vert\tilde{Q}_{1}-Q_{1}\Vert}{||Q_{1}\Vert}\geq\epsilon$

Thus the right-hand side of (11) has a lower bound,$|IQ\hat{\Pi}_{1}I1\Vert\tilde{Q}_{1}-Q_{1}\Vert\Vert\overline{P}_{1}\Vert\geq\Vert Q\hat{\Pi}_{1}\Vert\Vert Q_{1}\Vert\Vert\overline{P}_{1}\Vert\epsilon$

and soworks as a kind of condition number.
However, in many cases the above condition number can be very large even when the output
error
oris quite reasonable. This happens when the consistency inequality (10) isvery loose. Then the evaluation (11) does not give a good estimate.
A simple solution is proposed in [26, 27]. Let us introduce diagonal matrices
andandinsert them as:$\kappa=||Q\hat{\Pi}_{1}\Vert\Vert Q_{1}\Vert\Vert\overline{P}_{1}\Vert$$\kappa$$\Vert e_{1}\Vert$$\Vert\tilde{v}_{1}-v_{1}\Vert$$S_{B}$$S_{A}$$e_{1}=(Q\hat{\Pi}_{1}S_{A}^{-1})(S_{A}(\tilde{Q}_{1}-Q_{1})S_{B})(S_{B}^{-1}\overline{P}_{1})g1$

so that we have$\Vert e_{1}.\Vert\leq\Vert Q\hat{\Pi}_{1}S_{A}^{-1}\Vert\Vert S_{A}(\tilde{Q}_{1}-Q_{1})S_{B}\Vert\Vert S_{B}^{-1}\overline{P}_{1}\Vert\Vert g_{1}\Vert$

By letting $Z_{1}=S_{A}Q_{1}S_{B}$ and.

we have$\tilde{Z}_{1}=S_{A}\tilde{Q}_{1}S_{B}$$\Vert e_{1}\Vert\leq\Vert Q\hat{\Pi}_{1}S_{A}^{-1}\Vert\Vert\overline{Z}_{1}-Z_{1}\Vert\Vert S_{B}^{-1}\overline{P}_{1}\Vert\Vert g1\Vert$.

Now the “condition number” becomes$\kappa_{S}=\Vert Q\hat{\Pi}_{1}S_{A}^{-1}\Vert\Vert Z_{1}\Vert\Vert S_{B}^{-1}\overline{P}_{1}\Vert$,

which can be controlled by appropriately choosing. Our previous works[26, 27, 28]
andproposed to choose
andso that the columns ofand the rows ofhave theunit norm. This method is quite well known as a simple device to improve the condition numbers ofmatricesand. The effects of anddepend on the matrices, and we have no theoreticalguarantee of their goodness. We only know it works very well in the experiments.
Now if we want to limit the “relative error“$S_{A}$$S_{A}$$Q\Pi_{1}$$\overline{P}_{1}$$S_{B}$$S_{B}$$Q\hat{\Pi}_{1}S_{A}^{-1}$$S_{A}$$S_{B}^{-1}\overline{P}_{1}$$S_{B}$,$\frac{||e_{1}||}{||g_{1}||}\leq\delta$then it is enough to approximate $Z_{1}=S_{A}Q_{1}S_{B}$ as$\frac{\Vert\tilde{Z}_{1}-Z_{1}\Vert}{\Vert Z_{1}\Vert}\leq\frac{\delta}{\kappa_{S}}$.

Practically it is $observed\sim$ that this requirement is still too pessimistic. So in our implementation therequired precision foris mitigated as follows.
Before discussing it, note that we can bound some norms by using the proposed.
and
Let us use 2-norm which is bounded by the Frobenius norm as. Remember that the
Frobenius norm of a matrix is the positive square root of the sum of the 2-norms of the rows or thecolumns of the matrix. Reviewing the previous subsection, we can count the number of rows andcolumns of the matrices. We have$Z_{1}$$S_{A}$$\Vert A\Vert_{2}\leq\Vert A\Vert_{F}$$\Vert Q\hat{\Pi}_{1}S_{A}^{-1}\Vert\leq\sqrt{n_{0}-m}$,$||S_{B}^{-1}\overline{P}_{1}\Vert\leq\sqrt{T-n_{0}+1}$.

$S_{B}$

24
Thus $\kappa_{S}\leq\sqrt{(n_{0}-m)(T-n_{0}+1)}\Vert Z_{1}\Vert$ .
. We start computation (ofwith a constantIn our implementationis replaced byaccordingly.
preprocessing) with $\kappa_{0}=1$ , and if the resulting error is too large, then we reduce
After a few itcrations, the resulting error meets the requirement . If the required precision is too
small, then
may be smaller than the round-off error level. Then our program terminateswith a message that the precision cannot be achieved.
. Smaller
Note that the attainable precision is determined mostly by
andcomes ffom one common
is preferable. Remember that $Z_{1}=S_{A}Q_{1}S_{B}$ and, andcorresponds tocorresponds to the sampling points of the interpolation andmatrixpoints. Since
sampling
depends
points.
wethaton
selection
of
the
Thus
know
the
the otherfinite,pointssamplingby trying all
wecanminimize
of
the
set
is
of
choicesofthethe number
What wedirectly.
possibilities, but it is computationally impractical. So we do not minimizeli:
do is just an LQ (transposed QR) decomposition with pivoting for$\kappa_{0}\Vert Z_{1}\Vert$$\kappa_{S}$$\kappa_{0}$$\kappa_{0}$$\delta$$\delta$$\delta/(\kappa_{0}\Vert Z_{1}\Vert)$$\delta$$\Vert Z_{1}\Vert=\Vert S_{A}Q_{1}S_{B}\Vert$$\delta$$Q_{1}=\hat{P}_{1}\overline{P}_{1}^{-1}.\overline{P}_{1}$$P_{1}$$\Vert Z_{1}\Vert$$\hat{P}_{1}$$\hat{P}_{1}$$\overline{P}_{1}$$\Vert Z_{1}\Vert$$\Vert Z_{1}\Vert$$\Vert Z_{1}$$S_{A}P_{1}$$S_{A}P_{1}=(\begin{array}{l}\overline{L}_{1}\hat{L}_{1}\end{array})U_{1}=(\begin{array}{l}I\hat{L}_{1}\overline{L}_{1}^{-1}\end{array})\overline{L}_{1}U_{1}$.

(Note that we use $U$ instead of the usual $Q.$ ) The pivoted LQ decomposition approximately minimizes, rather, but it does not completely equal to$Z_{1}$$\hat{L}_{1}\overline{L}_{1}^{-1}$$Z_{1}=\hat{L}_{1}\overline{L}_{1}^{-1}S_{A}S_{B}$.

We don’t consider $S_{A}S_{B}$ , so it’s suboptimal. Also we ignore the effects of the recursion. But it works
well, as is reported in the next section.

Last, let us discuss about the numerical stability. The Driscoll-Healy algorithm [6], which isnumerically unstable, also uses exactly the same divide-and-conquer scheme as our algorithm. However, our algorithm is numerically stable, as the next section reports. The numerical instability of
the Driscoll-Healy algorithm comes from the use of Fourier-Chebyshev expansion using ultrasphericalpolynomial (4). It is successful for low values for $m$ , but for higher $m$ , we meet difficulties at representing the associated Legendre functions using trigonometric functions. Our algorithm uses thesame relation (4) to derive the fast interpolation (5), but it has enough degrees of freedom to attainnumerical stability, that is, choice of the sampling points. The use of Fourier-Chebyshev expansion
in the Driscoll-Healy algorithm effectively limits the sampling points to be the Chebyshev knots,large for large $m$ . In our new algorithm it is effectively solved by the pivoted
and that makes
LQ decomposition. In the old algorithm we have another method of sampling point selection for
numerical stability 18 .
As is discussed in [21], the split Legendre functions in our old algorithm considerably increase
numerical instability. Our new algorithm removes the split Legendre functions, so it may improve
the numerical stability compared to the old algorithm. The experimental results in the next section
confirm this.
$\overline{L}_{1}^{-1}$

$[$

3.5

$]$

Preprocessing algorithm

Our algorithm uses divide-and-conquer strategy. It recursively divide the problem into two subproblems. At some level of recursion, the problem becomes so small that the direct computation
is faster. At that point the recursion should be stopped. To minimize the computational cost, ourimplementation uses the branch-and-bound technique in the preprocessing phase. Figure 1 shows apseudo-code of our preprocessing algorithm. It is almost same as that given in [27].
The function fast-pp $(m, n_{0}, n_{1}, \mathcal{M}, c)$ computes the preprocessing of fast transform on the setof the evaluation points
for specified $m$ and $n_{0}\leq n<n_{1}$ . Parameter is the limit of thecomputational costs for the prescribed transform, and the return value of fast-pp is the obtained
computational costs. Function dir-pp takes similar parameters and returns the computational costsof the direct computation.
This algorithm compares three methods: (1) compute the prescribed transform by the directcomputation without interpolation, (2) compute the values on the set of sampling points (hereafter$\mathcal{M}$

$c$

25

1

function fast-pp $(m, n_{0}, n_{1}, M. c)$2$c_{D}=dir_{-}pp(m, n_{0}, n_{1}, \mathcal{M})$

3
4
5
6

7
8
9

10
11
12
13

14
15
16
17

18
19
20
21
22
23
24

new cost limit
compute interpolation matrix and its FMM approximation
letbe the set of the sampling points
letbe the computational costs of the fast interpolation
if $(c_{i}>c_{n})$prepare direct computation without interpolation (method 1)
returnelse//try use of interpolation$c_{n}= \min\{c,$ $c_{D}\}//$$\mathcal{X}\subset \mathcal{M}$$c_{i}$$c_{D}$$c_{d}=$

dir $-pp(m, n_{0}, n_{1}, \mathcal{X})$$c_{t}= \min\{c_{n}-c_{i},$ $c_{d}\}//$

new cost limit$\nu=(n_{0}+n_{1})/2//$ divide the problem into two halves$c_{0}=$

fast $-pp(m, n_{0}, \nu, \mathcal{X}, c_{t})$fast-pp $(m. \nu, n_{1}, \mathcal{X}, c_{t}-c_{0})$if $(c_{0}+c_{1}<c_{d}$ and $c_{i}+c_{0}+c_{1}<c_{n})$prepare interpolation with divide-and-conquer (method 3)
return $c_{1}+c_{0}+c_{1}$else if $(c_{i}+c_{d}<c_{n})$prepare interpolation with direct computation (method 2)
return $c_{i}+c_{d}$elseprepare direct computation without interpolation (method 1)
return$c_{1}=$$c_{D}$

Figure 1: Pseudo-code of preprocessing algorithm
by the direct computation and then interpolate them on $M$ , and (3) compute the values onusing divide-and-conquer and then interpolate them on $M$ .
At lines 4 to 6 it generates the interpolation matrix and its FMM approximation, and evaluates
its computational costs . If the costs of interpolation $q$ is larger than (given limit as parameter) or(costs of direct computation), then we must abandon methods (2) and (3) that use interpolation.
This check is done at lines 7 to 9. Otherwise, we will try methods (2) and (3). At line 11 the costsof the direct method for computing values on
is evaluated as . Lines 14 and 15 call fast-pp
recursively, with appropriate cost limits. The costs of the method (2) are now given as $c_{i}+c_{d}$ , andthe costs of the method (3) are $c_{i}+c_{0}+c_{1}$ . Comparing them with the costs of the direct method
without interpolation, the method of the lowest costs is chosen at lines 16 to 24. The function
does not abort even if the obtained computational costs are over the given limit , because the resultmust be discarded in the caller.
Let us consider the computational complexity of the preprocessing algorithm given in Figure 1.
Let the number of points $M=|\mathcal{M}|$ and $N=n_{1}-n_{0}+1$ . Note that usually we have $M\geq N$ . Using
Suda-Kuriyama algorithm[28] for generalized FMM, the preprocessing costs for the FMM in line 4
is $O(MN)$ . The computation of the interpolation matrix in line 4 uses pivoted LQ decomposition,
which requires computations of $O(MN^{2})$ , which is the major costs at this recursion level. So thecomputational costs are 1/8 of this level for each of the halves at the next recursion level. Therefore
the computational costs of the lower recursion levels are less than half of those at this level. Thus
the computational costs at the root of the recursion tree determine the complexity order, which is$O(K(T-m+1)^{2})$ . Summing up for every $m$ , and assuming that $K=O(T)$ , the total computational costs for the preprocessing is computed as $O(T^{4})$ . Those heavy costs solely come from the
LQ decomposition to obtain the interpolation matrix. If we had any method to compute theinterpolation matrix faster, then the computational complexity of the preprocessing phase would be
reduced.
$\mathcal{X})$$\mathcal{X}$$c$$c_{t}$$c_{D}$$\mathcal{X}$$c_{d}$$c_{D}$$c$$Q_{1}$

26
Evaluation of the new algorithm

4This section reports the performance of the new algorithm. The following data are taken from [23].

4.1

Evaluation in spherical harmonic transform

Table 1: Comparison of speedup rates of spherical harmonic transform by old and new algorithms

Table 1 compares the speedup rates of our old and new algorithms. The required precision is set
to $10^{-10}$ . Here speedup is computed with the number of floating point number operations relative tothat of the direct computation, and it is not of the CPU times.
Table 1 reports considerable speedup even for $T=255$ . This comes mostly from the dropping[27],and the effects of the FMM are small in this size. For larger sizes, the new algorithm performs slightly
better than the old one.

4.2

Evaluation in Legendre polynomial transform

Next the new algorithm is examined further by limiting to $m=0$ , that is, to the Legendre polynomialtransform. The computational costs are the higher for the smaller $m$ , and also the FMM has themore effects on the performance. Also the recursion is the deepest for $m=0$ . The number of points$M$ is now set as $M=T+1$ .

Figure 2: Relative computational costs of Legendre polynomial transforms by the old algorithm

Figures 2 and 3 show the floating point operation counts relative to the direct computation. The
x-axis shows the attained precision. $L12701d$ and $Ll27new$ ” give the performances for $T=127$ ,”etc.

Comparing those figures, we can see that $T=127$ gives small difference between the old and the
new algorithms, but for $T=2047$ the old algorithm takes about 1.5 times $mo$ re computations than
the new algorithm.

27

Figure 3: Relative computational costs of Legendre polynomial transforms by the new algorithm

For the old algorithm, the gradient is different for $T=1023$ with lower precision and for $T=2047$
than the other cases. This is because the split Legendre functions are used only for those conditions.
However for the new algorithm, we found no such phase change.
The plots for the old algorithm stop at precision of $10^{-15}$ to $10^{-11}$ . This is because higher precision
was not attained because of the numerical stability. For the new algorithm higher precisions wereattained. Those results imply that the numerical stability is much improved in the new algorithm.
Further discussion on the performance will be found in [23].

5Summary and future works

This paper discusses our fast algorithm for as sociated Legendre function transforms with generalized
FMM. This algorithm is formerly reported in our previous reports [22, 23], but the details of thealgorithm were not presented. In this paper, the transform algorithm and the preprocessing algorithm
are discussed. The performance is reported in English for the first time.
There are several remaining works to do.
The new algorithm performs better than the old algorithm, but its computational complexity is
not bounded theoretically. To bound it, the complexity bound for the interpolation with generalized
FMM must be obtained, and to that purpose, we must have some analytical expression for theinterpolation matrix entries.
Any specific property of the associated Legendre functions is utilized in the new algorithm.
The new algorithm just compresses the transform matrix. The only exception is that the forward
transform can be obtained by transposition of the inverse transform matrix, but this assumption canbe removed if we compute the forward transform matrix directly. Thus the new algorithm can beapplied to any transform at least conceptually.
Last, the implementation should be improved and released to the public use. Before doing that,
we should improve the implementation of generalized FMM, because its current performance issomewhat lower than the previous one (but with dramatically reduced preprocessing time).

Acknowledgments
This research is partly supported by Grant-in-Aid for Scientific Research on Priority Areas, MEXT
and CREST, JST.

28
References[1] B. K. Alpert and V. Rokhlin: A Fast Algorithm for the Evaluation of Legendre Expansions,J. Sci. Stat. Comput., Vol. 12, No. 1, pp. 158-179, 1991.
$SI\mathcal{A}M$[2] G. Beylkin, R. R. Coifman, and V. Rokhlin: Fast Wavelet Transforms and Numerical Algorithms
I, Comm. Pure Appl. Math., No. 44, pp. 141-183, 1991.
[3] J. P. Bovd: Multipole expansions and pseudospectral cardinal functions: A new generation of
the fast Fourier transform, J. Comput. Phys., Vol. 103, pp. 184-186, 1992.

[4] A. Dutt, M. Gu, and V. Rokhlin: Fast algorithms for polynomial interpolation, integration, anddifferentiation, SIAM J. Numer. Anal., Vol. 33, No. 5, pp. 1689-1711, 1996.
[5] L. Greengard and V. Rokhlin: A Fast algorithm for particle simulations, J. Comput. Phys.,Vol. 73, pp. 325-, 1987.

[6] J. R. Driscoll and D. Healy, Asymptotically Fast Algorithms for Spherical and Related Transforms, Proc. 30th IEEE FOCS, pp. 344-349 (1989).
[7] D. M. Healy Jr., D. Rockmore, P. J. Kostelec, and S. S. B. Moore, FFTs for 2-Sphere
Improvements and Variations, Tech. Rep. PCS-TR96-292, Dartmouth Univ., 1996.

–[8] R. Jakob-Chien and B. K. Alpert, A Fast Spherical Filter with Uniform Resolution, J. Comput.
Phys., Vol. 136, pp. 580-584, 1997[9] J. Keiner and D. Potts, Fast Evaluation of Quadrature Formulae on the Sphere, Math. Comp.,
Vol. 77, No. 261, pp. 397-419 (2008).
[10] S. Kunis and D. Potts, Fast spherical Fourier algorithms, J. Comp. Appl. Math., Vol. 161,
pp. 75-98 (2003).

[11] M. J. Mohlenkamp, A Fast Transform for Spherical Harmonics, J. Fourier Anal. Appl., Vol. 2,
pp. 159-184, 1999.
[12] A. Mori, R. Suda and M. Sugihara, An improvement on Orszag’s Fast Algorithm for Legendre
Polynomial hansform, Trans. Inform. Process. Soc. Japan, Vol. 40, No. 9, pp. 3612-3615 (1999).
[13] S. A. Orszag, Fast Eigenfunction T}ansforms, Science and Computers, Adv. Math. Supplementary Studies, Vol. 10, pp. 23-30, Academic Press, 1986.
[14] D. Potts, G. Steidl, and M. Tasche, Fast and stable algorithms for discrete spherical Fouriertransforms, Linear Algebra Appl., pp. 433 –450, 1998.
[15] V. Rokhlin and M. Tygert: Fast algorithms for spherical harmonic expansions, SIAM J. Sci.
Comp. Vol. 27, No. 6, pp. 1903-1928 (2006).
[16] G. Szego, Orthogonal Polynomials, p. 193. AMS, 1959.

[17] R. Suda, A Fast Spherical Harmonics Transform Algorithm, Infor. Proc. Soc. Japan SIG Notes,
No. 98-HPC-73, pp. 37-42, 1998.

[18] R. Suda: Stability Control of Fast Spherical Harmonics Transform, Info. Proc. Soc. Japan SIG
Notes, No. 2001-HPC-88, pp. 43-48 (2001).
[19] R. Suda: Fast spherical harmonics transform of FLTSS and its evaluation, The 2002 Workshop
on the Solution of Panial Differential Equations on the Sphere, Fields Inst., U. Toronto (2002).

[20] R. Suda: Fast spherical harmonic transform routine FLTSS applied to the shallow water test
set, Mon. Wea. Rev., Vol. 133, No. 3, pp. 634-648 (2005).

[21] R. Suda: Stability analysis of the fast Legendre transform algorithm based on the fast multipole
method, Proc. Bstonian Acad. Sci. Phys. Math., Vol. 53, No. 2, pp. 107-115 (2004).

29[22] R. Suda: Fast Spherical Harmonic Transform with generalized Fast Multipole Method, 2005
International Conference on Scientific Computing and Differential Equations (SciCADE05),
p. 86 (2005)[23] R. Suda: High Performance Implementation and Evaluation of the Spherical Harmonic Transforms in the Fast Orthogonal Function Transform Routine FXTPACK, Info. Proc. Soc. Japan
SIG Technical Report, Vol. 2005, No. 97 (2005-HPC-I04), pp. 31-36 (2005).
[24] R. Suda: FLTSS home page:

http: $//www$ . na. cse. nagoya-u. ac. $jp/\sim_{reiji}/fltss/$[25] R. Suda and M. Takami: A Fast Spherical Harmonics Transform Algorithm, Math. Comp.,
Vol. 71, No. 238, pp. 703-715 (2002).

[26] R. Suda and M. Takami: Error Analysis and Control of Fast Spherical Harmonics Transform,
Info. Proc. Soc. Japan SIG Notes, No. 2000-HPC-84, pp. 7-12 (2000).
[27] R. Suda and M. Takami: Error Analysis and Control of Fast Spherical Harmonics Transform,
Tmns. Info. Proc. Soc. Japan HPS, Vol. 42, No. SIG12 (HPS4), pp.49-59 (2001).

[28] R. Suda and S. Kuriyama: Another Preprocessing Algorithm for Generalized One-Dimensional
Fast Multipole Method, J. Comp. Phys., Vol.195, pp. 790-803 (2004).

[29] N. Yarvin and V. Rokhlin: A Generalized One-Dimensional Fast Multipole Method with Application to Filtering of Spherical Harmonics, J. Comp. Phys., Vol. 147, pp. 594-609 (1998).