Гильбертовы матрицы и их инверсии

В этом примере показано, как вычислить инверсию Гильбертовой матрицы с помощью 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 = 2.1065e+18

Поэтому инвертирование Гильбертовых матриц численно неустойчиво. Когда вы вычисляете обратную матрицу, H*inv(H) должен возвратить единичную матрицу или матрицу близко к единичной матрице в некотором допуске на погрешность.

Во-первых, вычислите инверсию H при помощи inv функция. Предупреждение выдано из-за числовой нестабильности.

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.0003    0.0044    0.0212   -0.2436    0.4383   -0.6683   -0.4409    0.6464   -0.9986   -0.4477    8.3067   -0.2302   -0.4680   -0.4770    1.4327   -1.4875    8.0000    2.0000
    0.4314    1.0000   -0.0002    0.0003    0.0528   -0.0284    0.9397    0.6603   -0.1158    1.3085   -0.5132    0.2835    0.5888   -1.7213    0.7195   -1.7957    0.6355    2.8214   -8.0000    2.0000
    0.7987   -0.8048    0.9997   -0.0041    0.1027    0.1222    0.5629   -0.1333    0.5837    1.0126   -0.9596    0.1601   -0.3606   -2.5592    0.5373   -2.9668   -0.2067    3.8693   -8.0000         0
    1.0486   -1.6648    0.4041    0.9885    0.0958   -0.1427    1.3149    0.2797   -0.1724    0.8018   -0.1066    0.3468    5.6058   -2.5304    2.5296   -3.6839   -0.4191    3.5598  -16.0000    1.0000
    1.2164   -2.4406    1.0946   -0.1853    1.1204   -0.1147    1.2598   -0.9847    0.2424   -0.0230    0.5197    0.8731    3.8277   -4.7249    0.5983   -2.3525   -0.3189    3.4062   -4.0000    3.0000
    1.3286   -3.1054    1.9447   -0.5911    0.1962    0.8550    1.2190   -1.0755    0.2927    1.1158    0.2254    0.9657    3.0420   -2.0190    1.2346   -1.7681   -0.4577    3.5038         0         0
    1.4027   -3.6617    2.8529   -1.1814    0.3594   -0.2091    2.4910   -0.4949    0.6416    0.4732    0.2943   -0.5944    4.9283   -4.1578    1.6593   -3.0688   -0.1768    4.6058         0    2.0000
    1.4501   -4.1206    3.7522   -1.8860    0.5367    0.0995    0.2101    1.4302   -1.0067    1.5837   -1.5936    3.3039   -0.4865    0.0223    0.9902   -1.7856   -0.3147    2.2500         0         0
    1.4784   -4.4954    4.6045   -2.6478    0.8130    0.1115    0.0798    0.1202    1.1006    0.0361   -0.2599    3.4626   -0.4334   -0.7285   -0.3979    0.2157    0.6016    0.2411         0   -1.0000
    1.4930   -4.7986    5.3893   -3.4322    1.1685   -0.0259    0.4050   -0.0399    1.2269    0.7081    0.0601    0.9479    3.7747   -2.1672    1.9212   -1.5514    0.2045    2.7134   -4.0000         0
      ⋮

Теперь используйте 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.4295    0.0537
   -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   -0.0268
    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.2147   -0.0268
    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   -0.0268
   -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.2147    0.0537
    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         0
   -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.2147         0
    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         0
   -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    0.0268
    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.2147         0
      ⋮

Чтобы избежать ошибок округления, используйте точные символьные расчеты. Для этого создайте символьную Гильбертову матрицу.

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})$

Документация

Гильбертовы матрицы и их инверсии

Документация Symbolic Math Toolbox

Поддержка