ヒルベルト行列と逆行列

この例では、Symbolic Math Toolbox™ を使ってヒルベルト行列の逆行列を計算する方法を説明します。

定義 :ヒルベルト行列は要素が単位分数になる正方行列です。 $H_{i j} = \frac{1}{i + j - 1}$ たとえば、3x3 のヒルベルト行列は、 $H = [\begin{array}{ccc} 1 & \frac{1}{2} & \frac{1}{3} \\ \frac{1}{2} & \frac{1}{3} & \frac{1}{4} \\ \frac{1}{3} & \frac{1}{4} & \frac{1}{5} \end{array}]$ となります。

純粋な数値法では失敗する悪条件の行列の場合でも、シンボリック計算では正確な結果が得られます。

20 行 20 列の数値ヒルベルト行列を作成します。

H = hilb(20);

この行列の条件数を求めます。ヒルベルト行列は悪条件の行列です。つまり、大きな条件数をもち、こういった行列がほぼ特異であることを示しています。条件数の計算でも数値誤差が発生しやすいことに注意してください。

cond(H)

ans = 3.3066e+18

したがって、ヒルベルト行列の逆行列は数値的に不安定です。逆行列を計算する場合、H*inv(H) は単位行列または許容誤差範囲内の単位行列に近い行列を返さなければなりません。

はじめに、関数 inv を使用して H の逆行列を計算します。数値的に不安定なため警告がスローされます。

H*inv(H)

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  9.542396e-20.

ans = 20×20

    1.0000   -0.0000   -0.0012   -0.0012    0.0118   -0.1624   -0.7867    0.9817   -0.5909   -1.1786   -0.1111   -3.7477    4.3067    2.3698   -0.2680   -2.0770    0.4077    2.9125    2.2648    0.1082
    0.0871    1.0000   -0.0006   -0.0009   -0.0090   -0.5388    1.3683   -0.2444   -0.1039   -0.1231   -0.9150   -2.8117    6.5888   -0.1975    1.0052   -2.4624   -0.2692    5.2024    1.2478    0.0755
    0.2054   -0.3020    0.9994    0.0024    0.0061   -0.4914    1.2560   -0.5424    0.0382   -1.1778   -1.3460   -1.3399    5.6394   -0.1046    0.9009   -2.6941   -0.3658    5.2330    1.3704    0.1685
    0.3158   -0.8256    0.4401    1.0012    0.0370   -0.5177    1.6844   -0.0735    0.1373   -0.4658   -0.6256   -3.0880    7.6058   -2.1391   -0.0357   -3.9665   -0.8322    5.9076    0.9175    0.4910
    0.4028   -1.3337    0.9924   -0.0390    1.0449   -0.4168    1.0098   -1.5264    1.1591   -1.7210    0.4572   -2.4602    7.8277   -0.3915    0.5983   -2.6859   -0.2981    5.7396   -1.8456    0.3786
    0.4670   -1.7512    1.4662   -0.0101    0.0223    0.6950    1.5640   -1.0055    1.0727   -0.7692   -0.1271   -2.4243    7.0420   -1.7390    1.2946   -2.6481   -0.6277    4.8238   -0.6456    0.7012
    0.5129   -2.0689    1.7984    0.1538   -0.1460   -0.3917    2.7314   -1.5910    0.8051   -0.8093   -0.2009   -3.9021    8.9283   -2.9270    1.7746   -3.9919   -0.6623    3.9135   -0.7386    0.7756
    0.5446   -2.2991    1.9870    0.4589   -0.4710   -0.0614    0.9184    0.8377   -0.3169    0.0831   -1.5936   -0.2424    3.5135    0.2075    0.6568   -2.7485   -0.3657    2.5463    0.0261    0.0726
    0.5656   -2.4583    2.0538    0.8936   -0.9286    0.2678    0.2762   -0.1833    0.9221   -0.8345   -0.9206    0.6412    3.5666   -0.5856   -0.1836   -1.2843    0.8606    1.5268    2.0609   -0.4121
    0.5788   -2.5622    2.0255    1.4119   -1.3842    0.1638    1.2865   -0.7210    1.0717    0.3686    0.1248   -2.1556    5.7747   -1.0638    1.6798   -2.0859   -0.1640    0.9892   -0.4810    0.3520
      ⋮

