[math-fun] Floating point matrix product anomalies
Spin a 14gon RotationMatrix[\[Pi]/7.] // InputForm {{0.9009688679024191, -0.4338837391175581}, {0.4338837391175581, 0.9009688679024191}} one million clicks, less one: Nest[%.# &, IdentityMatrix@2, 999999] // InputForm // tim 0.218347,1 {{-1.000000000000001, -5.551115123125783*^-16}, {5.551115123125783*^-16, -1.000000000000001}} An error this small after a million clicks means the error can not grow, even with a trillion clicks. Couldn't we do even better with just twenty-odd multiplies in a repeated squaring (addition chain) process? MatrixPower[%%, 999999] // InputForm // tim 0.000035,1 {{-1.000000000055938, -5.657918578094723*^-12}, {5.657918578094723*^-12, -1.000000000055938}} Horrific! And 35 musec wasn't enough time for Mathematica to go off and outsmart itself with eigenvalues, e.g. But just to make sure, I rolled my own repeated squarer: mpowr[%%%, 999999] // InputForm // tim 0.000176,1 {{-1.0000000000272093, -3.690142635903726*^-11}, {3.690142635903726*^-11, -1.0000000000272093}} Ten times worse! Here's mpowr: mpowr[M_, 0] := IdentityMatrix[Length[M]] mpowr[M_, n_?EvenQ] := mpowr[M.M, n/2] mpowr[M_, n_] := mpowr[M, n - 1].M Now repeat this whole experiment with 2 pi/7, whose rotation matrix is merely the square of the previous. RotationMatrix[2 \[Pi]/7.] // InputForm {{0.6234898018587336, -0.7818314824680298}, {0.7818314824680298, 0.6234898018587336}} Nest[%.# &, IdentityMatrix@2, 999999] // InputForm // tim 0.215342,1 {{1.0000000000023095, 3.3821279110668456*^-12}, {-3.3821279110668456*^-12, 1.0000000000023095}} This error will eventually balloon exponentially! Addition chains are no help: mpowr[%%, 999999] // InputForm // tim 0.000178,1 {{1.0000000000544191, 7.38031857849819*^-11}, {-7.38031857849819*^-11, 1.0000000000544191}} MatrixPower[%%%, 999999] // InputForm // tim 0.000069,1 {{1.0000000000302447, 5.991646068181922*^-11}, {-5.991646068181922*^-11, 1.0000000000302447}} But a rotation is three skews: In[271]:= Inactive[Dot][{{1,-Tan[t/2]},{0,1}},{{1,0},{Sin[t],1}},{{1,-Tan[t/2]},{0,1}}] Out[271]= {{1,-Tan[t/2]},{0,1}}.{{1,0},{Sin[t],1}}.{{1,-Tan[t/2]},{0,1}} In[272]:= Activate[%] Out[272]= {{1-Sin[t] Tan[t/2],-Tan[t/2]-Tan[t/2] (1-Sin[t] Tan[t/2])},{Sin[t],1-Sin[t] Tan[t/2]}} In[273]:= FullSimplify[%] Out[273]= {{Cos[t],-Sin[t]},{Sin[t],Cos[t]}} And skews are reversible, hence lossless. The three skews: In[448]:= %271/.t->2\[Pi]/7. Out[448]= {{1,-0.481575},{0,1}}.{{1,0},{0.781831,1}}.{{1,-0.481575},{0,1}} exactly equal the 2pi/7 matrix: In[449]:= Activate[%]//InputForm Out[449] {{0.6234898018587336, -0.7818314824680298}, {0.7818314824680298, 0.6234898018587336}} But if we multiply one-at-a-time, In[450]:= Nest[#.%%[[1]].%%[[2]].%%[[3]]&,IdentityMatrix@2 ,999999]//InputForm//tim During evaluation of In[450]:= 15.1245,1 Out[450] {{0.999999999999998, 5.5067062021407764*^-14}, {-5.3512749786932545*^-14, 0.9999999999999998}} No loss! Go for two million. In[452]:= Nest[#.%448[[1]].%448[[2]].%448[[3]]&,IdentityMatrix@2 ,2*999999]//InputForm//tim During evaluation of In[452]:= 20.0325,1 Out[452] {{0.999999999999998, 5.5067062021407764*^-14}, {-5.3512749786932545*^-14, 0.9999999999999998}} No loss. But merely squaring the one million matrix In[454]:= %450.%450//InputForm Out[454] {{0.999999999999996, 1.101341240428154*^-13}, {-1.0702549957386496*^-13, 0.9999999999999996}} doubled the error! Floating point matrix multiplication is *seriously* nonassociative. --rwg
participants (1)
-
Bill Gosper