実用段階に到達した精度保証付き数値計算

大石 進一 Shin’ichi Oishi

2000 年 9 月 8 日

概要

本稿では数値線形代数の分野に属する問題については、その数値解の精度保証が高速かつ容易に行えるようになりつつあることを、連立一次方程式や行列の固有値の精度保証付き数値計算法を例にして示す。

1 はじめに

あけましておめでとうございます。 21 世紀最初の年始に当たり、今後非線形理論とその応用の研究がどのように発展するか展望を述べよとのご下命を頂きました。しかし、非線形理論の広がりは実に大きく、とても全分野をサーベイできませんので、私の研究分野を中心にその夢を語らせて頂くことでお許し願いたいと思います。

1.1 科学技術の起源

21 世紀を迎えて、今後どのように学問が発展し、人類がどのように生活していくかを考えるには、過去を振り返るのがまずは基本であると思います。私事ですが、最近、微分積分学の本を執筆する機会がありました。その執筆の過程でアルキメデスからニュートンまでの 2000 年間を辿りながら、人類が微積分学をどのように発見し、科学文明の基礎を築いたかを少し調べることができました。

人類は 400 万年前に誕生し、 100 万年前に原人が現れ、 25 万年前に旧人が、 4 年前に新人が現れたといわれています。このような人類史の大部分は狩猟や植物の採取生活の時代です。約 1 万年前になって、食糧生産が開始され、農耕 · 畜産生活への移行という人類史の第一革命「新石器革命」が起きました。その結果、食料生産に余力が生じ、人類最古の文明である古代オリエント文明が紀元前 3200 年ごろ興りました。これは紀元前 4 世紀のアレキサンドロス大王の統一まで 3000 年間栄えました。この中で、メソポタミア文明では紀元前 3200 年ころ絵文字が発明され、楔形文字に発展し、エジプト文明では聖刻文字が発明されました。文字によって伝承される学問の興りはしたがって、この頃まで遡ると考えられます。実際、紀元前 1900 年から 1600 年頃には，古代バビロニア人がプリムトン 322 と呼ばれる粘土版に三角形の 3 辺の比が a : b : c で，これらが互いに素であるとき，次の命題が成り立つことを記しています。

「u と v を互いに素な自然数とする．この 2 数のどちらかが偶数で u > v であるとき， a = 2uv, b = u^{2} − v^{2}, c = u^{2} + v^{2} とすると， a, b, c の最大公約数は 1 で a^{2} + b^{2} = c^{2} が成立する．」

この命題は 3 辺の長さの比が 3 : 4 : 5, 5 : 12 : 13, 8 : 15 : 17 などの三角形は直角三角形であることを示しています。

c a b 図 1: 3 辺の長さの比が a : b : c の直角三角形

紀元前 2000 年頃には古代エジプトでは次のような方法で直角を作っていたといわれています．ある長さを単位として 12 の長さのロープを用意します．これを地面に長さ 4 の分だけまっすぐ這わせ，そこから長さ 3 の分を垂直成分として上に伸ばす．このとき，斜辺がピンと張るようにすれば， 3 の長さの辺と 4 の長さの辺に挟まれた角は直角となります．これは数学的な性質をもとにした建築技術の嚆矢の一つです．

こうして、古代メソポタミアや古代エジプトでは、数学はかなり高度に発展しました。しかし、証明という概念を出したのは古代ギリシア人でした。古代ギリシャ人の中には、自由を自覚し，物事を権威ではなく論理によって理解しようとした層の人がいました。これが哲学の始まりです。古代ギリシャ数学も自然哲学としてスタートしました。数学的な事実が、権威ではなく、論理によって正しさが示されるようになった訳です。これが証明の始まりです。また、この精神が現代の科学技術の思想の起源であると考えられます。古代ギリシャ数学の祖は紀元前 6 世紀に活躍したミレトス生まれのタレスであるといわれます。ミレトスは現在のトルコのエーゲ海沿いにある町です。タレスはエジプトで測量学を学び、それをギリシャに持ち帰ました。タレスの業績は、測量学のテクニックの成立理由を数学的に証明したことであるといわれています。タレスが証明を与えた数学的な命題としては、 2 辺とその挟角が等しい三角形は合同である、対頂角は等しい、 3 つの角度が等しい三角形は相似であるなど現代でも中学で幾何学として習うものが多くあります。

直角三角形の 3 辺の長さを a, b, c とすると，斜辺の長さを c として c^{2} = a^{2} + b^{2} が成り立つというのがピタゴラスの定理です。ピタゴラス (紀元前 572 年頃サモスに生まれた) がこの定理を証明したことは、記録は残されていませんが、広く古代の人々に信じられていました。ピタゴラスはピタゴラス学派を作り，自然現象の理解に整数が重要な役割を果たすとを認識していました。これは自然を理解するのに数学が有用となるという現代における視点の起こりの 1 つと言えます。

紀元前 323 年にアレキサンドロス大王が急死し、エジプトはプトレマイオス 1 世によって統治されました。彼は、アレキサンドリアを首都とし、そこに研究機関としてムセイオンを設立し、図書館を付属施設として建設しました。そして、学者を集めて研究を支えました。ユークリッドもそのような一人です。古代ギリシャ数学は幾何学であるといってよいのです、ユークリッドはそのような数学を体系化して紀元前 3 世紀に「ストイケイア」 13 巻をまとめました。いわゆる「ユークリッド原論」です。ユークリッド原論は幾つかの定義をもとに幾何学の命題が厳密な推論によって証明されるというスタイルを取っており，現代数学のスタイルのもとともなっています。ニュートンのプリンキピアも原論のスタイルを模して書かれています。

アルキメデスは古代ギリシャのシシリー島シラクサで紀元前 287 年から 212 年に活躍しました。少年時代，天文学者の父から教育を受け，エジプトのアレキサンドリアの大図書館でエウドクソス (紀元前 406-355?，後述) やユークリッドの著作を学んだものと思われています。アルキメデスは機械学に優れ，彼の発明したアルキメデススクリュウによって、当時のエジプトではナイル川の水を耕地にくみ上げるていたようです。また、数学者としても大変な天才でした。第 2 次ポエニ戦争によって、シラクサをローマ人が攻撃した際、アルキメデスは様々な武器を発明してその防衛に当たった話はプルタークをはじめとしていろいろな歴史書に書かれています。紀元前 212 年にシラクサはローマ軍の手によって陥落し，その時，数学の証明に没頭していたアルキメデスは気づかず，ローマ兵士の手によって命を落としたとされています。

エウドクソスはエジプト，ギリシャなどにおいて学問を学んだ後，アテネに滞在して数学者，天文学者，医者として活躍しました．数学における業績としては，極限の取り扱いとして「取り尽くし法」を考案したことです。アルキメデスはこの取り尽くし法によって様々な図形の面積や，立体の体積の公式に対する証明を与えた．この取り付くし法により、アルキメデスは現代の数学でも通じる厳密な論法によって図形の面積を計算する公式の証明法の枠組みを与えたといってよいと思います。実際、アルキメデスの著作はかなりの数が残されており，そこに現代数学に通じるような議論が展開されているのを見ると，非常に不思議な感覚に襲われます。そのような著作の中で特別の運命を辿ったものがあります。紀元後 1 世紀のへロンや 10 世紀のビザンチンの辞典「スダ」などにより，アルキメデスの著作として「方法」という本が存在したことが知られていましたが，発見されていませんでした。それが 1906 年デンマークの古典学者ハイベルグによって，イスタンブールの僧院において「方法」の 10 世紀の写本が発見されたのでした。この写本は羊皮紙に書かれており， 12,3 世紀頃その上に祈とう文が上書きされたため 1906 年までアルキメデスの「方法」であるとはわからなかったということです。幸いなことに，佐藤徹訳 · 解説「アルキメデス方法」 (東海大学出版会) によって，「方法」の日本語訳を優れた解説とともに読むことができます。

