１.　学会の概要

East Asia SIAMはSIAMの１セクションであり，東アジア地域における応用数理研究の振興を目的としている。これまで，東アジアの各国で毎年１回会議を開催しており，今年は札幌，アモイ（中国），デジョン（韓国）に続いて５回目となる。今回は約80名が参加した。発表は数理モデリングに関する話が多く，他に偏微分方程式や確率微分方程式の数値解法など，解析系の数値計算が中心であった。金融工学における数値計算に関連する発表も，私の発表を含めて5件あった。

２.　興味深かった講演

興味深かった発表について，以下にまとめる。

(1) Laplace Transformation Method for the Black-Scholes Equation (Student paper contest, 1st prize)
　　Hyoseop Lee (Seoul Nationa University)

ブラック-ショールズ方程式の数値解法として，時間方向にラプラス変換を行って解く方法を提案した。本解法は，利子率rやボラティリティσが資産価格に依存する場合にも適用可能である。この解法について，r，σに関するある条件の下で，ラプラス変換後の（資産価格を独立変数とする）偏微分方程式の解の存在と一意性を証明した。また，数値計算上はラプラス逆変換の計算が問題となるが，そのための複素積分における積分経路の最適な決め方を提案した。この経路を使うと，ラプラス変換におけるzの点数に対して，誤差が指数関数的に減少することが示せる。数値実験の結果，この傾向が実際に確認され，多くの場合にzは15点程度で十分な精度が達成できることがわかった。また，ラプラス変換後の各zに対する偏微分方程式は独立に解けるため，並列性が高いのもこの方法の特長である。数値実験では，15プロセッサ時に13.5倍程度の加速を達成していた。

(2) Sinc-Galerkin Method with DDM for the Option Pricing under Jump-Diffusion Model
　　Jun Liu and Hai-Wei Sun (South China normal University)

資産価格がジャンプ拡散モデルに従う場合，オプション価格V(x,t)の従う方程式はPIDE（Partial Integro-Differential Equations）となる。本発表では，代表的なジャンプ拡散モデルであるマートンモデルについて，効率的なアルゴリズムを提案した。PIDEの時間方向の離散化に当たっては，微分作用素の部分を後退差分，積分作用素の部分を前進差分で扱うIMEX型公式を用いる。ただし，これだけでは１次の精度であるため，リチャードソン加速を併用することにより高次化を行う。一方，空間方向については，Sinc関数を基底関数としたGalerkin法を用いる。これにより，格子点数をNとするとき，O(-c√N)という指数関数的な誤差の減少が達成できる。また，係数行列はToeplitz行列となる。積分作用素の計算で現れる数値積分も，Sinc格子点での値を使って計算する。ここでも誤差はO(-c√N)となる。この数値積分は，関数とガウス分布の畳み込み積となるため，改良型の高速ガウス変換を用いて，計算量をO(N)に削減する。

ヨーロピアン型のオプションでは，満期におけるペイオフが１階導関数に不連続を持つため，これが原因で本来のSinc-Galerkin法の収束次数が低下する。そこで，対数価格xの空間において，対数行使価格kの上下で領域を分割し，各領域で別々に計算を行う。境界では，オプション価格のV(x,t)の０階，１階の微分が一致することを要請する。これにより，係数行列はToeplitz行列を対角線に沿って並べたブロック対角に近い形になる。後退差分のためにこの行列を解く部分では，Toeplitz行列のための効率的なソルバを修正して利用できる。数値実験の結果，提案手法は理論通りに指数的な誤差の減少を示し，広く使われるクランク-ニコルソン法に比べて高精度であることがわかった。

感想： 高速ガウス変換を使い，ジャンプ拡散モデルの下でのオプション価格計算を高速化するという方法は，私も従来から色々と試みてきたので，本発表は興味深かった。本発表では改良型の高速ガウス変換を用いているが，ここで扱っている1資産の問題の場合には，Greengard & Strain が提案した従来型の高速ガウス変換のほうが効率が良いのではないかと思う。ただし，発表者らは，マートンモデル以外にKouのモデルへの適用も視野に入れており，そのための拡張に関しては，エルミート関数の性質を多用する従来型の高速ガウス変換より，テイラー展開に基づく改良型高速ガウス変換の方が素直に行くということであった。