ここで、ヒルベルト行列の精度を高める MATLAB® 関数 invhilb を使用します。この関数は 15 行 15 列までの厳密なヒルベルト行列の逆行列を見つけます。20 行 20 列のヒルベルト行列の場合、invhilb は逆行列の近似を求めます。

H*invhilb(20)

ans = 20×20
10¹⁰ ×

    0.0000   -0.0000    0.0000   -0.0000    0.0000   -0.0004    0.0013   -0.0037    0.0047   -0.2308    0.4019    0.2620   -0.4443   -7.5099    3.4753    3.2884   -1.1618    0.2301    0.0437   -0.0222
   -0.0000    0.0000   -0.0000    0.0000   -0.0000   -0.0001   -0.0009    0.0172   -0.0628    0.1251   -0.8975    2.7711   -2.8924   -1.4888   -2.5102    6.4897   -3.0581    0.7925   -0.0037    0.0037
    0.0000    0.0000    0.0000   -0.0000   -0.0000   -0.0001    0.0009    0.0042   -0.0303    0.0271   -0.1028    0.2862   -0.9267   -4.0597    1.8194    4.7727   -1.3255    0.4667   -0.0211   -0.0062
    0.0000    0.0000   -0.0000   -0.0000    0.0000   -0.0001    0.0004    0.0002   -0.0056    0.0526   -0.1136    0.3616   -1.6920   -3.7214    0.2709    4.4169   -1.1602    0.5541   -0.0024   -0.0014
   -0.0000    0.0000   -0.0000    0.0000   -0.0000   -0.0000   -0.0006    0.0068   -0.0394    0.1449   -0.7818    1.9314   -2.1273   -1.0745   -2.2988    4.6349   -1.8923    0.3793   -0.0435    0.0098
    0.0000   -0.0000    0.0000   -0.0000    0.0000   -0.0002    0.0014   -0.0028    0.0304   -0.1053   -0.0724   -0.2073   -0.3769   -5.0729    2.7552    2.3488   -0.5046    0.2240    0.0675   -0.0066
   -0.0000    0.0000   -0.0000    0.0000   -0.0000    0.0000   -0.0013    0.0101   -0.0520    0.1307   -0.7715    2.9297   -3.3374    1.4084   -3.8703    7.1053   -2.7578    0.8815   -0.1648    0.0037
    0.0000   -0.0000    0.0000   -0.0000   -0.0000   -0.0002    0.0018   -0.0040    0.0231   -0.0892    0.2000    0.2766   -0.8262   -4.7229    1.6513    2.1857   -0.4749   -0.1747    0.1257   -0.0082
   -0.0000    0.0000   -0.0000    0.0000   -0.0000    0.0000    0.0004    0.0144   -0.0513    0.1165   -0.7063    1.9776   -2.0602   -1.0763   -3.2081    3.8539   -2.4580    0.5675   -0.0189    0.0035
    0.0000    0.0000    0.0000   -0.0000    0.0000   -0.0002    0.0008   -0.0019    0.0027   -0.0167   -0.0497    0.8865   -0.6720   -1.2729    0.0571    2.3625   -1.2616    0.2443    0.0184   -0.0026
      ⋮

丸め誤差を回避するために、厳密なシンボリック計算を使用します。ここでは、シンボリックヒルベルト行列を作成します。

Hsym = sym(H)