シラクサ陥落によってアルキメデスを死なせたローマ帝国は，紀元前 146 年にはギリシャとカルタゴを支配下におきました．そして，紀元前 30 年にはエジプトを制圧しました。政治をするのがローマ人の役割で，科学や数学は卑しいもののするものであるというローマ帝国の思想により，ギリシャ時代に栄えた科学 · 技術 · 数学は支援を失い衰退していきました。紀元前 43 年にはシーザーがアレキサンドリアを占領しましたが，その際，大図書館はほとんど焼失してしまいました。また，ローマ時代にはキリスト教が普及し，異教徒の思想として科学技術が攻撃されました。異教徒の興した数学をすることは罪にさえ問われるような時代となったのでした。

こうして，紀元 4 世紀にギリシャが滅亡すると， 9 世紀頃のアラビアまでその数学は引き継がれることなく放置されたままでした。アラビアでは紀元 570 年マホメットが生まれ， 610 年イスラム教を唱えます。わずか 150 人程の信者でスタートしたイスラム教が 630 年にはメッカを占領し， 632 年アラビア半島を統一します。その後イスラム教国ウマイヤ朝ができました。イスラム教の国王は教主 (カリフ) を兼ねることが多く、カリフは商業で得た財力をもとに，学問，文芸を振興し，数学も発展しました。アラビアの数学者はギリシャ数学やインド数学の書物を翻訳したり，アルゴリズム (algorithm) や代数 (algebra) という言葉の発祥となるような創造的な研究もしました。 15 世紀に至りヨーロッパでもルネッサンスにより学問が振興されるようになり，アラビアの数学がヨーロッパに入ります。ギリシャ数学もアラビアの学者の翻訳を通じて入ってきました。また，活版印刷が始まり，ユークリッド原論やアルキメデスの書のラテン語訳が出版されて流布されていきました。こうして，アルキメデスから実に 2000 年経過して，ようやく科学 · 数学が大きく動き始め、それが現代につながっているわけです。これ以上の詳細を議論してる余裕はありませんが、興味をもたれた読者は拙著等を参照していただければ思います。また、古代ギリシャ数学の著作のラテン語訳本の表紙などを展示するホームページがあります。大変興味深いのでぜひ参照してください。

1.2 個人史

私は高校生の頃から、難しい方程式を解く道具を作るのが夢でした。高校時代、難しい積分を計算するためには、大学の数学を先読みすればよいと考えて、本屋に行き、ルベック積分の本を買ってきて読んだことがありました。しかしルベック積分論を呼んでも、難しい関数の積分の値が計算できる訳ではないことに気づき、大変な違和感を覚えた記憶があります。これに続いて関数解析の勉強を大学生になっても続けていましたが、吉田耕作先生の位相解析の本がよくわからず、 3 年生になって、ある夜すらすらと筋が見えるような気がしたのを今でも覚えている。要するに吉田先生は自分の好きなところから程度を落とさず書いていたのだ。また、同じ頃、ファインマンの量子電磁気学に大変興味をもつ。図書館に入りびたりで、 Dirac の初期の論文や渡辺慧先生の Progress の論文なども読みふけった。卒論は量子電磁気学で書きたいと思ったが、引き受けての研究室がなさそうであると勝手にあきらめた。そして、岩波講座の現代物理学の基礎でソリトンに出会った。非線形偏微分方程式を量子力学の逆散乱理論 (線形) で解く方法があるという、何だか夢のような話が紹介されていた。夏休みに大きなサーベイ論文を読んで、その中に逆散乱法で解かれていない方程式として紹介されていた広田方程式を逆散乱法で解いて卒論にした。これを英語の論文にしようと考えたが、なんとなく自信がなかった。ただの演習問題に過ぎないという予感があった。修士に入り、日本人研究者の広田良吾先生の発明された広田の方法によるソリトンの解法が、逆散乱法よりしっくりするのでこれに夢中になる。簡単に広田の双線形方程式で書かれる方程式のある一般系には必ず 2 ソリトン解が存在することが証明できた。戸田盛和先生のセミナーで戸田先生、和達三樹先生、小寺武康先生の前で話をしたら、面白いと喜んだり、こんなことがこんなに簡単に示せるのかと驚かれたりした。それで、英語の論文にしようと思ったが、英語の論文がなかなか書けなかった。これが修士 1 年であったが、結局、博士 1 年の秋になるまで英語の論文が書けなかった。後で考えると 2 年間、英語の論文の書き方にこだわったことが、後に論文を書くことの基礎になったのでよかったと思う。随分、英文論文を書くためにという系統の本を読んだ。また、教師として、非常に短期間によい英語の論文を書く学生もいることに舌をまくことがある。こういうときは本気で褒めることにしている。この後、ソリトンの周期解を求める問題やパンレベ性など難しい問題を勉強したりしたが、ある日、広田の方法がある意味でフーリエ解析の非線形への拡張であると見なせることがわかり、一般解を構成できることが示せた。結論的には、ソリトン方程式の解がある積分方程式のフレッドホルム行列式となることを示した。これを名古屋大学のプラズマ研究所のソリトンの研究会で報告すると市川芳彦先生が今回の発表の中で一番面白かった講演であったと講評して下さった。これは、どうあっても英語の論文にしなくてはという思いで、やっと英語の論文を書き上げた。戸田先生が真っ赤になるほど原稿を修正してくださった。後で、東大天文台の吉田春夫さんに論文のリプリントを差し上げたら、これほど無駄なく書いてある論文は珍しいと言って褒めてくださった。戸田先生の添削の賜物であると思った。これが博士 1 年の秋であった。この論文では 2 次元 KdV 方程式にこだわった。この方程式は Kadomtsev-Petviashivili 方程式と呼ばれていたので、 Korteweg-de Vries 方程式が KdV 方程式に簡約されたのと同じで、きっと簡約されるだろうということで KP 方程式と短く呼ぶことにした。一つ突破口ができて、英語の論文の書き方もマスターすると、後は次から次へと英語の論文を書くことができ、結局 8 編英語の論文を書いて学位を取得することができた。

ここで転機が訪れる。私の勝手に敬愛していた数学者に京都大学の佐藤幹夫先生がいた。天才である思っていた。天才なら、なんでこんなに面白いソリトンに興味を示さないのだろうと、これも勝手に思っていた。ところがやはり興味を示し始めたのである。京都大学の数理研でソリトンの研究会が開かれるようになり、私は、研究会の後、佐藤先生のグループに呼び出されて、長時間、話を聞いて下さった。そのとき、私は最新の論文のプレプリントをお渡しした。これには、 KP 方程式のソリトン解を行列形に書き直す話も書いてあった。プレプリントをお渡しして半年位たったとき、突然、数理解析研究所でソリトンの研究会を開くと、佐藤先生のグループから連絡が入った。出席してみると、無限次元グラスマン多様体のプリュッカー関係式が広田の双線形方程式であるという、発表であった。いわゆる、ソリトンの佐藤理論の誕生の瞬間であった。その意義はここでは述べるスペースがないが、この理論の誕生によって私は、計画していた連作 5 本分の論文が没になったし、多くの物理数学的な研究者がソリトンの分野から撤退したように感じる。それほどソリトンの基礎の枠組みを明らかにした理論であった。私はすでに学位も取得していたので、佐藤先生とまったく違う道を行こうと決心した。佐藤先生の理論を追及する先にオリジナリティを出すのは、論文は生産できたとしても、真の意味の展開を行うことは極めて難しいと判断したからであった。

