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 = 
5.1944e+19

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 =  1.309732e-19.
ans = 20×20

    1.0000   -0.0000   -0.0004   -0.0133   -0.0078    0.1494    0.4224   -0.1053   -1.5628    1.5000    0.4620    0.9110    2.8195   -0.3235   -0.6250   -5.6448   -0.9671   -4.5672    3.7431   -0.3125
   -0.0558    1.0000    0.0000   -0.0050    0.0078    0.1464    0.2495    0.3592   -1.6037   -0.1875   -0.0668    0.9002    2.9993   -1.2306   -1.6875   -6.1419   -1.9599   -2.3883    0.7089    0.3750
   -0.2370    0.5524    1.0000   -0.0034    0.0195    0.2554    0.0272    0.5033   -1.1266    0.1875   -0.3678    0.9255    0.9066    0.7549   -0.8125   -6.0203   -1.3259    0.5149   -0.5031    0.1875
   -0.3796    1.1272   -0.2467    0.9963    0.0273   -0.0371    0.4644   -0.5493   -0.1693   -0.9375   -0.2952    0.1800    2.9372   -1.2208   -2.0000   -1.5512   -1.3187    0.1597    0.2638         0
   -0.4789    1.5782   -0.4757   -0.0590    0.9922    0.1828   -0.0083    0.3378   -0.1781   -0.5000   -0.1041    0.6158    2.4213   -1.8971   -1.2500   -1.5460   -0.2774   -2.0273    1.0890   -0.1875
   -0.5455    1.8979   -0.5939   -0.2425    0.0664    1.0698    0.2080    0.4684   -0.9442   -0.9375   -0.5483    0.7483    1.0083    0.8748    1.0625   -4.4184   -0.7901    1.6298   -0.7021    0.1250
   -0.5889    2.1108   -0.5976   -0.5488    0.2051    0.1245    1.0389   -0.3502   -0.7717   -0.3750   -0.2460    0.5758    1.9188   -0.3853   -0.3125   -3.8710   -1.0324   -0.5902    0.1486    0.0625
   -0.6160    2.2432   -0.5093   -0.9648    0.4648   -0.1813    0.7114   -0.1442    0.2887   -1.0000   -0.2212   -1.5915    1.7356   -1.7679   -1.4375   -1.7909   -4.9411    5.2257   -2.2542    0.6250
   -0.6317    2.3172   -0.3569   -1.4416    0.6758    0.0942   -0.2973    0.4413   -0.9635    1.0000   -0.9617    1.9564   -3.2818    2.5536    1.9062   -7.7002    2.1550   -1.9163   -0.7808    0.2500
   -0.6393    2.3492   -0.1620   -1.9724    1.0527   -0.1114    0.2407    0.0095   -0.2193    1.5000   -0.7571    1.0671   -0.7490   -0.6288    0.1250   -4.0656    0.5029   -0.3419    0.5529    0.0938
      ⋮

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.1057    0.4019    0.2620   -0.4443   -7.5099    1.5435    3.2884   -1.1618    0.2301    0.0437    0.0029
   -0.0000    0.0000   -0.0000    0.0000   -0.0000   -0.0001   -0.0009    0.0172   -0.0628    0.2441   -0.8975    2.7711   -2.8924   -1.4888   -4.3218    6.4897   -3.0581    0.7925   -0.0037    0.0029
    0.0000    0.0000    0.0000   -0.0000   -0.0000   -0.0001    0.0009    0.0042   -0.0303   -0.1460   -0.1028    0.2862   -0.9267   -4.0597   -4.8050    4.7727   -1.3255    0.4667   -0.0211   -0.0033
    0.0000    0.0000   -0.0000   -0.0000   -0.0000   -0.0001    0.0004    0.0002   -0.0056   -0.1653   -0.1136    0.3616   -1.6920   -3.7214    0.9328    4.4169   -1.1602    0.5541   -0.0024   -0.0050
   -0.0000    0.0000   -0.0000    0.0000   -0.0000   -0.0000   -0.0006    0.0068   -0.0394    0.1946   -0.7818    1.9314   -2.1273   -1.0745    1.2885    4.6349   -1.8923    0.3793   -0.0435    0.0045
    0.0000   -0.0000    0.0000   -0.0000    0.0000   -0.0002    0.0014   -0.0028    0.0304    0.0562   -0.0724   -0.2073   -0.3769   -5.0729    4.9325    2.3488   -0.5046    0.2240    0.0675   -0.0138
   -0.0000    0.0000   -0.0000    0.0000   -0.0000    0.0000   -0.0013    0.0101   -0.0520    0.0721   -0.7715    2.9297   -3.3374    1.4084    0.5570    7.1053   -2.7578    0.8815   -0.1648   -0.0038
    0.0000   -0.0000    0.0000   -0.0000    0.0000   -0.0002    0.0018   -0.0040    0.0231   -0.2278    0.2000    0.2766   -0.8262   -4.7229    3.8319    2.1857   -0.4749   -0.1747    0.1257   -0.0141
   -0.0000    0.0000   -0.0000    0.0000   -0.0000    0.0000    0.0004    0.0144   -0.0513    0.1850   -0.7063    1.9776   -2.0602   -1.0763   -0.0201    3.8539   -2.4580    0.5675   -0.0189    0.0036
    0.0000    0.0000    0.0000   -0.0000   -0.0000   -0.0002    0.0008   -0.0019    0.0027   -0.0679   -0.0497    0.8865   -0.6720   -1.2729    0.8925    2.3625   -1.2616    0.2443    0.0184   -0.0091
      ⋮

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)