相原 研輔

スライドショーには JavaScript が必要です。

講演者：相原研輔（東京理科大学）
題　目：クリロフ部分空間法における丸め誤差の収束性への影響
日　時：2017年1月27日(金) 18:00-
場　所：早稲田大学西早稲田キャンパス63号館 1階数学応用数理会議室

オープニング

司会：では、時間になりましたので始めさせていただきます。第20回目の数理人セミナー、講演者は世話人の相原さんにお願いさせていただきました。では紹介は短いですが、どうぞよろしくお願いします。

講演

相原：どうもありがとうございます。東京理科大学の相原です。よろしくお願いいたします。この記念すべき第20回の数理人セミナーで講演の機会をいただきましてありがとうございます。発表タイトルは「クリロフ部分空間法における丸め誤差の収束性への影響」です。

　一番初めに、もうご存じの方が多いかもしれませんけど、少し自己紹介をさせてください。名前は相原研輔といいまして、今現在東京理科大学の理学部第一部数理情報科学科で助教をしています。ちょっと宣伝をさせていただくと、この数理情報科学科は来年度から応用数学科へ名称変更します。実はもともと15年くらい前に応用数学科から数理情報科学科に名称変更したのですが、ここでまた改めて応用数学科に名前を変えて再出発ということになりました。自分はこの数理情報科学科の出身でして、博士は2014年に博士（理学）を取得しました。この後にお話しする研究分野でいうと、近い分野の方は結構、工学とか理工学の学位を取られている方が多いですが私自身は理学になります。ちょっとめずらしいかなと思います。

司会：それは応用数学科になっても同じですか？
相原：同じですね。理学研究科であることは変わらないので、理学研究科と理学部は変わらないので同じです。
司会：ごめんなさい、そういう意味じゃなくて、今のところだといろいろ取れるということなんですか、そういうわけじゃない？
相原：いや、これしか取れないです。自動的に理学になります。

　研究分野は数値解析、主に数値線形代数の中で研究をいろいろ行っています。下に書いてありますけれどざっくり言うと線形計算。線形計算というのは、行列に関する数値計算手法です。それにおいて高速、高精度なアルゴリズムの開発や改良に取り組んでいます。現在は共同研究等も含めて４つくらい取り組んでいる内容があるんですけれども、一番古いのは学生のころから大規模な連立一次方程式に対する反復法というものについて研究を行っていまして、今日もこれの話をさせていただきたいと思います。
　今、こちらにいらっしゃる方には大変申し訳ないのですが、前半と後半に分かれていて、話はつながっているんですが、前半は本当に基本的なところからお話をさせていただいて、後半になったら温度差がありますが、最近自分が研究している内容に絡む専門的な話も少しさせていただきたいと思います。最初は退屈かもしれませんがお付き合いください。またもちろん途中でいろいろ質問を入れていただければと思います。

A：学生さんとも共同研究をしているんですか？
相原：そうですね。私が所属する研究室に修士の学生さんがいて、この辺の研究をしているので、その修士の学生を指導しつつ共同研究という形で進めています。いろいろ停滞もしていますが、頑張ってやっています。

　発表の流れがこちらです。先ほど言いましたが、最初は連立一次方程式とそれに対するタイトルにあるクリロフ部分空間法というところについて、本当に基本的な話をさせていただいてから、今日の本題として、浮動小数点演算を行う上では避けられない丸め誤差の影響についていろいろとお話をしたいと思います。特にキーワードとしては、収束速度と近似解精度、この２つに絞ってお話できればと思います。
　その次に、反復法の収束性を改善するにはどうしたらいいか、ということで誤差解析やそれに基づく漸化式の修正などについて考えていきます。
　最後に簡単ですが数値実験をお見せしてまとめになります。

　では初めに、本研究で考えているのは、ごくごくスタンダードな正則行列を係数に持つ連立一次方程式 $Ax=b$ という形をしたものです。 $A$ というのは与えられた $n\times n$ の正則な行列で、右辺項 $b$ というのが $n$ 次元のベクトルとしたとき、 $x$ を求めるという問題です。
　厳密な解というのは、 $x^*$ で表しますが、当然 $b$ に $A^{-1}$ を掛けたものです。しかし、ここで対象としているのは行列の特徴として非常に大規模なものです。大規模といっているのは、 $n$ が数万から大体数千万くらいのものを対象としています。さらに、成分のほとんどが $0$ であるような疎行列です。こうなってくると、まず大規模であることから当然手では計算できないのでコンピュータによる数値計算が不可欠になってきます。こんなところから…？と思われるかもしれませんが、先ほど理学部ですと言いましたが、理学部で純粋系の先生方とお話をしていると、今でも連立一次方程式の研究をしていますというと、そんなものは $A$ の逆を掛けるだけではないのかとか、線形代数で習ういわゆる掃き出し法を行えばよいでしょうというふうに言われることがあります。なので、なかなかそう上手くはいかない、というところから説明をしないといけません。結構そこは異分野の先生とお話しするときには気を使っている部分です。
　あとは、行列が疎であるという性質を使って後々出てくる反復法を使うということです。こういう行列が現れる応用分野というのは、本当にさまざまで、科学技術計算では頻繁に出てきまして、主には流体の問題とか電磁場の解析、熱の問題とか画像処理とか、いろいろなところで出てきます。私自身は、アルゴリズムの改良に取り組んでいるので、何か特定のアプリケーションがあってそれを解こうというわけではないので、こういう応用分野については知識が乏しい部分があるんですが、逆に言うとどんな分野であれこういう方程式が出てきたら適用できるような「アルゴリズム」を考えていこうということです。

　ただ、何も応用の話をしないのもあれなので、一般論になりますけれども、連立一次方程式がどうやって出てくるのかという話をざっくりさせていただきます。多くは自然現象などに対する数値シミュレーションで頻繁に出てくると思います。何か現象を解析したい、シミュレーションをしたいという場合には、それを何か数学的なモデル、数理モデルで表しまして、これはなかなか解析的には解けない微分方程式となるものが多いので、コンピュータで計算できるように近似方程式を立てます。その方程式をコンピュータで解いて得られた近似解を現象にフィードバックする。これがいわゆるシミュレーションの大まかな流れかと思います。周りに書いたのはちょっとした例なんですが、何かこういうプレートの熱問題で、温度の定常分布を求めようなんていうものは、 $xy$ 平面でモデルを立てて、例えばこういうラプラス方程式ですけれども、現実の問題だとなかなか厳密解を求められないことが多いので近似して数値計算をします。連続の問題をこのようにメッシュを切って離散化し、各格子点を未知数として方程式を立てて解いてあげます。ここで、その未知数というのが解になっているような連立一次方程式が立てられるので、それを解くことで現象にフィードバックできるわけです。
　どういうふうにして大規模な問題が出てくるのかというと、このメッシュを切るときにこれだと余りにも小さいんですが、連続の問題に近づけようと思ったらメッシュをできるだけ細かく切ればそれだけ連続の問題に近くなってきますので、そうすると途端に問題の規模が大きくなってくるということです。ざっくりした話ですが、こういうシミュレーションなどでは大規模な連立一次方程式がばんばん出てきます。

