Rounding can make a numerically singular matrix regular so that is probably what we are seeing here. You initial covariance matrix is actually singular, but when copied from the mail or read from a file it becomes non-singular. You try the following
julia> A = randn(15, 14);B = A*A'; # B is singular except for floating point errors julia> norm(inv(B)*B - I) 219.10806398211932 julia> C = round(B, 9); julia> norm(inv(C)*C - I) 3.5885069748862715e-6 2015-01-15 21:54 GMT-05:00 ShaoWei Teo <[email protected]>: > Thanks for helping with the troubleshoot. > > Oh my, this is really weird! I simplified the problem and post it here, > which I realised could not reproduce the result. More on this later. > > The original problem is to calculate the inverse of a covariance matrix. > Let's say my input is:- > h = > 0.001508 0.013789 -0.00104 0.005551 0.004493 -0.00884 0.013144 0.017682 > -0.00487 0.000629 -0.05168 -0.00023 0.002338 0.007286 0.007799 0.000526 > -0.00334 -0.00314 -0.00392 -0.00752 -0.00797 -0.00741 -0.00271 -0.00521 > -0.02595 -0.00841 -0.00058 0.007486 -0.01012 -0.0026 0.002072 0.003723 > 0.002328 0.007764 0.007441 -0.00261 0.018319 0.023409 0.000208 0.016593 > -0.02026 0.000187 0.004628 0.004931 0.015028 -0.00523 0.002056 -0.00962 > 0.031875 -0.0027 -0.00029 0.005102 -0.00589 -0.00548 0.004088 -0.01538 > -0.00087 -0.02139 -0.01081 0.006092 -0.00728 0.019137 -0.00297 0.001239 > 0.010175 -0.00556 0.019964 0.014191 -0.00666 0.008469 -0.05717 -0.00424 > -0.03361 -0.00972 0.011719 0.002002 0.006843 -0.00058 0.023126 0.014145 > 0.008645 0.008639 0.013656 0.009018 0.010999 0.058236 -0.00115 -0.01331 > 0.009285 0.014607 0.007854 -0.02343 0.011951 -0.01901 -0.022 0.004033 > -0.03805 -0.03497 0.010243 -0.02748 -0.05395 0.009717 -0.01425 0.009262 > -0.0257 0.003396 0.001709 -0.00489 0.024266 0.018609 -0.01335 0.013088 > 0.012458 -0.00599 0.014979 0.049237 -0.00056 0.01015 -0.00691 0.019328 > 0.002141 -0.01978 0.009303 -0.0404 -0.02996 0.005715 -0.02302 -0.02434 > 0.008846 -0.02494 -0.00902 0.011864 -0.00023 0.014309 -0.027 0.010168 > -0.03179 0.006432 0.010828 -0.05265 0.006049 -0.05946 -0.04611 0.004026 > -0.02481 -0.01016 0.010103 -0.00369 0.007764 -0.01631 -0.00382 0.031399 > -0.00182 0.016132 0.039641 -0.00495 0.040496 0.046879 -0.00186 0.032187 > 0.012492 -0.0045 0.033889 -0.00991 0.027912 -0.00523 0.028394 -0.00907 > 0.004339 0.006011 -0.00727 0.01518 0.004874 -0.00712 0.005307 0.053483 > -0.00014 0.036373 -0.00403 0.012531 > My cond(cov(h)) is about 6e17 while norm(cov(h) * inv(cov(h)) - > eye(cov(h))) is 1408.688. > > I calculated inv(cov(h)) in Float64 to be:- > 4.75E+19 3.89E+18 -1.11E+20 -1.89E+19 -1.35E+19 -3.23E+19 -2.19E+19 > 1.48E+19 8.26E+19 3.33E+19 -3E+18 -5.33E+19 -6.8E+17 -1.21E+19 -2.26E+19 > 1.14E+19 3.97E+17 -2.55E+19 -4.83E+18 -1.5E+18 -3.85E+18 -5.46E+18 3.2E+18 > 1.42E+19 6.92E+18 -7.8E+17 -8.86E+18 -1.5E+17 -2.23E+18 -3.66E+18 > -1.09E+20 -2.19E+19 7.71E+19 1.32E+19 1.69E+19 -5.48E+18 -7.86E+18 > -9.20E+18 -4.46E+19 -7.46E+18 1.62E+18 -5.59E+19 8.33E+18 3.78E+19 > 2.44E+18 -1.40E+19 -3.22E+18 1.13E+19 1.56E+18 3.46E+18 2.36E+18 -1E+18 > -1.5E+18 -1.04E+19 -1.7E+18 2.02E+17 -4.88E+18 1.1E+18 5.44E+18 1.32E+18 > -2.70E+19 -3.3E+18 2.21E+19 5.03E+18 4.52E+17 -8.11E+18 8.83E+17 -2.2E+18 > -9.9E+17 -2.9E+18 5.99E+17 -1.08E+19 1.52E+18 6.55E+18 4.46E+17 -5.92E+19 > -1.05E+19 3.87E+18 3.32E+18 -2.1E+18 -3.68E+19 -1.27E+19 1.24E+18 3.30E+19 > 1.06E+19 -4.3E+16 -6.71E+19 5.03E+18 1.70E+19 -1.08E+19 -1.41E+19 > -4.91E+18 -1.17E+19 -2.41E+18 1.96E+18 -7.24E+18 -8.54E+18 1.64E+18 > 7.68E+18 6.67E+18 -5.3E+17 -2.55E+19 1.92E+18 6.66E+18 -4.06E+18 1.61E+19 > 3.10E+18 -9.71E+18 -1.8E+18 -1.9E+18 2.49E+18 1.46E+18 1.07E+18 3.73E+18 > 4.8E+17 -2E+17 9.73E+18 -1.2E+18 -5.32E+18 1.63E+17 1.26E+20 2.18E+19 > -6.20E+19 -1.43E+19 -5.43E+18 4.30E+19 1.29E+19 5.47E+18 -7.31E+18 > -3.5E+18 -1.5E+18 9.79E+19 -9.62E+18 -3.80E+19 9.95E+18 3.03E+19 6.18E+18 > -4.4E+18 -1.3E+18 -2.28E+18 1.10E+19 6.51E+18 9.52E+16 -6.15E+18 -2.46E+18 > 1.8E+17 2.28E+19 -2.45E+18 -9.54E+18 7.3E+17 -2.79E+18 -7.1E+17 1.68E+18 > 2.15E+17 5.85E+17 7.67E+16 -4.3E+17 -2.1E+17 -1.56E+18 6.91E+16 3.89E+16 > -2.42E+18 2.7E+17 1.18E+18 -2.3E+17 -3.91E+19 -1.00E+19 -7.10E+19 > -9.82E+18 -6.21E+18 -5.39E+19 -2.96E+19 1.08E+19 8.47E+19 2.61E+19 > -2.94E+18 -7.40E+19 4.43E+18 1.06E+19 -1.11E+19 -2.7E+18 1.79E+16 > 9.56E+18 1.67E+18 7.85E+17 2.75E+18 2.29E+18 -1.3E+18 -6.85E+18 -2.8E+18 > 3.03E+17 4.47E+18 -3.7E+16 5.68E+17 1.46E+18 -2.27E+19 -1.9E+18 4.34E+19 > 7.95E+18 3.61E+18 6.70E+18 7.85E+18 -5.46E+18 -2.58E+19 -1.09E+19 1.30E+18 > 1.02E+19 7E+17 5.64E+18 5.54E+18 -2.08E+19 -3.04E+18 -1.5E+18 9.1E+17 > -6.6E+17 -1.31E+19 -4.3E+18 7.45E+17 1.52E+19 1.35E+18 -3.6E+17 -1.03E+19 > 1.26E+18 4.36E+18 1.42E+18 > > In BigFloat, I get:- > 4.77E+80 1.49E+80 -4.60E+80 -5.23E+79 -1.22E+80 -8.97E+79 6.50E+79 > 5.96E+79 5.12E+80 -6.93E+79 -2.07E+79 8.19E+80 -6.40E+79 -2.69E+80 > 2.03E+80 1.72E+80 4.27E+79 -6.15E+79 -5.82E+78 -4.18E+79 -8.44E+78 > 3.44E+79 8.72E+78 7.25E+79 8.13E+78 1.17E+78 5.89E+79 -1.39E+79 -6.50E+79 > -2.47E+79 -4.73E+80 -2.15E+79 4.60E+80 1.27E+80 -8.47E+79 -3.20E+80 > 6.46E+79 -3.72E+79 1.97E+80 -1.88E+79 1.97E+79 -3.87E+80 2.45E+79 7.90E+79 > -9.53E+79 -2.44E+79 3.39E+78 9.66E+79 2.21E+79 -1.87E+79 -2.99E+79 > 2.38E+79 -9.25E+78 -8.70E+78 7.14E+78 6.36E+78 -1.14E+80 3.13E+78 7.19E+78 > -5.09E+79 -2.58E+80 -4.53E+79 1.28E+79 8.93E+78 2.08E+79 -9.30E+79 > -5.22E+79 1.75E+78 8.52E+79 -5.76E+78 -5.00E+78 -5.10E+79 1.53E+79 > 6.41E+79 4.95E+79 -3.92E+80 -1.78E+79 -1.10E+80 2.26E+79 -4.94E+79 > -3.33E+80 -5.71E+79 2.76E+79 5.16E+80 -4.69E+79 -1.61E+79 2.10E+80 > 3.82E+78 6.82E+78 1.76E+80 7.27E+79 4.29E+79 2.80E+79 2.27E+79 -6.88E+79 > -9.59E+79 5.09E+79 2.83E+78 1.49E+80 4.00E+78 5.11E+78 -1.48E+79 -9.98E+78 > -5.57E+79 -4.36E+79 7.90E+79 4.57E+78 -4.90E+79 -1.60E+79 1.16E+79 > 5.45E+79 -3.50E+78 3.05E+78 -4.78E+79 3.88E+78 -1.49E+78 3.95E+79 > -3.50E+78 -1.14E+79 2.98E+78 8.58E+80 6.84E+79 -1.08E+80 -8.58E+79 > 4.58E+79 5.18E+80 7.28E+79 -9.69E+78 -6.53E+80 1.09E+80 1.63E+79 -2.08E+80 > -2.31E+79 -9.20E+79 -2.69E+80 8.26E+79 2.92E+78 -3.90E+79 -2.15E+79 > 4.92E+79 1.43E+80 1.35E+78 -3.35E+78 -1.42E+80 -5.17E+79 -5.84E+78 > 2.86E+80 -9.75E+78 -1.74E+79 1.01E+80 -6.88E+78 2.87E+78 1.59E+79 > 4.86E+78 -4.32E+78 -8.06E+78 6.54E+78 -1.42E+78 8.87E+78 -6.05E+78 > 4.99E+77 1.42E+79 -8.27E+77 -3.20E+78 4.14E+78 1.40E+80 4.62E+79 > -3.01E+80 -1.46E+79 -1.94E+80 -4.24E+80 -3.52E+79 5.99E+79 6.23E+80 > 2.38E+80 6.86E+78 -8.09E+80 7.95E+78 -5.71E+79 -3.47E+80 -5.08E+79 > -1.30E+79 3.07E+79 3.07E+78 1.58E+79 1.31E+79 -7.27E+78 -4.60E+78 > -3.95E+79 -7.16E+78 -2.89E+76 -5.69E+78 3.92E+78 2.03E+79 1.11E+79 > -2.76E+80 -6.14E+79 1.31E+80 1.85E+79 5.71E+79 1.28E+78 -4.03E+79 > -1.69E+79 -1.05E+80 -2.15E+79 -3.40E+77 -7.95E+79 2.04E+79 9.64E+79 > 3.95E+79 -7.91E+79 -2.91E+79 -4.31E+79 -4.12E+78 -2.65E+79 -1.21E+80 > -4.92E+79 1.16E+79 1.00E+80 1.06E+80 4.31E+78 -4.73E+80 2.00E+79 5.34E+79 > -1.56E+80 > When I write the output to a .csv and read it in again. The trailing > numbers are removed, and then the calculation of matrix inverse may become > correct for Float64 and/or BigFloat... This is very weird. > > On Friday, January 16, 2015 at 3:40:23 AM UTC+8, Steven G. Johnson wrote: >> >> I can't reproduce your problem. With >> >> m = [2.69E-05 -2.25E-05 1.25E-05 -2.30E-05 >> -1.09E-05 9.02E-06 -1.42E-05 -4.29E-06 >> 1.47E-05 2.11E-06 -2.27E-06 1.61E-06 2.15E-05 >> 1.42E-05 -4.05E-05 >> -2.25E-05 0.000247072 -7.98E-06 2.24E-05 >> 6.45E-05 5.80E-06 8.79E-05 9.02E-05 -3.00E-06 >> 0.000103897 7.80E-05 1.98E-06 0.000156923 >> 8.07E-06 4.72E-05 >> 1.25E-05 -7.98E-06 4.21E-05 -3.97E-05 >> -4.36E-05 3.42E-05 -8.94E-05 -5.93E-05 >> 3.43E-05 -3.79E-05 6.44E-06 6.72E-06 2.05E-05 >> 4.11E-05 -7.26E-05 >> -2.30E-05 2.24E-05 -3.97E-05 0.000106283 >> 5.56E-05 -2.98E-05 9.56E-05 6.94E-05 -3.99E-05 >> 5.84E-05 -7.42E-05 2.62E-06 -7.48E-05 >> -4.52E-05 0.000142413 >> -1.09E-05 6.45E-05 -4.36E-05 5.56E-05 >> 0.000225041 -2.71E-05 0.000313224 0.000229185 >> -3.48E-05 0.000149339 -1.68E-05 -1.15E-05 >> 0.00013528 -4.35E-05 0.000184476 >> 9.02E-06 5.80E-06 3.42E-05 -2.98E-05 >> -2.71E-05 3.13E-05 -6.18E-05 -3.87E-05 >> 2.98E-05 -2.03E-05 7.29E-06 7.13E-06 2.97E-05 >> 3.44E-05 -5.49E-05 >> -1.42E-05 8.79E-05 -8.94E-05 9.56E-05 >> 0.000313224 -6.18E-05 0.000598084 0.000395901 >> -7.03E-05 0.000264016 -9.28E-05 -1.04E-05 >> 0.000167989 -7.99E-05 0.000271952 >> -4.29E-06 9.02E-05 -5.93E-05 6.94E-05 >> 0.000229185 -3.87E-05 0.000395901 0.000279212 >> -4.60E-05 0.000190289 -3.90E-05 -8.30E-06 >> 0.000143872 -5.35E-05 0.000199811 >> 1.47E-05 -3.00E-06 3.43E-05 -3.99E-05 >> -3.48E-05 2.98E-05 -7.03E-05 -4.60E-05 >> 3.38E-05 -1.98E-05 1.32E-05 5.89E-06 3.18E-05 >> 3.79E-05 -7.77E-05 >> 2.11E-06 0.000103897 -3.79E-05 5.84E-05 >> 0.000149339 -2.03E-05 0.000264016 0.000190289 >> -1.98E-05 0.000211512 -5.33E-05 -4.02E-06 >> 8.48E-05 -1.75E-05 0.000103224 >> -2.27E-06 7.80E-05 6.44E-06 -7.42E-05 >> -1.68E-05 7.29E-06 -9.28E-05 -3.90E-05 >> 1.32E-05 -5.33E-05 0.000201486 -1.31E-05 >> 0.000143599 1.20E-05 -6.10E-05 >> 1.61E-06 1.98E-06 6.72E-06 2.62E-06 -1.15E-05 >> 7.13E-06 -1.04E-05 -8.30E-06 5.89E-06 >> -4.02E-06 -1.31E-05 4.74E-06 -4.32E-06 >> 6.33E-06 -8.79E-06 >> 2.15E-05 0.000156923 2.05E-05 -7.48E-05 >> 0.00013528 2.97E-05 0.000167989 0.000143872 >> 3.18E-05 8.48E-05 0.000143599 -4.32E-06 >> 0.000328454 3.34E-05 -1.27E-05 >> 1.42E-05 8.07E-06 4.11E-05 -4.52E-05 >> -4.35E-05 3.44E-05 -7.99E-05 -5.35E-05 >> 3.79E-05 -1.75E-05 1.20E-05 6.33E-06 3.34E-05 >> 4.63E-05 -8.88E-05 >> -4.05E-05 4.72E-05 -7.26E-05 0.000142413 >> 0.000184476 -5.49E-05 0.000271952 0.000199811 >> -7.77E-05 0.000103224 -6.10E-05 -8.79E-06 >> -1.27E-05 -8.88E-05 0.000293684] >> >> I get cond(m) about 1e5, and norm(m * inv(m) - eye(m)) is 2.7e-11 (i.e. m >> * inv(m) is the identity to about 11 digits) >> >
