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. . For example, the 3x3 Hilbert matrix is
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
-0.6412 2.3515 0.0572 -2.5231 1.4707 -0.2228 -0.0079 0.4551 -1.0618 1.0000 0.8374 -0.2768 -0.4907 0.0096 0.6562 -3.7467 0.2445 1.3445 -0.0822 -0.1250
-0.6391 2.3329 0.2879 -3.0792 1.9023 -0.4960 0.0371 -0.0346 -0.2657 0.1250 -0.3391 1.1698 -0.0063 -1.4632 -0.2188 -4.1549 0.0700 1.6185 -0.9817 0.0938
-0.6340 2.2998 0.5206 -3.6240 2.3359 -0.5938 0 0.3750 -0.9863 -0.2500 -0.3125 1.5000 -1.0000 0.2500 0.7500 -4.6250 1.5000 0 1.1250 -0.2188
-0.6269 2.2568 0.7498 -4.1541 2.8242 -0.8203 -0.0625 -0.1250 0.1895 -1.1250 0.3359 -1.1250 0.5000 -0.7344 -0.4688 0.6250 0.2500 -2.0000 1.3750 -0.4062
-0.6184 2.2073 0.9709 -4.6567 3.2715 -1.0000 -0.1875 0.5312 -0.3340 0.3438 -1.1055 -0.1250 -1.1250 1.3906 0.4688 -5.1875 0 -0.5000 -0.2500 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
-0.0000 0.0000 -0.0000 0.0000 -0.0000 -0.0001 0.0004 0.0050 -0.0400 -0.0382 -0.4875 1.2090 -1.1395 -0.8712 0.2819 3.4449 -1.6401 0.3319 0.0022 0.0055
0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0001 0.0009 -0.0015 -0.0172 -0.0155 -0.1909 0.1515 -0.2069 -2.4186 1.2751 3.6202 -0.9900 0.3459 0.0372 -0.0033
0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0008 -0.0053 -0.0272 -0.0822 -0.4647 0.8892 0.1879 -1.2012 -0.1879 4.7513 -1.4663 0.2265 -0.0168 0.0013
0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0002 0.0008 0.0001 -0.0568 -0.0960 -0.3590 0.5939 -1.6408 -1.7985 0.5570 2.4159 -0.5771 0.4521 0.0346 -0.0028
0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0000 0.0004 0.0010 -0.0072 0.0075 -0.0210 0.9664 -3.4762 -1.4294 -0.1812 2.0569 -1.8690 0.3523 0.0289 0.0052
⋮
To avoid round-off errors, use exact symbolic computations. For this, create the symbolic Hilbert matrix.
Hsym = sym(H)
Hsym =
Get the value of the condition number. It has been derived by symbolic methods and is free of numerical errors.
vpa(cond(Hsym))
ans =
Although its condition number is large, you can compute the exact inverse of the matrix.
Hsym*inv(Hsym)
ans =