　話を戻しまして、この講演では大規模な疎行列を係数に持つ連立一次方程式を考えていきましょう。解法としては、先ほども述べたような掃き出し法に基づくガウスの消去法などの直接法というものと、反復法に分かれます。こういう大規模で疎な問題ですと多くの場合反復法が有効になってきます。説明は不要かもしれませんが一応させていただきますと、反復法というのは適当な初期近似解 $x_0$ というのを与えて、それを更新しながら徐々に厳密解に近づけていくような数値計算法です。これはイメージ図ですけれども、厳密解 $x^*$ がここにあったとしたら、適当な解を与えて、これを更新していって近づけばいいなということです。後で詳しくお話ししますが、多くの場合、近似解などの更新には何かしらの漸化式が用いられます。これが後々メインの話になってきます。
　反復法の研究で何を目的にしているのか、というのがこちらです。速くて安くて旨いものです。下に牛丼屋（？）というふうに書いてあるんですけれども、速いというのは反復法の収束性、数学的に早く収束してくれるものが嬉しいし、安いというのは計算コストです。一反復当たりの計算量やメモリの使用量、そういったものが少ない方が好ましいです。旨いというのは、これは後々も出てきますが、数値計算になるので、厳密な解を得るというのはなかなか難しいですので、できるだけ近似精度の高いものを得ましょうということです。安くて速くて旨いものを求められるようなアルゴリズムを考えましょう。ちなみに、牛丼屋というのは学生時代にお世話になった先生の受け売りなんですが、こういうものが達成できるといい論文が書けるから、良い牛丼屋を目指してください、というような話をよくしていただきました。

　さて、その反復法なんですが、ご存知かもしれませんが反復法の中にも定常反復法と非定常反復法という２つに大別されています。本研究ではタイトルにもあるクリロフ部分空間法というものを使いますが、これは非定常な反復法です。この次の数理人セミナーで名古屋大学の宮武先生が定常反復法に関する講演をしていだたけるということのようなので、それも楽しみにしておりますが、そちらで定常反復法の話はたくさんあると思いますので、ここでは非定常の話のみになります。
　クリロフ部分空間法の基本的な発想として、こちらにGrowing subspaceと書きましたが、基本的には空間探索です。線形空間を広げながら、その中でいい解を探していきましょうという発想になります。クリロフ部分空間とは何かというと、このような記号で表していますが、一本のベクトルとそこに $A$ のべき乗を掛けていったもので張られるような空間、これをクリロフ部分空間と呼びます。この場合は右辺項を使っているので、 $A$ と $b$ とで張られるクリロフ部分空間です。
　クリロフ部分空間法では、空間条件というのがあって、 $k$ 反復目の近似解はこの $k$ 次のクリロフ部分空間の中から選びましょうというものです。クリロフ部分空間法は、こういうふうな枠組みで捉えることができます。そうすると空間は広がっていきますので、 $n$ 次元のユークリッド空間の中で考えると、どこかに解はあるはずなので、近似解を更新して空間を広げていけば、いつかはここの解が見つかるんじゃないですか、というような発想です。
　あと、関連する分野で研究をされている方は、最近IDR法という言葉を耳にしたことがあると思います。自分は博士課程時代にこちらの研究をずっとやっていたんですけれども、IDRというのはInduced Dimension Reductionの略で、日本語だと帰納的次元縮小というふうに訳されています。これもいい日本語かどうかはわからないんですけれども、最近はこういう言葉が使われています。IDRも一つの反復法群でして、発想がクリロフ部分空間法とは逆になります。考えているのはShrinking subspace、縮小するような空間を考えます。こちらの場合には、空間条件として残差ベクトルというものを考えます。残差というのは近似解 $x_k$ に対して $b-Ax_k$ で定義されるもので、residual になりますので $r_k$ で表現しています。この残差ベクトルが縮小する空間の中に入るというような条件を課すと、空間が縮小するにつれて、残差も $0$ に近づいていきます。そして、対応するような解が得られれば、それは良い近似解になるでしょうというような発想です。
　今日は、このIDR法については紹介程度で、一般的にクリロフ部分空間法について話をさせていただきますが、ちなみに現在はこのIDR法もクリロフ部分空間法の一種である、そういう解釈もできる、というような論文も出ていまして、枠組みとしてはクリロフ部分空間法のほうが大きいものになっていますが、発想がもともと全然違う、ということは興味深いです。

　では、具体的な解法の説明になっていくんですけれども、クリロフ部分空間法の最も代表的なものは共役勾配法です。Conjugate Gradient法でCG法と呼ばれます。1952年にHestenesとStiefelによって出された方法です。これは有名過ぎますけれども、係数行列が正定値対称の場合に有効な方法です。そしてクリロフ部分空間法は全般そうなんですけれども、先ほどの空間条件だけですと解は一意に定まりませんので、何かしらここに条件を加えることで解を構成し、アルゴリズムが一つ定まるというものです。
　CG法の場合には、空間条件に加えてこのような残差の直交条件を課します。 $k$ 番目の残差というのが $k$ 次のクリロフ部分空間法と直交する。CG法は、最適化などに詳しい方はご存知かと思いますが、線形のCG法というといわゆる二次関数の最小化問題に基づくもので、壺型の二次関数があって、そこの最小点を求めることと $Ax=b$ を解くことが等価なので、勾配が $0$ になるところを求めればよいわけです。勾配というのがこの残差に一致するので、それが $0$ になるようなところを求めるという発想に基づいて解法が導出されています。よく最適化分野ではそのように習いますし、もともとの発想はそうなんですが、クリロフ部分空間法として考えると、こういう直交条件を課すというふうに解釈することができます。
　直交性があると何が嬉しいかというと、先ほどの反復の中で最終的には $x^*$ が欲しいわけなんですけれども、もしこの反復を進めて空間が広がっていって最大で $n$ 次元のユークリッド空間が張れたとしたらそれと直交するようなベクトルというのは $0$ しかありませんので、残差が必ず最後は $0$ になるという原理です。だから、反復法とは言いつつも有限回の操作で解が求まる、理論的にはそういう性質を持った解法です。
　具体的に、では解をどうやって更新していくんですかというと、漸化式の形式がこのようになっています。近似解というのは $x_k$ に探索方向 $p_k$ というものを計算して、そこに係数 $\alpha_k$ 、これはステップ幅とも言いますが、それを掛けて足してあげて次の解とします。それに対応するような残差というのも、こういう漸化式で更新します。これは $p_k$ に $A$ が掛かって符号がプラスからマイナスに変わっただけの形式なんですが、これが残差ベクトルに相当します。これは簡単な計算で出すことができまして、近似解の漸化式に両辺 $A$ を掛けてあげて、右辺の $b$ から引いてあげるとこの残差の漸化式が出てきます。探索方向 $p_k$ に関しても、同じような漸化式で計算しますが、係数 $\alpha_k, \beta_k$ に関しては、直交条件などから自動的に決まるものになっています。まずこれがCG法です。

