大石進一 （早稲田大学理工学術院教授）
精度保証付き数値計算の理論とツールについて論じていきます。 理論としては数値線形代数の高速精度保証、線形計画問題、非線形方程式の解法、関数方程式、計算機援用証明などから話題を選んで解説する。ツールについては既存の商用ツール(MATLAB)、フリーのツール(Scilab)、自作のツール(SLAB)について議論する。
精度保証付き数値計算が実用化されつつある。その基礎として最も重要なものは、倍精度浮動小数点数の規格IEEE754である。これは1985年にKahanという優れた数値計算の学者が中心となって制定したものである。IEEEは米国の電子工学会であるが、規格の制定も行っている。まず、IEEE754の規格を読むことから始めよう。
浮動小数点数の規格IEEE754はつぎのサイトでその規格書を読むことができる。
IEEE Standard for Binary Floating-Point Arithmetic
この規格書を読んでみましょう。PentiumIII,IVやSunのCPUなどほとんどすべての計算機はこの規格を満たしています。また、佐賀大皆本さんのIEEE754に関する資料(pdf)が公開されました。その中での注意に従って計算してみると、次のようなことがWindows上のMATLABでも起きることがわかります。
>> a=0/0
警告: ゼロ割です.
a = NaN
>> b=a^0
b = 1
数値計算を行う環境として、各種の数値計算ツールが利用可能です。これを紹介するのも今日の目的です。
ここでは、Scilab(サイラボと読む)とMATLAB(マットラボと読む)について紹介します。また、筆者が作ったSlabもあります。今日はこれらのツールの使い方も勉強しましょう。
皆さんのパソコンにはScilabが既にインストールされています。Scilabはインタプリタと呼ばれる種類のプログラム言語です。まず、その使い方を勉強してみましょう。
ここにScilabの使い方があります。これを読んでみましょう。ktermを立ち上げて
> scilab &
と打ち込むとScilabの入力画面が立ち上がります。
参考文献
ポアソン分布のグラフ\ref{fig:fig-poi.eps}を描くことを考えよう．
まず，ポアソン分布
P_k(t)=\frac{(\lambda t)^k}{k!}e^{-\lambda t},~(k=0,1,2,\cdots)
を計算するための関数を定義する方法を示そう．
エディタで次を内容とするファイルpoisson.sciを作る．
function [p]=poisson(l,n)
p=zeros(n+1,1);
p(1)=exp(-l);
for i=1:n,
p(i+1)=l*p(i)/i;
end;
このファイルがScilabのパスの通っているフォルダにあるものとしてScilabに呼び出す．
--> getf('poisson.sci')
その後，次のようにすると，図が描かれる．
-->p1=poisson(2,10);
-->p1=poisson(1,10);
-->p2=poisson(2,10);
-->p3=poisson(3,10);
-->n=0:1:10
n =
! 0. 1. 2. 3. 4. 5. 6. 7.
8. 9. 10. !
-->plot2d([n' n' n'],[p1 p2 p3])
レポート、論文、書籍などを書くのにtexが有用である。texの使い方については、つぎのメモがある。
さて、以上のようにScilabにグラフを書かせることは簡単であるが、これを文書の形に残すためのリテラシを簡単にメモしておこう。
Scilabでグラフを描かせた後、ScilabGraphの中のFIleメニューを開く。この中にExportがあるので、これを選択する。Export Typeを
Postscipt-Latexに選びOKとすると、保存するデイレクトリとファイル名を聞かれるので、例えばexp2.epsというファイル名
にして保存する。platex2etexのファイルでは、
\documentclass{jbook}
\usepackage{graphicx}
とスタイルファイルとしてgraphicx.styを利用する。図を描かせる位置で以下のように引用する。
\begin{figure}[htbp]
\begin{center}
\includegraphics[keepaspectratio=true,height
=60mm]{exp2.eps}
\end{center}
\caption{保留時間}
\label{fig:exp2.eps}
\end{figure}
こうして、図がtexファイル中に引用できる。図を描かせるには様々な方法があるので、これはあくまで一つのスタイルである。
Scilabはフリーのツールである。MATLABと呼ばれる商用のソフトウエアもある。これはグラフィックスもより高度で、ディジタル信号処理ツールボックスなどもある。また、数値計算ツールのデファクトスタンダードとなる文法を提供してきたことも見逃せない。
参考文献
Calcの使い方について述べて行く。最初は電卓のように感じるが、後にその奥の深さがわかると思う。
calcをインストールすると、コマンドラインから、
$ calc
と入力するとCalcが起動する。helpと入力するとhelpにどのようなトピックスがあるか表示される。そこで、例えば、help fullとすると、すべてのhelp情報を見ることができる。
こうして、立ち上がった状態で
> 5*(2+9)
> 55
のように、数式をC言語のように記述すると、電卓のように答えを返す。ここでさらに、
> .*10.5
> 577.5
とすると、最後に計算された値が'.'というシンボルによって呼び出され(今の場合は55)計算を続行することができる。Calcでは表示上は数は小数点で表され
るが内部ではすべて有理数で記憶されている。変数を導入することもできる。
> a=.
とすると。現在の値が変数aに記憶される。
ここで、階乗の計算を行ってみよう。
> a=100
とする。a!を計算するには次のようにする。
> fact(a)
93326215443944152681699238856266700490715968264381621468592963895
21759999322991560894146397615651828625369792082722375825118521091
6864000000000000000000000000
この計算は正確であり、浮動小数点演算ではできない種類のものである(桁数が多すぎる!)。
小数点以下何桁まで精度を保つかも指定できる。
> config("display", 50)
> epsilon(1e-50)
としておけば、次のようにしてcos(1)の値が小数点以下50桁まで正しく計算される。
> cos(1)
~.54030230586813971740093660744297660373231042061792
初等関数を計算する精度が制御できることがわかる。
参考文献
連立一次方程式
1.1x + 2.2y - 3.3z = 1
2.1x - 3.2y + 3.3z = 2
-3.1x + 4.2y + 4.3z = 3
をScilabとMATLABとCalcで解け。
まず、次のファイルを次の手順でダウンロードしよう。
では、scilabにおいて丸めのモードの切り替えができることを見てみよう。
と打ってみよう。丸めのモードが正しく設定されていることがわかる。
尚、丸めのDLLはVC++(ver. 6)で作成した。(Down.c, Near.c, Up.c)
では、少しscilabを動かしてみよう。ここにscilabの動かし方があるので見てみよう。
次の解説（ｈｔｍｌ）を読んで精度保証付き数値計算の原理を理解しよう。次の解説（電子情報通信学会2001年1月号掲載)
大石進一:実用段階に到達した精度保証付き数値計算(dviファイル、約100k)
も参考にされたい。
参考文献
精度保証付き数値計算 大石 進一 著 コロナ社 2,200円 1999/12発行 ISBN4-339-02605-0
数値計算 大石 進一
著 裳華房 2,000円 1999/03発行 ISBN4-7853-1514-8
連立一次方程式の解の精度保証をしてみよう。まず、次のファイルをダウンロードしてみよう。やはりSciproに保存する。
このプログラムを動かしてみよう。
これで精度保証ができている。その原理をみてみよう。
Matlab6についての丸めの制御法を示す。丸めの制御の命令は
>> system_dependent('setround',mode);
である。ただし，
mode=-Infのときround downへ
mode=+Infのときround upへ
mode=0のときround to zeroへ(チョッピング)
mode='nearest'のときround to nearestへ(default mode)
である。
例えば丸めを制御するためのm-file(setround.m)として
function setround(x)
if x==1
system_dependent('setround', Inf)
elseif x==-1
system_dependent('setround',-Inf)
else
system_dependent('setround', 'nearest')
end
をつくっておけばsetround(1)で上への丸め、setround(-1)で下への丸め、setround(0)で'nearest'への丸めという風に丸めを制御することができる。
実行例
>> setround(1)
>> a=1/3;
>> sprintf('%30.25f',a)
ans =
0.3333333333333333700000000
>> setround(-1)
>> a=1/10;
>> sprintf('%30.25f',a)
ans =
0.0999999999999999920000000
Computational Gridに関する話題を集めておいたので読んでみて欲しい。
昨日の復習。
64台のクラスタを利用して10000次元までの蜜係数行列をもつ連立一次方程式を解いた結果を述べた。
PCクラスタ上での連立一次方程式の精度保証付き数値計算(大石進一、荻田武史)
その結果、Oishi-Rumpの提案した高速解法は5000次元の連立一次方程式(良条件)ぐらいまでは有効であることがわかる。また、10000次元でも高速化していないアルゴリズムなら解けることもわかる。
などについて論じた。
用語 MPI (Message Passing Interface) 、BLAS(Basic Linear Algebraic Subroutines)
Bauer-Fikeの定理を利用して固有値の精度保証を行う方法について解説する。次の解説（電子情報通信学会2001年1月号掲載)
大石進一:実用段階に到達した精度保証付き数値計算(dviファイル、約100k)
も参考にされたい。また、講義では、成分ごとのBauer-Fikeの定理を示し、それを利用する方法についても解説する。この方法では局所的に固有値が一意であることの十分条件も示される。
まず、次の小島政和氏(東京工業大学教授)の資料（数理計画法の最近の発展 --- 内点法を中心にして ---）(pdf,100k)を見てみよう(またはより詳しい解説)。小島氏は内点法の研究で大変優れた成果をあげておられるが、その主双対内点法はまさに精度保証の理論を利用しているといっても過言ではない。ここでは双対定理を利用した計画問題の精度保証法について議論する。
例証 Slab
Scalabilityの確保
規模と条件数
(1) 関数の評価は規格がない。
例 超越関数の精度保証付き数値計算
(2) データ型が複雑である。
自動微分、区間、区間定数を含む微分係数など
演算子多重定義の時
精度保証に利用する離散化方程式の規模が大きくなる
特異性の問題
精度保証理論はほぼ存在定理の(構成的な)証明
微分方程式の解の存在証明(常微分、偏微分)
ケプラー問題など証明の最後のつめに使う。
カオス(Conley Index(Mischaikowのサーベイpdf)など構造安定な量の計算) 談話会終了後、石井(九大)さんからLorenzアトラクタの存在検証がW.Tuckerによってなされていることを教えていただいた。これはSmaleの第14問題(Lorenz微分方程式のアトラクタはWilliams, Guckenheimer, とYorkeが提案したgeometric Lorenz attractorであるか?)と関連する。そして、 Tucker (2002) はこれを肯定的に解いた(pdf,180k)。（中尾先生のご指摘、関連のページ)
計算幾何学などへの応用
最近の精度保証関係の国際シンポジウム(2002 SIAM Workshop on Validated Computing)のWorkshopプログラムを見ることができる。
精度保証がつくMATLAB-likeな数値計算用インタプリタを著者は作成した。そのKnow Howを議論する。インタプリタを作ることがオブジェクト指向プログラミングそのものとなったともいえる。
プログラミング言語の売りは「何がすごいところか」である。したがって，数値計算用インタプリタを作るとしたら，その売りは「すごい数値計算用のプログラムが含まれる」というところになる。その意味で，数値計算用のプログラムパッケージとして，プロが何を開発してきて，何が成功しているかを知らなくてはならない。この知識に関しては精度保証付き数値計算に詳しく解説している。
さて，数値計算パッケージの研究，実用化においては，もちろん色々な成果があるが，その筆頭の一つは
LAPACK+最適化BLAS
であろう。このパッケージが組み込まれた使いやすいインタプリタというのが，ここでの作成の目標となるインタプリタである。実はこのようなインタプリタは世の中に存在する。商用のインタプリタMATLAB6である。したがって，これはよい手本となる。
コンピュータサイエンスの進展は著しい。実は，目玉のパッケージがある場合には，その目玉のよさを発揮させるインタプリタやコンパイラの作成理論とツールは非常によく発達しているのである。この理論とツールを生かして，「LAPACK+最適化BLASを売りとする数値計算ツール」を作ることを目標にしてみよう。
以下，Linux,Windows+cygwin，Mac OSX(Developing Toolはインストール済みとする)などUnix系OSを利用していると想定する。これらの環境では，通常bison(yacc),flex(lex)などのコンパイラ作成ツールがもともとインストール済みである。また，メモリのガーべージコレクションのためのgcもGPLライセンスで利用できる(のでダウンロードして利用できる)。ここでは，これらのツールを使うことを前提に数値計算用インタープリタを作っていく方法を説明する。
C言語を用いて数値計算用のインタプリタを作るとき，flexは字句解析用のCプログラムを自動生成するツールで，bisonは構文解析用のCプログラムを自動生成するツールである。それぞれ，lexおよびyaccの上位コンパチである（使い方はbisonとyaccで若干異なる)。flexとbisonについてはweb上に解説があるのでそれらをまず読んでみて欲しい。
これらには簡単な電卓用のプログラムが附いている。すなわち，これらの説明を読めば直ちに
(1) bisonの機能を用いて，かけ算わり算を足し算引き算より優先させて計算する倍精度浮動小数点数を用いた電卓を作ることができる。少し頑張れば
(1+） 実部，虚部がともにdoubleの複素数が扱える，電卓
(1++) 行列(double)が扱えるようにする。ただし，加減乗算にはBLASの関数を用いる。
(1+++) 行列(複素数)が扱えるようにする。ただし，加減乗算にはBLASの関数を用いる。
(1++++) 変数を扱えるようにする。
(1+++++) sin,cosなどやべき乗もできるようにする。
などできる。
(2) 制御文(if文やwhile文やfor文)を扱えるようにする。
(3) ユーザが関数を定義できるようにする。
(2),(3)の問題を解くためには，解説１，２では不十分となる。これらができるようになるための説明のあるwebをリンクしよう。
ずばり，面白いページは次である。
前橋和弥氏のページ ここの「電卓を作ってみよう」がとても面白い。
これをベースに
(M+) 前橋氏の電卓にsin,cos,べき乗などの関数を組み込んでみよう。また，ユーザ定義関数の返す値を指定できるようにする。
(M++) 前橋氏の電卓に行列ができるようにしてみよう。もちろん，加減乗算はBLASの関数を使う。また，ついでだから，
Ax=b
を
> x=A\b
でできるようにしてみる。複素行列も扱えるようになったら，固有値問題も解けるようにしてみたい。
> [P D] = eig(A)
で行列Aの固有値（を対角要素に並べた行列）Dと対応する固有ベクトルを並べた行列Pが計算できるようにしたい。
(M+++) 前橋氏の電卓でgcでメモリの確保をしたらどうなるだろう。
(4) ファイルも読み込めるようにするというのは前橋氏の電卓はすでにできる。スクリプトと関数定義を読み込めるようにしてみよう。
(5) 絵が書けるようになりたい。gnuplotを呼び出す仕組みにするのは簡単だろう。ただ，これだけでは寂しい。Javaで絵を描くのはいいかもしれない。
尚，前橋和弥氏は氏のwebの参照やプログラムの利用を快諾して下さいました。ここに深く感謝の意を表します。ありがとうございました。また，bisonやflexの参考webを作成されている方にも謝意を捧げる。
次いで，現代応用解析Bで出題した課題の概要を示そう。
以上に基づき、実際にインタプリタを作成した。現在のところ，ユーザがオブジェクトを定義できるようにする機能以外は実装されており，総ステップ数が１万を優に超えながらも快適に動作している。これをSlabと名づけている。Slabをどのように公表するかは未定であるが，これを今日は詳細に紹介したいと考えている。
まず、Scilabによる非線形方程式の解法について述べる。
Newton法の解法ルーチンの例(NEWTON2.sci)
を示す。これはScilabのプログラムである。これには自動微分が使われている。自動微分のプログラムは次のようになる。
さらに、要素としてなんでも型が入る行列クラスが必要となる。それは、つぎで定義される。
これらをダウンロードしてフォルダSciproに入れておこう。そして、非線形方程式としてつぎの例を考える。
非線形方程式をNewton法で解くにはつぎのようにする。
-->getf('E:\SCIPRO\NEWTON2.SCI');
-->getf('E:\SCIPRO\MAT.SCI');
-->getf('E\SCIPRO\AUTODIF.SCI');
-->getf('E:\SCIPRO\FUNC.SCI');
-->x=newton(2,'func')
x =
x(1)
!mat dim !
x(2)
! 2. 1. !
x(3)
0.9995920
x(4)
0.0285618
このようにScilabにおけるクラス（オブジェクト）はt-listを用いて定義される。
もう少し複雑な例を解いてみよう。つぎの非線形関数のゼロ点を求めてみる。
function y=func(x)
y=x;
y(3)=x(3)*x(3)-x(4)*x(4)-3*x(3)+2;
y(4)=2*x(3)*x(4)-3*x(4)-1;
y(5)=3*sin(x(3)+x(4)+2*x(5))+x(5)*x(5)-3;
y(6)=x(6)*x(6)+3*sin(x(3))-3;
y(7)=cos(x(7))*x(7)+x(3)*x(3)-x(5)*x(5)+x(4)-2;
y(8)=x(8)+2*x(3)-x(7)+x(6)+x(3)*x(4)-5;
endfunction
この関数定義をつぎのファイルに入れた。これくらいになると、手で微係数を求めるのは少々めんどくさくなる。自動微分のありがたみがわかる。
さて、Scilabで解いてみよう。
-->getf('E:\SCIPRO\FUNC3.SCI');
Warning :redefining function: func
-->x=newton(6,'func')
x =
x(1)
!mat dim !
x(2)
! 6. 1. !
x(3)
2.3002426
x(4)
0.6248105
x(5)
- 0.4799244
x(6)
0.8737101
x(7)
14.396064
x(8)
12.484653
いよいよKrawczyk作用素を定義してみよう。それはつぎのようになる。
著者が考案した幾つかの方法を議論する。
つぎの論文に示した方法である。
Shin’ichi Oishi: Numerical verification of existence and inclusion of solutions for nonlinear operator equations, Journal of Computational and Applied Mathematics, Vol.60, 1995, pp.171-185
非線形楕円型方程式等を含む非線形作用素方程式の解の存在・唯一性・存在範囲を精度保証付き数値計算にもとづく計算機援用証明によって証明する方法を示した。また、そのためのソフトウェアライブラリを作製し、検証例を示した。
この論文を書いた頃(1993)は著者は有理数演算を用いていたが、現在では倍精度浮動小数点数を用いている。この方法でDufiing方程式を解くには次のようにする。以下、Scilabのプログラムである。
Duffing方程式の決定方程式と精度保証プログラム(DUFFING.sci)
これにはフーリエクラス(FOURIER.sci)が用いられている。
常微分方程式のホモクリニックオービットやヘテロクリニックオービットの数値的存在検証のために開発した方法。
Shin'ichi
OISHI: Two Topics in Nonlinear System Analysis through Fixed Point Theorems,
IEICI Trans. Fundamentals Vol.E77-A,
No.7, (1994) pp.1144-115
に一般の境界値問題の場合に対する検証法を書いた。また、次の論文に常微分方程式のホモクリニックオービットやヘテロクリニックオービットの数値的存在検証法を示した。
Shin'ichi Oishi: "A Numerical Method of Proving the Existence and Inclusion of Connecting Orbits for Continuous Dynamical Systems", Journal of Universal Computer Science, vol.4 no.2, (1998), 193-201
参考文献 非線形解析入門 大石 進一 著 コロナ社 2,800円 1997/04発行 ISBN4-339-02600-X
これについてはつぎの文献を参照しながら議論する。
参考文献 精度保証付き数値計算 大石 進一 著 コロナ社 2,200円 1999/12発行 ISBN4-339-02605-0
最近は村重氏(東京大学)と特異積分方程式の解の数値的検証法について研究を進めている。