講演者:相原 研輔(東京理科大学)
題 目:クリロフ部分空間法における丸め誤差の収束性への影響
日 時:2017年1月27日(金) 18:00-
場 所:早稲田大学西早稲田キャンパス63号館 1階 数学応用数理会議室
オープニング
司会:では、時間になりましたので始めさせていただきます。第20回目の数理人セミナー、講演者は世話人の相原さんにお願いさせていただきました。では紹介は短いですが、どうぞよろしくお願いします。
講演
相原:どうもありがとうございます。東京理科大学の相原です。よろしくお願いいたします。この記念すべき第20回の数理人セミナーで講演の機会をいただきましてありがとうございます。発表タイトルは「クリロフ部分空間法における丸め誤差の収束性への影響」です。
一番初めに、もうご存じの方が多いかもしれませんけど、少し自己紹介をさせてください。名前は相原研輔といいまして、今現在東京理科大学の理学部第一部数理情報科学科で助教をしています。ちょっと宣伝をさせていただくと、この数理情報科学科は来年度から応用数学科へ名称変更します。実はもともと15年くらい前に応用数学科から数理情報科学科に名称変更したのですが、ここでまた改めて応用数学科に名前を変えて再出発ということになりました。自分はこの数理情報科学科の出身でして、博士は2014年に博士(理学)を取得しました。この後にお話しする研究分野でいうと、近い分野の方は結構、工学とか理工学の学位を取られている方が多いですが私自身は理学になります。ちょっとめずらしいかなと思います。
司会:それは応用数学科になっても同じですか?
相原:同じですね。理学研究科であることは変わらないので、理学研究科と理学部は変わらないので同じです。
司会:ごめんなさい、そういう意味じゃなくて、今のところだといろいろ取れるということなんですか、そういうわけじゃない?
相原:いや、これしか取れないです。自動的に理学になります。
研究分野は数値解析、主に数値線形代数の中で研究をいろいろ行っています。下に書いてありますけれどざっくり言うと線形計算。線形計算というのは、行列に関する数値計算手法です。それにおいて高速、高精度なアルゴリズムの開発や改良に取り組んでいます。現在は共同研究等も含めて4つくらい取り組んでいる内容があるんですけれども、一番古いのは学生のころから大規模な連立一次方程式に対する反復法というものについて研究を行っていまして、今日もこれの話をさせていただきたいと思います。
今、こちらにいらっしゃる方には大変申し訳ないのですが、前半と後半に分かれていて、話はつながっているんですが、前半は本当に基本的なところからお話をさせていただいて、後半になったら温度差がありますが、最近自分が研究している内容に絡む専門的な話も少しさせていただきたいと思います。最初は退屈かもしれませんがお付き合いください。またもちろん途中でいろいろ質問を入れていただければと思います。
A:学生さんとも共同研究をしているんですか?
相原:そうですね。私が所属する研究室に修士の学生さんがいて、この辺の研究をしているので、その修士の学生を指導しつつ共同研究という形で進めています。いろいろ停滞もしていますが、頑張ってやっています。
発表の流れがこちらです。先ほど言いましたが、最初は連立一次方程式とそれに対するタイトルにあるクリロフ部分空間法というところについて、本当に基本的な話をさせていただいてから、今日の本題として、浮動小数点演算を行う上では避けられない丸め誤差の影響についていろいろとお話をしたいと思います。特にキーワードとしては、収束速度と近似解精度、この2つに絞ってお話できればと思います。
その次に、反復法の収束性を改善するにはどうしたらいいか、ということで誤差解析やそれに基づく漸化式の修正などについて考えていきます。
最後に簡単ですが数値実験をお見せしてまとめになります。
では初めに、本研究で考えているのは、ごくごくスタンダードな正則行列を係数に持つ連立一次方程式
厳密な解というのは、
あとは、行列が疎であるという性質を使って後々出てくる反復法を使うということです。こういう行列が現れる応用分野というのは、本当にさまざまで、科学技術計算では頻繁に出てきまして、主には流体の問題とか電磁場の解析、熱の問題とか画像処理とか、いろいろなところで出てきます。私自身は、アルゴリズムの改良に取り組んでいるので、何か特定のアプリケーションがあってそれを解こうというわけではないので、こういう応用分野については知識が乏しい部分があるんですが、逆に言うとどんな分野であれこういう方程式が出てきたら適用できるような「アルゴリズム」を考えていこうということです。
ただ、何も応用の話をしないのもあれなので、一般論になりますけれども、連立一次方程式がどうやって出てくるのかという話をざっくりさせていただきます。多くは自然現象などに対する数値シミュレーションで頻繁に出てくると思います。何か現象を解析したい、シミュレーションをしたいという場合には、それを何か数学的なモデル、数理モデルで表しまして、これはなかなか解析的には解けない微分方程式となるものが多いので、コンピュータで計算できるように近似方程式を立てます。その方程式をコンピュータで解いて得られた近似解を現象にフィードバックする。これがいわゆるシミュレーションの大まかな流れかと思います。周りに書いたのはちょっとした例なんですが、何かこういうプレートの熱問題で、温度の定常分布を求めようなんていうものは、
どういうふうにして大規模な問題が出てくるのかというと、このメッシュを切るときにこれだと余りにも小さいんですが、連続の問題に近づけようと思ったらメッシュをできるだけ細かく切ればそれだけ連続の問題に近くなってきますので、そうすると途端に問題の規模が大きくなってくるということです。ざっくりした話ですが、こういうシミュレーションなどでは大規模な連立一次方程式がばんばん出てきます。
話を戻しまして、この講演では大規模な疎行列を係数に持つ連立一次方程式を考えていきましょう。解法としては、先ほども述べたような掃き出し法に基づくガウスの消去法などの直接法というものと、反復法に分かれます。こういう大規模で疎な問題ですと多くの場合反復法が有効になってきます。説明は不要かもしれませんが一応させていただきますと、反復法というのは適当な初期近似解
反復法の研究で何を目的にしているのか、というのがこちらです。速くて安くて旨いものです。下に牛丼屋(?)というふうに書いてあるんですけれども、速いというのは反復法の収束性、数学的に早く収束してくれるものが嬉しいし、安いというのは計算コストです。一反復当たりの計算量やメモリの使用量、そういったものが少ない方が好ましいです。旨いというのは、これは後々も出てきますが、数値計算になるので、厳密な解を得るというのはなかなか難しいですので、できるだけ近似精度の高いものを得ましょうということです。安くて速くて旨いものを求められるようなアルゴリズムを考えましょう。ちなみに、牛丼屋というのは学生時代にお世話になった先生の受け売りなんですが、こういうものが達成できるといい論文が書けるから、良い牛丼屋を目指してください、というような話をよくしていただきました。
さて、その反復法なんですが、ご存知かもしれませんが反復法の中にも定常反復法と非定常反復法という2つに大別されています。本研究ではタイトルにもあるクリロフ部分空間法というものを使いますが、これは非定常な反復法です。この次の数理人セミナーで名古屋大学の宮武先生が定常反復法に関する講演をしていだたけるということのようなので、それも楽しみにしておりますが、そちらで定常反復法の話はたくさんあると思いますので、ここでは非定常の話のみになります。
クリロフ部分空間法の基本的な発想として、こちらにGrowing subspaceと書きましたが、基本的には空間探索です。線形空間を広げながら、その中でいい解を探していきましょうという発想になります。クリロフ部分空間とは何かというと、このような記号で表していますが、一本のベクトルとそこに
クリロフ部分空間法では、空間条件というのがあって、
あと、関連する分野で研究をされている方は、最近IDR法という言葉を耳にしたことがあると思います。自分は博士課程時代にこちらの研究をずっとやっていたんですけれども、IDRというのはInduced Dimension Reductionの略で、日本語だと帰納的次元縮小というふうに訳されています。これもいい日本語かどうかはわからないんですけれども、最近はこういう言葉が使われています。IDRも一つの反復法群でして、発想がクリロフ部分空間法とは逆になります。考えているのはShrinking subspace、縮小するような空間を考えます。こちらの場合には、空間条件として残差ベクトルというものを考えます。残差というのは近似解
今日は、このIDR法については紹介程度で、一般的にクリロフ部分空間法について話をさせていただきますが、ちなみに現在はこのIDR法もクリロフ部分空間法の一種である、そういう解釈もできる、というような論文も出ていまして、枠組みとしてはクリロフ部分空間法のほうが大きいものになっていますが、発想がもともと全然違う、ということは興味深いです。
では、具体的な解法の説明になっていくんですけれども、クリロフ部分空間法の最も代表的なものは共役勾配法です。Conjugate Gradient法でCG法と呼ばれます。1952年にHestenesとStiefelによって出された方法です。これは有名過ぎますけれども、係数行列が正定値対称の場合に有効な方法です。そしてクリロフ部分空間法は全般そうなんですけれども、先ほどの空間条件だけですと解は一意に定まりませんので、何かしらここに条件を加えることで解を構成し、アルゴリズムが一つ定まるというものです。
CG法の場合には、空間条件に加えてこのような残差の直交条件を課します。
直交性があると何が嬉しいかというと、先ほどの反復の中で最終的には
具体的に、では解をどうやって更新していくんですかというと、漸化式の形式がこのようになっています。近似解というのは
では次に、この空間条件に加える条件を変えてみましょう。他にどんな条件が考えられるかというと、残差の最小条件というのがあります。これはいわゆる共役残差法、Conjugate Residual法でCR法と呼ばれるものです。この場合には、直交条件ではなくて
CR法は、基本的にはCG法と同じような漸化式を使いまして、係数の
残差の最小化に基づく方法でも、最終的に
アルゴリズムを少し紹介しますと、教科書などに載っているようなものですが、こちらがCG法です。係数
CR法もほぼほぼ同じような形で書くことができます。違いとしては
では一般のクリロフ部分空間法の話に戻りまして、今紹介したのは短い漸化式を用いるような算法です。短い漸化式といっているのは例えば、
CG法とかCR法というのが対称向けで、では非対称の場合はどうなのかというと、これらを拡張したBi-CG法とかBi-CR法というものがよく知られています。ちなみに、CR法と数学的に等価なMINRES法というのもよく知られていますが、ここでは割愛させていただきます。
非対称向けのBi-CG法やBi-CR法がどうやって出てくるのかというと、アルゴリズムの導出の仕方は本当にさまざまなものがあるんですが、次のようなものを考えると、比較的わかりやすいかなと思います。
もともとの
これは、方程式が
A:その対称化した方程式に前処理をつけると非対称になるんですか?
相原:そうですね。
A:非対称な系に対して、CGもしくはCRを適用していると考えるんですか?
相原:そうとも考えられるかもしれませんが、両側の前処理を考えれば対称な系に適用していることになります。ただ、厳密には正定値な対称行列にCG法を適用しているわけではないです。これは形式的なもので、あくまでアルゴリズムを導出する際にこういう捉え方ができるということです。
少し補足をすると、対称化されたものに前処理付きのCG法、CR法を適用するという発想に対して、逆に対称化されていない、
アルゴリズムですけれども、先ほどCG法とCR法をお見せしましたが、ほぼ同じような形で書けます。左側がBi-CG法、右側がBi-CR法です。チルダがついていて少し漸化式が増えてはいるんですが、近似解と残差の漸化式は全く同じです。Bi-CR法について、補助ベクトルを導入しているのは対称の場合と同じになります。なので、CG法が組めたら、Bi-CG、Bi-CRまで同じような形で自然と書けますので、簡単に実装することができます。
クリロフ部分空間法というのは反復法群で本当にいろいろなものが現在提案されています。特に対称向けだったらそんなに多くはないんですが、非対称向けはこの年号が1952年から始まって、Bi-CRはちょっと新しいですが、今はBi-CGやBi-CRを改良したものがたくさん提案されています。ご存知かもしれませんが、ざっくり紹介させていただくと、このBi-CG法に安定化多項式というものを使って残差ノルムの収束性を改善したHybrid Bi-CG法と呼ばれる枠組みの解法群とか、あとQMR、これはquasi-minimal residualの略で、Bi-CG系統だと残差の直交化というものを使いますので、そこに何かしら残差ノルムの擬似的な最小化という性質を組み入れたのがQMR系統の解法です。同じ発想でBi-CRに安定化多項式を掛けたHybrid Bi-CR法やQMR化をしたものがあります。あと先ほど紹介した帰納的次元縮小に基づくIDR系統はまだ新しい算法です。もともとは1989年でちょっと古いんですが、最近リメイクされて注目されていた解法群です。ここに挙げただけでもごく一部なんですが、沢山の解法とそのいろいろなvariantsが提案されています。クリロフ部分空間法は何か一つ万能な解法があれば一番いいんですけれども、なかなかそうはいかいないので、さまざまな方程式の特徴に合わせて得意不得意があるようないろいろな算法が提案されているということです。
今、ここには「短い漸化式」と書きましたが、だとすると当然「長い漸化式」の算法もあります。それはいわゆるGMRES法とかGCR法といった算法が長い漸化式と表現されます。その中でもrestart版とかtruncated版という効果的なものがありますけれども、この講演では割愛させていただきます。
ここまで、連立一次方程式とクリロフ部分空間法の本当に基本的なところをお話しさせていただきました。ここから本題の丸め誤差というものについて考えていきたいと思います。
計算する際には、今現在の計算機ですと倍精度浮動小数点数というのが基本になります。これも今さら言うこともないかもしれませんが、符号部1bit、指数部11bit、仮数部が52bitで計64bitで表現するような数になります。こういうものを使わなければいけないので、実数の範囲で考える0.1という数も、2進数に変換すると循環小数になってしまうので、倍精度浮動小数点数で表そうとすると仮数部を52bitに制限しなければいけなくて、それがいわゆる「丸める」という操作になります。近似値しか扱えないので、ここで丸め誤差が生じます。
こういう浮動小数点数を使うと、当然近似計算で、どれくらいの精度を持っているのか?という話になります。計算精度としては、
あとは、
反復法の話に戻りまして、その丸め誤差という観点から理想的な反復法とはどういうものか、という話をさせていただきます。キーワードとして、収束速度と近似解精度という二つを取り上げます。
まず収束速度ですけれども、これは数値解析の世界でいういわゆる何次収束するか、という話ではありません。残差は高々
近似解精度ですが、これは後でお話ししますので、先に数値例を一つ紹介したいと思います。
よく学生とクリロフ部分空間法の話をしていると、手計算で簡単に解を求められるのに、クリロフ部分空間法を適用すると解けないような問題はあるか…という話題が出てきます。そんなのは、専門家の方はいくらでもご存じかとは思いますが、その最たる例として、こういうcyclicな行列というのが挙げられます。これは右上に成分が一個あって、下副対角にも何か成分が入っていて、ほかは全部0です。これの固有値分布というのは、複素平面上でこんなふうに円の形になります。これは簡単な計算で確認することができます。
こういう問題に、ここではBi-CG法を適用してみたときの収束履歴を見ていきます。このグラフは相対残差ノルムの収束履歴をあらわしていまして、横軸が反復回数、縦軸が相対残差の2-ノルムを対数でプロットしたものです。なお、この計算はまだ倍精度ではなくて、約8倍精度の高精度な計算を行っています。だから、これはほとんど誤差が入らない、理論的な性質をある程度維持していると考えられます。実際に、
この赤いラインが倍精度で計算したものです。まず
この行列は各行に一つしか成分がないので、簡単に割っていけば厳密な解というのがわかるんですが、こういうものにクリロフ部分空間法を適用すると、あっという間に収束しないような場合が出てくるということです。
B:これは
相原:そうです。
B:それでもう一回掛けるとまた下がるという。
相原:そうです。くるくる回るということです。
B:クリロフ部分空間法が収束しようと思ったら、解を
相原:ここの残差が減っていないということですか?
B:最後までいかないと、非ゼロ成分の埋まり方だけで理解できますか?これは
相原:
B:そうかそういうふうに理解するのか。
相原:全てのゼロ点がそういう何か固有値になるまで反復をしないと残差が結局減らないというふうに解釈ができると思います。
B:初期値によらずに、やはりどうしてもこういう形になってしまいますか?
相原:そうですね。右辺には結構寄ります。今右辺項を乱数で取って、結局各固有ベクトルの成分が均等にばらばらに入るように組み立てているので。右辺が特定の固有ベクトルだけで張られているような場合だと、もちろんそこで落ちます。
極端な例ではありましたが、収束速度の話をさせていただきました。通常はある程度滑らかに減少していくのに、それが滑らかに落ちずに反復回数が増えてしまうとか、そういうものも収束速度の低下ということになります。
一方で、近似解精度に関してですけれども、精度というと当然誤差をはかりたくなるわけでして、誤差は真の解
一方、この後にお話しさせていただきたいのは、残差は見ているんですが、この
この後residual gapについて、もう少し詳しくお話ししたいと思いますが、理想的なという意味では、きちんと要求精度を満たす解法が望ましいということになります。
ではresidual gapは何が問題なのかという話なんですが、これは収束判定について考えるとよくわかります。一般的な収束判定はCG法などを思い浮かべていただくといいんですが、漸化式から求まる残差
これも数値例を紹介したいと思います。行列の話は後ほどしますので、あるテスト行列に対して今度はBi-CR法を適用してみます。すると、こんなふうに残差ノルムの振る舞いが出てきます。横軸が反復回数、縦軸が相対残差ノルムで収束判定が相対残差の意味で
収束性の改善については、いろいろな手法があります。あまり一つ一つ細かくはお話しできないんですが、素朴な発想を幾つかお話しさせていただくと、まず一つは多倍長計算の利用です。先ほど擬似的な8倍精度で計算したものは、うまくいっていましたので、多倍長を使えばいいでしょうという発想があります。実際に多倍長計算を使うと収束速度、residual gapともに改善されます。ただし、普通の倍精度の演算に比べると、かなり計算は重たくなります。また、コスト面を考慮して部分的に利用する混合精度演算ではどうかというと、経験的に部分的な利用では改善は難しいということがわかっています。ただし、我々は多倍長というのは収束性を改善するために直接的には使わないんですが、誤差解析の際に有効活用することができます。というのは、先ほどのように多倍長で計算したものと倍精度で計算したものを比較してみると、丸め誤差の影響というのを可視化することができるので、どのくらいの誤差が混入しているか、あるいは混合精度演算とかも使って、どういうところで誤差が大きく拡大しているかとか、そういうのを調べるのにはすごく役立てることができます。そういう意味で多倍長を研究の中で使ったりしています。
あとリスタートです。リスタートは、適当なところで反復を打ち切ってそこで得られた近似解を新しい初期値として再出発する手法になります。リスタート手法に関しては、丸め誤差の影響を排除することができます。そこまでに積もった丸め誤差というものを完全に捨てる事になりますので有効ですが、一般的にはリスタートするまでにつくったクリロフ部分空間の基底というものも捨ててしまうことになるので、そこから再出発した場合には、反復を続けた場合に比べて反復回数は増加してしまうということがあります。
あと真の残差の計算です。これはresidual gapというのを改善するのだったら、漸化式を使わないでダイレクトに真の残差を計算すればいいじゃないか、という考え方です。これは少し古い教科書などを見るとそういう手法が実際に書かれたりしていますが、現在はいろいろ解析されていまして、真の残差を利用するとresidual gapは改善できるんだけれども、先にお話ししたような収束速度の低下に繋がってしまうということが知られています。これは簡単に実験的にも確認することができるんですが、真の残差をそのまま使うのは良くないというふうに今は考えられています。
では、これらはそれぞれメリット・デメリットがあるんですけれども、今日はもう少し思考を凝らした「ベクトル更新のグループ化」というものについてお話しできればと思います。ここから込み入った話にはなるんですが、丸め誤差の影響を抑制するように漸化式を修正して、さらに低コストで収束速度を維持しつつresidual gapも改善したい、そういうものを目標につくられた手法です。
では、ここから収束性の改善の話をさせていただきます。普段、学会等で講演する場合にはここからお話をするので、ここにいらっしゃる方にとって前半はかなり前振りが長いなと思われたかもしれませんが、ここからがより専門的な内容になります。
まずは、収束性の改善のために丸め誤差の影響というものを考えなければいけないので、誤差評価について紹介させていただきます。ここでは、浮動小数点演算で行う行列ベクトル積について考えていきたいと思います。
行列ベクトル積
考えるのは近似解と残差を更新するこういった漸化式です。本来はここに係数
この不等式の意味するところは何かというと、この漸化式で
また少し見方を変えてみると、こういう残差の急激な減少というのは、
A:局所的な誤差の上限のところの
相原:すみません。ここは抜け落ちています。
A:何か不思議ですね。ここで出ている結論が、本当は残差ノルムの急激な減少の更新が、それをすると収束速度が全体としては落ちる。
相原:そうですね。その話はこの後でもさせていただこうと思っていたんですが、一般に反復法を習うと残差を減らしましょうということで、残差が減るというのはすごく嬉しいことなんですが、局所的に見て急激に、特に残差が一旦大きくなってから急激に減少してしまったりするのは、実は丸め誤差の観点からすると望ましくないということがあります。そのあたりは反復法を専門に研究している方でないと、多分そんなに知られていないことなので、ぜひそういう認識を使っている方々に持ってもらえたらと思っています。また後でその話もします。
では二つ目の項目なんですけれども、今の局所的な誤差というものから、今度はresidual gapを評価しましょう。Residual gapの上限としては、こんな形で与えることができます。細かい部分は置いておいて、大事な部分は何かというと残差ノルムの和です。なので、残差ノルムが振動して初期残差ノルムというのを一つの基準として考えると、それよりも非常に大きくなるような中間的な残差ノルムというものがあると、この上限というのが相対的に大きくなり得ます。これによってresidual gapが大きくなるので、近似解精度の劣化に繋がってしまうということになります。
結論として何が言いたいのかというと、相対的に大きな局所誤差の蓄積を避けて、さらにresidual gapを小さくするには、残差ノルムというのは滑らかに減少していってくれることが理想的な振る舞いということになります。
あまりこういう話というのは普段学会とかでもほとんどないと思うんですけれども、振動のパターンというのを四通りで書いてみて、それと収束性との関係というのを、現実にはこんなに単純ではないんですが、あくまでイメージとして持っていただくといいかなということで書かせていただきました。左上が一番理想的なパターンで右下が最悪なパターンです。振動の大きさとノルムが初期残差に比べてどれくらい大きくなるかというのを表現しています。
左上の理想的なパターンは振動が小さくて、
一方で右上は、初期残差に比べて大きなノルムはないけれども、一回一回の振動が激しい場合です。これは局所的に大きな誤差が積み重なってしまうので、速度に結構影響します。ただ初期残差よりは小さいので、residual gapに関してはきっと大丈夫でしょう。
では左下はというと、振動自体は小さいです。振動は小さいけれども初期残差に比べて小さな振動をしながらノルムが一度大きくなってから収束していくものです。これは速度に関しては問題ないだろうと思われますが、residual gapが大きく出てしまうと考えられます。
最後に、これは2桁で基準を取っていますが、振動も大きいしノルムも大きい、という最悪なパターンです。どうしようもないということになります。
もちろん現実はこんな単純には分類できないんですが、イメージとしてこういうようなものを考えてみると、収束振る舞いと収束性との関係が見やすいかと思います。ようするに左上のものが理想的なわけです。残差が減れば嬉しいのだけれども、ただ急激に減ればいいというものではないという、先ほどのお話になります。
B:左下で言うと、最後のresidual gapは3回目から4回目の急激な減少が最後residual gapにあらわれると?
相原:Residual gapに影響するのは最大のピークです。これで言うとピークが4乗くらいですので。
B:その差がresidual gapになるんですね。
相原:はい。だから4桁くらい、一番上のピークが
B:さっきの解説だと、
相原:Residual gapに関しては、グローバルなものですので、上限を相対評価するときには初期残差で割って評価をするんですね。だから初期残差に比べて、0番から
時間もちょっと押してきました。あと、最後のほうは数値実験をざっとお見せするだけなので、改善手法について少しお話をさせていただければと思います。
残差更新のグループ化というアイデアがあります。これは私のオリジナルではなくて、1994年のNeumaier先生と96年のSleijpen先生、van der Vorst先生の論文で示されています。
近似解と残差の漸化式を考えたときに、ある単調増加数列
元々の列として、
B:ここで飛ばした
相原:1回は求めています。一度は求めて、何かここで大きな誤差の混入といった状況を検知したら、そこで再計算するということです。
B:アルゴリズムの順番としてはどうなっているんですか?上の列と下の列の。
相原:順番としては、
それでは数列
まず、同時にはなかなか難しいので1番から考えてみますと、
B:再計算するノルムは
相原:
B:本当の正解は
相原:そうですね。
B:ただ、そうじゃなくて、ある意味、解ベクトルとの関係をうまく持たせて、あえて摂動させて何か全体としてうまく収束させようとしているのか。
相原:これも一応
B:丸め誤差が入ると違うんじゃないんですか。質問の仕方を変えると、それを
相原:丸め誤差が入るので異なりますが、この場合は実は
B:じゃ、一緒になる?
相原:はい。
B:最初のほうで
相原:その話との関係は、今は
B:むしろそっちのほうが誤差を含んでいる。
相原:はい、含んでいます。なので、後半で真の残差を計算するとまずいというのが正確な表現になります。
B:わかりました。ありがとうございます。
相原:今日お話した先ほどのSleijpen先生の論文とは違うんですが、van der Vorst先生がそういう最初の段階だけ真の残差をうまく計算してあげて、置き換えをうまくやればresidual gapを防げますよといった手法も別の論文として出されていて、後半さえ使わなければうまくはいきます。パラメータを用いた適当な条件というのも、実は最初の段階だけ部分的にこういう計算をしましょうという発想になっています。
A:今、
相原:はい、
A:
相原:ここの数ですか。
B:はい。それはその置き換えをするまでは、1回置き換えをして次の置き換えをするまでの間には、ずっとその
相原:とっておかないといけないですね。ただ、1反復するごとにこの
2番目が残差ノルムの大きな振動、急激な減少を回避するにはどうしたらいいかという話なんですが、これは何かというと、先ほどは大きな中間残差を避けるように
そこで、大きな振動がある場合には、例えば
ただし、先ほどと同じようにここでの再計算にもやはり行列ベクトル積が余分にかかってしまうので、頻繁にはできないわけです。だからこれも適当なパラメータを入れてあげて、本当に振動が大きい場合だけとか、部分的にこういう処理を入れてあげることで、strategicalな手法ではあるんですが、効率的に丸め誤差の影響を抑制できるということになります。
では、ここまでがメインの話になりますけれども、以上の話はCG法などの漸化式をメインに考えてきましたが、冒頭でCR法あるいはBi-CR法は少し漸化式の形式が異なりますという話をしました。そういう場合にも適用できるかというと、それは適用できます。実際に効果があることを示した論文が最近公開されまして、私がやった内容になっています。
発想自体は全く同じなんですけれども、概要だけお話させていただきます。Bi-CG系統とBi-CR系統で漸化式はどう違うかというと、Bi-CR系統の場合には補助ベクトル
漸化式による違いなんですが、比較のために並べてみますと、Bi-CG系統の漸化式では先ほどお見せしたようなこんな感じです。局所誤差は隣り合う残差ノルムの和で書けるし、residual gapも残差ノルムで書けます。Bi-CR系統の漸化式の場合には、局所誤差は新たに導入された補助ベクトル
A:Bi-CGのresidual gapに比べてBi-CRのresidual gapは随分大きくなりそうな…。
相原:そうですね。
A:
相原:そうなんです。これが実際にシャープかどうかというのはなかなか議論が難しくて、結構over estimateかもとは思っています。ただし、これを元に同じようにグループ化というアイデアを入れてみると、実験的には結構効果があります。この評価自体はそれほどシャープではないかもしれないですが。
では、最後に数値実験なんですが、数値実験のほうは私自身がBi-CR系でそういうことをやったので、実験はBi-CR法でお見せしたいと思います。実際に効果があるということを示します。
計算条件はこのようになっていて、初期値は0で収束判定条件は
テスト行列は、これは最初にお見せした行列でして、フロリダのSparse Matrix Collectionからこの行列をテストとして持ってきました。詳細はここに書いてあるとおりで、右辺は乱数ベクトルで与えています。
最初にお見せしたように、まず従来のBi-CR法の収束履歴ですが、大きく振動してしまっていて、最終的にここでギャップが生じていて4桁くらいは差が出ているという状態のものです。
これに対して、中でグループ化というものを行って計算する提案手法がこちらです。今、横軸が反復回数の
A:上にそのノルムが跳ね上がった分、下の限界精度の出方みたいな、例えばこの場合で言えば
相原:ここですか。そうですね。これよりもっと下げようと思うと…。
A:どこかの頭を叩かなければいけないですか。
相原:頭を叩かないといけないですね。おっしゃる通りです。そういう意味では、そもそも漸化式から求まる残差は、この先は計算していないんですけれども、これがどこまで下がるかというのもあるので、これは多分下がっていくので、真の残差の計算、まさに
今、これはギャップが改善できましたという話なんですが、もうちょっと難しい問題にします。
これはよく反復法のテストに使われる問題で、こういう偏微分方程式の離散化から出てくるような問題を対象に解いてみました。解法は同じBi-CR法です。
これは従来の方法とそのグルーピングを行った方法を一緒に重ねてしまっていますが、赤いラインが従来法で青いラインが提案法です。振動も大きくてノルムも大きいので、さっきの4つのパターンでいうと最悪なパターンになります。
従来の方法は
この表なんですけれども、反復回数は4282回で、ここのAdd. MVsというのは、先ほどグループ化をして残差を再計算するときに余分な行列ベクトル積がかかりますと言いましたが、それにどれくらいかかったかという回数です。175回かかっています。この場合はそもそも収束していないので、収束するようになっただけで十分かと思うんですが、Bi-CR法だと1反復当たり2回行列ベクトル積がかかりますので、全体として行列ベクトル積は8500回くらいなんです。行列ベクトル積8500回に対して、その余分なコストという意味では大体175回の行列ベクトル積なので、ほとんど大したことないです。微々たる追加の計算で、これくらいの収束改善が行えるということなので、これはやらないよりは標準的に組み込んでおいたほうがきっと良いでしょうと思います。
最後にまとめになりますけれども、本講演ではクリロフ部分空間法の中でも短い漸化式を用いる解法において、その丸め誤差の影響を議論してきました。基本的には残差ノルムの振動、Bi-CR系だったらその補助ベクトルの振動になりますけれども、そういったものは収束速度や近似解精度に大きく影響しますということです。その丸め誤差の蓄積、拡大を回避するには滑らかな減少というのが大事です。途中でも議論が出ましたが、残差というのはただ減ればいいというものではなくて、急激に減っても困るし、少しずつ滑らかに減少していってくれるのが丸め誤差の観点からすると理想的な振る舞いになります。それを達成するために残差更新のグループ化というものについて紹介をさせていただきました。
今後の課題というか展望なんですけれども、今現在、修士課程の学生と共同研究をやっています。残差ノルムの振る舞いを滑らかにするには、古典的なスムージングという技法があります。今日は紹介できなかったんですが、残差ノルムを滑らかにするという意味では有用なんですが、今までは丸め誤差の影響を抑えるというふうには使われていなかったんです。このスムージングを少し工夫してあげると、数値的に安定な、その丸め誤差の影響を抑制するような実装法というのができそうだということで検討しています。既に講演を聞いている方もいらっしゃるかもしれませんが、学生と一緒にやっています。これができると何が嬉しいのかというと、同じような手法ではあるんですが、こういうグループ化といった手法に比べると、残差ノルムは単調減少するし、実装は簡単だし、パラメータなどの設定も不要だし、余分な計算コスト、行列ベクトル積もかかりません。そういう利点がいろいろあるので、スムージングでできるんだったらそのほうが嬉しいということで、このグループ化を超えるような手法ができないかなと考えています。
以上です。ありがとうございました。
質疑
司会:どうもありがとうございました。既にいろいろお話しいただいていると思いますが、何か議論とかコメントがあれば。
A:こういった研究でつくられたコードというのは公開されているんですか。
相原:よくそういった質問を、応用分野の方に言われるんですが、すみません公開はしていません。公開できればと思うんですけれども、公開できるほどの腕がないというか、いろいろ研究がまだまだ途中というのもあるんですけれども、いずれそういう公開とかもできたらと思っています。
司会:ほかにありませんか。いいですか、質問しなくて。
クロージング
司会:では一旦終わりにして、また何か話すことがあれば、また後でお話ししましょう。どうもありがとうございました。