(3) Shift-Invert Arnoldi Approximation to the Toeplitz Matrix Exponential
　　Spike T. Lee (University of Macau)

資産価格がジャンプ拡散モデルに従う場合のオプション価格の満たすPIDEは，空間方向をSinc-Galerkin法により離散化すると，Toeplitz行列を係数行列とする大規模連立常微分方程式となる。この方程式を解く方法として，行列の指数関数exp(-At)を用いる方法がある。そこで，AがToeplitz行列の場合に，exp(-At)y（yはベクトル）を効率的に計算する方法を提案した。

exp(-At)yの数値計算法としては，アーノルディ法を用いてクリロフ部分空間に射影する方法がよく知られている。しかし，この方法では，tが大きくなるにつれて，クリロフ部分空間の次元を増やさなければならないという問題点がある。一方，行列の有理関数まで用いる方法[1]では，項の数がより少なくて済むことが知られている。特に，有理式の分母の次数が１の場合，この方法はAの代わりに(A+σI)^-1に対してアーノルディ方法を適用するshift-and-invert型のアーノルディ法と見なすことができる。この方法では，アーノルディ法の各ステップで(A+σI)^-1とベクトルの積が必要となり，これが計算量の大部分を占める。そこで，Toeplitz行列Aの逆行列を陽的に求める Gohberg-Semencul 公式[2]を利用した。この公式によれば，逆行列A^-1の第1列，第n列のみから作られる4つのToeplitz行列X，Y，X'，Y'を用いて，A^-1はA^-1=c(XY^T-Y'X'^T)と表せる。そこで，アーノルディ法の反復前に高速Toeplitzソルバを用いてA^-1の第1列，第n列を求め，反復中ではGS公式を用いて，しかもX，Y，X'，Y'に関する行列ベクトル積をFFTを用いて計算することにより，1反復当たりの計算量をO(n log n)に抑えた。

さらに，反復回数についても理論的評価を行った。Toeplitz行列Aの第(i,i+k)成分をa_kとするとき，a_kをフーリエ級数展開の第k成分とするような関数f(x)をAの母関数と呼ぶ。本研究では，f(x)が Re(f)≧0 かつ ‖Im(f)/Re(f)‖_∞=O(1) という2つの条件を満たすならば，必要なアーノルディ法の反復回数がtによらないことを証明した。この性質は，数値実験でも確認できた。

参考文献
[1] I. Moret and P. Novati: "RD-rational approximates of the matrix exponential", BIT, Vol. 44, pp. 595-615 (2004).

[2] I. Gohberg and A.A. Semencul: "On inversion of finite-section Toeplitz matrices and their continuous analogues" Mat. Issled. Kishinev, Vol. 7, No. 2, 201-224 (1972) (In Russian).

感想： 大変高度な内容で，後で発表者が修士課程の学生だと聞いて驚いた。英語も完璧だった。アーノルディ法によるexp(-At)yの計算は，私も最近別の分野の問題で使っているが，本問題と同様に，tが大きくなると，必要なクリロフ部分空間の次元が急速に大きくなるという問題がある。本研究の方法は，直接にはToeplitz行列の場合にのみ適用可能であるが，部分空間の次元がtによらないという性質は非常に興味深い。上記の参考文献などを調べてみたい。

(4) Error-Free Transformations and Its Applications
　　Shin'ichi Oishi (Waseda University)

浮動小数点演算でx=a+bを計算すると，丸め誤差のためxはa+bと正確には一致しない。そこで，x+y=a+bが厳密に成り立つようなyを同時に求める方法がある。必要な計算量は6flopsである。これは(a,b)から(x,y)へのerror-free transformation（EFT）と呼ばれ，Knuthにより1969年に提案された。一方，積についても，a*b=x+y（xは浮動小数点演算で計算した積）が厳密に成り立つyを求める方法が，Dekkerにより提案された。これらのEFTを繰り返し用いると，加算，減算，乗算で構成されたアルゴリズムは，複数の浮動小数点数の和を求める問題へと，誤差なしに変換できる。一方，複数の浮動小数点数の和を求める問題については，指定した任意の精度で高速に計算できるAccSumというアルゴリズムが，Rump，Oishi，Ogitaらにより提案されている。これらを組み合わせることで，悪条件の問題に対しても，指定した精度で結果を返すアルゴリズムを構成できる。