　では次に、この空間条件に加える条件を変えてみましょう。他にどんな条件が考えられるかというと、残差の最小条件というのがあります。これはいわゆる共役残差法、Conjugate Residual法でCR法と呼ばれるものです。この場合には、直交条件ではなくて $k$ 反復目の残差はそのノルムが最小になるように選びましょうという発想です。
　CR法は、基本的にはCG法と同じような漸化式を使いまして、係数の $\alpha_k, \beta_k$ が異なる解法というふうに解釈されます。ただし、実装する際には漸化式にも違いがあります。ここでお話しをさせていただくと、ここに行列 $A$ と探索方向 $p_k$ の積というものがありますが、この部分に補助ベクトル $q_k$ というのを導入しまして、 $q_k$ というのはベクトル $Ap_k$ に等価なものなんですが、これ自身も漸化式でつくります。この $p_k$ の漸化式の両辺に $A$ を掛けただけなんですが、こういう形でアルゴリズムが実装されます。なので、CR法とCG法は理論的には $\alpha_k, \beta_k$ の計算が違うだけなんですが、実装するときには漸化式の違いもあります。これは後々、後半の話に響いてくるので頭の隅に置いておいていただければと思います。
　残差の最小化に基づく方法でも、最終的に $n$ 反復したら $n$ 次元空間の中で最小値は $0$ ですので、結果として有限回の操作で解が求まるという性質は変わりません。

　アルゴリズムを少し紹介しますと、教科書などに載っているようなものですが、こちらがCG法です。係数 $\alpha_k$ を計算して、近似解と残差を漸化式で更新して、 $\beta_k$ を計算して探索方向を計算する。あとはこれを繰り返すだけです。
　CR法もほぼほぼ同じような形で書くことができます。違いとしては $\alpha_k, \beta_k$ の計算が少し違うというのと、補助ベクトルの $q_k$ というのが導入されているという点だけです。

　では一般のクリロフ部分空間法の話に戻りまして、今紹介したのは短い漸化式を用いるような算法です。短い漸化式といっているのは例えば、 $x_{k+1}$ が $x_k$ にスカラー倍のベクトルを足すといった形で書けているものです。ここまでは、短い漸化式で書ける対称行列向けの算法を紹介しました。言いそびれましたが、CG法は正定値向けであるのに対して、CR法は正定値でないような場合にCG法に比べると良い収束性を示す場合が多いです。
　CG法とかCR法というのが対称向けで、では非対称の場合はどうなのかというと、これらを拡張したBi-CG法とかBi-CR法というものがよく知られています。ちなみに、CR法と数学的に等価なMINRES法というのもよく知られていますが、ここでは割愛させていただきます。
　非対称向けのBi-CG法やBi-CR法がどうやって出てくるのかというと、アルゴリズムの導出の仕方は本当にさまざまなものがあるんですが、次のようなものを考えると、比較的わかりやすいかなと思います。
　もともとの $Ax=b$ に対してこういう拡大された行列を考えまして、ここに $A$ と $A^\top$ を置きます。 $\tilde x$ と $\tilde b$ というのを並べてあげて、ブロックの大きい方程式を考えます。これは対称行列になりますので、こういう対称化された方程式にCG法とCR法を形式的に適用するとBi-CG法とBi-CR法の算法を導出することができます。ただし、ここに前処理付きと書いています。前処理についてはほとんど説明していませんけれども、何か方程式を事前に解きやすい形に変換をしてから解法を適用するのを前処理と言います。この場合でいうと、このような前処理行列を使って方程式を変換してから、CG法やCR法を適用すると、Bi-CG法とBi-CR法が出てくるということになります。
　これは、方程式が $Ax=b$ と $A^\top \tilde x =\tilde b$ という二つが並んでいますので、これら二つの方程式を同時に解いているようなイメージになります。欲しい解はこちらの一個だけなんですが、算法自体が両方を同時に解いているような形式になります。

A：その対称化した方程式に前処理をつけると非対称になるんですか？
相原：そうですね。
A：非対称な系に対して、CGもしくはCRを適用していると考えるんですか？
相原：そうとも考えられるかもしれませんが、両側の前処理を考えれば対称な系に適用していることになります。ただ、厳密には正定値な対称行列にCG法を適用しているわけではないです。これは形式的なもので、あくまでアルゴリズムを導出する際にこういう捉え方ができるということです。
　少し補足をすると、対称化されたものに前処理付きのCG法、CR法を適用するという発想に対して、逆に対称化されていない、 $A$ と $A^\top$ を左上と右下に並べたような非対称な方程式に対して内積の入れ方を変えてあげて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に制限しなければいけなくて、それがいわゆる「丸める」という操作になります。近似値しか扱えないので、ここで丸め誤差が生じます。
　こういう浮動小数点数を使うと、当然近似計算で、どれくらいの精度を持っているのか？という話になります。計算精度としては、 $z$ を実数とすると、このflというのはそれを浮動小数点数に丸めた値でして、そのときの精度としては、 $\delta_1$ というのが相対誤差になっていて、その大きさというのが ${\bf u}$ で抑えられる。この ${\bf u}$ というのが相対精度(unit roundoff)と呼ばれるものです。ここの仮数部の値に依存しますが、倍精度だったら $2^{-53}$ で、10進でいうと16桁くらいの有効桁を持っています。
　あとは、 $x,y$ というのを倍精度浮動小数点数とすると、その二つの演算についてです。この op はオペレーションの意味で、四則演算等を行うと、そこでも誤差が発生して、それもunit roundoffで抑えられるぐらいの相対誤差 $\delta_2$ が混入します。なので、実数が与えられたらそれを浮動小数点数に丸めたときにも誤差が入るし、それらを演算したときにも誤差が入るので、結局は理論通りの収束性は維持できないということなります。

　反復法の話に戻りまして、その丸め誤差という観点から理想的な反復法とはどういうものか、という話をさせていただきます。キーワードとして、収束速度と近似解精度という二つを取り上げます。
　まず収束速度ですけれども、これは数値解析の世界でいういわゆる何次収束するか、という話ではありません。残差は高々 $n$ 反復で本来 $0$ になるはずのものが、丸め誤差の影響によって、残差が減少する速さが低下してしまったり、ひどいときには残差が停滞したり、あるいは発散して収束しないということが起きますので、そういう意味で収束速度という言葉を使っています。つまり、あまり慣れないかもしれませんが、理想的という意味では、残差が丸め誤差の無い世界での理論的な減少の速さを維持できるアルゴリズムが望ましいということになります。
　近似解精度ですが、これは後でお話ししますので、先に数値例を一つ紹介したいと思います。

　よく学生とクリロフ部分空間法の話をしていると、手計算で簡単に解を求められるのに、クリロフ部分空間法を適用すると解けないような問題はあるか…という話題が出てきます。そんなのは、専門家の方はいくらでもご存じかとは思いますが、その最たる例として、こういうcyclicな行列というのが挙げられます。これは右上に成分が一個あって、下副対角にも何か成分が入っていて、ほかは全部0です。これの固有値分布というのは、複素平面上でこんなふうに円の形になります。これは簡単な計算で確認することができます。
　こういう問題に、ここではBi-CG法を適用してみたときの収束履歴を見ていきます。このグラフは相対残差ノルムの収束履歴をあらわしていまして、横軸が反復回数、縦軸が相対残差の2-ノルムを対数でプロットしたものです。なお、この計算はまだ倍精度ではなくて、約８倍精度の高精度な計算を行っています。だから、これはほとんど誤差が入らない、理論的な性質をある程度維持していると考えられます。実際に、 $n=100$ の問題で残差ノルムは全然減らないんですけれども、100反復目ですとんと落ちる。つまり、先ほど言いました理論がしっかり成り立っている。こちらは $n=1000$ のときです。1000反復目までは全然残差は減少しなくて、すとんと落ちる、こういう性質が見えます。これは多倍長を使っているので数学的な性質が確認できるんですが、では倍精度で計算するとどうなるか、ということです。重ねてみます。
　この赤いラインが倍精度で計算したものです。まず $n=100$ では、初めのうちは重なっていて、 $n$ 反復目に落ちるには落ちるんですけれども、大体 $10^{-6}$ ぐらいで止まっているわけです。先ほど有効桁は大体16桁と言いましたけれども、これは解の精度として残差の意味で6桁ぐらいしか出ていないので、これが丸め誤差の影響ということになります。また、 $n=1000$ のときには、もはや残差が全然減りません。1000反復を超えても全然減らないということで、かなり丸め誤差の影響を受けていることが分かります。
　この行列は各行に一つしか成分がないので、簡単に割っていけば厳密な解というのがわかるんですが、こういうものにクリロフ部分空間法を適用すると、あっという間に収束しないような場合が出てくるということです。