ソリトン理論は代数的な恒等式の上になりたつ理論体系である。それなら、代数的な関係の成立しない、非線形問題を解くためのツールを作ろうと考えた。例えば、カオス系を厳密に解析できる道具を作りたい。それには、数値計算を利用するしかない。非線形問題を、線形問題で近似してまずその解を求め、次にこれを少しずつ与えられた非線形問題に近づけていきながら解の変化を追跡すると言うホモトピー法などを使って、非線形回路方程式の直流解を求める研究などをした。しかし、通常の数値計算は誤差を把握していないため、数値計算結果から、解の存在を含め、数学的にはなんら主張ができない。 1990 年になって、情報処理学会の会誌で精度保証付き数値計算の特集が出て、数値計算の誤差が丸め誤差も含めて把握する方法があることが書かれていた。私はこれだと思った。しかし、最初は、丸め誤差の把握までできるとは到底信じられなかったので、とにかく自分で実際にプログラムを書いて確かめたいと思った。そこで、その年の連休に C で精度保証用のプログラムを作ってみようとチャレンジしてみた。もともと私はすべて手計算で数式を扱うことが研究の主方法でプログラム作りをしないほうであったから、 180 度の方向転換であった。このときは、実質的には多倍長精度のライブラリを書こうとした。精度保証には桁数が必要と考えたからである。しかし、多倍長精度のプログラムは初心者には難しいものの代表である。ものの見事に失敗した。 C のメモリの確保と開放の概念に不慣れであるのだから当たり前である。このような時に、柏木君が Calc という有理数演算数値計算ツールがあると教えてくれた。そして、 Calc はオブジェクト指向プログラミングもでき、 C のように変数の型にうるさくないので非常に気持ちよくプログラミングできることも教えてくれた。そこで、 Calc の上に区間演算を実装して、区間解析を試して見たり、フーリエ多項式で表される関数を一つのオブジェクトとして実装して、非線形常微分方程式の解の存在検証ができることを確認した。カオス解をもつダフィング方程式の周期解や周期倍分岐現象の存在検証が精度保証付き数値計算を使って実際にできるのだから精度保証付き数値計算は本物である、と確信することができた。とにかく、工学や理学に現れる実際の非線形微分方程式という難しくかつ応用がある本物の方程式を厳密に解くことができるのだから、精度保証付き数値計算は本物であるという確信がもてたので、今度はこれを高速化し磨きをかけることに力を注ぐことにした。道具として精度保証付き数値計算を確立したいと言うわけである。確かに私はホモルーベンス (道具をつくる動物) であることを示そう。有理数演算はソフト的に実現しえいるので、実行時間が極めて遅い。また演算子多重定義も遅くなる理由である。当時、有理数演算で、 100 元の連立方程式を解くと、一晩では終わらない位であった。これを例えば 10 時間とすると、 36000 秒である。浮動小数点演算では 100 元の連立方程式は 0.01 秒位では楽に解けるので有理数演算で厳密に解くと 3 百 60 万倍は浮動小数点演算に比べて時間がかかる。これを何とかしなくてはいけない。そこで、浮動小数点演算を使うことを必然とし、 C++ などのオブジェクト指向言語を利用して区間クラスを定義して、どのくらい高速化されるか試してみた。その結果、普通の浮動小数点演算の 1000 倍程度までには高速化できることがわかった。しかし、区間演算を利用して、例えば区間ガウスの消去法を使うので、解ける連立一次方程式のクラスまで制限されてしまう。区間演算は爆発するというのが正直なところであるから。

これらの問題を攻略するにあったって、何が遅くなる原因か調べてみた。その結果、区間演算を演算子多重定義することに多くの時間がかかることがわかった。こうして、目標が定まった。区間演算を演算子多重定義することなく、うまく実装せよ。この目標が定まった後も多くの紆余曲折があった。とにかく、こうやって精度保証付き数値計算の研究を始めて 10 年たったころ、非常に素直なかつ高速な精度保証方式を思いついた。適用範囲はまだ線形数値計算代数の分野に限られるが、それは簡単なアイディアだけにどんな問題にも適用できそうな方法である。標語で言うと、線形連立方程式の数値解は、もう一回数値解を計算する手間で精度保証できる。著者のもっているノートパソコンは PentiumIII/600MHz であるが、これに最速の数値計算パッケージを導入すると 1000 元連立方程式の数値解が 2 秒ほどで求められる。その精度保証には 3 秒弱で成功する。ついに精度保証が数値解を求めるのとほぼ同じ計算時間でできるようになった。著者にとっては 3 百 60 万倍の高速化が達成されたことになる。ちょっと信じられないというのが実感である。以下、節を改めて、その概要を論じよう。

2 丸めのモード制御による精度保証付き数値計算法

現代の数値計算は浮動小数点演算用の高速なハードウエアを利用して行われている。以下、本稿では浮動小数点演算の結果、問題の解が得られたとき、その解を数値解と呼ぶことにする。不思議なことに、現在の所、得られた数値解が真の解の存在を意味するかとか、真の解と数値解の誤差を厳密に評価することは普通行われていない。これは、このような真の解の存在や厳密でシャープな誤差評価が大変難しいものであるという認識によるものと思われる。しかし、振り返って、離散問題の計算量理論を考えてみよう。そこでは、真の解を導くのにかかる計算時間を問題にするのであって、真の解との関係が明らかでない数値が得られたとしても評価されないことである。これは、数値計算の扱う問題が連続問題であるからといって全く同じことであると思われる。数値計算においても、やはり、真の解に関わる数学的に正しいことが保証された情報が引き出されて始めて、計算が実行されたと考えるべきであろう。数値解の精度保証とは、数値解がどの位真の解に関して情報をもつかを引き出す数値計算である。従来は、数値解の精度保証は極めて困難である、あるいは実用的でないという認識により、これを行わないことが一般的であった。

本稿は、最近の数値解の精度保証の技術の進歩により、少なくとも、多くの線形問題に対しては、その数値解を計算するのの、数倍の手間で精度保証が可能となりつつあることを報告する。特に、線形連立方程式の数値解の精度保証は、数値解を得るのと同じ手間で精度保証できる。また、行列の固有や特異値も高速に精度保証できることを示そう。数値解析の分野で、数値線形代数という言葉がある。これは、線形連立方程式の数値解を求めたり、行列の固有値や特異値を数値計算する、数値解析の一分野を指す言葉である。数値線形代数はすべての数値計算の基礎となる、数値計算分野において極めて重要な位置を占める。本稿で報告することは、数値線形代数の分野に関しては、精度保証は数値解を求めるのとさして違わない手間でできるようになったということである。別の言葉でいえば、数値線形代数分野の精度保証技術は実用の域に到達したということである。