応用として，ある点が線分の右側にあるか左側にあるか，あるいは線分上にあるかを判定する問題が挙げられる。この問題は3×3の行列式の計算に帰着でき，さらにEFTを繰り返し使うことで，複数の浮動小数点数の和を正確に求める問題へと変換される。そこでAccSumを使うことで，必要な精度，すなわち結果の符号が正しいことを保証できる精度で，和を求めることができ，点と線分の位置関係が厳密に判定される。この計算法を用いて凸包の計算を行ったところ，点が密集しているような困難な問題に対しても，正しい凸包を構成できた。

EFTを組み合わせて計算を行うと誤差が全く生じないため，計算順序を気にせずに，計算量を減らすためのアルゴリズムの変形，所要メモリ量を最小化する変形などが自由にできる。本発表に先立つRump教授の発表では，この性質を利用してAccSumの計算量を最小化する方法が紹介されていた。

感想： 精度保証付き数値計算というと難しいイメージがあったが，EFTとAccSumを組み合わせて行う方式は，わかりやすく，かつ応用範囲の広い方法であると感じた。機会があれば使ってみたい。

(5) Solving Very Large Sparse Equations from Three-Dimensional Modelling (Keynote Speech)
　　Iain S. Duff (Rutherford Appleton Laboratory)

3次元のシミュレーションにおいて生じる連立1次方程式Ax=bの解法について，直接法を中心に現状のサーベイを行った。3次元のモデルをk×k×kの格子で離散化する場合を考えると，行列の次元はk³となる。この方程式を直接法で解こうとすると，行列の疎性を利用し，うまくオーダリングを行っても，計算量はk⁶，メモリ量はk⁴以下にはできないことが理論的に知られている。これは，2次元でk×kの格子でモデリングを行った場合，計算量がk³，メモリ量はk² log kにできるのに比べると，大きな差である。そのため，直接法のみを使って解こうとすると，最新の計算機を使っても，k=100程度が限界である。

そこで，反復法とのハイブリッドを行うことにより，頑健さという直接法の利点と，計算量・メモリ量の少なさという反復法の利点の両方を取り入れることが重要となる。ハイブリッド手法の例としては，(i) マルチグリッド法において粗い格子上での方程式を解くのに直接法を使うこと，(ii) 領域分割法において，部分領域での問題を直接法で解くとともに，境界領域での問題を解くのに直接法を前処理として使うこと，(iii) ブロック型の反復法において，各ブロックに対する方程式を直接法で解くこと，(iv) 不完全分解を前処理として使うこと，(v) 元の行列に対する近似行列の正確な分解を前処理として使うこと，などが挙げられる。

このうち，手法(v)の例として，静的ピボッティングと呼ばれる方法がある。いま，不定値対称疎行列の修正コレスキー分解を考える。この計算では，通常，analysisと呼ばれるステップでオーダリングとデータ構造を決定し，それに基づいて数値的分解（修正コレスキー分解）を行う。しかし，analysisの過程ではピボットの値がわからないため，そこで決めたオーダリングに基づいて分解を行っていくと，ピボットが0に近い値になり，分解を続けられなくなることが起こりうる。このような場合，analysisで決めたオーダリングを修正して分解を続けることが1つの方法である。しかしこの方法では，フィルインを最小化するように定めた元々のオーダリングの性質を崩してしまう。さらに，並列化を行っている場合は，オーダリングに基づいてデータの分散などを決めているため，それを途中で変更することは弊害が大きい。そこで，オーダリングは変えずに，ピボットに√ε （εはマシンイプシロン）程度の数を足し込むことで分解を続ける方法がある。これを静的ピボッティングと呼ぶ。静的ピボッティングでは，結果的に，A+E=LDL^T（Eは全要素が√ε以下の対角行列 ）という分解をしていることになる。

静的ピボッティングでは，誤差項Eから生じる解の誤差を減らすため，反復改良を行う。しかし，この方法では，誤差が十分に減らない場合があった。これに対して，発表者らは最近，Ax=bを解くのにFGMRES法を用い，その前処理として上記の静的ピボッティングにより得られるLDL^T分解を使う方法を提案した。この方法では，多くの問題に対し，FGMRES法を5，6反復程度実行することで，十分な精度の解が得られることが確認されている。