B：これは $A$ をどんどん掛けていったら非ゼロ成分が動いていく形のようですが。 $A$ を2乗したら並んで下に下がる。
相原：そうです。
B：それでもう一回掛けるとまた下がるという。
相原：そうです。くるくる回るということです。
B:クリロフ部分空間法が収束しようと思ったら、解を $A^{-1}$ を加工して表現しようとなると、最後の対角成分が埋まるころには…、最後の1反復までいかないと線形和で書き表せないからという解釈が成り立つんですか？
相原：ここの残差が減っていないということですか？
B：最後までいかないと、非ゼロ成分の埋まり方だけで理解できますか？これは $A^{-1}$ というのはどういう形なんですか、簡単に書ける？
相原： $A^{-1}$ は簡単に書けると思いますけれども…、残差が減らない要因としては、固有値分布が円になって、今は右辺項を書いていないですけど乱数でつくっているんですが、右辺ベクトルを固有ベクトル展開したときに、全ての固有値の固有ベクトル成分が均等に入ってしまっています。そういうときに、クリロフ部分空間法は右辺項に対して何か多項式を掛けて結局それを0にしたいわけなんですけれども、いわゆる残差多項式と呼ばれるものですね。その残差多項式のゼロ点が固有値になるようにつくりたいわけですが、その固有値のどこを取ったとしても、どこかしらには大きい成分が残ってしまっていて…。
B:そうかそういうふうに理解するのか。
相原：全てのゼロ点がそういう何か固有値になるまで反復をしないと残差が結局減らないというふうに解釈ができると思います。
B：初期値によらずに、やはりどうしてもこういう形になってしまいますか？
相原：そうですね。右辺には結構寄ります。今右辺項を乱数で取って、結局各固有ベクトルの成分が均等にばらばらに入るように組み立てているので。右辺が特定の固有ベクトルだけで張られているような場合だと、もちろんそこで落ちます。

　極端な例ではありましたが、収束速度の話をさせていただきました。通常はある程度滑らかに減少していくのに、それが滑らかに落ちずに反復回数が増えてしまうとか、そういうものも収束速度の低下ということになります。

　一方で、近似解精度に関してですけれども、精度というと当然誤差をはかりたくなるわけでして、誤差は真の解 $x^*$ と $k$ 回目の反復解 $x_k$ との差です。ただ、一般にはこの真の解が分からないというのが前提ですので、誤差をダイレクトに求めるというのは当然できないわけです。なので、反復計算をするときには、先ほどから出ている残差というのを使うしかないのです。精度保証とかそういったことも考えられるかもしれませんが、ここではもう少し一般的な話をさせていただきます。この $b-Ax_k$ は評価できますので、これが小さくなればそれだけ解に近づいている、というふうに考えられますが、数値解析の教科書などにもよく書いてあるように、残差と誤差の間には条件数といったものがありますので、条件数が大きくなってくると、ここの間には差異が生じるので気をつけましょうということになります。
　一方、この後にお話しさせていただきたいのは、残差は見ているんですが、この $b-Ax_k$ の値というのは、一般的には反復の中では計算されません。どういうことかというと、これは「真の残差」というふうに呼んだりするんですが、計算の中では先ほどの漸化式を使いますので、「漸化式から求まる残差」しか計算されないわけです。そうすると、この残差 $r_k$ と $b-Ax_k$ というのは、本来は等価になるようにアルゴリズムが組み立てられているわけなんですけれども、丸め誤差のある世界では残差 $r_k$ と $b-Ax_k$ の間に乖離が生じてしまいます。そのようなものをresidual gapというふうに我々は呼んでいます。なので、ここにも差が生じてしまうので残差 $r_k$ と誤差はさらに離れてしまうということになります。
　この後residual gapについて、もう少し詳しくお話ししたいと思いますが、理想的なという意味では、きちんと要求精度を満たす解法が望ましいということになります。

　ではresidual gapは何が問題なのかという話なんですが、これは収束判定について考えるとよくわかります。一般的な収束判定はCG法などを思い浮かべていただくといいんですが、漸化式から求まる残差 $r_k$ のノルムが何か適当なtoleranceより小さくなったら止める、これが標準的なものです。ここで疑問が生じます。反復が終了したとき、つまりこの関係が満たされたときに、真の残差 $b-Ax_k$ に対してこの同じ関係が成り立ちますか、ということです。これは $r_k$ とは本来は等価なものなので、無限精度、つまり丸め誤差の無い世界だったら当然Yesということになるんですが、丸め誤差のある世界ではNoとなり、これを満たさないことがあります。どういうことかというと、簡単な不等式で書くことができまして、この真の残差ノルムの値を評価しようと思ったら、そこに漸化式から求まる残差 $r_k$ を引いて足してあげて、上から抑えるとこんな形になります。ここの項がいわゆるresidual gapと呼ばれるもので、漸化式から求まる残差 $r_k$ と真の残差 $b-Ax_k$ との差分です。本来は、この項は0なんですが、丸め誤差があると0ではなくなります。なので、漸化式から求まる残差 $r_k$ がいくら小さくなったとしても、この差分が大きくなってくると真の残差ノルムの上限も大きくなってしまい、真の残差ノルムが適当なレベルで停滞してしまうことの要因の一つです。なので、結論としては、真の残差の意味で良い精度の解が欲しいと思ったら、residual gapをまずは小さくすることが大事で、収束判定などにも絡んできます。

　これも数値例を紹介したいと思います。行列の話は後ほどしますので、あるテスト行列に対して今度はBi-CR法を適用してみます。すると、こんなふうに残差ノルムの振る舞いが出てきます。横軸が反復回数、縦軸が相対残差ノルムで収束判定が相対残差の意味で $10^{-10}$ 未満としています。この青いラインは漸化式から求まる残差で、振動していますけれども、減少していって4403回目で $10^{-10}$ に到達したのでこのくらいの値で収束、反復が停止したわけです。ここに真の残差ノルムを毎反復計算したものを重ねると次のようになります。この赤いラインが真の残差ノルムの振る舞いで、最初は一致しているんですが、最後の方を見ていただくと、ここで停滞してしまっています。結果的に、この反復が停止したところでは、大体 $10^{-7}$ くらいで止まっているので、4桁くらいここには乖離が生じています。なので、反復が停止して良い精度の近似解が得られたなと思うかもしれませんが、真の残差を評価すると実際には解の精度が劣化している、そういうことがよく起きます。この後の話はこれをどうやって改善していったらいいか、ということになります。