このような進展のブレークスルーは著者とドイツの研究者 Rump の共同研究で開発された丸めのモードの制御方式による精度保証技術である。これは、浮動小数点数のよく普及した規格である IEEE754 を満たす計算機環境があること及び IEEE754 の丸めのモードの切り替えができることを仮定している。そして、浮動小数点演算の丸めのモードを適切に切り替えて数値計算するという新しいアイディアによって、数値線形代数分野の問題の数値解の精度保証を可能とする技法である。これにより、一から精度保証のための計算機環境を構築することなく、スーパーコンピュータの環境でも、並列計算機の環境でも、パソコンの環境でも、今までに数値解を得るのに必要だった計算時間とほぼ同じ計算時間または数倍の計算時間で精度保証が可能となることが示されつつある。

本稿では、数値解を得るための最適な環境を構築している読者が、丸めのモードを切り替える命令さえ実装できれば (これは IEEE754 に従う CPU を使っている計算機であれば簡単である)、その環境の中で精度保証できる方法を示すことを目的とする。また、まだ、そのような数値計算環境を構築されていない読者も多いと考えて、共通のモデル環境として、まず、数値解を得る為の環境を構築する方法を 2 通り示す。具体的には、 1 番目の環境は Fortran や C を使う環境で、 LAPACK と高速 BLAS の導入により、高速で安定な環境を構築するモデルである。 2 番目は Matlab や Matlab 互換な Octave や Scilab により数値計算環境を構築する例である。 1 番目の環境に対して、 2 番目の環境 (Matlab は Matlab5 のバージョンとして) では 5 倍程計算時間が遅くなるが、インタプリタで行列型が組み込み型として演算子多重定義されていることを利用するとプログラム開発はそれに見合うだけ簡単となる。また、例え 2 番目の環境でも、 C や Fortran で一から書いた CPU のメモリの特性を考えていない普通のプログラムよりはずっと高速である。また、最近開発版が出された、 Windows/Penmtium III 用の MATLAB6 は様々な線形問題の数値解を LAPACK と高速 BLAS の組み合わせとほぼ同じ時間で計算する。すなわち、 Matlab6 では、プログラム開発が簡単な上に、実行速度も組み込みの命令を使えば最高速となっている。

尚、日本語における精度保証付き数値計算の専門書 [1],[2],[3] が出版され、深くこの分野が学べるようになった。数値計算とは精度保証付き数値計算であるという立場からの数値計算の入門書も出版されている [4]。本稿では、これらの資料を補足して、数値線形代数分野においては、著者と Rump が開発した丸めモードの制御による精度保証数値計算方式によって、従来の数値計算パッケージをそのまま利用して、高速に精度保証ができることを明らかにし、数値線形代数分野での精度保証がほぼ実用段階に到達したことを報告する。

3 精度保証付き数値計算環境の構築

精度保証がどれくらいのコストでできるかを論じるときに、通常の数値解を求める場合の、何倍の時間で計算できるかを論じるべきであるというのが、ここ数年来の著者の立場である。本稿でも、読者が現在使われている環境の中で、丸めの制御が可能であれば、精度保証が高速にできることを示すことを目的としている。

しかし、数値計算環境を構築されていない読者も多いと考え、通常の数値計算環境がどのようなものかをモデル的に示すために、通常の数値計算のための計算機環境をどう構築するかという点から議論しよう。そして、このような数値計算環境の中で計算した数値解を精度保証することを考える。この場合、精度保証にどのくらい時間がかかるかは、その環境での数値解の計算時間との比較で行うことにする。以下、まず、 2 種類の数値計算環境のモデルを示す。

3.1 LAPACK での計算機環境

