球面調和関数は球面上の直交関数系で、球面上や球座標を用いた場合の関数の表現・近似・解析に向いた関数です。これを 1 次元下げると単位円上の三角関数系となりますので、球面調和関数変換はフーリエ変換の 2 次元版ということができます。
球面調和関数変換の最大のユーザは気象・気候シミュレーションです。日々行われている天気予報のうち、地球全体を対象とする部分は球面調和関数変換によっています。また、地磁気の解析や、画像処理などにも用いられています。そのほか、球座標を用いるさまざまな物理シミュレーションや数値計算にも用いられます。
球面調和関数変換は、緯度経度格子点上の関数値と、球面調和関数展開の係数との間の相互変換の計算です。
この変換は、東西方向のフーリエ変換と、南北方向のルジャンドル陪関数変換とに分解することができます。
緯度経度格子が N × 2N の場合、フーリエ変換には FFT (Fast Fourier Transform) が使えて計算量が O(N2 log N) になるのですが、ルジャンドル陪関数変換は通常定義式を直接計算するため計算量が O(N3) になります。
この O(N3) という計算量は通常のシミュレーションで最も大きな計算量となり、精度を上げようとして大規模シミュレーションにすると、全体の計算時間の大部分を占めるようになります。これは昔から知られていて、ルジャンドル陪関数変換の高速化は長年の懸案でした。
我々は、高速多重極子展開法 (FMM: Fast Multipole Method) を用いた高速変換アルゴリズムを提案しています。計算量は O(N2 log N) です。
我々のアルゴリズム以外にもいくつか高速ルジャンドル陪関数変換アルゴリズムの提案があります。もっとも優れたものは Mohlenkamp が提案した FWT (高速ウェーブレット変換)を利用したアルゴリズムで、計算量は O(N2 log2 N) ですが、実際の速度は我々のアルゴリズムとほぼ同じです。
このほかに FFT に基づく D-H 法、FCT (高速コサイン変換) を用いた Potts らの方法、FMM と FFT を組み合わせた Alpert-Rokhlin の方法、FFT と直接計算を用いた Orszag の法を森らが改良したものなどがありますが、これらはいずれも数値的に不安定であるか、球面調和関数変換まではできないかどちらかになっています。これらはいずれも FFT を使っており、どうやらルジャンドル陪関数変換を FFT に帰着させるのには無理があるようです。
我々のアルゴリズムは高速多重極子展開法 (FMM) にその高速性を負っています。FMM が高性能であればあるほど、単純に我々の高速球面調和関数変換法も速くなります。
そこで、高性能の FMM について研究を進めてきました。最初は Taylor 展開に基づいていましたが、その後 Chebyshev 補間に基づく方法、さらには低階数近似に基づく方法に改良されました。
低階数近似に基づく方法は、特定の関数系を対象にしていた従来の FMM とは異なり、数値的に行列を圧縮します。このため、高精度の近似展開が未知の関数系に対しても効率的な計算を達成することができます。これは「一般化 FMM」と呼ばれています。
実際に、この一般化 FMM を用いて、従来は分離ルジャンドル関数を用いて定義されていた我々のアルゴリズムを、分離を用いずに直接ルジャンドル陪関数変換を近似することに成功しました。
我々のアルゴリズムは FMM とルジャンドル陪関数変換で二重の分割統治法になっており、数値アルゴリズムとしてはかなり複雑なものといえます。しかし上述のとおり従来手法のかなりの部分が数値的な不安定性を問題として持っていることから、誤差と安定性の解析・最適化がどうしても必要でした。
これはかなり難しい問題でしたが、我々は行列積のノルムを対角スケーリングした一貫性不等式で評価する手法、および数値的な安定性の指標としてαβ積という概念を提案し、誤差と安定性の解析および最適化を実現しました。
これらの手法は、ルジャンドル陪関数変換、一般化高速多重極子展開法、さらには低階数近似にも一貫して用いられて、全体としての誤差制御と安定性の最適化を達成しています。
我々の高速変換アルゴリズムは汎用性が高く、直交関数変換であれば原理的に適用可能です(数値的安定性はその都度検討する必要があります)。
また、一般化 FMM および低階数近似も幅広い応用が期待されます。