感想： 本講演は1時間程度の基調講演であり，ここで紹介した以外にも，マルチグリッド法を利用したヘルムホルツ方程式の解法など，Rutherford Appleton Laboratory と CERFACS における最近の研究成果が数多く紹介されていた。しかし，自分の発表の直前だったため，あまり集中して聞けずに残念だった。

(6) Unified Tight Frame Approach for Missing Data Recovery in Images (Invited Talk)
　　Raymond Chan (The Chinesse University of Hong Kong)

ランダムノイズが加わったり文字を上書きされた画像から元の画像を復元する方法として，tight frame と呼ばれる数学的概念を使った新しい手法が紹介された。画素を並べてできる2次元のベクトルをyとするとき，画像の復元とは，元の画像fにノイズnが加わってできた画像 y=f+n からfを求める問題と定義される。fを求めるための基準としては，(i) fができるだけyに近いこと，(ii) fができるだけ滑らかであること（元の画像は，エッジの部分を除いて滑らかであると想定されるから）の2つがある。そのため，βをパラメータ，D を何らかの微分作用素として，min{‖f-y‖+β‖Df‖} という最小化問題を解くことが行われてきた。たとえば ‖Df‖=∫|∇f|とした場合，この最小化問題は，変分法により，f-y+β∇・(∇f/|∇f|)=0 という偏微分方程式と等価となる。これは非線形の偏微分方程式であるため，仮想的な時間項 ∂f/∂t を導入して時間依存の問題に直し，t → ∞ の極限として解を求めること（time marching method）が行われてきた。

一方，全く別の原理に基づく方法として，tight frame を利用する方法がある。Tight frame とは，冗長性を持つ正規直交基底のことである。たとえば，Rⁿにおいて，連続する3つの要素が(1,2,1)で残りが0のベクトル，連続する3つの要素が(1,0,-1)で残りが0のベクトル，連続する3つの要素が(1,-2,1)で残りが0のベクトルをすべて合わせた集合を考えると，これは冗長性を持つRⁿの基底となる。これらを用いて，冗長性を持つRⁿの正規直交基底，すなわち tight frame を系統的に構成できる。この場合，tight frame はある n×n 行列 H₀，H₁，H₂ を縦に並べた行列 A として表現でき，正規直交性より A^*A=I となる（ただし AA^*=I は成り立たない）。Tight frame を用いて画像復元を行うには，まず復元対象の画素のみを取り出し，行列Aをかけて tight frame の各基底方向の成分を取り出し，ローパスフィルタをかけ，A^*により元の基底に戻す。式で書くと，f^(r+1)=(I-P_Λ)A^*T_λAf^(r)+P_Λy となる。ただし，P_Λ は復元対象でない画素のみを取り出す射影演算子である。f⁽⁰⁾=y としてこの操作を繰り返すことにより，画像が復元できる。実験の結果，salt-and-pepper noise（画素の一部がランダムに白あるいは黒のドットで置き換えられるノイズ）の除去，画像に上書きされた文字の除去の問題（inprinting）で，tight frame を用いた方法は偏微分方程式に基づく方法よりも良い結果を与えることが示された。

Tight frame を用いた手法は，他にも様々な画像処理に適用できる。たとえば，粗い画像から，2倍の解像度を持つ画像を求める問題。この場合，粗い画像は，求めたい高解像度の画像に，(1,2,1)というフィルタを掛けたものと見なせる。したがって，tight frame で展開して考えると，(1,2,1)に相当する成分の値がわかっているときに，他の成分を求めるという問題となる。これは，画像復元において，元の空間と tight frame で表現した空間との役割を入れ替えた問題となっており，transformed domain inprinting と呼ばれる。この問題も，行列 A^*，A を使って空間を移ることで解くことができる。応用としては，動画として記録された多数の連続する粗い画像から，高い解像度の静止画を求める問題が挙げられ，tight frame による手法は，単に双2次補間を行う手法に比べてより良い結果を得られることが示された。

(7) A Dynamic Programming Approach to Optimizing the Blocking Strategy for the Householder QR Decomposition
　　(Student paper contest, 2nd prize)
　　Takeshi Fukaya, Yusaku Yamamoto and Shao-Liang Zhang (Nagoya University)