最初は、 Fortran あるいは C で数値計算する環境である。 Fortran や C コンパイラについは議論を省こう。このような言語を用いた数値計算ライブラリは高度に発展しており、安定性、高速性の観点から、このようなライブラリを利用して数値計算することが絶対と言っていいほど必要である。このようなライブラリの代表の一つはLA-PACK(レイパックと発音する、 phase(http://phase.etl.go.jp/mirrors/ netlib/) からダウンロードできる) である。 LAPACK は現代コンピュータのメモリ階層 (キャッシュヒット率) を考えた高速で信頼性のある数値線形代数の直接解法パッケージである (図2)。 LAPACKのユーザマニュアル(第3版)はhttp://www.netlib.org/lapack/lug/

LAPACK BLAS CPU1  CPU2  並列計算機ポータブルにするため BLAS の CPU,OS ごとに最適化されたの関数を用いて書かれている。 BLAS を作る。

図 2: LAPACK と BLAS

lapack lug.html から見ることができる。メモリ階層を考えると言うことは、現代では一番高速なメモリであるキャッシュにデータを送り込んだら、いかにその中にあるデータだけで計算が進むかを工夫することである。このような観点から、 LAPACK では行列の分割解法が利用されている。どの位の行列に分割するかは、与えられたハードウエアによって異なる。また、アセンブラを利用すれば高速化も図れる。しかし、このような方向では、ハードウエアごとに異なったコードを考えなくてはならない。 LAPACK の開発者は数値線形代数分野のライブラリを作成するに当たり、基本となる線形代数演算を関数仕様としてまとめ標準化することにより、高速化とコンピュータ依存を避けるポータブル化を実現している。このような基本線形代数演算関数のライブリを BLAS(\underline{B}asic \underline{L}inear \underline{A}Lgebraic \underline{S}ubprograms) という。 BLAS にはベクトルの和、差、スカラー倍などのレベル 1 の関数、行列とベクトルの積などのレベル 2 の関数、行列と行列の積などのレベル 3 の関数がある。 LAPACK は BLAS の関数によって記述される。その BLAS はメーカなどがハードウエアごとに最適なものを作成するという仕組みである。したがって、ユーザは LAPACK のポータブルな関数でプログラムを書き、自分の環境にあった BLAS を捜してきて、その上でプログラムを実行することになる。 BLAS には、メモリ階層を考えず、 Fortran あるいは C で書かれたプログラム作成時用のものがある。これを標準 BLAS という。自分の計算機用に最適化された BLAS を用いた場合、標準 BLAS に対して 10 倍以上は高速となるのが普通である。 Fortran や C で書かれた普通の数値計算プログラムは安定性を得るのに大変であるが、それをおいても、標準 BLAS を用いた LAPACK のライブラリに速度の面でも勝つことは大変である。これに対して、最適化 BLAS を使って LAPACK を走らせると標準 BLAS を用いた場合に比べて 10 倍以上は高速なのであるから、ユーザが LAPACK にある解法ライブラリを使うのではなく、自前のプログラムを書いて使う意義は、プログラムや数値計算の学習という観点以外には全くと言っていい程存在しない。

LAPACK は phase(http://phase.etl.go.jp/mirrors/ netlib/) から手に入る。 Pentium II,III, Celeron で Linux を使う場合は最適化 BLAS は http://www.cs.utk.edu/ ˜ ghenry/distrib/archive.htm を通じて入手することができる (Pentium III 用の BLAS はまでできていないようなので PentiumII 用のもので代用させる)。また、 C のコンパイラの特性を考えてループ展開などの技法を駆使して、 C で書かれているにもかかわらず、コンパイル時にループ展開の数などのパラメータを最適化することによって最適化 BLAS に匹敵するオブジェクトコードを生成する試みとして ATLAS(http://phase.etl.go.jp/mirrors/ netlib/atlas/) がある。 ATLAS をコンパイルすると非常に時間がかかるが、これは一回のみである。 ATLAS を最適化 BLAS の代わりに使用することもできる。したがって、 Fortran あるいは C コンパイラがあれば、 ATLAS と LAPACK の組み合わせで、最適化 BLAS と LAPACK の組み合わせに匹敵する環境を構築することができる。

以下、簡単に、 Pentium II, III 上で Linux を用いている場合のインストール法を説明しよう。 phaseなどからLAPACKのターボールをダウンロードしてきて、 makeする。できあがったライブラリlapack alpha.aやlibF77.aを例えば/usr/local/lib/に置く。そして、 Linux 用の最適化 BLAS である lsblaspii1.2c 11.99.a も/usr/local/lib/ に置けば準備は完成である。 LAPACK のユーザマニュアル (第 3 版) などを参照しながら、プログラムを完成させコンパイルすればよい。 CLAPACK を利用した場合の、プログラム (test.c という名前であるとする) のコンパイル用 Makefile の一例を以下に示そう。

     LIB = /usr/local/lib/lapack_alpha.a
           /usr/local/lib/lsblaspii1.2c_11.99.a
           /usr/local/lib/libF77.a
     test: test.o $(LIB)
     cc -o test test.o $(LIB)
     test.o: test.c
     cc -c test.c

こうして、 LAPACK と最適化 BLAS、または、 LAPACK と ATLAS を組み合わせた数値計算環境が構築された。この環境の下では、 Pentium III/500MHz の CPU の場合、 1000 × 1000 のフルな一般係数行列に対する連立一次方程式は 3 秒を超えない程で解ける。

3.2 Matlab、 Octave あるいは Scilab による計算機環境

LAPACKをFortranやCLAPACKをCから呼び出して使う方法は、最適化BLAS との併用により、最高速での線形計算環境の一つを提供する。しかし、 Fortran や C のプログラミングは、場面により、煩雑に過ぎる場合がある。このような場合、使い易い GUI をもたせた数値計算用インタプリタ言語の利用が勧められる。計算時間は Matlab5、 Octave、 Scilab の場合 LAPACK と最適化 BLAS での計算時間の 5 倍程となるが、 Fortran や C での (普通の数値計算の教科書の例題に示されているような) メモリ階層を考えていない単純なコードによる数値計算より 2 倍前後は高速となる。また、先に書いたように Windows/Penmtium III 用の MATLAB6 では、少なくとも、連立一次方程式の数値解を LAPACK と高速 BLAS の組み合わせとほぼ同じ時間で計算する。

このような数値計算言語としては商用の Matlab が代表的で、 GUI やグラフ作成能力を含めると高い性能を誇っている。欧米では、大学ユーザーへの思い切った割引などで、 Matlab の普及には目を見張るものがあり、教育、研究にかなり使われている。 Matlab は商用であり、 Windows 版がある。 Windows でのインストールは標準的である。そこで、 Linux 上でフリーな Matlab クローンでオープンソースである Octave (http://www.che.wisc.edu/octave) あるいは Scilab(http://wwwrocq.inria.fr/scilab) を用いて数値計算環境を整える方法を示そう。 Octave と Scilab は LAPACK などのパッケージを基にしており、やはり極めて高い数値計算性能と優れたユーザインタファイスを持っている。 Octave は Linux において、ターボー

ルを展開して、 ./configure と make と make install によってインストール

する。 Scilab も同様であるが、少し Makefile が不備な面があり、若干修正が必要である。 Scilab のインストールの詳細は著者による Matalab, Octave, Scilab の使い方の html ファイル (http://www.oishi.info.waseda. ac.jp) などを参照されたい (そのテキストも出版される [5])。

Matlab5、 Octave,Scilab では、 Pentium III/500MHz の CPU の場合、 1000×1000 のフルな一般係数行列に対する連立一次方程式は 15 秒程で解ける。また、 Matlab6 を同じ計算機の Windows 上で実行した場合、 3 秒程で解ける。

  Demmel[6] の Matlab コードが (http://www.cs.berkeley.edu/~demmel/) にあ

るのでよいプログラム例として参考にされたい。

4 浮動小数点演算

現在の浮動小数点演算は IEEE754 規格に基ずくものが標準的である。これは優れた数値計算学者である Kahan らの努力により 1985 年に制定された (http://www.cs. berkeley.edu/ ˜ wkahan/ ieee754status/)。現在では Intel の Pentium CPU を含むほとんどのマシンが IEEE754 に従っているが、どのような CPU が例外かなどは Interval FAQ の記録 What machines support IEEE754?(http://studsys.mscs. mu.edu/˜ georgec/ IFAQ/casares1.html) が参考になる。また、 IEEE754 にしたがってどのように浮動小数点演算が実装されているかは、例えば、 Intel の PentiumIII のマニュアル (http://developer.intel.com/design/intarch/pentiumiii/ Manual.htm) などを見るとよい。

IEEE754 によって標準化されているのは 2 進の倍精度浮動小数点数である。浮動小数点数を利用した精度保証付き計算には IEEE754 の丸め (以下に示すものの中で上と下への丸め) を利用する。 c を実数 (c ∈ ℝ) とする。

(1) 上向きの丸め (round upward) c 以上の浮動小数点数の中で最も小さい数に丸る。これを △ : ℝ → 𝔽 と表す。プログラム中では setround(Inf) で上への丸めのモードに切り替える命令を表すことにする。

(2) 下向きの丸め (round downward) c 以下の浮動小数点数の中で最も大きい数に丸める。これを ▽ : ℝ → 𝔽 と表す。プログラム中では setround(-Inf) で下への丸めのモードへの切り替えの命令を表すことにする。

これ以外にも近似計算のデフォールトの丸めである最近点への丸め (round to nearest) とチョッピング (切り捨て)(round toward 0) の丸めがある。 IEEE754 では浮動小数点数演算 (𝔽 上での四則演算) は丸めとの関係によりつぎのように定義されている。 · ∈ {+, −, ×, /} ， ◯ ∈ {△, ▽} のとき

x ⨀ y = ◯(x · y) (任意の x，y ∈ 𝔽 について)  (1)

この式は，左辺の浮動小数点数の四則演算の結果 x ⨀ y は，右辺の数学的に正しい (実数としての) 四則演算の結果 x · y を指定された丸めを行って得られた数 ◯(x · y) に一致するように計算する規格を表している。また，平方根も

(\sqrt{x})_{fp} = ◯(\sqrt{x}) (任意の x ∈ 𝔽 について)  (2)

と，浮動小数点数演算によって計算された平方根 (\sqrt{x})_{fp} は，正確な実数演算で計算された平方根\sqrt{x} を指定された丸めの方向へ丸めた数となる規格である。注意すべきことは，指数関数や三角関数などはこのような規格をみたすようには規定されていないことである。

Octave や Scilab は C++を用いて書かれており、丸めの向きの変更は C の関数をリンクすることによって，実現される。その実現法は著者の本 [3],[4],[5] やホームページを参照されたい。また、 Matlab6 では丸めのモードの切り替えの命令が装備されている。

5 数値線形代数の精度保証

著者はドイツの研究者 Rump との共同研究で、丸めのモード制御計算による精度保証法を開発し、線形連立方程式の求解問題に対しては数値解を求めるのと同じ計算量で数値解の精度保証ができることを示した [7]。人間の計算でも、答えの検算は、もう一度やり直すのが普通であるから、数値解を求めるのと同じ計算時間で、近似解と真の解の誤差がシャープに求まるというのは理論的にもこれ以上高速化が望めないのではないかと推察される。

基礎は内積計算である。以下、説明を単純にするために、 Matlab6 で計算しているとしよう。 n 次元ベクトル u = (u_{1}, u_{2}, ⋯, u_{n})^{t} と v = (v_{1}, v_{2}, ⋯, v_{n})^{t} を考える。ただし、 u_{i}, v_{i} は i = 1, 2, ⋯, n に対して倍精度浮動小数点数であるとする。このとき、 u と v の内積 z を計算することを考えよう。よく考えると

   > setround(-Inf);
   > zd=u'*v;
   > setround(Inf)
   > zu=u'*v;

とすると z ∈ [zd, zu] となることがわかる。 Matlab ではベクトル u を縦ベクトルとすると、 u′ は u を転置した横ベクトルとなる。したがって、 u, v に対して u′ ∗ v はベクトルの内積演算となる。内部では LAPACK のベクトルの内積を計算するライブラリ関数が呼び出されて計算される。この計算は丸めの制御計算方式の典型となる。すなわち、従来の高速ライブラリの関数を呼び出して精度保証計算を進める方式で、従来使用されていなかった丸めの制御命令が、関数呼び出しの前と後に位置し、関数の呼び出し (ひいては高速化) の妨げにならないような計算方式となっている。

5.1 線形連立方程式の数値解の高速精度保証

以上の内積の精度保証計算が、通常の近似内積計算の 2 倍の計算速度で実行できることはプログラムを見てすぐにわかる。内積計算が高速に精度保証できると、いわゆる数値線形代数のアルゴリズムは高速に精度保証することができる。

まず、行列方程式

Ax = b  (3)

を解く問題を考えよう。ただし、 A を n × n 行列、 b を n ベクトルとする。行列 A が入力してあるとすると、式 (3) は

   > c=A\b;

という一行の命令で解くことができる。このとき、 LAPACK の線形方程式の解法関数が呼び出されて、式 (3) の数値解が求められる。

この数値解を精度保証することを考える。次の定理はよく知られており、これが基本となる。

定理 1 R を A の逆行列の近似とする。もし、

∥RA − I∥_{∞} < 1  (4)

が成立するなら、 A^{−1} が存在して、次の評価が成立する:

∥x^{∗} − c∥_{∞} ≤ \frac{∥R(Ac − b)∥_{∞}}{1 − ∥RA − I∥_{∞}} ≤ \frac{∥R∥_{∞}∥(Ac − b)∥_{∞}}{1 − ∥RA − I∥_{∞}}  (5)

ただし、 x^{∗} は (3) の真の解とする。 □ この定理に基づいて、数値解 c の精度保証をしてみよう。 A の数値的な逆行列は

   > R=inv(A)

とこれも一行で求められる。 RA − I の包み込みを計算するプログラムは次のようになる。

   > setround(-Inf);
   > Gd=R*A-eye(n,n);
   > setround(Inf);
   > Gu=R*A-eye(n,n);

ただし、 eye(n,n) は n × n 単位行列を作る命令である。このプログラムを実行すると RA − I ∈ [Gd, Gu] となる。これから ∥RA − I∥_{∞} の上界を計算するプログラムは次のようになる。

   > Gmax=max(abs(Gd),abs(Gu));
   > setround(Inf);
   > Gnorm=norm(Gmax,inf);

abs も max もそれぞれ、行列の成分毎に絶対値及び最大値を取る演算として定義されている。また、 norm(Gmax,inf) は行列 Gmax の最大値ノルムを計算する関数である。 ∥R∥_{∞} も ∥(Ac − b)∥_{∞} も同様に厳密な上界が計算できる。したがって、定理 1 より ∥RA − I∥_{∞} < 1 という評価が得られる場合には ∥x^{∗} − c∥ の上界が計算できることがわかる。

このプログラムの計算量を評価してみる。逆行列 R の数値計算に n^{3}flops かかり、 RA − I ∈ [Gd, Gu] という包み込みの計算に 2n^{3}flops かかるので、総合して 3n^{3}flops の計算量がかかることがわかる。これは、ガウス消去による数値解の計算量 n^{3}/3 の 9 倍の手間である。以上の計算法からわかるように、 Matlab などの命令の間に

setround(Inf); あるいは setround(-Inf); という丸めの制御命令が入っている

ため、 Matlab の関数がそのまま呼び出せることに注意すると、実際の精度保証計算でも数値解を得る手間の 9 倍の実行時間で精度保証できることがわかるであろう。従来は、区間演算型などを導入して、演算子多重定義のオーバーヘッドとキャッシュヒット率を落とした効果、並びにパイプライン処理が途切れるなどの効果で、精度保証に 100 倍ぐらいはかかっていたことに比べれると、大幅な高速化がされたことになる。しかし、 9 倍という計算時間は、まだまだ精度保証のコストが高いともいえる。そこで、精度保証の高速化を図ることを考えよう [7]。次のような方針を取ることとする。逆行列 R を計算する代わりに L̃ と Ũ の近似逆行列 X_{L} と X_{U} を計算する。ただし、 L̃ と Ũ はガウス消去によって A を LU 分解して得られた L と U とする。そして、 R の代わりに X_{U} · X_{L} を用いる。ただし、行列の積は計算しないものとする。 X_{L} の計算量は n^{3}/6 flops で X_{U} の計算量も n^{3}/6 flops となる。以下、これ以外の精度保証の計算には O(n^{2}) の計算量の計算しか現れないことを示そう。

定理 1 から

∥x^{∗} − x̃∥_{∞} ≤ \frac{∥(X_{U} · X_{L})(Ax̃ − b)∥_{∞}}{1 − ∥(X_{U} · X_{L})A − I∥_{∞}}  (6)

が成り立つ。ここで、 ∥(X_{U} · X_{L})(Ax̃ − b)∥ という項は行列とベクトルの積となるので O(n^{2}) flops で計算できる。高速化の基本的なアイディアは次である。 L, U, X_{L} 及び X_{U} をガウスの消去アルゴリズムで最近点への丸めのモードで計算したとする。このとき、最近得られた Higham による補題 (ftp://ma.man.ac.uk/oub/) を用いる。

補題 1 (Higham[8]) 最近点への丸めのモードの下で、 A ∈ ℝ^{n×n} に対してガウスの消去法で LU 分解 L̃ と Ũ を計算したとする。このとき、

L̃Ũ = A + ΔA, |ΔA| ≤ γ_{n}|L̃| |Ũ|,  (7)

が成立する。ただし、

γ_{n} = \frac{nu}{1 − nu}.  (8)

ここに、 u は単位の丸めで IEEE 754 の倍精度浮動小数点数に対しては 2^{−53} で与えられる。 □

補題 2 T を三角正方行列とする。 T の左逆 X_{T} をガウスの消去法で計算したとする。このとき、

|X_{T} T − I| ≤ γ_{n}|X_{T} | |T |.  (9)

が成立する。 □ この二つの補題を使うと、

∥X_{U}X_{L}A − I∥_{∞} ≤ ∥X_{U}X_{L}(A − L̃Ũ)∥_{∞} + ∥X_{U}X_{L}L̃Ũ − I∥_{∞}
≤ γ_{n}∥ |X_{U}| |X_{L}||L̃| |Ũ| ∥_{∞} + γ_{n}∥ |X_{U}| |Ũ| ∥_{∞}
+γ_{n}∥ |X_{U}| |X_{L}| |L̃| |Ũ| ∥_{∞}
= 2γ_{n}∥ |X_{U}| |X_{L}| |L̃| |Ũ| ∥_{∞} + γ_{n}∥ |X_{U}| |Ũ| ∥_{∞} (10)

という公式が導出できる。ここで e = (1, 1, ⋯, 1)^{t} ∈ ℝ^{n} とすると明らかに

∥ |X_{U}| |X_{L}| |L̃| |Ũ| ∥_{∞} = ∥(|X_{U}|(|X_{L}|(|L̃|(|Ũ|e))))∥_{∞}  (11)

及び

∥ |X_{U}| |Ũ| ∥_{∞} = ∥(|X_{U}|(|Ũ|e))∥_{∞}.  (12)

が成立する。これは ∥ |X_{U}| |X_{L}| |L̃| |Ũ| ∥_{∞} と ∥ |X_{U}| |Ũ| ∥_{∞} が O(n^{2}) で計算できることを示している。よって、今回の方法では X_{L} と X_{U} の計算に合計で n^{3}/3 flops の計算が必要で、他の計算は O(n^{2}) 以下であることがわかった。こうして、 ∥x^{∗} − x̃∥_{∞} の厳密な上界が\frac{1}{3}n^{3} + O(n^{2}) flops で計算できることがわかった。

例を示そう。上述の方法により，近似解とその精度保証をする C 言語プログラム (工夫していない普通のもの)、 LAPACK を利用して作成したプログラム (Intel の Pentium 用に最適化された BLAS を用いて実行する) による計算時間を測定した (表 5.1)。 CPU は Intel PentiumIII 600MHz である。

表 1: 高速精度保証の計算時間 (sec)

\begin{matrix}n & C & LAPACK \\ 200 & 0.3 & 0.05 \\ 400 & 2 & 0.24 \\ 600 & 8 & 0.7 \\ 800 & 19 & 1.5 \\ 1000 & 37 & 2.8 \end{matrix}

5.2 行列の固有値の精度保証

固有値や特異値の数値計算の精度保証の方式も丸めの制御方式で高速に実行できることも著者は示している [9]。

B, C, D と E を n × n 実行列で要素が 𝔽 に属するとする。準備として、まず、 ∥BCD − E∥_{∞} の厳密な上界を数値計算する方法を示そう。そのために最初に BC の包み込みを計算する関数を定義しよう。

function [\underline{P} , \overline{P}] = prod(B, C); setround(−Inf); \underline{P} = B ∗ C; setround(Inf); \overline{P} = B ∗ C; F, R, S を n × n 実行列とする。 R のすべての要素が非負ならば (このことを以下 R ≥ 0 と表す)、次が成立する:

[F − R, F + R]S = [F S − R|S|, F S + R|S|]  (13)

次のプログラムは [\underline{P} , \overline{P}] ⊆ [F − R, F + R] を満たす n × n 行列 F と R を計算する (R ≥ 0):

function [F, R] = cr(\underline{P} , \overline{P}); setround(Inf); F = \underline{P} + 0.5 ∗ (\overline{P} − \underline{P}); R = F − \underline{P};

A, F, R を n × n 行列で R ≥ 0 を満たすとする。次のプログラムは n × n 行列\underline{D} と \overline{D} で [\underline{D}, \overline{D}] ⊇ A ∗ [F − R, F + R] を満たすものを計算する:

function [\underline{D}, \overline{D}] = iprod(A, F, R); G = abs(A); setround(−Inf); \underline{D} = A ∗ F + (−G) ∗ R; setround(Inf); \overline{D} = A ∗ F + G ∗ R;

次のプログラムは ∥BCD − E∥_{∞} の厳密な上界を計算する:

function t = errt(B, C, D, E); [\underline{P} , \overline{P}] = prod(C, D); [F, R] = cr(\underline{P} , \overline{P}); [\underline{H}, \overline{H}] = iprod(B, F, R); setround(−Inf); \underline{H} = \underline{H} − E; setround(Inf); \overline{H} = \overline{H} − E; Hmax = max\{abs(\underline{H}), abs(\overline{H})\}; t = norm(Hmax, inf);

さて、行列の固有値の精度保証のための理論を示そう [9]。ここでは、 Bauer-Fike の定理 [10] を精度保証のために拡張した結果を示す [9]。 Theorem 1 (Oishi[9]) A を n × n 一般複素行列とする (これを A ∈ M_{n}(ℂ) と書く)。数値計算の結果、対角行列 T ∈ M_{n}(ℂ) と Q, P ∈ M_{n}(ℂ) が得られて

QAP = T + E, QP = I + F  (14)

となったとする。ただし、 E, F ∈ M_{n}(ℂ) で I は n × n 単位行列とする。 μ を A の固有値とする。すると、すべての ∥ · ∥_{p}, (p = 1, 2, ⋯) と ∥ · ∥_{∞} に対して次が成立する。

min_{j} |T_{jj} − μ| ≤ ∥E∥ + ∥A∥ · ∥F ∥  (15)

証明 λ_{j}, (j = 1, 2, ⋯, n) を T の固有値とする。 min_{1≤j≤n} |μ − λ_{j}| を評価しよう。 μ が T の固有値と一致するときは、 min_{1≤j≤n} |μ − λ_{j}| はゼロとなり定理の結論は明らかである。よって、そこで μ が T の固有値でない場合を考えよう。次が成り立つ。

(T − μI)^{−1}Q(A − μI)P = I − (T − μI)^{−1}(T − QAP + μ(QP − I))  (16)

A − μI が特異なことから、

1 ≤ ∥(T − μI)^{−1}(T − QAP + μ(QP − I))∥  (17)

が従う。これより、

∥(T − μI)^{−1}∥ = \frac{1}{min_{1≤j≤n} |μ − λ_{j}|}  (18)

を得る。これと

|μ| ≤ ∥A∥,  (19)

より式 (15) の成立がわかる。 □ 注意 1 不等式 (15) の右辺が上から ε で評価されたとする。このとき、この定理から A の固有値は ⋃_{j=1}^{n}{z ∈ ℂ : |λ_{j} − z| ≤ ε}, (j = 1, 2, ⋯) の中に存在することがわかる。 □

A ∈ M_{n}(ℂ) とする。

   > [P,T]=eig(A);

とすると、その対角成分が A の固有値の近似となっている対角行列 T と対応する固有値を並べた行列 P が得られる。したがって、 AP ≈ P となることが期待される。近似固有ベクトルを並べた P は近似なので、数値計算誤差が入るために、 P は一般に正則となる。例えば、

A = ( \begin{matrix}1 & 1 & 0 \\ 0 & 1 & 1 \\ 0 & 0 & 1 \end{matrix} )  (20)

とすると、 Octave は

T = ( \begin{matrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix} ) , P = ( \begin{matrix}1 & −1 & 1 \\ 0 & 4.44 ⋯ × 10^{−16} & −4.44 ⋯ × 10^{−16} \\ 0 & 0 & 1.97 ⋯ × 10^{−31} \end{matrix} ) . (21)

という結果を与える。すなわち、数値計算誤差が P を正則にする。 Q を P^{−1} を数値計算に求めたものとする。このとき、対角行列 T ∈ M_{n}(ℂ) と Q, P ∈ M_{n}(ℂ) で

QAP = T + E, QP = I + F,  (22)

としたとき、数値計算誤差 E, F ∈ M_{n}(ℂ) が十分小さくなるものが計算できることが期待できる。このとき、定理 1 より ∥ · ∥_{p}, (1 ≤ p ≤ ∞) に対して A の固有値の包み込みが

min_{j} |T_{jj} − μ| ≤ ∥E∥ + ∥A∥ · ∥F ∥.  (23)

と計算できる。以下、 ∥E∥_{∞} + ∥A∥_{∞} · ∥F ∥_{∞}. の厳密な上界を計算するアルゴリズムを示そう。そのために、 2 つの複素行列 A, B ∈ M_{n}(ℂ) の積の包み込みを計算するためのアルゴリズムを示そう。

function [\underline{R_{E}}, \overline{R_{E}}, \underline{I_{M}}, \overline{I_{M}}] = cprod(A, B)

A_{r} = real(A); A_{i} = imag(A); B_{r} = real(B); B_{i} = imag(B); setround(−Inf); \underline{R_{E}} = A_{r} ∗ B_{r} + (−A_{i}) ∗ B_{i}; \underline{I_{M}} = A_{r} ∗ B_{i} + A_{i} ∗ B_{r}; setround(Inf); \overline{R_{E}} = A_{r} ∗ B_{r} + (−A_{i}) ∗ B_{i}; \overline{I_{M}} = A_{r} ∗ B_{i} + A_{i} ∗ B_{r};

このアルゴリズムは n × n 行列\underline{R_{E}}, \overline{R_{E}}, \underline{I_{M}}, \overline{I_{M}} で

A ∗ B ⊆ [\underline{R_{E}}, \overline{R_{E}}] + \sqrt{−1}[\underline{I_{M}}, \overline{I_{M}}].  (24)

となるものを計算できる。ここで real(A) と imag(A) は、それぞれ、複素行列 A の実部と虚部を計算する命令である。

次の MATLAB 的なプログラムは複素行列 B, C ∈ M_{n}(ℂ) に対して ∥BC − I∥_{∞} の厳密な上限の評価を与える。

function d = cerrb(B, C, eye(n)); [\underline{R_{E}}, \overline{R_{E}}, \underline{I_{M}}, \overline{I_{M}}] = cprod(B, C); setround(−Inf); \underline{R_{1}} = \underline{R_{E}} − eye(n); setround(Inf); \overline{R_{1}} = \overline{R_{E}} − eye(n); R_{1} = max(abs(\underline{R_{1}}), abs(\overline{R_{1}})); I_{1} = max(abs(\underline{I_{M}}), abs(\overline{I_{M}})); E = R_{1} + i ∗ I_{1}; d = norm(E, inf);

更に、次の MATLAB 的なプログラムは複素行列 B, C, D, E ∈ M_{n}(ℂ) に対して ∥BCD − E∥_{∞} の厳密な上界を計算する。

function c = cerrt(B, C, D, E); B_{r} = real(B); B_{i} = imag(B); E_{r} = real(E); E_{i} = imag(E); [\underline{R_{E}}, \overline{R_{E}}, \underline{I_{M}}, \overline{I_{M}}] = cprod(C, D); [R_{c}, R_{r}] = cr(\underline{R_{E}}, \overline{R_{E}}); [I_{c}, I_{r}] = cr(\underline{I_{M}}, \overline{I_{M}}); [\underline{B_{r}R}, \overline{B_{r}R}] = iprod(B_{r}, R_{c}, R_{r}); [\underline{(−B_{i})I}, \overline{(−B_{i})I}] = iprod(−B_{i}, I_{c}, I_{r}); [\underline{B_{r}I}, \overline{B_{r}I}] = iprod(B_{r}, I_{c}, I_{r}); [\underline{B_{i}R}, \overline{B_{i}R}] = iprod(B_{i}, R_{c}, R_{r}); setround(−Inf); \underline{R_{1}} = \underline{B_{r}R} + \underline{(−B_{i})I} − E_{r}; \underline{I_{1}} = \underline{B_{r}I} + \underline{B_{i}R} − E_{i}; setround(Inf); \overline{R_{1}} = \overline{B_{r}R} + \overline{(−B_{i})I} − E_{r}; \overline{I_{1}} = \overline{B_{r}I} + \overline{B_{i}R} − E_{i}; R_{1} = max(abs(\underline{R_{1}}), abs(\overline{R_{1}})); I_{1} = max(abs(\underline{I_{1}}), abs(\overline{I_{1}})); T = R_{1} + i ∗ I_{1}; c = norm(T, inf);

次の MATLAB 的なプログラムは与えられた A ∈ M_{n}(ℂ) に対して ∥E∥_{∞}+∥A∥_{∞}· ∥F ∥_{∞} の厳密な上界を計算する。

function e = vgeig(A);

[P, T] = eig(A);
Q = inv(P);
e_{1} = cerrt(Q, A, P, T);
f = cerrb(Q, P, eye(n));
e = e_{1} + norm(A, inf) ∗ f;

つぎの表はこのアルゴリズムによってすべての固有値を計算した結果を示している。ただし、 CPU は Pentium III/500 MHz で Matlab6 を使っている。

表 2: Verification of all eigenvalues

\begin{matrix}次元 n & 計算時間 (sec) & 保証された数値計算誤差 \\ 50 & 0.6 & 3.5 × 10^{−12} \\ 100 & 1.2 & 2.2 × 10^{−11} \\ 150 & 3.9 & 6.8 × 10^{−11} \\ 200 & 9.1 & 1.2 × 10^{−10} \\ 300 & 29.7 & 4.6 × 10^{−10} \end{matrix}

6 むすび

本稿では、スペースの関係で省略したが、行列の特異値の計算も同様にできる [9]。楕円型偏微分方程式の有限要素解析に現れるような、数万元の線形連立方程式もそのスパース性を生かした反復解法を利用した精度保証法により、前節のパソコンで数秒でその数値解の精度保証が可能となっている (著者の研究室の仕事)。このように、数値線形代数の分野では、精度保証が高速かつ簡単に行えるようになっており、その実用化の段階を迎えた。多くの読者にこの事情を理解していただき、日常的に、数値解の精度保証が行われるようになることを期待して筆を置くことにしたい。

謝辞

著者の日頃の研究において、色々ご検討賜る愛媛大学山本哲朗先生ならびに、ハンブルグ工科大学の Rump 教授に厚く感謝いたします。また、本稿の執筆に当たり、成蹊大学工学部弓削康平先生には大変お世話になったことを深謝いたします。

参考文献

[1] 大石進一: 非線形解析入門、コロナ社 (1997)

[2] 中尾充宏、山本野人: 精度保証付き数値計算、日本評論社 (1998)

[3] 大石進一: 精度保証付き数値計算、コロナ社 (2000)

[4] 大石進一: 数値計算、裳華房 (1999)

[5] 大石進一: “Linux 数値計算ツール”, コロナ社 (2000.9)

[6] J. Demmel: “Applied Numerical Linear Algebra”, SIAM (1997)

[7] Shin’ichi Oishi and Siegfried M. Rump: “Fast verification of solutions of matrix equations”, submitted to Numer. Math..

[8] Nicholas J. Higham: “Accuracy and Stability of Numerical Algorithms”, SIAM (1996).

[9] Shin’ichi OISHI: “Fast Enclosure of Matrix Eigenvalues and Singular Values via Rounding Mode Controlled Computation”, to appear in Linear Algebra and its Applications

[10] F. Bauer and C. Fike: “Norms and exclusion theorems.”, Numer. Math., 2, pp.137-141 (1960).