　収束性の改善については、いろいろな手法があります。あまり一つ一つ細かくはお話しできないんですが、素朴な発想を幾つかお話しさせていただくと、まず一つは多倍長計算の利用です。先ほど擬似的な8倍精度で計算したものは、うまくいっていましたので、多倍長を使えばいいでしょうという発想があります。実際に多倍長計算を使うと収束速度、residual gapともに改善されます。ただし、普通の倍精度の演算に比べると、かなり計算は重たくなります。また、コスト面を考慮して部分的に利用する混合精度演算ではどうかというと、経験的に部分的な利用では改善は難しいということがわかっています。ただし、我々は多倍長というのは収束性を改善するために直接的には使わないんですが、誤差解析の際に有効活用することができます。というのは、先ほどのように多倍長で計算したものと倍精度で計算したものを比較してみると、丸め誤差の影響というのを可視化することができるので、どのくらいの誤差が混入しているか、あるいは混合精度演算とかも使って、どういうところで誤差が大きく拡大しているかとか、そういうのを調べるのにはすごく役立てることができます。そういう意味で多倍長を研究の中で使ったりしています。
　あとリスタートです。リスタートは、適当なところで反復を打ち切ってそこで得られた近似解を新しい初期値として再出発する手法になります。リスタート手法に関しては、丸め誤差の影響を排除することができます。そこまでに積もった丸め誤差というものを完全に捨てる事になりますので有効ですが、一般的にはリスタートするまでにつくったクリロフ部分空間の基底というものも捨ててしまうことになるので、そこから再出発した場合には、反復を続けた場合に比べて反復回数は増加してしまうということがあります。
　あと真の残差の計算です。これはresidual gapというのを改善するのだったら、漸化式を使わないでダイレクトに真の残差を計算すればいいじゃないか、という考え方です。これは少し古い教科書などを見るとそういう手法が実際に書かれたりしていますが、現在はいろいろ解析されていまして、真の残差を利用するとresidual gapは改善できるんだけれども、先にお話ししたような収束速度の低下に繋がってしまうということが知られています。これは簡単に実験的にも確認することができるんですが、真の残差をそのまま使うのは良くないというふうに今は考えられています。
　では、これらはそれぞれメリット・デメリットがあるんですけれども、今日はもう少し思考を凝らした「ベクトル更新のグループ化」というものについてお話しできればと思います。ここから込み入った話にはなるんですが、丸め誤差の影響を抑制するように漸化式を修正して、さらに低コストで収束速度を維持しつつresidual gapも改善したい、そういうものを目標につくられた手法です。

　では、ここから収束性の改善の話をさせていただきます。普段、学会等で講演する場合にはここからお話をするので、ここにいらっしゃる方にとって前半はかなり前振りが長いなと思われたかもしれませんが、ここからがより専門的な内容になります。

　まずは、収束性の改善のために丸め誤差の影響というものを考えなければいけないので、誤差評価について紹介させていただきます。ここでは、浮動小数点演算で行う行列ベクトル積について考えていきたいと思います。
　行列ベクトル積 $Ap$ というものを行うと、それは真の行列ベクトル積の結果 $Ap$ と誤差項 $\Delta p$ で表現することができます。ただし、この $\Delta$ というのは $n\times n$ の行列でして、このような不等式を満たすものです。あとは書かれているとおりです。これ以降は、こちらのSleijpen先生、van der Vorst先生の論文に倣って、ほぼここに書かれている内容を紹介させていただきたいと思います。また、行列ベクトル積以外の演算については、丸め誤差の影響は無視して考えていきます。内積演算とかスカラー積、そういうものは一切無視して考えていますのでご了承ください。

　考えるのは近似解と残差を更新するこういった漸化式です。本来はここに係数 $\alpha_k$ が入っていたんですけれども、入っても入らなくてもそれほど議論が変わらないので簡単のためにここでは省略させていただきます。ここに行列ベクトル積 $Ap_k$ というのがあるので、ここで混入する誤差を評価します。 $Ap_k$ を浮動小数点演算で計算すると誤差が入ってくるので、この $\Delta _k p_k$ というのがどれくらいのサイズなのか、この一回の局所的な行列ベクトル積の誤差の上限というのを評価すると、こんなふうに書けます。上限は残差ノルムで書けていまして、 $r_k$ と $r_{k+1}$ のノルムの和で表せて、あとは定数倍です。ただし、 $\kappa$ はこういう定義で、あと ${\bf u}^2$ の項は無視しています。
　この不等式の意味するところは何かというと、この漸化式で $r_k - Ap_k$ というのを計算して結果として出てくるのは $r_{k+1}$ なので、結果として出てくる残差 $r_{k+1}$ にどれくらいの誤差が入るのかということを考えたいわけです。もし $r_k$ から $r_{k+1}$ に急激に減少して $r_{k+1}$ のノルムが非常に小さくなるといったことが起きると、ここの誤差項というのは相対的に大きい $r_k$ の値に依存して大きくなる場合がありますので、結果として出てくる $r_{k+1}$ は相対的に大きな誤差を含んでしまうことが起こり得ます。わかりにくいかもしれませんが、 $r_{k+1}$ は小さくて、 $r_k$ が大きい場合、この $Ap_k$ の計算で発生する誤差項はこの大きい $r_k$ に依存するので、 $r_{k+1}$ から見ると相対的に大きな誤差になるわけです。
　また少し見方を変えてみると、こういう残差の急激な減少というのは、 $r_k -Ap_k$ を計算して非常に小さなものができるということなので、ここの引き算で何桁か、大幅な桁落ち等が起きるというふうにも考えることができます。残差に相対的に大きな誤差が入ってくると、直交性が崩れたりとかそういうことが起きますので、収束速度の低下にもつながると考えられます。

A：局所的な誤差の上限のところの $\Delta_k$ は $k$ に依存しないというのは省略して書いてあるんですか？一個前の上の式は $\Delta_k$ と書いてあるので。
相原：すみません。ここは抜け落ちています。 $\Delta_k$ です。ただ、 $\Delta_k$ は基本的には $k$ に依存するんですが、上限を評価するときには依存しない形になるので、あまり気にしなくても大丈夫です。
A:何か不思議ですね。ここで出ている結論が、本当は残差ノルムの急激な減少の更新が、それをすると収束速度が全体としては落ちる。
相原：そうですね。その話はこの後でもさせていただこうと思っていたんですが、一般に反復法を習うと残差を減らしましょうということで、残差が減るというのはすごく嬉しいことなんですが、局所的に見て急激に、特に残差が一旦大きくなってから急激に減少してしまったりするのは、実は丸め誤差の観点からすると望ましくないということがあります。そのあたりは反復法を専門に研究している方でないと、多分そんなに知られていないことなので、ぜひそういう認識を使っている方々に持ってもらえたらと思っています。また後でその話もします。

