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

Этот пример показов, как вычислить обратную матрицу Гильбертовой матрицы с помощью 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

Поддержка