線形代数演算
シンボリック ヒルベルト行列
次の例は、シンボリック バージョンの 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
x
と y
の差分を計算すると、2 つの解に微小な差があることがわかります。
x-y
ans = 8.0779356694631608874e-27 -6.4623485355705287099e-26 7.1085833891275815809e-26
vpa
と digits(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]
N
が Hs
のヌル空間にあることを確認します。
Hs*N
ans = 0 0 0