　では二つ目の項目なんですけれども、今の局所的な誤差というものから、今度はresidual gapを評価しましょう。Residual gapの上限としては、こんな形で与えることができます。細かい部分は置いておいて、大事な部分は何かというと残差ノルムの和です。なので、残差ノルムが振動して初期残差ノルムというのを一つの基準として考えると、それよりも非常に大きくなるような中間的な残差ノルムというものがあると、この上限というのが相対的に大きくなり得ます。これによってresidual gapが大きくなるので、近似解精度の劣化に繋がってしまうということになります。
　結論として何が言いたいのかというと、相対的に大きな局所誤差の蓄積を避けて、さらにresidual gapを小さくするには、残差ノルムというのは滑らかに減少していってくれることが理想的な振る舞いということになります。

　あまりこういう話というのは普段学会とかでもほとんどないと思うんですけれども、振動のパターンというのを四通りで書いてみて、それと収束性との関係というのを、現実にはこんなに単純ではないんですが、あくまでイメージとして持っていただくといいかなということで書かせていただきました。左上が一番理想的なパターンで右下が最悪なパターンです。振動の大きさとノルムが初期残差に比べてどれくらい大きくなるかというのを表現しています。
　左上の理想的なパターンは振動が小さくて、 $10^{0}$ が初期相対残差ノルムの値だとすると、初期残差に比べて大きくならない、つまり振動はしつつも全体としては小さな振動で滑らかに減少してくれるものです。これが収束速度の意味でもresidual gapの意味でも良好という結果になります。
　一方で右上は、初期残差に比べて大きなノルムはないけれども、一回一回の振動が激しい場合です。これは局所的に大きな誤差が積み重なってしまうので、速度に結構影響します。ただ初期残差よりは小さいので、residual gapに関してはきっと大丈夫でしょう。
　では左下はというと、振動自体は小さいです。振動は小さいけれども初期残差に比べて小さな振動をしながらノルムが一度大きくなってから収束していくものです。これは速度に関しては問題ないだろうと思われますが、residual gapが大きく出てしまうと考えられます。
　最後に、これは2桁で基準を取っていますが、振動も大きいしノルムも大きい、という最悪なパターンです。どうしようもないということになります。
　もちろん現実はこんな単純には分類できないんですが、イメージとしてこういうようなものを考えてみると、収束振る舞いと収束性との関係が見やすいかと思います。ようするに左上のものが理想的なわけです。残差が減れば嬉しいのだけれども、ただ急激に減ればいいというものではないという、先ほどのお話になります。

B：左下で言うと、最後のresidual gapは3回目から4回目の急激な減少が最後residual gapにあらわれると？
相原：Residual gapに影響するのは最大のピークです。これで言うとピークが4乗くらいですので。
B：その差がresidual gapになるんですね。
相原：はい。だから４桁くらい、一番上のピークが $10^{4}$ くらいまでいくと、最後は例えば極端な話 $10^{-12}$ とか、4桁上がったら4桁くらいは劣化してしまうようなイメージを持っていただくとよいと思います。
B：さっきの解説だと、 $r_{k+1}$ と $r_k$ の差ではなかったか、初期値との差？
相原：Residual gapに関しては、グローバルなものですので、上限を相対評価するときには初期残差で割って評価をするんですね。だから初期残差に比べて、0番から $k$ 番目の間で残差ノルムがどれくらい大きくなるかというところを見ます。一方、局所的な方は、これは全体ではなくて本当にごく一部分で、 $k$ 番目から $k+1$ 番目でどれくらい減少がありますかという話です。例えば、3～4桁残差ノルムが減って嬉しいじゃないかと思うかもしれませんが、3～4桁減っているというのは、その分だけ有効桁が減っているようなイメージを持っていただくとよいかと思います。

　時間もちょっと押してきました。あと、最後のほうは数値実験をざっとお見せするだけなので、改善手法について少しお話をさせていただければと思います。
　残差更新のグループ化というアイデアがあります。これは私のオリジナルではなくて、1994年のNeumaier先生と96年のSleijpen先生、van der Vorst先生の論文で示されています。
　近似解と残差の漸化式を考えたときに、ある単調増加数列 $\pi (j)$ というものを考えて、ある $k$ 番目、適当な番号のときに近似解や残差というのをこういう式で置き換えるという手法です。置き換えるというのは、再計算するという意味になります。ただし、ベクトル $y_j$ というのは、探索方向ベクトルのいくつかを足し合わせて、グルーピングしたものになります。この式で言うとイメージしにくいかもしれませんが、こちらの図をご覧ください。
　元々の列として、 $x_0, x_1, x_2,...$ とこんなふうに反復列があったとすると、 $\pi (j)$ というのは、0番は一致させて、置き換えるタイミングを $\pi (1)$ は2番、 $\pi (2)$ は5番というふうに適当に選んであげるわけです。そして、そこまでの探索方向ベクトルの和というのをグルーピングします。このグループ化したものに対して $Ay_j$ という行列ベクトル積を行って、残差を飛び飛びに計算するような式を立ててあげます。これで言うと2番目で計算し直して5番目で計算をし直して、ということを行います。そうすると何が変わるのかというと、もともとこの $Ap_k$ という行列ベクトル積を行っていたのに対して、こちらの置き換えたほうの $r_{\pi (j)}$ の列というのは、行列ベクトル積を行うのはこの $Ay_j$ ということになるので、ここで発生する丸め誤差の影響が $Ap_k$ で発生する誤差の影響よりも少なくなるように、この数列 $\pi (j)$ を選ぶことができれば、誤差の影響が少なくなるのではないか、という発想です。

B：ここで飛ばした $x_1$ は求めているんですね。
相原：1回は求めています。一度は求めて、何かここで大きな誤差の混入といった状況を検知したら、そこで再計算するということです。
B：アルゴリズムの順番としてはどうなっているんですか？上の列と下の列の。
相原：順番としては、 $x_0$ を与えて、 $x_{\pi (0)}$ は等価なのでただ置くだけです。 $x_1$ を計算して、 $x_2$ を計算して、何かここで検知されたら、 $x_2$ を $x_{\pi (1)}$ を計算して置き換える。それで、 $x_3$ を計算して…、という順です。

　それでは数列 $\pi (j)$ というのをどうやって選ぶかという発想についてお話しします。理想的な収束振る舞いというのは先ほどの話から、まずは大きな中間残差ノルムを回避したいという点で、これはresidual gapを小さくするためのものです。もう一つは、残差ノルムの大きな振動について、大きいところから小さいところへの急激な減少を回避したいという点です。これら二つの項目があります。
　まず、同時にはなかなか難しいので1番から考えてみますと、 $\pi (j)$ をどうやって選ぶのかというと、この青いラインがもともとの収束振る舞いだとします。 $r_0$ からスタートして $r_1, r_2, r_3, ...$ と、こんなふうになっているとすると、 $\pi (j)$ はできるだけ滑らかになるように、この大きい中間残差を回避するように選んであげればよいわけです。なので、ここで言うと2番目と5番目で、ここのベクトルを計算し直すということです。ただし、残差の再計算には余分な行列ベクトル積がかかってしまうので、あまり必要以上には計算したくないわけです。計算コストが多くかかってしまうので。なので、詳細は割愛させていただきますが、実際には適当なパラメータを導入してあげて、何か一定の条件を満たしたときだけこういう再計算をしてあげましょうという手法が開発されています。

