Here is matrix inversion implemented in pure Haskell and a demonstration of its use with Complex Rational exact arithmetic:
A= [ 9 8*i 9*i 4 8*i 3 5*i 3*i ; 6*i 6 7*i 5 0 3*i 5*i 7 ; 4 4*i 4*i 8 1*i 3 1 9*i ; 4 6*i 7 5 4*i 7 4 4 ; 0 3 6 0 4 2*i 0 7 ; 7 2 0 4*i 3 3*i 1 4*i ; 4*i 3 8 1 2*i 3 1 3*i ; 3 3 4*i 3*i 6 3 7 1*i ]
Ainv= [ -651341+282669*i 70449+105643*i 470260+543736*i -90558-322394*i 940582+401976*i -1471210-32706*i 304731+202312*i 428573+497242*i ; 544367-80697*i -1154274+187794*i 1139302-87593*i -21921+1224087*i 1488234-1168263*i -1627965+1312344*i -1248400-535541*i 337176-362289*i ; -31970+991102*i 93690-139399*i 123282-341373*i -421343-757909*i 50198+495419*i 8593-654052*i -468577+527063*i 897227+455914*i ; 695254+159129*i -629809+1226054*i -2027820-1557744*i -161289+1062653*i -1299979-787107*i 915864+11456*i 1378121-78438*i 1043558-421873*i ; -925872-186803*i 1371020-146791*i -2291428-34310*i 2099076-57282*i -3405480+2662395*i 1616586-1020331*i 313380-265939*i -1386116-96576*i ; -1420678-1211999*i 1206717-844239*i 93772+1350941*i 1290846-32987*i -511+2133602*i 1126031+1379618*i -2659652-75501*i -2000675-207615*i ; 2039672+316590*i -813732+665441*i 417543-103784*i -2312041+106452*i 1746852-2305394*i -871774-685499*i 1641748-24893*i -261188+84637*i ; -23113-302279*i -610268-221894*i 1101428+322959*i -838351-211053*i -239412-540356*i 160747+259506*i 736020+689615*i -180808+391291*i ] / (-14799118+6333791*i)
Aesthetically striking is how simple input matrices become complicated when inverted: complexity seemingly appears out of nowhere. The complexity arises from a convolution of several effects: numerators and denominators intermingle as fractions get added; real and imaginary parts intermingle as complex numbers get multiplied; matrix elements intermingle in row operations.
(Perhaps this mixing could be put to cryptographic use someday? Row operations do vaguely look like Salsa20. AES has a field inversion to define its S box.)
The code is kind of a mess, not really trying to achieve performance: we compute the LU decomposition to get the determinant, which is used just to provide the scalar factor in the final displayed result of the inverse to avoid fractions inside the matrix. (It is neat that multiplying through by the determinant makes fractions disappear.) We do a separately implemented Gauss-Jordan elimination on the original matrix to compute the inverse, then multiply by the determinant computed via LU.
We made modifications to two library packages.
We forked the base module Data.Complex to a new module Data.ComplexNum to provide a Complex type that does not require RealFloat (so only requires Num). (Complex normally requires RealFloat because abs and signum of a Complex number requires sqrt.) Avoiding RealFloat allowed us to form Complex Rational, i.e., Complex (Ratio Integer), a powerful data type. We defined abs and signum as error.
Our starting point for matrix operations was Data.Matrix in the matrix package (by Daniel Díaz). The original luDecomp function calls abs to select the largest pivot during LU decomposition. The output type of abs is constrained in the Num class to be the same as the input type (abs :: a -> a), so Complex numbers weirdly provide a Complex output type for abs even through mathematically the absolute value of a complex number is real. Unfortunately Complex numbers do not provide Ord, so we cannot select the "largest" absolute value as the pivot. Therefore, we added a new entry point to LU decomposition, luDecompWithMag, to allow passing in a custom function to compute the magnitude. Since we are avoiding square root, we provided magnitude-squared (magnitude2) as the pivot-choosing function.
We fixed some (but not all) space leaks in LU decomposition so processing matrices of size 100x100 no longer requires gigabytes of memory.
We added to Data.Matrix two different implementations of matrix inversion both which use monadic imperative style for in-place modification of the augmented matrix with row operations: one (invertS) uses Control.Monad.Trans.State.Strict and the other (invert) uses MVector and Control.Monad.ST. The latter is a little bit faster. Both implementations required care to avoid space leaks, which were tracked down using ghc heap profiling.
The modifications to Data.Matrix can be tracked in this github repository.
Runtime (in seconds) and maxresident memory usage (k) of inversion of random Complex Rational matrices, where each entry in the original matrix is a random single-digit pure real or pure imaginary number (as in the example above with size 8), of sizes 10 through 320: [ 10 0.01 2888 ; 20 0.26 4280 ; 30 1.40 5324 ; 40 4.84 7132 ; 50 12.32 10204 ; 60 27.05 14292 ; 70 52.67 19416 ; 80 94.36 26588 ; 90 157.18 33752 ; 100 250.12 44000 ; 110 373.71 56524 ; 120 546.32 69600 ; 130 768.92 85988 ; 140 1056.05 103632 ; 150 1419.47 124900 ; 160 1877.97 148696 ; 170 2448.16 176364 ; 180 3130.31 208104 ; 190 3963.96 238824 ; 200 4958.93 277504 ; 210 6138.61 318704 ; 220 7543.16 364544 ; 230 9252.55 413700 ; 240 11115.32 467196 ; 250 13217.16 524300 ; 260 15690.80 594956 ; 270 18520.73 658448 ; 280 21761.36 734460 ; 290 25509.77 812288 ; 300 30048.34 1022028 ; 310 34641.25 1218636 ; 320 39985.63 1419340 ]
Runtime and memory for Hilbert matrices (Rational) of sizes 10 through 510: [ 10 0.00 2736 ; 20 0.03 3676 ; 30 0.14 4164 ; 40 0.51 4608 ; 50 1.45 5804 ; 60 3.65 8196 ; 70 7.86 9224 ; 80 15.59 11096 ; 90 27.85 13444 ; 100 45.03 17320 ; 110 67.97 21408 ; 120 98.61 26876 ; 130 137.02 27132 ; 140 187.48 29664 ; 150 248.99 37792 ; 160 325.17 38816 ; 170 420.00 43920 ; 180 529.02 55796 ; 190 659.03 56304 ; 200 811.39 65520 ; 210 980.33 70912 ; 220 1202.42 72608 ; 230 1426.34 82852 ; 240 1696.02 89996 ; 250 2014.39 97192 ; 260 2300.64 109496 ; 270 2693.48 118692 ; 280 3124.52 130980 ; 290 3589.06 138164 ; 300 4158.67 210900 ; 310 4711.61 223188 ; 320 5323.13 243692 ; 330 6047.59 262104 ; 340 6789.74 241640 ; 350 7621.69 269268 ; 360 8433.56 338924 ; 370 9447.28 360428 ; 380 10496.28 358356 ; 390 11662.98 422872 ; 400 12935.35 375772 ; 410 14203.69 477164 ; 420 15655.88 420840 ; 430 17143.12 547820 ; 440 19003.82 517104 ; 450 20592.28 550880 ; 460 22471.48 574428 ; 470 24484.14 604120 ; 480 26672.65 662492 ; 490 28958.52 766936 ; 500 31302.45 646132 ; 510 33932.76 731120 ]
Slope of runtime on a log-log plot is about 4.3.
According to Wikipedia, the worst case run time (bit complexity) of Gaussian elimination is exponential, but the Bareiss algorithm can guarantee O(n^5): Bareiss, Erwin H. (1968), "Sylvester's Identity and multistep integer-preserving Gaussian elimination", Mathematics of Computation 22 (102): 565–578, doi:10.2307/2004533. Future work: try the code presented here on the exponentially bad matrices described in the paper; implement Bareiss algorithm.
Also someday, it may be a fun exercise to create a data type encoding algebraic expressions which is an instance of Num, then use the polymorphism of the matrix code presented here to invert symbolic matrices.
No comments :
Post a Comment