Hsym = 
 $(\begin{array}{c} 1 & \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9} & \frac{1}{10} & \frac{1}{11} & \frac{1}{12} & \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} \\ \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9} & \frac{1}{10} & \frac{1}{11} & \frac{1}{12} & \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} \\ \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9} & \frac{1}{10} & \frac{1}{11} & \frac{1}{12} & \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} \\ \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9} & \frac{1}{10} & \frac{1}{11} & \frac{1}{12} & \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} \\ \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9} & \frac{1}{10} & \frac{1}{11} & \frac{1}{12} & \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} \\ \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9} & \frac{1}{10} & \frac{1}{11} & \frac{1}{12} & \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} \\ \frac{1}{7} & \frac{1}{8} & \frac{1}{9} & \frac{1}{10} & \frac{1}{11} & \frac{1}{12} & \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} & \frac{1}{26} \\ \frac{1}{8} & \frac{1}{9} & \frac{1}{10} & \frac{1}{11} & \frac{1}{12} & \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} & \frac{1}{26} & \frac{1}{27} \\ \frac{1}{9} & \frac{1}{10} & \frac{1}{11} & \frac{1}{12} & \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} & \frac{1}{26} & \frac{1}{27} & \frac{1}{28} \\ \frac{1}{10} & \frac{1}{11} & \frac{1}{12} & \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} & \frac{1}{26} & \frac{1}{27} & \frac{1}{28} & \frac{1}{29} \\ \frac{1}{11} & \frac{1}{12} & \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} & \frac{1}{26} & \frac{1}{27} & \frac{1}{28} & \frac{1}{29} & \frac{1}{30} \\ \frac{1}{12} & \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} & \frac{1}{26} & \frac{1}{27} & \frac{1}{28} & \frac{1}{29} & \frac{1}{30} & \frac{1}{31} \\ \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} & \frac{1}{26} & \frac{1}{27} & \frac{1}{28} & \frac{1}{29} & \frac{1}{30} & \frac{1}{31} & \frac{1}{32} \\ \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} & \frac{1}{26} & \frac{1}{27} & \frac{1}{28} & \frac{1}{29} & \frac{1}{30} & \frac{1}{31} & \frac{1}{32} & \frac{1}{33} \\ \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} & \frac{1}{26} & \frac{1}{27} & \frac{1}{28} & \frac{1}{29} & \frac{1}{30} & \frac{1}{31} & \frac{1}{32} & \frac{1}{33} & \frac{1}{34} \\ \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} & \frac{1}{26} & \frac{1}{27} & \frac{1}{28} & \frac{1}{29} & \frac{1}{30} & \frac{1}{31} & \frac{1}{32} & \frac{1}{33} & \frac{1}{34} & \frac{1}{35} \\ \frac{1}{17} & \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} & \frac{1}{26} & \frac{1}{27} & \frac{1}{28} & \frac{1}{29} & \frac{1}{30} & \frac{1}{31} & \frac{1}{32} & \frac{1}{33} & \frac{1}{34} & \frac{1}{35} & \frac{1}{36} \\ \frac{1}{18} & \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} & \frac{1}{26} & \frac{1}{27} & \frac{1}{28} & \frac{1}{29} & \frac{1}{30} & \frac{1}{31} & \frac{1}{32} & \frac{1}{33} & \frac{1}{34} & \frac{1}{35} & \frac{1}{36} & \frac{1}{37} \\ \frac{1}{19} & \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} & \frac{1}{26} & \frac{1}{27} & \frac{1}{28} & \frac{1}{29} & \frac{1}{30} & \frac{1}{31} & \frac{1}{32} & \frac{1}{33} & \frac{1}{34} & \frac{1}{35} & \frac{1}{36} & \frac{1}{37} & \frac{1}{38} \\ \frac{1}{20} & \frac{1}{21} & \frac{1}{22} & \frac{1}{23} & \frac{1}{24} & \frac{1}{25} & \frac{1}{26} & \frac{1}{27} & \frac{1}{28} & \frac{1}{29} & \frac{1}{30} & \frac{1}{31} & \frac{1}{32} & \frac{1}{33} & \frac{1}{34} & \frac{1}{35} & \frac{1}{36} & \frac{1}{37} & \frac{1}{38} & \frac{1}{39} \end{array})$

条件数の値を取得します。シンボリックメソッドから導出されたので数値誤差がありません。

vpa(cond(Hsym))

ans =  $24521565858153031724608315432.509$

条件数が大きくても厳密な逆行列を計算できます。

Hsym*inv(Hsym)

ans = 
 $(\begin{array}{c} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{array})$