Main Content

線形代数演算

シンボリック ヒルベルト行列

次の例は、シンボリック バージョンの 3 行 3 列ヒルベルト行列に基づく、基本的な線形代数演算を実行する方法を示します。

3 行 3 列のヒルベルト行列を作成します。format short を指定した場合、MATLAB® の出力は以下のとおりです。

H = hilb(3)
H =
    1.0000    0.5000    0.3333
    0.5000    0.3333    0.2500
    0.3333    0.2500    0.2000

計算された H の要素は小さな整数の比で表される浮動小数点数です。H はクラス double の MATLAB 配列です。

H をシンボリック行列に変換します。

H = sym(H)
H =
[   1, 1/2, 1/3]
[ 1/2, 1/3, 1/4]
[ 1/3, 1/4, 1/5]

シンボリック線形代数演算

H に関するシンボリック演算では、浮動小数点近似 hilb(3) ではなく、無限精度のヒルベルト行列 sym(hilb(3)) に対応する結果が得られます。

H の逆行列を求めます。

inv(H)
ans =
[   9,  -36,   30]
[ -36,  192, -180]
[  30, -180,  180]

H の行列式を求めます。

det(H)
ans =
1/2160

線形方程式系を解くために、バックスラッシュ演算子を使えます。たとえば、H*x = b の解を求めます。

b = [1; 1; 1];
x = H\b
 x =
  3
 -24
  30

これら 3 つの結果 (逆行列、行列式、線形方程式の解) は、無限精度の有理数ヒルベルト行列に対する厳密な結果です。

可変精度の演算

前の演算と、20 桁の精度の可変精度算術演算を対比します。

digits(20)
V = vpa(H)
V =
 
[                    1.0,                    0.5, 0.33333333333333333333]
[                    0.5, 0.33333333333333333333,                   0.25]
[ 0.33333333333333333333,                   0.25,                    0.2]

個々の要素の小数点表現は、MATLAB が可変精度の算術演算を使用していることを示します。各演算の結果は有効小数桁数 20 で丸められます。

行列の逆行列を求めますが、誤差は行列の条件数 (hilb(3) の場合では約 500) だけ拡大されることに注意してください。

cond(V)
ans =
 
524.0567775860608

無限精度と可変精度の逆行列の差を計算します。

ih = inv(H)
ih =
 
[   9,  -36,   30]
[ -36,  192, -180]
[  30, -180,  180]
iv = inv(V)
iv =
 
[   9.0,  -36.0,   30.0]
[ -36.0,  192.0, -180.0]
[  30.0, -180.0,  180.0]

これらの行列は同一であるように見えますが、差分を計算すると同じではないことがわかります。

dhv = ih - iv
dhv =
 
[ -5.4929962552349494034e-26,  2.4556924435168009098e-25, -2.1971985020939797614e-25]
[  2.4556924435168009098e-25, -1.2666203129718236271e-24,  1.1373733422604130529e-24]
[ -2.1971985020939797614e-25,  1.1373733422604130529e-24, -1.0856745539758488233e-24]

方程式 V*y = b を解きます。解は H*x = b の解と同じように見えます。

y = V\b
y =
 
   3.0
 -24.0
  30.0

xy の差分を計算すると、2 つの解に微小な差があることがわかります。

x-y
ans =
 
  8.0779356694631608874e-27
 -6.4623485355705287099e-26
  7.1085833891275815809e-26

vpadigits(16) を使用すると、標準倍精度の MATLAB ルーチンを使用した場合と同等の精度が得られます。

特異値のシンボリックな調査

H が特異行列になる H(1,1) の値 s を求めます。

syms s
Hs = H;
Hs(1,1) = s
Z = det(Hs)
sol = solve(Z)
Hs =
[   s, 1/2, 1/3]
[ 1/2, 1/3, 1/4]
[ 1/3, 1/4, 1/5]

Z =
s/240 - 1/270

sol = 
8/9

s の解を Hs に代入します。

Hs = subs(Hs, s, sol)
Hs =
[ 8/9, 1/2, 1/3]
[ 1/2, 1/3, 1/4]
[ 1/3, 1/4, 1/5]

Hs の行列式がゼロであることを確認します。

det(Hs)
ans =
0

Hs のヌル空間と列空間を求めます。両方の空間は自明ではありません。

N = null(Hs)
C = colspace(Hs)
N=
3/10
 -6/5
    1

C =
[     1,     0]
[     0,     1]
[ -3/10,   6/5]

NHs のヌル空間にあることを確認します。

Hs*N
ans =
 
 0
 0
 0