B：再計算するノルムは $b-Ax_k$ を計算しているわけじゃないんですよね。最初のほうで $b-Ax_k$ を計算してはいけないという話があったので、その値とその $r_{\pi (1)}$ は違うはずですよね？丸め誤差が入っていると。
相原： $r_{\pi (1)}$ は初期値が0だと $r_0$ は $b$ と一致するので、ここは $b-Ax_k$ を計算しているようなものです。2番目以降は $r_{\pi (1)}$ を右辺項とするようなものに対して $b-Ax_k$ を計算しているような感じなんです。実は、これは英語で表現するときには、真の残差はtrue residualと呼ばれますけれども、この残差ベクトルはlocal true residualと呼ばれるんですね。そのローカルな意味で $b-Ax_k$ のような形になっている。
B：本当の正解は $b-Ax_k$ のはずですよね？
相原：そうですね。
B:ただ、そうじゃなくて、ある意味、解ベクトルとの関係をうまく持たせて、あえて摂動させて何か全体としてうまく収束させようとしているのか。
相原：これも一応 $b-Ax_k$ にはなっています。
B:丸め誤差が入ると違うんじゃないんですか。質問の仕方を変えると、それを $b-Ax_k$ で計算しても同じ効果があるんですか？
相原：丸め誤差が入るので異なりますが、この場合は実は $b-Ax_k$ で計算しても効果はあります。
B:じゃ、一緒になる？
相原：はい。
B：最初のほうで $b-Ax_k$ を計算すると収束速度が劣化するという話との関係が…。
相原：その話との関係は、今は $b-Ax_k$ を計算したときにどれくらいの誤差が入るのかという議論をしていないんですが、 $b-Ax_k$ を計算するとそこに入ってくる誤差項といのうは、 $b$ の大きさに依存してくるんです。なので、実はその真の残差の計算は残差ノルムが大きい間は計算しても大丈夫なんです。それは大きさ的に残差ノルムがまだ大きい間に真の残差を計算しても、そこで入ってくる誤差というのは相対的には大したことがない。でも残差ノルムが小さくなってきたときに、真の残差を計算すると、そこの大きさが $b$ の大きさに依ってくるので。
B:むしろそっちのほうが誤差を含んでいる。
相原：はい、含んでいます。なので、後半で真の残差を計算するとまずいというのが正確な表現になります。
B:わかりました。ありがとうございます。
相原：今日お話した先ほどのSleijpen先生の論文とは違うんですが、van der Vorst先生がそういう最初の段階だけ真の残差をうまく計算してあげて、置き換えをうまくやればresidual gapを防げますよといった手法も別の論文として出されていて、後半さえ使わなければうまくはいきます。パラメータを用いた適当な条件というのも、実は最初の段階だけ部分的にこういう計算をしましょうという発想になっています。

A：今、 $\alpha_k$ は省略しているというのはまだ続いていて、 $p_k$ の線形結合をただ単に足し算で書かれているのは、 $\alpha_k$ がもしかしたら隠れているんですか？
相原：はい、 $\alpha_k$ が隠れています。 $\alpha_0, \alpha_1$ がここに隠れています。
A: $\pi (j)$ をうまく使って適当な頻度でうまい置き換えをしてやるというのは、そんなにあまりインターバルは大きくないように取っているので、あまりたくさん $p_k$ を取っておく必要はない？
相原：ここの数ですか。
B:はい。それはその置き換えをするまでは、1回置き換えをして次の置き換えをするまでの間には、ずっとその $p_k$ をとっておかなければいけないんですか？
相原：とっておかないといけないですね。ただ、1反復するごとにこの $p_k$ に足し込んでいくので、メモリ自体はそんなに使わないです。全部保持いくわけでなくて1反復する度に足し込んでいくということなので、数本分だけ余分なメモリが必要になります。

　2番目が残差ノルムの大きな振動、急激な減少を回避するにはどうしたらいいかという話なんですが、これは何かというと、先ほどは大きな中間残差を避けるように $\pi (j)$ を選びましょうと言いましたが、ご質問があったように、この間は計算していないのかというと、そんなことはなくて実際には計算しているわけです、一度は。計算しているので、何かここで、例えば $r_3$ から $r_4$ に向かってこのように大きな減少があると、ここで発生する局所誤差は相対的に大きく、これが混入してしまっていると収束速度に影響するのではないかと思われます。実際その通りで、こういう場合にも何か対処してあげないといけません。
　そこで、大きな振動がある場合には、例えば $r_{\pi (1)}$ から $r_{\pi (2)}$ の間で急激な減少が起きた場合は、残差を再計算してあげましょうという発想が出てきます。このイメージとしては、これから $r_3, r_4$ を計算して、ここで振動が大きいということを検知したら、この $r_4$ というのを捨てて、この形式で計算し直してあげるのです。これは $r_{\pi (1)}$ をベースにしてグルーピンクしたものを使って $A$ を掛けて再計算するわけです。そうすると、 $r_3$ から $r_4$ を作った経路ではなくて、 $r_{\pi (1)}$ から $r_4$ に飛んで計算したようなイメージになっている。そうすると、小さいものから大きいものをつくるときには局所的な誤差は大きくないので、こっちからの経路だったら問題ないだろうということになります。
　ただし、先ほどと同じようにここでの再計算にもやはり行列ベクトル積が余分にかかってしまうので、頻繁にはできないわけです。だからこれも適当なパラメータを入れてあげて、本当に振動が大きい場合だけとか、部分的にこういう処理を入れてあげることで、strategicalな手法ではあるんですが、効率的に丸め誤差の影響を抑制できるということになります。

　では、ここまでがメインの話になりますけれども、以上の話はCG法などの漸化式をメインに考えてきましたが、冒頭でCR法あるいはBi-CR法は少し漸化式の形式が異なりますという話をしました。そういう場合にも適用できるかというと、それは適用できます。実際に効果があることを示した論文が最近公開されまして、私がやった内容になっています。

発想自体は全く同じなんですけれども、概要だけお話させていただきます。Bi-CG系統とBi-CR系統で漸化式はどう違うかというと、Bi-CR系統の場合には補助ベクトル $q_k$ というのを導入してあげて、これを漸化式で作りましょうという考え方を示しました。この場合には、先ほど $Ap_k$ という部分に行列ベクトル積を使っていたんですが、今度は残差に $A$ を掛けるような計算になります。この行列ベクトル積を行う場所が異なれば、当然入ってくる丸め誤差も変わってくるので、あらためてこの形式に対して局所的な誤差を解析したりresidual gapの大きさを解析したりして考え直してあげないといけないということです。ただ、手法自体は同様の手順で考えることができます。

　漸化式による違いなんですが、比較のために並べてみますと、Bi-CG系統の漸化式では先ほどお見せしたようなこんな感じです。局所誤差は隣り合う残差ノルムの和で書けるし、residual gapも残差ノルムで書けます。Bi-CR系統の漸化式の場合には、局所誤差は新たに導入された補助ベクトル $q_k$ のノルムの和で書くことができて、residual gapに関してはちょっと異なるんですが、残差に $A$ を掛けたもので評価することができます。なので、先ほどのグループ化というのを行うに当たって、何をすればいいかというと、今度は基準を残差ではなくてこの補助ベクトルで考えてあげればいいわけです。つまり、補助ベクトルの振動を回避するようにうまくグルーピングをしてあげれば、Bi-CR系統でも効果が得られるはず、というのが私の論文の中でやった内容です。

