は～い。ともちゃんです。半年ぶりの登場ですよ。まずは、第97回の『行列式（1）』の5.で取り上げた行列式の基本的な性質を復習を兼ねて、列挙しておくわね。1.はじめに
2.連立一次方程式
（1）行列式の基本的性質
（2）Cramer(クラメル）の公式
① クラメルまたはクラーメル
② 2元連立一次方程式
③ クラメルの公式の証明
④ 一般の連立一次方程式
⑤ クラメルの公式利用上の注意
（3）クラメルの公式に基づくユーザー定義関数
ⓞ DERIVEでの連立一次方程式の解法
① Dkを求めるユーザー定義関数
② Fdk_user関数の注意点
③ Fdk_user関数の実行結果
④ クラメルの公式に基づくユーザー定義関数
⑤ Fcramer_user関数の注意点
⑥ Fcramer_user関数の実行結果
⑦ Random_Matrix関数
（4）値がゼロとなる行列式
① f(k,j)=p(k)、f(k,j)=p(j)
② f(k,j)=p(k)+q(j)
③ f(k,j)=(k+j)^m
④ f(k,j)=k、jの多項式
（5）Hilbert(ヒルベルト)行列と連立一次方程式
① Hilbert(ヒルベルト)行列
② ヒルベルト行列が係数行列の場合
（6）Row_Reduce関数
① ヒルベルト行列が5次の場合
② ヒルベルト行列が10次の場合
（7）ヒルベルト行列の場合にクラメルの公式を使う
① ヒルベルト行列が10次の場合
今回は、第97回の『行列式（1）』の続きです。連立一次方程式のCramer（クラメルまたはクラーメル）の公式を中心とした基本的な内容です。では、ともちゃん、どうぞ！
目次へ戻る
は～い。ともちゃんです。半年ぶりの登場ですよ。まずは、第97回の『行列式（1）』の5.で取り上げた行列式の基本的な性質を復習を兼ねて、列挙しておくわね。
それが、こちらの7つ。
（1）転置行列の行列式の値（誤解が生じない場合、以下『値』は、略す）は、変わらない、
（2）行(列)の交換で行列式の符号が反転する、
（3）二つの行(列)が等しい行列式はゼロとなる、
（4）一つの行(列)の要素を定数倍で行列式も定数倍となる、
（5）複数の行列式の和に分解できることがある、
これは、たとえば、
あるいは、
のような分解のことね。
さらには、第97回では、取り上げていないけど、
なども分解です。
aで始まる要素をA、bで始まる要素をBで模式的に表すと、
と表現すると分かりやすいかしら。
下図のような、3次行列式では、
aで始まる要素をA、bで始まる要素をB、cで始まる要素をCで模式的に表すと、
と表現できるでしょう。3×3×3×=27個の行列式の和に分解されます。念のため、DERIVEで検算して、確認しました。
（6）行(列)×定数を他の行(列)に加えても行列式の値は、変わらない、
（7）行(列)の1つ以外の要素が0の場合は、行列式の次数を1つ下げられる、
具体的には、3次行列式、|A|が下図であれば、
のように2次行列式に変形できるということ。
目次へ戻る
ウィキペディアによれば、『クラメルの公式の名は、ガブリエル・クラーメル (1704～1752) に因むもので、クラーメルは任意個の未知数に関する法則を1750年に記している。なお特別の場合に限れば、コリン・マクローリンが1748年に公表している（また、恐らくはそれを1729年ごろにはすでに知っていたと思われる。）』ということで、18世紀のスイスの数学者の名前にちなむのね。なお、Cramerのカタカナ表記は、様々なサイトで揺れており、クラメル、クラーメル、クラメール、とばらけちゃっている。
『ああ！Cramer、あなたを何て、呼んだらいいの？』
本稿では、『クラメル』と呼ぶことにしますよ。
目次へ戻る
未知数が2個の場合、その間の関係式が次のような一次方程式であれば、
a11 x1+a12 x2＝b1、-（1）
a21 x1+a22 x2＝b2、-（2）
となって、係数行列Aを下図のように表して、
また、未知数ベクトルX、定数ベクトルBを下図のように定めると、
方程式は、AX=B、と簡潔に表される。（以下、X、Bをx、b、と書く場合もあります。）
さてと、中学校に戻って考えると、x1を求めるとき、x2を消去するために、
（1）式×a22－（2）式×a12を作ると、
(a11×a22－a21×a12)×x1=(b1×a22－b2×a12)、
これから、x1は、その係数がゼロでなければ、求められるわ。
だけど、よっく見ると、x1の係数、(a11×a22－a21×a12)は、Aの行列式、|A|、そのものなのよ。
また、右辺もAの1列目を列ベクトルBで置き換えた行列の行列式に等しいことに気がつけば、x1は、
と書き表せることが分かるでしょう。
x2も同様に、
（1）式×a21－（2）式×a11を作れば、
(a12×a21－a22×a11)×x2=(b1×a21－b2×a11)、
ここで、x2の係数が|A|となるように項の順序を入れ替えるため、－1を両辺に掛けると
(a22×a11－a12×a21)×x2=(b2×a11－b1×a21)、となるので、
x1のときと同様に、x2の係数は、Aの行列式|A|、そのものであり、右辺は、Aの2列目を列ベクトルBで置き換えた行列の行列式に等しいことが分かるので、x2は、
と書き表される。これが、2元連立一次方程式の解を求める『クラメルの公式』ね。
目次へ戻る
未知数が3個の場合、その間の関係式が次のような一次方程式であれば、
a11 x1+a12 x2+a13 x3＝b1
a21 x1+a22 x2+a23 x3＝b2
a31 x1+a32 x2+a33 x3＝b3
2元の場合と同様の記号を使って、係数行列Aを下図のように記載する。
また、3次の列ベクトルとして、
X＝［x1，x2，x3］’ を未知数ベクトル、B＝［b1，b2，b3］’ を定数ベクトルとおいたとき、上の3元連立一次方程式は、
A X＝B、と簡潔に書くことができるのは、まったく同じね。（行（横）ベクトルの右側のアポストロフィーは、転置を表す）
そこで、2元のクラメルの公式を拡張して、3元の場合も、x1=|D1|/|A|とクラメルの公式により、解を求められることを証明してみるわよ。ここで、D1は、Aの1列目をBで置き換えた行列とします。すなわち、
従って、|D1|は、
となる。
このx1=|D1|/|A|が成り立つことを証明するために、|D1|を1列目について、展開して考えるんだけど、その際、b1～b3は、元の方程式、
a11 x1+a12 x2+a13 x3＝b1
a21 x1+a22 x2+a23 x3＝b2
a31 x1+a32 x2+a33 x3＝b3
の左辺を代入すると、
こうすると、明らかに、2.（1）で取り上げた行列式の基本的性質の（5）が使えるため、
と分解できる。
ここで、基本的性質（4）の『一つの行(列)の要素を定数倍で行列式も定数倍』を使えば、
さらに、x2とx3の係数の行列式は、□で囲んだ列同士が等しいので、基本的性質（3）の『二つの行(列)が等しい行列式はゼロ』を使えば、これらの2つの項は、いずれもゼロとなる。よって、|D1|は、
すなわち、|D1|=x1×|A|、である。
これから、|A|<>0 のときは、x1=|D1|/|A|であることが証明された。
|D2|、|D3|についても、同様に証明できることは、容易に分かるわ。もちろん、『行列と行列式』（古屋茂 著：培風館：1970年増補第2刷）に記載されているように、一般のn元の場合を証明するのが本当だけど、原理は、上とまったく同じ。だけど、Σ記号が多数現れるので、かえって、分かりにくいでしょう。
目次へ戻る
今回の行列式を取り上げた動機が『非対称道中双六』の問題に端を発しているので、n次の正方行列（|A|<>0、の場合）に興味がある訳ね。だけど、一般の連立一次方程式では、『不定』（解が複数個のパラメーターを含む。方程式の数が少ない場合）や『不能』（方程式を満たす解がない場合）も現れる。このようなことをハッキリと説明するためには、行列の『ランク（位）』を考える必要があります。ランクは、m行n列の係数行列の『小行列式』のうち、ゼロでない最大の次数と定義されます。ここでは、これ以上、立ち入らないことにしましょう。
目次へ戻る
クラメルの公式は、連立一次方程式（A X=B）の解、未知数ベクトル、Xが係数行列 Aと定数ベクトル Bにどのように依存しているかを陽に示す重要な公式ね。これを使えば、Xk=|Dk|/|A|、k=1～n。と明確に解を表示できます。ただし、Dkは、Aのk列目をBで置き換えた行列の行列式、|A|<>0。
とは言うものの、計算上は、|A|や|Dk|という行列式の値を求めるのに時間とメモリーを浪費する訳よ。つまり、係数行列の要素が相異なる記号変数の場合、DERIVEで計算できる行列式の次数は、7～8次元が限界なの。（詳しくは、後述）
一方、係数行列の要素がすべて概数（DERIVEの標準的な精度では10桁）であれば、制約は、相当、緩やかになるのね。それでも、100次の行列式の近似計算には、約30秒が必要。これは、普通、待てる限界に近い時間でしょ。
さらに、100元の連立一次方程式の解をクラメルの公式で求めるためには、|A|は、1回計算すれば、良いけれど、|Dk|は、100回計算しないといけないので、結局、1個の行列式の計算に30秒を要するとして、全部で100+1回で、約30秒×101≒約1時間かかることになります。これでは、時間がかかりすぎるでしょう。
※ 1度だけで済めば、これでも良いケースもあるでしょう。手計算の頃であれば、100元の連立一次方程式を解くのに1時間程度とは、すごく速いと驚かれたでしょう。
目次へ戻る
DERIVEには、A X=Bの形の連立一次方程式を解く方法がすでに用意されています。それは、主として、以下の3つです。
・ 逆行列を使う方法（係数行列が正方行列で行列式が非ゼロ（正則）、かつ、次数の小さな場合）、
X=A^-1×Bとして求める。ただし、DERIVEでは、要素が記号変数では、5次元程度まで。
要素が数値、特に概数であれば、大きい次数も可能。
・ Solve関数を使う方法（連立方程式以外でも、おなじみの方法）、
DERIVEのメニューの解くから、『連立方程式』を選択して進む、あるいは、同等な方法として、
Solve(［式1，式2，・・］，［未知数1，未知数2，・・］)を厳密に又は近似的に計算、
すでに定義済みの係数行列A、定数列ベクトルB、未知数列ベクトルXがあれば、Solve(A X=B，X)、とする。
※ 逆行列を使う方法のように正方行列である必要は無い。方程式の数と未知数の数が異なる場合は、利用者が表示された解をもとに判断する必要が出てきます。
・ Row_Reduce 関数を使う方法（後述）、
目次へ戻る
ここでは、クラメルの公式に基づき計算を行うDERIVEのユーザー定義関数を作ってみましょう。これは、クラメルの公式の理解を深めるためのものです。
まず、DERIVE上で、Dk（Aのk列を列ベクトル、Bで置き換えた行列をDkと書く）をどのように求めるかを考えてみましょう。シンプルに考えて、関数内で使うローカル変数として、行列、D_:=A_と置いた後に、D_のk行j列の要素、D_↓k↓j:=B↓k、j=1～n、と置き換えるのが簡単でしょう。この考えのもと、Dkを与えるユーザー定義関数、Fdk_user(A_,v_,k_)は、A_を係数行列、v_を定数ベクトル（行ベクトル）としたとき、次のように表されます。
Fdk_user(A_, v_, k_)
:=PROG(
d_ ≔ A_,
n_ ≔ DIMENSION(v_),
VECTOR(d_↓j_↓k_ ≔ v_↓j_, j_, 1, n_),
d_
)
※ 関数名の末尾に付ける『_user』は、ユーザー定義関数を表すマイルールです。
（DERIVEの画面コピーは下図）
目次へ戻る
・Prog( )関数は、第69回『方程式の数値解法（4）（ハレー法、DERIVEのPROGとRETURN）』で説明しているけど、( )内のカンマで区切られた手順を順次実行する関数ですよ。（Prog関数内にProg関数を使うなどの入れ子も可能なので、チョー便利ね）
・アンダーバー付き変数は、第97回でも触れているけど、関数内部でのみ利用するローカル変数名に、A_やn_のようにアンダーバーを付けて区別するのがマイルールなのよ。（DERIVEは、ローカル変数とグローバル変数の区別がないので、誤って、外部でローカル変数の値を変えないための用心）
・変数に値をセットする際は、変数名＝値、ではなくて、変数名：＝値、と書く必要があります。
・行列やベクトルの要素の値を読み出す（参照する）際、たとえば、行列の場合、関数を使うと、ELEMENT(行列名，行番号，列番号)であり、添え字を使うと、行列名↓行番号↓列番号、です。この両者に本質的な違いはありません。
しかし、行列やベクトルの要素の値を更新する際は、添え字を使う方法のみが許されます。例えば、行列であれば、ELEMENT(行列名，行番号，列番号)：＝値、と書くと文法エラーとなります。このエラーは、たとえば、三角関数で説明すると、変数：＝Sin(π/6)、としたとき、左辺の変数の値は、0.5となりますが、逆に、Sin(θ)：＝0.5、としても、θの値をπ/6に更新できないことを考えれば、当然のことです。ただ、案外とうっかりしがちなので気をつけましょう。
目次へ戻る
係数行列Aが下図であり、定数ベクトルを［b1，b2，b3］としたとき、
D1を求めるためには、Fdk_user(A，［b1，b2，b3］，1)として、
が得られました。確かに、Aの1列目が［b1，b2，b3］’によって、置き換えられています。
k=2の場合の実行結果、D2は、
同様に、D3は、
となります。
目次へ戻る
Fcramer_user(A_，v_)は、クラメルの公式に基づいて、連立一次方程式の解を求めるユーザー定義関数です。ここで、A_は、係数行列、v_は、行ベクトルで定数ベクトルを表す。
Fcramer_user(A_, v_)
≔ PROG(
n_ ≔ DIMENSION(v_),
w_ ≔ DET(A_),
IF(w_ = 0, RETURN 0),
VECTOR(DET(Fdk_user(A_, v_, k_))/w_, k_, 1, n_)
)
DERIVEの画面コピーは、下図の通り。
目次へ戻る
・第69回で既出のReturn手続きがIf関数内で使われています。これは、|A_|がゼロの場合は、0を返してProg関数を飛び出し、以降のVector関数を実行しないためです。
上の関数内部で引用されている前出のFdk_user関数は、関数の外部で定義されているとしています。なお、Fcramer_user関数内で定義したい場合は、第69回で説明しているように関数内部で定義する関数の記法と引数に対して、DERIVE特有の記法が必要ですので、ご注意ください。ここでは、そのようなわかりにくさを避けて、Fdk_user関数がFcramer_user関数の外で定義されているとしました。
目次へ戻る
3元連立一次方程式に対するFcramer_user関数の実行結果は、下図の通りです。
計算結果は、2行目以降が［x1，x2，x3］と行ベクトル形式でXを返していることが分かります。
ただし、2.（3）の冒頭に書いたように、要素がすべて記号変数の行列式は、次元が5，6までは、それぞれ、0.02秒程度で求められます。7となると、|A|の計算だけで、約2秒を要し、8では、現実的な時間内で結果が返ってきませんでした。
しかし、8次元の行列式の値も、1行目で余因子関数を利用して展開すれば、
|A|=Σ(k=1～8)(ELEMENT(A,1,k)×Cofactor(A,1,k)、
あるいは、同値の方法として、Σ(k=1～8)(A↓1↓k)×Cofactor(A,1,k)、と計算することにより、約80秒で|A|を求められました。これらから、要素が相異なる記号変数からなる行列の行列式は、7次元程度がDERIVEの計算限界であることが分かります。
※ Cofactor関数は、第97回『行列式（1）』で紹介したDERIVEの余因子関数（行列，i行，j列）：行列からi行j列を削除した行列の行列式×(－1)^(i+j)である。（スカラー量。符号が付いている点に注意）
目次へ戻る
これまでは、行列式の要素がすべて記号変数の場合ね。ここでは、要素がすべて数値である行列式の（DERIVEでの）計算限界を調べてみましょう。でも、DERIVE画面から行列の要素をいちいち手入力するのは、面倒だわ。そこで、第67回「方程式の数値解法（2）（Newton法、逐次近似法）」で紹介した、Random_Matrix関数を使ってみましょう。
RANDOM_MATRIX(最大行数，最大列数，整数)、であり、最大行数×最大列数の行列（今回は、正方行列を考えるので第1引数と第2引数は同一。正方行列以外の行列や縦ベクトル、横ベクトルも可能）が作成され、その要素は、第3引数を1としたときは、－1＜要素＜1の有理数（近似計算を指示するとDERIVEの精度＝通常、有効数字10桁の数値となる）、
一方、2以上の整数を指定した場合は、整数であり、その範囲は、－整数＜要素＜整数となる。なお、第3引数に1.5などの小数点を含んだ数字を与えると意味のない結果が現れるので注意。言葉で説明すると、上のように、ちょっと、ややこしいので、以下で例を示しますよ。
例1：第1，第2引数を5、第3引数を1とするとき、RANDOM_MATRIX(5，5，1)は、厳密な計算では、下図となります。
すべての要素が、－1＜要素＜1、の有理数となっています。
ここで、近似計算を指示（あるいは、最初から近似計算を指示）すると、
と10桁の有効数字を持つ小数で表示され、その行列式の値は、(0.1076416623)となりました。
例2：第1，第2引数を10、第3引数を10とするとき、RANDOM_MATRIX(10，10，10)は、下図の結果となりました。
第1，第2引数が10なので、10次行列であり、その要素は、第3引数が2以上の整数（10）のため、－10＜要素＜10、の整数となっているのが分かるかしら。ちなみに、この行列の行列式=3610597193、とほぼ瞬時に求められました。
例3：第1，第2引数を100、第3引数を10とすると、第1，第2引数が100なので、100次行列であり、その要素は、第3引数が10のため、－10＜要素＜10、の整数となる。大きすぎるので、貼り付けませんが、0.1秒程度で生成されました。しかし、その行列式は、
－13310623491144754619940801589（次行に続く）
5176734344984505000170751441012338161302731688806（次行に続く）
1847416620921378165598528776491231543007725185919（次行に続く）
4460747560821579927008852
≒－1.331062349×10^151、と途方もない大きさの数が約11秒かかって求められたわ。
※ Random_matrix関数を利用しているので、その値は、計算の度に変わってしまう。いずれも、一例。
目次へ戻る
要素が規則的な行列の行列式は、意外にゼロとなるので、驚くわよ。そのいくつかを挙げてみましょう。なお、以下では、行列Aのk,j 要素がf(k,j)であるとしますね。どのような関数だとDET(A)がゼロとなるのでしょうか？
f(k,j)=p(k)、f(k,j)=p(j)、ここで、pは、任意関数。これらは、いずれもゼロとなります。
たとえば、f(k,j)=sin(k)のようなもの。
このDET(A)=0、
これは、kのみに依存する要素だと、同一の値を持つ複数の行が現れるから、行列式の値は、ゼロとなる。同様に、jのみに依存する要素だと、同一の値を持つ列が現れるので、ゼロとなる。
なお、上の三角関数などでは、行列式の次数が5程度となると、厳密計算でゼロとなることを検証することが難しくなる。その際は、sin(1)やsin(2)などを、仮の変数、c1、c2などで置き換えると7次、8次程度までは、検証可能。
目次へ戻る
f(k,j)=p(k)+q(j)、ここで、p、qは、任意関数。行列式の次数が3次以上でゼロとなります。
3次の場合にゼロとなる理由は、下図の通りです。
2行目－1行目、3行目－1行目で、それぞれの行の共通因数をくくり出すと、同一の値を持つ2つの行ができるので、行列式の値は、ゼロとなります。4次以上でも、証明は、同じですね。
目次へ戻る
f(k,j)=(k+j)^m、のとき。mを正整数として、行列式の次数がm+2以上であれば、行列式は、ゼロとなります。
m=3の場合で、5次行列式は、下図。
この行列式は、ゼロとなります。
m=5では、行列式の値のみを挙げるとして、行列式の次数を1～10まで変化させて、
([32, -26281, -35402880, 18231008400, 1183818240000, -2985984000000, 0, 0, 0, 0])
と行列式の値は、変化します。最初にゼロとなったのは、数が7次の時ね。すなわち、5+2=7。
m=10では、今度は、1～15次元まで計算すれば、
最初に行列式の値がゼロとなるのは、12次元の時。やはり、10+2ということ。
目次へ戻る
この場合も、行列式の次数が大きくなれば、ゼロとなります。以下では、行列式の変形をするのだけど、DERIVEでは、行列式の演算を行うと、最終的な計算結果が返ってくるので、行列に対する操作として、説明していますよ。（2021/10/5に修正しました）
n次行列の要素を、a(k,j)として、n行目－n-1行目で、a’(n,j)=a(n,j)－a(n－1,j)=Δa(n－1,j)、
ここで、a'は、新しい行列の要素を示す。また、以下では、jは、記述に無関係なので略して、a(1)、a(2)などと記すわよ。
ついでに、Δは、（前進）差分演算子で、Δa(k)=a(k+1)－a(k)、と行番号に働く。（Δの定義を第82回に合わせました）
同様に、n－1行目－n－2行目で、a'(n－1)=a(n－1)－a(n－2)=Δa(n－1)、となる。
以下、・・・、最後に、2行目－1行目で、a'(2)=a(2)－a(1)=Δa(1)、となり、
一般に、a'(k)=Δa(k－1)、k=n～2、ただし、a'(1)=a(1)
今度は、新行列のn行目～3行目について、同様のことを行うと、新新行列のk行目は、
a''(k)=a'(k)－a'(k－1)=Δa(k－1)－Δa(k－2)=Δ（a(k－1)－a(k－2)）=Δ^2 a(k－2)、k=n～3、
ただし、a''(2)=a'(2)=Δa(1)、a''(1)=a(1)、である。
同様に、新新行列のn行目～4行目について、同様のことを行うと、新新新行列のk行目は、
a'''(k)=a''(k)－a''(k－1)=Δ^2 a(k－2)－Δ^2 a(k－2)=Δ^3 a(k－3)、k=n～4、
ただし、a''''(3)=Δ^2 a(1)、a'''(2)=Δa(1)、a'''(1)=a(1)、である。
このような操作を繰り返すと、最終的に、作られる行列のk行目の要素を、大文字で、A(k)と書けば、
A(k)=Δ^(k－1) a(1)、とすることができる。k=n～2，また、A(1)=a(1)である。
丁寧に書くと、4次の場合、下図のようにすることができます。（Δは、行のみに働く前進差分演算子）
このとき、a(k)がkのn－2次式であれば、A(n)=Δ^(n－1) a(1)は、ゼロとなるので、その行列の行列式は、ゼロとなります。（なぜならば、その行で行列式を展開すると、係数がすべてゼロだからなのね）、逆に言えば、要素が、kのn次式の場合、n+2次以上の行列式では、その値がゼロとなることが証明できたと思うわ。
例：(k+j)^2、k=1～4、j=1～4、の4次行列は、
4行－3行、3行－2行、2行－1行後は、
4行－3行、3行－2行後は、
ここまででも、よいのだけど、ダメ押しで、4行－3行を作ると、
4行目がすべてゼロとなり、この行列の行列式は、ゼロだと分かりました。
目次へ戻る
前節に見て来たように、規則的な要素を持つ行列の行列式は、ゼロとなるケースが案外と多い。一方、要素をk、jの分数式、1/(αk+βj+γ)、などや無理式、√(k+j)、で考えると、行列式は、一般には、ゼロとはならなくなります。
具体的な分数式の例として、α=β=1、γ=－1、とした場合が、有名な『Hilbert(ヒルベルト)行列』と言われるものとなるのね。すなわち、ウィキペディアの記法にならって、ヒルベルト行列を、H
と記載すれば、その（k，j）要素は、1/(k+j－1)、である。
※ ヒルベルトは、ドイツの著名な数学者。『【David Hilbert】 ドイツの数学者。数学の各分野にわたって業績が多く、基礎論では形式主義を唱えた。理論物理学の基礎に関する論文もある。（1862～1943）』（広辞苑第7版）
ヒルベルト行列のn=3と5の具体的な形を下図に示すわよ。
n=3
n=5
n_次のヒルベルト行列を与えるユーザー定義関数を作ったので、下に示すわね。
Fhi_user(n_) ≔ VECTOR(VECTOR(1/(k_ + j_ - 1), k_, 1, n_), j_, 1, n_)
ちなみに、n次のヒルベルト行列の行列式は、ウィキペディアの『ヒルベルト行列』の項によると、数列CnをCn=Π(k=1～n－1)(k！)、としたとき、n次のHの行列式、DET(H)=（Cn）^4／C2n、と表される。
上の例では、n=3では、DET(H)=1/2160、n=5では、DET(H)=1/266716800000、となり、いずれも、実際に計算してものと一致する。
n_次のヒルベルト行列の行列式を表すユーザー定義関数は、次のようね。
Fh_user(n_) ≔ Fc_user(n_)^4/Fc_user(2·n_)、
Fc_user(n_) ≔ ∏(k_!, k_, 1, n_ - 1)、
Fcは、先の数列Cnを表すユーザー定義関数、
※ ここでは、深く立ち入りませんが、Hの逆行列は、ウィキペディアによれば、次となる。
H^-1の（k，j）要素は、(－1)(k+j)×(k+j－1)×Comb(n+k－1，n－j)Comb(n+j－1，n－k)（Comb(k+j－2，k－1)）^2、
ここで、Comb(n，r)は、組み合わせを示すDERIVEの関数で、数学の nCrと同じ。
上の例では、n=3のHの逆行列は、
となり、先のHとの積は、下図のように単位行列となりました。
n＝5でも同様に、
上図の左側の行列が逆行列、右側がヒルベルト行列。
乗算の結果は、
と単位行列となります。
n_次のヒルベルト行列の逆行列を与えるユーザー定義関数は、次のようね。
Fhr_user(n_) ≔ VECTOR(VECTOR((-1)^(k_ + j_)·(k_ + j_ - 1)·COMB(n_ + k_ - 1, n_ - j_)·COMB(n_ + j_ - 1, n_ - k_)·COMB(k_ + j_ - 2, k_ - 1)^2, k_, 1, n_), j_, 1, n_)
目次へ戻る
『計算機による数値計算法』（一松信・戸川隼人 著：新曜社）（1975年2月）によれば、係数行列がヒルベルト行列における近似計算では、桁落ちが起こりやすいという。そこで、前節で定義した関数を使って、試してみましょう。なお、以下では、便宜的に、未知数ベクトルを（Xではなくて）Yで、また、計算をクラメルの公式ではなく、Solve関数を使っています。
同書の例にならって、n個の未知数がすべて1であるとき、定数ベクトルBの要素は、Bk=(k + n) - (k)、と表される。
（k=1～n）、ここで、ψは、ディガンマ関数であり、ψ(z+n)=ψ(z)+Σ（k=1～n）（1/(z+k－1)）、という漸化式を満たす。
ただし、ψ(1)=－γ（オイラー定数=0.57772・・）、（この表式は、DERIVEが作り出すもので、ここでは、ディガンマ関数やオイラー定数の知識は、必要ありません）
まずは、上記の書籍に出てくるn=3の場合を見てみましょう。
Aは、前に挙げた下図の通り。
定数ベクトルBは、
となるので、未知数ベクトルYを
Y ≔ [y1, y2, y3]、として、SOLVE(A·Y` = B`, Y)、を厳密に計算すると、
(y1 = 1 ∧ y2 = 1 ∧ y3 = 1)が得られる。
近似計算でも、この場合は、厳密解と同一の結果となる。
では、次に、10次のヒルベルト行列の場合をやってみるわよ。
A:=Fhi_user(10)
B:= VECTOR((k + 10) － (k), k, 1, 10)
Y:=[y1, y2, y3, y4, y5, y6, y7, y8, y9, y10]
このとき、SOLVE(A·Y` = B`, Y)、の厳密解は、
(y1 = 1 ∧ y2 = 1 ∧ y3 = 1 ∧ y4 = 1 ∧ y5 = 1 ∧ y6 = 1 ∧ y7 = 1 ∧ y8 = 1
∧ y9 = 1 ∧ y10 = 1)
となり、当然ながら、真の値と一致する。
ところが、計算精度は、10桁で、近似解を計算すると、
(y1 = 0.9999996778 ∧ y2 = 1.000029303 ∧ y3 = 0.9993475278 ∧ y4 = 1.006156436
∧ y5 = 0.9697051614 ∧ y6 = 1.085500923 ∧ y7 = 0.8565380387 ∧ y8 = 1.141339867 ∧ y9 = 0.9245444917 ∧ y10 = 1.016838935)、
となり、いずれも真の値と大きな違いが生じており、特に、赤字で示した、y7やy8では、0.14程度の大きさの絶対誤差が生じていることが分かる。10桁の精度で計算しているのだから、少なくとも、5～6桁程度は、合っていると思うと大きな間違いであることがわかった。だが、この場合でも、計算精度を20桁にすれば、近似計算でも真の値と一致する結果が得られる。
なお、上のように分数を小数にして計算すると小数に丸める際の誤差も含まれることになる。そこで、同書の注にあるように、AとBを、360360×2×17×19倍（※）して、すべての要素を整数の行列、ベクトルにしてから、10桁の精度で計算すると、以下のようになった。
(y1 = 0.9999998216 ∧ y2 = 1.000014547 ∧ y3 = 0.9997007865 ∧ y4 = 1.002657605
∧ y5 = 0.9875305203 ∧ y6 = 1.033864438 ∧ y7 = 0.9449526614 ∧ y8 = 1.052810026 ∧ y9 = 0.9724372161 ∧ y10 = 1.006032477)
これを見ると、最大の絶対誤差は、0.06程度で収まっている。決して良くはないが、上述の場合の0.14と比較すると絶対誤差が半分以下に収まった。
※ 360360=2^3·3^2·5·7·11·13、は、8次までのヒルベルト行列では、この整数をかけてから計算するとよいとある。今回の10次では、360360倍だけでは、整数化できないので、さらに、17×19×2、を掛けてすべて整数になるようにしてみた。
同書の脚注によれば、『ヒルベルト行列も3次くらいであるならば、10進3桁といった極端に短い桁数にしない限り、それほど悪条件ではない。しかし、n次のヒルベルト行列を10進n桁で機械的に計算すると、ほぼ、同様の現象を生ずることが知られている。見方によっては、これは、敏感な問題であって、係数の1/3を0.333333、などと小数に丸めたために、大きな誤差が生ずるとも解釈できる。また、一般に残差が小さいことと、解が真の値に近いこととは必ずしも一致せず、どちらが「よりよい」解であるかは、目的によって何ともいえないことにも注意しておく』とある。
目次へ戻る
DERIVEでは、連立一次方程式の解法として、2.（3）ⓞで簡単に説明したように、逆行列による方法、Solve関数を使う方法、と並んで、Row_Reduce関数を使う方法があります。連立一次方程式の進んだ数値解法については、別の機会に取り上げることにして、ここでは、Row_Reduce関数の使い方を説明しておくわね。
係数行列が5次のヒルベルト行列、A、未知数がすべて1の場合を取り上げます。下図。
この場合の定数ベクトルBは、次となる。
B:=[137/60, 29/20, 153/140, 743/840, 1879/2520]
Row_reduce関数を使って、連立一次方程式、AX=Bを解く場合は、
Row_Reduce(A，B’)とする。（’は、転置。Bが上では、横ベクトルなので）
今の例では、
とすればよい。
結果は、下図のように返ってくる。
この5行6列の行列の6列目が未知数Xの値、6列目を除いた左側の5行5列が単位行列であれば、計算自体は、完了している。
上では、厳密解を求めたが、近似計算でも同一の結果となった。
目次へ戻る
では、10次のヒルベルト行列で近似計算を指示した場合、どうなるかというと、精度10桁計算では、
Solve関数を利用した場合より、絶対誤差が著しく大きいことに驚くわ。真の値が1に対して、最大6近い解が出ている。
精度20桁計算では、
となり、10桁と比較して、20桁の場合、その絶対誤差は、小数点10桁以下である。
一方、前に述べたように、A、Bをそれぞれ、360360*17*19*2倍してから、Row_Reduce関数を使って、計算させると、10桁の精度でも、
となって、
分数を整数化しない場合と比較して、かなりの改善となっている。絶対誤差は、最大でも0.03程度である。
ヒルベルト行列に限らず、係数行列や定数ベクトルにある分数は、できるだけ、整数化するようにしておくと、近似計算の精度が向上するという教訓がハッキリと分かる一つの例でしょ。
目次へ戻る
Solve関数を使った例は、すでに挙げている。Row_Reduce関数を使った場合、近似計算では、同じ計算精度では、Solve関数よりも、見劣りした計算結果が出ている。では、クラメルの公式を使った場合をあらためて検算して見ましょう。5次では、差がなかったので、ここでは、10次のヒルベルト行列で試みてみましょう。
なんども、記載するのは、気が引けるので、Aは、10次のヒルベルト行列、Bは、
B ≔ [7381/2520, 55991/27720, 44441/27720, 485333/360360, 420983/360360,
74587/72072, 134159/144144, 2074783/2450448, 634871/816816, 33464927/46558512]、
のままとする。厳密解は、未知数がすべて1となる。
このとき、クラメルの公式に基づく、ユーザー定義関数では、
Fcramer_user(A, B)として、厳密に計算すると、([1, 1, 1, 1, 1, 1, 1, 1, 1, 1])が返ります。
さて、問題は、近似計算の場合ですね。
計算結果は、
[-0.1008874372, 1.774981110, 0.9777582164, 1.239524089, -0.1634513357,
4.178135058, -4.187356000, 5.990579519, -1.609859812, 1.571999365]とRow_Reduce関数を使った場合と同程度の誤差が出て来る。
20桁の計算精度では、
[1.0000000000073500543, 1, 1.0000000000011492299, 1, 1.0000000000023309099,
0.99999999999107283725, 1.0000000000145579701, 0.99999999998589962139,
1.0000000000074176087, 0.99999999999836598176]、
となり、Row_reduce関数の場合と同様に著しく改善される。
A、Bを360360*17*19*2倍して、すべて整数としていた場合は、10桁の計算精度で、
([1.102105180, 1.098951350, 1.017689899, 1.015879184, 1.007051159, 0.9825853284, 1.029389561, 0.9711222145, 1.015385443, 0.9965727956])
やはり、大幅に改善できることが分かりますね。誤差の大きさは、0.03程度に収まる。誤差の大きさなどは、Row_Reduce関数と同程度となることが分かったかしら。
目次へ戻る
作成日：2021/10/1
若干の修正：2021/10/2
前進差分の定義を変更：2021/10/3
2.（4）⑤の証明を修正：2021/10/5
句読点の誤り等を微修正：2021/10/6
