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 = 2.3622e+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.0012 -0.0012 0.0118 -0.1624 -0.7867 0.9817 -0.5909 -1.1786 -0.1111 -3.7477 4.3067 2.3698 -0.2680 -2.0770 0.4077 2.9125 2.2648 0.1082
0.0871 1.0000 -0.0006 -0.0009 -0.0090 -0.5388 1.3683 -0.2444 -0.1039 -0.1231 -0.9150 -2.8117 6.5888 -0.1975 1.0052 -2.4624 -0.2692 5.2024 1.2478 0.0755
0.2054 -0.3020 0.9994 0.0024 0.0061 -0.4914 1.2560 -0.5424 0.0382 -1.1778 -1.3460 -1.3399 5.6394 -0.1046 0.9009 -2.6941 -0.3658 5.2330 1.3704 0.1685
0.3158 -0.8256 0.4401 1.0012 0.0370 -0.5177 1.6844 -0.0735 0.1373 -0.4658 -0.6256 -3.0880 7.6058 -2.1391 -0.0357 -3.9665 -0.8322 5.9076 0.9175 0.4910
0.4028 -1.3337 0.9924 -0.0390 1.0449 -0.4168 1.0098 -1.5264 1.1591 -1.7210 0.4572 -2.4602 7.8277 -0.3915 0.5983 -2.6859 -0.2981 5.7396 -1.8456 0.3786
0.4670 -1.7512 1.4662 -0.0101 0.0223 0.6950 1.5640 -1.0055 1.0727 -0.7692 -0.1271 -2.4243 7.0420 -1.7390 1.2946 -2.6481 -0.6277 4.8238 -0.6456 0.7012
0.5129 -2.0689 1.7984 0.1538 -0.1460 -0.3917 2.7314 -1.5910 0.8051 -0.8093 -0.2009 -3.9021 8.9283 -2.9270 1.7746 -3.9919 -0.6623 3.9135 -0.7386 0.7756
0.5446 -2.2991 1.9870 0.4589 -0.4710 -0.0614 0.9184 0.8377 -0.3169 0.0831 -1.5936 -0.2424 3.5135 0.2075 0.6568 -2.7485 -0.3657 2.5463 0.0261 0.0726
0.5656 -2.4583 2.0538 0.8936 -0.9286 0.2678 0.2762 -0.1833 0.9221 -0.8345 -0.9206 0.6412 3.5666 -0.5856 -0.1836 -1.2843 0.8606 1.5268 2.0609 -0.4121
0.5788 -2.5622 2.0255 1.4119 -1.3842 0.1638 1.2865 -0.7210 1.0717 0.3686 0.1248 -2.1556 5.7747 -1.0638 1.6798 -2.0859 -0.1640 0.9892 -0.4810 0.3520
0.5859 -2.6238 1.9256 2.0017 -1.9942 0.6104 -0.3480 -0.0803 0.1818 -1.1004 0.6539 -0.2920 1.3655 -0.4783 -0.5803 -0.7372 1.0466 -0.6521 2.4559 -0.5805
0.5887 -2.6534 1.7763 2.6275 -2.6306 0.6814 0.3989 -0.1594 -0.1037 -0.3355 -0.1718 0.7493 2.1269 0.2033 0.8929 -1.1647 0.3534 1.0040 0.2170 -0.0984
0.5882 -2.6590 1.5934 3.2715 -3.2617 0.9219 0.8125 -0.8750 0.8750 -0.5469 0.1875 -2.2500 5.5000 -2.5000 0.5000 -2.1250 -0.6562 1.5000 -0.5000 0.2500
0.5852 -2.6467 1.3896 3.9117 -3.9126 1.1515 1.0874 -1.4450 1.1922 -0.4310 0.0610 -3.5599 6.4263 -1.4031 1.6006 -3.2961 -0.8064 4.3220 -1.9198 0.7373
0.5804 -2.6214 1.1733 4.5546 -4.6713 1.7987 -0.3512 0.2821 -0.2734 -0.8129 0.3491 0.5070 0.6902 0.5044 0.0835 0.0187 0.1956 -0.2224 1.9464 -0.1961
⋮
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
-0.0000 0.0000 -0.0000 0.0000 -0.0000 -0.0001 0.0004 0.0050 -0.0400 0.0368 -0.4875 1.2090 -1.1395 -0.8712 -0.3376 3.4449 -1.6401 0.3319 0.0022 0.0017
0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0001 0.0009 -0.0015 -0.0172 0.0267 -0.1909 0.1515 -0.2069 -2.4186 -0.1889 3.6202 -0.9900 0.3459 0.0372 -0.0030
-0.0000 0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0007 0.0097 -0.0384 0.0201 -0.3053 1.3925 -2.2012 -1.1811 -0.7248 3.8588 -1.6106 0.5679 -0.0543 0.0025
0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0001 0.0007 0.0020 -0.0014 -0.0255 -0.1178 0.7277 -0.7834 -2.9257 1.2532 2.2557 -0.7382 0.3044 -0.0141 -0.0050
-0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0001 0.0000 0.0043 -0.0304 0.0353 -0.3375 0.7563 -1.3325 -2.0856 -1.4510 3.1561 -1.1128 0.4057 -0.0310 0.0002
⋮
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 =