A:Bi-CGのresidual gapに比べてBi-CRのresidual gapは随分大きくなりそうな…。
相原：そうですね。
A: $k-j$ があるので、大きな値になりそう。
相原：そうなんです。これが実際にシャープかどうかというのはなかなか議論が難しくて、結構over estimateかもとは思っています。ただし、これを元に同じようにグループ化というアイデアを入れてみると、実験的には結構効果があります。この評価自体はそれほどシャープではないかもしれないですが。

　では、最後に数値実験なんですが、数値実験のほうは私自身がBi-CR系でそういうことをやったので、実験はBi-CR法でお見せしたいと思います。実際に効果があるということを示します。

　計算条件はこのようになっていて、初期値は0で収束判定条件は $10^{-10}$ というのを基準にしています。従来のBi-CR法とその中でベクトル更新のグループ化を行ったBi-CR法の収束性を比較しています。

　テスト行列は、これは最初にお見せした行列でして、フロリダのSparse Matrix Collectionからこの行列をテストとして持ってきました。詳細はここに書いてあるとおりで、右辺は乱数ベクトルで与えています。

　最初にお見せしたように、まず従来のBi-CR法の収束履歴ですが、大きく振動してしまっていて、最終的にここでギャップが生じていて4桁くらいは差が出ているという状態のものです。

　これに対して、中でグループ化というものを行って計算する提案手法がこちらです。今、横軸が反復回数の $k$ で毎反復プロットしているので、この残差ノルムの振る舞いとしては滑らかになっているようには見えないんですけれども、この内部ではグループ化して、ある程度滑らかな収束になるように調整されています。赤いラインが真の残差ノルムで青いラインが漸化式から求まる残差ノルムなので、先ほどと見比べてみると、きれいに一致していて最終的に10桁まではちゃんと出ています。また、収束の速さはほとんど変わっていないです。数学的には何も変えていないので、数学的には等価なものですが、実装するときに先ほどのグループ化を行うことで、丸め誤差の影響が変わっているということです。

A:上にそのノルムが跳ね上がった分、下の限界精度の出方みたいな、例えばこの場合で言えば $10^{6}$ だけ大きくなったとしたら、今出ている $10^{-10}$ 程度しか精度が出ないことが期待されるんですけど、もっと精度を出すということはできてはいないんですか？今ここでは $10^{-10}$ しか出ていないですけれども。
相原：ここですか。そうですね。これよりもっと下げようと思うと…。
A：どこかの頭を叩かなければいけないですか。
相原：頭を叩かないといけないですね。おっしゃる通りです。そういう意味では、そもそも漸化式から求まる残差は、この先は計算していないんですけれども、これがどこまで下がるかというのもあるので、これは多分下がっていくので、真の残差の計算、まさに $b-Ax_k$ の計算というのをstrategicalにここの中に入れてあげるとひょっとしたら改善するかもしれないです。ただ、そういうのは、どのタイミングでとか、結構ややこしい話にはなってくるので、汎用的かどうかはわからないですけれども。

　今、これはギャップが改善できましたという話なんですが、もうちょっと難しい問題にします。

　これはよく反復法のテストに使われる問題で、こういう偏微分方程式の離散化から出てくるような問題を対象に解いてみました。解法は同じBi-CR法です。

　これは従来の方法とそのグルーピングを行った方法を一緒に重ねてしまっていますが、赤いラインが従来法で青いラインが提案法です。振動も大きくてノルムも大きいので、さっきの4つのパターンでいうと最悪なパターンになります。
　従来の方法は $n$ 反復しても全然収束していません。少し減っているようにも見えるんですが、この後も収束しないというような振る舞いです。それに対して提案手法は、これは内部でグループ化を行っているんですが、きれいに収束していて、反復終了時点での真の残差ノルムの値が $10^{-11}$ くらいまでしっかり出ているということになります。
　この表なんですけれども、反復回数は4282回で、ここのAdd. MVsというのは、先ほどグループ化をして残差を再計算するときに余分な行列ベクトル積がかかりますと言いましたが、それにどれくらいかかったかという回数です。175回かかっています。この場合はそもそも収束していないので、収束するようになっただけで十分かと思うんですが、Bi-CR法だと1反復当たり2回行列ベクトル積がかかりますので、全体として行列ベクトル積は8500回くらいなんです。行列ベクトル積8500回に対して、その余分なコストという意味では大体175回の行列ベクトル積なので、ほとんど大したことないです。微々たる追加の計算で、これくらいの収束改善が行えるということなので、これはやらないよりは標準的に組み込んでおいたほうがきっと良いでしょうと思います。

　最後にまとめになりますけれども、本講演ではクリロフ部分空間法の中でも短い漸化式を用いる解法において、その丸め誤差の影響を議論してきました。基本的には残差ノルムの振動、Bi-CR系だったらその補助ベクトルの振動になりますけれども、そういったものは収束速度や近似解精度に大きく影響しますということです。その丸め誤差の蓄積、拡大を回避するには滑らかな減少というのが大事です。途中でも議論が出ましたが、残差というのはただ減ればいいというものではなくて、急激に減っても困るし、少しずつ滑らかに減少していってくれるのが丸め誤差の観点からすると理想的な振る舞いになります。それを達成するために残差更新のグループ化というものについて紹介をさせていただきました。
　今後の課題というか展望なんですけれども、今現在、修士課程の学生と共同研究をやっています。残差ノルムの振る舞いを滑らかにするには、古典的なスムージングという技法があります。今日は紹介できなかったんですが、残差ノルムを滑らかにするという意味では有用なんですが、今までは丸め誤差の影響を抑えるというふうには使われていなかったんです。このスムージングを少し工夫してあげると、数値的に安定な、その丸め誤差の影響を抑制するような実装法というのができそうだということで検討しています。既に講演を聞いている方もいらっしゃるかもしれませんが、学生と一緒にやっています。これができると何が嬉しいのかというと、同じような手法ではあるんですが、こういうグループ化といった手法に比べると、残差ノルムは単調減少するし、実装は簡単だし、パラメータなどの設定も不要だし、余分な計算コスト、行列ベクトル積もかかりません。そういう利点がいろいろあるので、スムージングでできるんだったらそのほうが嬉しいということで、このグループ化を超えるような手法ができないかなと考えています。
　以上です。ありがとうございました。

質疑

司会：どうもありがとうございました。既にいろいろお話しいただいていると思いますが、何か議論とかコメントがあれば。
A：こういった研究でつくられたコードというのは公開されているんですか。
相原：よくそういった質問を、応用分野の方に言われるんですが、すみません公開はしていません。公開できればと思うんですけれども、公開できるほどの腕がないというか、いろいろ研究がまだまだ途中というのもあるんですけれども、いずれそういう公開とかもできたらと思っています。
司会：ほかにありませんか。いいですか、質問しなくて。

クロージング

司会：では一旦終わりにして、また何か話すことがあれば、また後でお話ししましょう。どうもありがとうございました。

プロフィール

2014年3月東京理科大学大学院理学研究科博士後期課程修了．博士(理学）．東京理科大学理学部第一部数理情報科学科助教を経て，現在，東京都市大学知識工学部情報科学科教育講師．数値線形代数，特に大規模な線形方程式に対する反復法の研究に従事．日本応用数理学会，日本オペレーションズ・リサーチ学会，各会員．