7.8 Intlab_V9/intval/verifylss.m の冒頭のコメント

ヴァージョンが上がるごとにファイルサイズは単調増加している。

-rw-r--r--  1 mk  lab00   8041  1  3  2003 Intlab_Unix_V4.1.2/intval/verifylss.m
-rw-r--r--  1 mk  lab00  38913 12 12  2017 Intlab_V10.2/intval/verifylss.m
-rw-r--r--  1 mk  lab00  17050  3 16  2010 Intlab_V6/intval/verifylss.m
-rw-r--r--  1 mk  staff  37919  5 23  2014 Intlab_V9/intval/verifylss.m

しかしコメントは本質的には違っていないようである。 V4.1.2 の段階のプログラムは読みやすそうである。 解読する場合はそれから読んだらどうだろうか？

ハイブリッドな(密行列係数)連立1次方程式ソルバー

最初の段階は、近似解の誤差を包み込む修正を施した Krawcyzk 作用素に基づく。 Krawcyzk 作用素については

R. Krawczyk: Newton-Algorithmen zur Bestimmung von Nullstellen mit Fehlerschranken, Computing 4, 187-201, 1969.

R. E. Moore: A Test for Existence of Solutions for Non-Linear Systems, SIAM J. Numer. Anal. 4, 611-615, 1977

を見よ (Krawczyk [2], Moore [18])。

反復スキームと、イプシロン・インフレーションについては、Rump による次 の論文で与えられている (Rump [29], [22])。

S. M. Rump: Kleine Fehlerschranken bei Matrixproblemen, Ph.D. Thesis, University of Karlsruhe, 1980. (これは入手できてない。)

S. M. Rump: Verification methods for dense and sparse systems of equations, in : J. Herzberger (ed.), Topics in Validated Computations - Studies in Computational Mathematics, Elsevier, Amsterdam, 63-136, 1994.

疎行列係数は s.p.d (対称正定値) であるとする。そうでない場合は、通常の 方程式として扱う (疎行列は、係数行列が対称正定値の場合はそれを生かした 計算をするが、そうでない場合は特別扱いはしない、という意味かな？「通常 の方程式」は「正規方程式」かな？)。対称な場合は

  X = verifylss(A,b,1);

と呼び出さねばならない。これは条件数を２乗にすることに注意 ³。

7回反復しても成功しない場合は、第2段階。それは

A. Neumaier: A simple derivation of the Hansen-Bliek-Rohn-Ning-Kearfott enclosure for linear interval equations, Reliable Computing 5 (1999), 131-136, and Erratum, Reliable Computing 6 (2000), 227. (入手済み)

による (Neumaier [30], [31])。

Krawcyzk 作用素に基づく第1段階を飛ばして、直接第2段階から始めること も出来る。そうするには intvalinit('LssSecondStage') とする。 詳しいことは intvalinit を見よ⁴。

第2段階の結果は、非常に悪条件の連立1次方程式については(第1段階)より 良い品質となる可能性があり、通常は同程度の品質である。

第2段階は、第1段階よりも少し遅い。800MHzのラップトップでの時間計測 (timing) は次のようになる⁵。

t1 = 第1段階のみ, t2 = 第2段階のみ

	A and b point data			A and b point data
n	t1	t2	t2/t1	t1	t2	t2/t1
100	0.04	0.05	1.33	0.04	0.05	1.31
200	0.17	0.25	1.49	0.23	0.31	1.35
500	2.1	3.2	1.52	2.8	3.9	1.40
1000	13.7	21.6	1.58	18.0	26.3	1.46

疎行列係数に対するソルバーは

S. M. Rump: Validated Solution of Large Linear Systems, Computing Supplementum 9, 191-212, 1993

(Rump [32]) に基づく。 現時点での実装は正値対称行列のみ⁶。

intvalinit('ImprovedResidual') とすると、 点行列について残差改良をする⁷。

intvalinit('QuadrupleResidual') とすると、 4倍精度演算を用いた残差改良をする。

解釈(インタープリテーション)のオーバーヘッドのため、これは、密な点行列 に対してのみ行われる⁸。

追加で必要な計算時間は以下の表を見よ。やはり 800MHz ラップトップ、点行 列と点ベクトル、相対誤差は...⁹

time[sec] t1 = doubleresidual, t2 = improvedresidual, t3 = quadrupleresidual max.rel.error r1 = doubleresidual, r2 = improvedresidual, r3 = quadrupleresidual

行列の条件数 1e10:

n	t1	t2	t3	t2/t1	t3/t1	t1	r2	r3
100	0.03	0.03	0.07	1.20	2.63	1.6e-05	6.0e-11	2.2e-16
200	0.13	0.17	0.37	1.34	2.85	4.6e-05	1.5e-10	2.2e-16
500	1.7	2.2	4.0	1.25	2.28	1.4e-04	6.7e-10	2.2e-16

行列の条件数 1e15:

n	t1	t2	t3	t2/t1	t3/t1	t1	r2	r3
100	0.04	0.05	0.17	1.21	3.71	2.8e-00	1.3e-05	1.4e-15
200	0.27	0.30	0.45	1.11	1.66	2.8e-00	1.3e-05	1.4e-15
500	3.2	3.5	4.4	1.11	1.39	2.8e-00	1.3e-05	1.4e-15

(時間のことはさておき、非常にきれいな結果と思う。 オプション1つで切り替えられるのは素晴らしい。)

第3段階の、過剰な、あるいは不定な連立1次方程式については

S. M. Rump: Solving Algebraic Systems with High Accuracy, in "A New Approach to Scientific Computation", eds. U. Kulisch and W. Miranker, Academic Press, 51-120, 1983.

を見よ (Rump [33])。