ハウスホルダーQR分解におけるブロック化手法を，対象マシンと問題サイズに応じて自動的に最適化する方法を提案した。従来，ハウスホルダーQR分解のブロック化手法としては，固定幅ブロック化と再帰的ブロック化が主に使われていた。しかし，これらは両方とも，1個のパラメータ（前者はブロック幅，後者は再帰段数）しか含まない比較的自由度の小さい手法である。ブロック化の手法はこれ以外にも可変幅ブロック化，固定幅ブロック化と再帰的ブロック化の組み合わせなど，様々に考えられ，対象とするアーキテクチャによっては，従来使われてきた以外のブロック化が効率的な場合もありうる。

そこで本研究では，従来知られているブロック化手法を全て含み，かつ，それ以外の様々なブロック化手法も含む広いクラスとして，一般化再帰的ブロック化というクラスを定義した。本手法では，列数nの行列のQR分解において，行列を幅n₁の左側ブロックと幅n₂の右側ブロックに分割し，(i) 左側ブロックのQR分解，(ii) そこで使ったハウスホルダー変換の，右側ブロックへの作用，(iii) 右側ブロックのQR分解，(iv) 両方のブロックに対するWY表現の合成，をこの順で行う。このうち，ステップ(i)，(iii)はさらに一般化再帰的ブロック化を用いて再帰的に行われ，一方，(ii)，(iv)は行列乗算を用いて実行される。一般化再帰的ブロック化に含まれる各ブロック化手法は，葉の数がnの2分木を用いて表現できる。

ここで，ステップ(ii)，(iv)の行列乗算の実行時間を，対象マシンでの実測値に基づく性能予測モデルにより算出し，(i)，(iii)の時間については再帰的に求めることにすると，2分木が与えられたとき，対応するブロック化手法の実行時間を予測することができる。予測実行時間が最小となる2分木を求める問題は，動的計画法によりO(n²)の計算量で近似的に解けるため，これにより，準最適なブロック化手法を実用的な問題で求めることができる。本手法をいくつかの計算機で実行したところ，多くの場合にメーカー提供のLAPACKルーチンと同等の性能を自動で達成することができた。特に，Intel Xeonに対しては，MKLを超える性能を達成できた。

感想： 本会議では，毎年，東アジア地区の学生が過去に発表した論文の中から自薦・他薦をもとに学生論文賞を選ぶことになっており，今年は本論文を含めて4件が受賞した。本論文は，自動チューニングに関する国際会議 iWAPT2008 において発表したものである。本会議は応用数学に関する発表が中心であり，ハイパフォーマンスコンピューティングをテーマとした発表は，（数値実験で並列化を行っている例はあるものの）ほとんどなかった。そのため，本会議にはあまり相応しくないテーマかとも思ったが，賞を頂いたということは，自動チューニングというテーマの重要性を理解して頂けたということであり，よかったと思う。

３.　全体的な感想

本会議では，東アジア地区における応用数学の研究の現状を概観でき，興味深かった。上記で紹介したのは，金融工学関係を中心とする一部の発表のみであるが，学生の発表など，レベルがかなり高いと感じた。

今回は小規模な会議だったため，参加者の多くと話ができたのもよかった。現地の実行委員会では，学会発表以外のプログラムにもかなり力を入れており，2日目の夕方に行われたモスクの見学やブルネイ川でのクルージング，実行委員の Victor Didenko 教授宅でのバーベキューパーティー，懇親会などを通じて，多くの参加者と交流することができた。充実したプログラムを用意して下さった実行委員会に感謝したい。

４.　写真

学会が行われた Universiti Brunei Darussalam の Chancellor Hall（左）と講演風景（右）

左より，市内にある旧モスク，その池，その夜景，市の外れにある新モスクの夜景

2日目の夕方に行われたブルネイ川でのクルージング。左より，船，マングローブ，クルージング風景

2日目の夕方に行われたブルネイ川でのクルージング（続き）。左より，川岸の木の上にいた猿，川岸で見かけたワニ，夕焼けの風景

左より，懇親会の風景，懇親会で行われた学生論文賞の授賞式（2枚）