Main Content

Hilbert Matrices and Their Inverses

This example shows how to compute the inverse of a Hilbert matrix using Symbolic Math Toolbox™.

Definition : A Hilbert matrix is a square matrix with entries being the unit fraction. Hij=1i+j-1. For example, the 3x3 Hilbert matrix is H=[11213121314131415]

Symbolic computations give accurate results for these ill-conditioned matrices, while purely numerical methods fail.

Create a 20-by-20 numeric Hilbert matrix.

H = hilb(20);

Find the condition number of this matrix. Hilbert matrices are ill-conditioned, meaning that they have large condition numbers indicating that such matrices are nearly singular. Note that computing condition numbers is also prone to numeric errors.

cond(H)
ans = 3.4374e+18

Therefore, inverting Hilbert matrices is numerically unstable. When you compute a matrix inverse, H*inv(H) must return an identity matrix or a matrix close to the identity matrix within some error margin.

First, compute the inverse of H by using the inv function. A warning is thrown due to the numerical instability.

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.0001   -0.0022    0.0587    0.8314   -2.3117   -0.8683    1.2591   -1.6349   -0.9861    1.4523    4.3067    0.3698    1.7320   -0.0770    0.4077    2.9125    2.2648    0.1082
    0.1802    1.0000    0.0002   -0.0010    0.0290    0.2469   -0.4234   -3.4230    1.4735   -0.3983   -0.0370   -2.0498    6.5888   -0.1975    1.0052   -2.4624   -0.7692    3.2024   -0.7522    0.0755
    0.3735   -0.4735    1.0004   -0.0025    0.0196    0.1023   -1.0735   -2.2242    1.1291    0.0097    0.0518   -0.7944    5.6394   -0.1046    1.9009   -0.6941   -0.8658    3.2330   -0.6296    0.6685
    0.4855   -0.7484   -0.0892    0.9990    0.0288    0.0054   -0.5003   -2.5844    0.0721    0.3358   -0.1093   -0.5662    7.6058   -2.1391    0.9643   -1.9665   -1.3322    3.9076   -1.0825    0.4910
    0.5432   -0.8111   -0.4820    0.2616    1.0136   -0.0105   -0.1986   -2.2764    1.3674   -1.0126    1.1863   -2.7936    7.8277   -0.3915    1.5983   -0.6859   -0.7981    3.7396   -1.8456    0.3786
    0.5691   -0.7468   -1.0618    0.7180    0.0106    0.8025    0.2065   -2.4355    1.7777    0.1033    0.1979   -2.4643    7.0420   -1.7390    2.2946   -1.6481   -1.1277    2.8238   -2.6456    0.2012
    0.5765   -0.6188   -1.7164    1.3002   -0.0543    0.0001    0.8564   -2.0045    1.6032   -0.3129    0.6886   -5.5175    8.9283   -2.9270    1.7746   -1.9919   -1.1623    3.9135   -2.7386    0.7756
    0.5732   -0.4647   -2.3749    1.9346   -0.2175    0.1400   -0.8640    0.3099   -0.7429    0.7839   -1.0358   -0.3165    3.5135    0.2075    0.6568   -1.7485   -0.8657    2.5463   -0.9739    0.3226
    0.5638   -0.3050   -2.9997    2.5794   -0.4041    0.2633   -0.5988   -0.5048    1.1364   -0.4863    0.5347   -0.7874    3.5666    0.4144   -0.1836   -0.2843    0.3606   -0.4732    1.0609   -0.1621
    0.5510   -0.1504   -3.5719    3.1963   -0.5577    0.0054    0.3813   -1.0399    1.3217    1.0949    0.2605   -0.8452    5.7747   -1.0638    1.6798   -1.0859   -0.6640    0.9892   -0.4810    0.3520
      ⋮

Now, use the MATLAB® invhilb function that offers better accuracy for Hilbert matrices. This function finds exact inverses of Hilbert matrices up to 15-by-15. For a 20-by-20 Hilbert matrix, invhilb finds the approximation of the matrix inverse.

H*invhilb(20)
ans = 20×20
1010 ×

    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
      ⋮

To avoid round-off errors, use exact symbolic computations. For this, create the symbolic Hilbert matrix.

Hsym = sym(H)
Hsym = 

(11213141516171819110111112113114115116117118119120121314151617181911011111211311411511611711811912012113141516171819110111112113114115116117118119120121122141516171819110111112113114115116117118119120121122123151617181911011111211311411511611711811912012112212312416171819110111112113114115116117118119120121122123124125171819110111112113114115116117118119120121122123124125126181911011111211311411511611711811912012112212312412512612719110111112113114115116117118119120121122123124125126127128110111112113114115116117118119120121122123124125126127128129111112113114115116117118119120121122123124125126127128129130112113114115116117118119120121122123124125126127128129130131113114115116117118119120121122123124125126127128129130131132114115116117118119120121122123124125126127128129130131132133115116117118119120121122123124125126127128129130131132133134116117118119120121122123124125126127128129130131132133134135117118119120121122123124125126127128129130131132133134135136118119120121122123124125126127128129130131132133134135136137119120121122123124125126127128129130131132133134135136137138120121122123124125126127128129130131132133134135136137138139)

Get the value of the condition number. It has been derived by symbolic methods and is free of numerical errors.

vpa(cond(Hsym))
ans = 24521565858153031724608315432.509

Although its condition number is large, you can compute the exact inverse of the matrix.

Hsym*inv(Hsym)
ans = 

(1000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001000000000000000000001)