2

I traced a bug in a larger program back to this problem. I am getting larger than expected errors for changing the order of multiplicative operations using numpy. Here is a self-contained example (the actual code at the end is short, but I copy in explicit matrices to make it run):

import numpy

M = numpy.array(
[[ 1.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00],
 [ 0.00000000e+00,  4.24142041e-01,  1.49421321e-08,  7.87081062e-09,  2.58494178e-08, -5.09534880e-02, -1.09351587e-08,  5.95629388e-10,  1.34707444e-08,  0.00000000e+00, -4.88069109e-01,  1.84060546e-08,  2.76013885e-08, -1.40665727e-08,  5.86332400e-02,  1.87260208e-08,  1.59941364e-08,  1.04757263e-08],
 [ 0.00000000e+00,  1.49421321e-08,  5.26397502e-16,  2.77281383e-16,  9.10651097e-16, -1.79504429e-09, -3.85235532e-16,  2.09834728e-17,  4.74561875e-16,  0.00000000e+00, -1.71942236e-08,  6.48428290e-16,  9.72371404e-16, -4.95552353e-16,  2.06559485e-09,  6.59700404e-16,  5.63458643e-16,  3.69050154e-16],
 [ 0.00000000e+00,  7.87081062e-09,  2.77281383e-16,  1.46058758e-16,  4.79688058e-16, -9.45544689e-10, -2.02923914e-16,  1.10531040e-17,  2.49976819e-16,  0.00000000e+00, -9.05710624e-09,  3.41561448e-16,  5.12199406e-16, -2.61033613e-16,  1.08805797e-09,  3.47499066e-16,  2.96803444e-16,  1.94398219e-16],
 [ 0.00000000e+00,  2.58494178e-08,  9.10651097e-16,  4.79688058e-16,  1.57539771e-15, -3.10537007e-09, -6.66445336e-16,  3.63007470e-17,  8.20977096e-16,  0.00000000e+00, -2.97454652e-08,  1.12176052e-15,  1.68217190e-15, -8.57289960e-16,  3.57341403e-09,  1.14126092e-15,  9.74765703e-16,  6.38445141e-16],
 [ 0.00000000e+00, -5.09534880e-02, -1.79504429e-09, -9.45544689e-10, -3.10537007e-09,  6.12119924e-03,  1.31367425e-09, -7.15547905e-11, -1.61828196e-09,  0.00000000e+00,  5.86332433e-02, -2.21117595e-09, -3.31583970e-09,  1.68986064e-09, -7.04379147e-03, -2.24961448e-09, -1.92142480e-09, -1.25848122e-09],
 [ 0.00000000e+00, -1.09351587e-08, -3.85235532e-16, -2.02923914e-16, -6.66445336e-16,  1.31367425e-09,  2.81928419e-16, -1.53564166e-17, -3.47300464e-16,  0.00000000e+00,  1.25833156e-08, -4.74541799e-16, -7.11614349e-16,  3.62662007e-16, -1.51167232e-09, -4.82791114e-16, -4.12358131e-16, -2.70083410e-16],
 [ 0.00000000e+00,  5.95629388e-10,  2.09834728e-17,  1.10531040e-17,  3.63007470e-17, -7.15547905e-11, -1.53564166e-17,  8.36451786e-19,  1.89171798e-17,  0.00000000e+00, -6.85403182e-10,  2.58479141e-17,  3.87610671e-17, -1.97539108e-17,  8.23395879e-11,  2.62972477e-17,  2.24608192e-17,  1.47112284e-17],
 [ 0.00000000e+00,  1.34707444e-08,  4.74561875e-16,  2.49976819e-16,  8.20977096e-16, -1.61828196e-09, -3.47300464e-16,  1.89171798e-17,  4.27830628e-16,  0.00000000e+00, -1.55010671e-08,  5.84575999e-16,  8.76619656e-16, -4.46754122e-16,  1.86219076e-09,  5.94738120e-16,  5.07973517e-16,  3.32708899e-16],
 [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  1.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00],
 [ 0.00000000e+00, -4.88069109e-01, -1.71942236e-08, -9.05710624e-09, -2.97454652e-08,  5.86332433e-02,  1.25833156e-08, -6.85403182e-10, -1.55010671e-08,  0.00000000e+00,  5.61631322e-01, -2.11802316e-08, -3.17614944e-08,  1.61866991e-08, -6.74704944e-02, -2.15484234e-08, -1.84047870e-08, -1.20546371e-08],
 [ 0.00000000e+00,  1.84060546e-08,  6.48428290e-16,  3.41561448e-16,  1.12176052e-15, -2.21117595e-09, -4.74541799e-16,  2.58479141e-17,  5.84575999e-16,  0.00000000e+00, -2.11802316e-08,  7.98748561e-16,  1.19778898e-15, -6.10432541e-16,  2.54444623e-09,  8.12633805e-16,  6.94081037e-16,  4.54604285e-16],
 [ 0.00000000e+00,  2.76013885e-08,  9.72371404e-16,  5.12199406e-16,  1.68217190e-15, -3.31583970e-09, -7.11614349e-16,  3.87610671e-17,  8.76619656e-16,  0.00000000e+00, -3.17614944e-08,  1.19778898e-15,  1.79618282e-15, -9.15393662e-16,  3.81560581e-09,  1.21861104e-15,  1.04083144e-15,  6.81716412e-16],
 [ 0.00000000e+00, -1.40665727e-08, -4.95552353e-16, -2.61033613e-16, -8.57289960e-16,  1.68986064e-09,  3.62662007e-16, -1.97539108e-17, -4.46754122e-16,  0.00000000e+00,  1.61866991e-08, -6.10432541e-16, -9.15393662e-16,  4.66514627e-16, -1.94455784e-09, -6.21044146e-16, -5.30441833e-16, -3.47425038e-16],
 [ 0.00000000e+00,  5.86332400e-02,  2.06559485e-09,  1.08805797e-09,  3.57341403e-09, -7.04379147e-03, -1.51167232e-09,  8.23395879e-11,  1.86219076e-09,  0.00000000e+00, -6.74704944e-02,  2.54444623e-09,  3.81560581e-09, -1.94455784e-09,  8.10543757e-03,  2.58867824e-09,  2.21102355e-09,  1.44816055e-09],
 [ 0.00000000e+00,  1.87260208e-08,  6.59700404e-16,  3.47499066e-16,  1.14126092e-15, -2.24961448e-09, -4.82791114e-16,  2.62972477e-17,  5.94738120e-16,  0.00000000e+00, -2.15484234e-08,  8.12633805e-16,  1.21861104e-15, -6.21044146e-16,  2.58867824e-09,  8.26760426e-16,  7.06146767e-16,  4.62507011e-16],
 [ 0.00000000e+00,  1.59941364e-08,  5.63458643e-16,  2.96803444e-16,  9.74765703e-16, -1.92142480e-09, -4.12358131e-16,  2.24608192e-17,  5.07973517e-16,  0.00000000e+00, -1.84047870e-08,  6.94081037e-16,  1.04083144e-15, -5.30441833e-16,  2.21102355e-09,  7.06146767e-16,  6.03129082e-16,  3.95033217e-16],
 [ 0.00000000e+00,  1.04757263e-08,  3.69050154e-16,  1.94398219e-16,  6.38445141e-16, -1.25848122e-09, -2.70083410e-16,  1.47112284e-17,  3.32708899e-16,  0.00000000e+00, -1.20546371e-08,  4.54604285e-16,  6.81716412e-16, -3.47425038e-16,  1.44816055e-09,  4.62507011e-16,  3.95033217e-16,  2.58736061e-16]]
)

A = numpy.array(
[[-0.00000000e+00, -7.49420654e-01,  2.36645133e-08,  2.83214009e-09, -6.52157725e-08,  9.00302357e-02,  2.07762711e-08, -3.93769333e-08, -3.96766769e-08,  0.00000000e+00, -6.51261884e-01,  2.35425627e-08,  3.46140156e-08, -1.89780598e-08,  7.82380887e-02,  2.51868179e-08,  2.13865316e-08,  1.44744812e-08],
 [ 3.05311332e-16,  3.63603065e-02,  1.77351557e-01,  1.52401227e-02, -1.31437146e-06,  3.02666853e-01, -5.56572048e-07, -1.60534678e-06, -4.65237579e-07,  4.42851120e-16,  1.11596525e-01, -8.16894013e-07,  9.17216520e-07, -2.96530450e-07,  9.28940840e-01,  3.66607356e-07, -9.91357211e-09, -2.27619767e-07],
 [-1.35525272e-20, -3.78275577e-05,  2.02944409e-04, -1.63077293e-05,  3.96391244e-01, -3.14375043e-04,  2.19758080e-01,  4.39817640e-01,  5.14025579e-01,  5.28070693e-20,  7.95260108e-06,  2.07648708e-01,  2.80669103e-01, -2.70743918e-02,  6.62922300e-05, -3.97663420e-01, -9.72844144e-02,  2.16082768e-01],
 [ 3.17637355e-22,  3.02225698e-07, -2.20579495e-06,  1.36032858e-05,  1.89121855e-01,  2.46318863e-06,  3.29601248e-01,  5.12358908e-01, -6.11878058e-01, -4.09855638e-21,  5.33328462e-08, -1.13839294e-01, -3.51810943e-01,  1.99126758e-01,  8.07593444e-07, -1.95429714e-01,  6.34840080e-02, -9.20601029e-03],
 [ 2.11758237e-22,  2.93821317e-07, -1.25434775e-06,  5.01041209e-07, -2.41225155e-01,  2.43924254e-06, -2.10410277e-01,  5.85509545e-01, -1.53207165e-01,  2.17127990e-22, -6.54394382e-08,  4.61741591e-02,  3.94747797e-01, -1.31509913e-01, -6.95848293e-07,  4.33611439e-01, -3.91114733e-01, -1.22481544e-01],
 [-3.17637355e-22,  1.30494444e-07, -6.28897449e-07, -5.69416759e-07, -4.89409426e-01,  8.16591863e-07,  2.25530164e-01,  1.36440113e-01,  1.56457139e-01,  6.09145504e-22, -1.46265910e-07, -7.06195107e-01,  1.76837517e-01, -1.84706615e-01, -1.10677502e-06, -2.65190566e-01,  1.75625222e-01,  3.53552474e-02],
 [-1.32348898e-23,  9.76074136e-09, -6.64706365e-08,  8.97102007e-08,  1.04453697e-02,  3.77554871e-08,  1.06780632e-01, -5.51530507e-02, -2.67072454e-02,  8.34309184e-23,  9.76785235e-09, -1.66380632e-01,  8.89817377e-02,  3.52794263e-01, -1.01989239e-07,  3.79383206e-01, -2.21311259e-02,  8.24771602e-01],
 [-2.64697796e-23, -3.81071470e-09, -1.52470457e-08,  2.31867964e-07,  1.61167689e-01,  1.21476062e-08,  4.62920198e-01, -6.30833114e-02,  1.55154734e-01, -6.81907803e-23,  3.36355606e-09, -1.04616530e-01,  3.04210903e-01,  4.74962450e-01, -3.26436018e-08,  3.66884819e-01,  1.89001475e-01, -4.81947218e-01],
 [ 3.30872245e-23, -3.82486279e-09,  4.57705305e-08, -1.53730823e-07, -1.09868035e-01,  9.49535030e-08, -6.11367953e-01,  2.09616196e-01,  6.92500565e-03, -1.51367999e-23, -1.37726681e-08,  2.39736500e-02,  1.68012777e-01,  5.33544789e-01, -7.45547417e-08, -2.35647324e-01,  4.47745356e-01, -2.63191932e-02],
 [ 0.00000000e+00,  2.95812127e-07, -4.47630753e-07, -6.75579326e-08, -5.05647373e-01,  2.12597159e-06,  3.29042472e-01,  1.47754834e-01,  4.59445954e-02, -4.42071975e-22,  7.64987025e-09,  5.88363526e-01, -7.18829924e-02, -1.25810121e-01, -3.22069856e-07,  1.14059812e-01,  4.67535857e-01,  1.15511821e-01],
 [-1.58818678e-22, -9.75341754e-08,  8.33486809e-07, -7.58797257e-07, -4.09124432e-01, -1.22408647e-06,  1.92397506e-01, -2.35640312e-01, -1.54747228e-01,  4.33075192e-22, -5.97477948e-08,  2.27856071e-01,  1.87299705e-01,  4.18606298e-01, -3.80248666e-07, -4.18737252e-01, -5.27560777e-01, -1.53387616e-02],
 [-3.17637355e-22, -2.73493680e-07,  1.20952799e-06, -1.15938669e-07,  2.19496103e-01, -2.38765220e-06,  5.63854677e-02, -2.33009173e-01, -5.12386378e-01, -1.42993918e-21,  1.90409524e-08,  7.60732865e-02,  6.64392889e-01, -2.87091516e-01, -3.51089068e-07, -1.51359887e-01,  2.62855799e-01,  1.00893426e-01],
 [ 0.00000000e+00, -9.82865059e-02,  5.33893628e-01, -8.64663760e-02, -1.52855583e-04, -8.18146794e-01, -8.26226436e-05, -1.62722671e-04, -1.97549957e-04,  1.91248595e-16,  2.01206803e-02, -7.98176370e-05, -1.09471810e-04,  1.05532690e-05,  1.67486593e-01,  1.53257887e-04,  3.75957532e-05, -8.28268558e-05],
 [ 3.46944695e-18,  9.86836968e-03, -3.21503728e-02, -9.96041972e-01,  9.89973139e-06,  8.21452559e-02,  7.87619622e-06,  1.42732304e-05,  4.19866736e-07,  1.08567995e-16, -5.53252953e-04,  2.12382866e-06,  5.52438245e-08,  2.06607768e-06, -4.60524623e-03, -8.61743850e-06, -7.93397037e-07,  3.40561040e-06],
 [ 0.00000000e+00,  5.60976585e-02,  8.26117626e-01,  1.38453289e-02,  1.66897114e-06,  4.66962542e-01,  5.21437516e-07,  1.01948728e-06,  6.95013869e-07, -3.76249560e-17, -3.69824778e-02,  1.99982217e-08,  2.01540087e-07,  1.03736187e-07, -3.07846242e-01, -1.07869609e-06, -3.39760103e-07,  4.32423200e-07]]
)

B = numpy.array(
[[ 0.00000000e+00, -1.00000000e+00,  1.06205003e-05],
 [ 6.51261884e-01,  2.02210886e-17,  2.58415811e-17],
 [ 2.29433542e-08,  1.34184985e-16, -2.74111890e-17],
 [ 1.20854772e-08,  7.60639281e-18,  2.96257867e-16],
 [ 3.96912801e-08,  5.98385402e-23, -1.93603571e-17],
 [-7.82380933e-02,  1.68322376e-16,  2.15108010e-16],
 [-1.67907242e-08, -1.39105759e-21, -1.31411424e-16],
 [ 9.14577381e-10, -9.94043901e-22, -1.02202642e-16],
 [ 2.06840670e-08,  2.09868713e-21,  2.04206049e-16],
 [ 0.00000000e+00,  1.06205003e-05,  1.00000000e+00],
 [-7.49420658e-01, -1.59224140e-17,  4.48478893e-17],
 [ 2.82621401e-08,  3.37439644e-21,  3.18670444e-16],
 [ 4.23813971e-08, -8.80703978e-22, -8.05157812e-17],
 [-2.15989497e-08, -1.22983410e-22, -1.35186018e-17],
 [ 9.00302037e-02, -1.32539872e-16,  3.73318253e-16],
 [ 2.87534420e-08,  1.37971615e-21,  1.40919484e-16],
 [ 2.45586865e-08,  5.08971450e-22,  5.55657960e-17],
 [ 1.60852746e-08,  1.11154797e-23, -1.54979207e-18]]
)

Z = A @ M @ B
print("original matrix\n", Z, '\n')

norm = max(max(row) for row in abs(Z))
print("largest element\n", norm, '\n')

Z1 = Z / norm
print("matrix normalized one way\n", Z1, '\n')

Z2 = A @ (M/norm) @ B
print("matrix normalized another way\n", Z1, '\n')

print("unexpected deviation \n", Z1-Z2, '\n')

Below is the output. I expect the matrices normalized in the two different "ways" to be the same to within machine tolerance (roughly 1e-16), relative to the largest element, which is ~1. So why am I getting errors on the order of 1e-7? I checked that all the dtypes are float64. Thanks for any help!

original matrix
 [[-2.70976453e-10  2.77813727e-27 -2.66468642e-27]
 [ 3.14803157e-10 -3.05306629e-16  4.42854363e-16]
 [ 2.40140485e-14  1.35530880e-20  5.28069254e-20]
 [ 4.02675253e-16 -3.17680884e-22 -4.09855301e-21]
 [-5.36313469e-16 -2.11755931e-22  2.17130239e-22]
 [ 3.12972364e-16  3.17643824e-22  6.09142131e-22]
 [-2.99648801e-17  1.32357759e-23  8.34307778e-23]
 [-1.00370207e-16  2.64690554e-23 -6.81910614e-23]
 [ 4.76605198e-18 -3.30873853e-23 -1.51364485e-23]
 [-4.00415286e-16 -4.69502499e-27 -4.42071975e-22]
 [-1.67417098e-17  1.58823277e-22  4.33073505e-22]
 [ 5.36721887e-16  3.17622168e-22 -1.42994255e-21]
 [ 2.42670193e-11  2.03115564e-21  1.91248595e-16]
 [-6.86052564e-12 -3.46829390e-18  1.08568032e-16]
 [-8.03322460e-11 -3.99595601e-22 -3.76249560e-17]] 

largest element
 3.148031574328576e-10 

matrix normalized one way
 [[-8.60780609e-01  8.82499813e-18 -8.46461149e-18]
 [ 1.00000000e+00 -9.69833439e-07  1.40676595e-06]
 [ 7.62827434e-05  4.30525797e-11  1.67745857e-10]
 [ 1.27913346e-06 -1.00914135e-12 -1.30194152e-11]
 [-1.70364704e-06 -6.72661395e-13  6.89733358e-13]
 [ 9.94184324e-07  1.00902363e-12  1.93499371e-12]
 [-9.51860850e-08  4.20446097e-14  2.65025226e-13]
 [-3.18834816e-07  8.40812894e-14 -2.16614922e-13]
 [ 1.51397846e-08 -1.05104998e-13 -4.80822639e-14]
 [-1.27195448e-06 -1.49141610e-17 -1.40428063e-12]
 [-5.31815181e-08  5.04516152e-13  1.37569619e-12]
 [ 1.70494442e-06  1.00895484e-12 -4.54233866e-12]
 [ 7.70863275e-02  6.45214508e-12  6.07518033e-07]
 [-2.17930649e-02 -1.10173415e-08  3.44875931e-07]
 [-2.55182466e-01 -1.26935068e-12 -1.19518992e-07]] 

matrix normalized another way
 [[-8.60780609e-01  8.82499813e-18 -8.46461149e-18]
 [ 1.00000000e+00 -9.69833439e-07  1.40676595e-06]
 [ 7.62827434e-05  4.30525797e-11  1.67745857e-10]
 [ 1.27913346e-06 -1.00914135e-12 -1.30194152e-11]
 [-1.70364704e-06 -6.72661395e-13  6.89733358e-13]
 [ 9.94184324e-07  1.00902363e-12  1.93499371e-12]
 [-9.51860850e-08  4.20446097e-14  2.65025226e-13]
 [-3.18834816e-07  8.40812894e-14 -2.16614922e-13]
 [ 1.51397846e-08 -1.05104998e-13 -4.80822639e-14]
 [-1.27195448e-06 -1.49141610e-17 -1.40428063e-12]
 [-5.31815181e-08  5.04516152e-13  1.37569619e-12]
 [ 1.70494442e-06  1.00895484e-12 -4.54233866e-12]
 [ 7.70863275e-02  6.45214508e-12  6.07518033e-07]
 [-2.17930649e-02 -1.10173415e-08  3.44875931e-07]
 [-2.55182466e-01 -1.26935068e-12 -1.19518992e-07]] 

unexpected deviation 
 [[-1.85496826e-07 -3.37450100e-24 -4.46082542e-24]
 [-2.90512590e-08  0.00000000e+00 -2.11758237e-22]
 [ 3.36542580e-12  6.46234854e-27  0.00000000e+00]
 [-1.03382940e-13  2.01948392e-28 -1.61558713e-27]
 [-6.32965290e-14  0.00000000e+00 -2.01948392e-28]
 [ 7.05435233e-15  2.01948392e-28 -4.03896783e-28]
 [-6.90238453e-16  6.31088724e-30  0.00000000e+00]
 [ 6.56248400e-15 -1.26217745e-29  0.00000000e+00]
 [ 1.33816164e-16 -1.26217745e-29 -6.31088724e-30]
 [-1.09952637e-13 -2.98596179e-30 -4.03896783e-28]
 [-2.49491876e-14  0.00000000e+00 -2.01948392e-28]
 [ 1.41929230e-14  0.00000000e+00 -8.07793567e-28]
 [ 1.07272085e-08  2.81919955e-25  0.00000000e+00]
 [ 2.16769630e-10  1.65436123e-24  5.29395592e-23]
 [ 4.15586621e-09 -5.77572400e-26  0.00000000e+00]]
  • 1
    Matrix multiplication is the sum of products, If the arrays are large that sum can accumulate precision errors. Recently I explored a SO where I found that the use of `math.fsum` made more of a difference than the precision of indiviual numbers (using `mpmath`). – hpaulj Jun 28 '23 at 16:25
  • @hpaulj Share the link? – jared Jun 28 '23 at 16:29
  • 2
    https://stackoverflow.com/questions/76551058/comparing-accurrate-dot-product-with-mpmath – hpaulj Jun 28 '23 at 16:37
  • @hpaulj Thanks. I will read that carefully, but I don't think it is the root of the problem. For example, n x n matrices in a product A.B.C would accumulate <~n^2 times the error of scalars a x b x c, assuming that a is on the order of elements of A, and so forth, because each element of A.B.C is a sum of n^2 terms. The <~ is because the errors of individual flops have random sign. n here is on the order of 1 to 10, so I should be getting 1e-14 precision at least. – JackOfAllTrades Jun 28 '23 at 17:15
  • 1
    Relative error `((Z1-Z2)/Z1` is on the order of 1e-7 for most terms, or 1e-16 for the small ones. `M/norm` has a wide range of values from 3e9 to 1e-6. – hpaulj Jun 28 '23 at 17:46

1 Answers1

2

Thanks to hint from @hpaulj, I figured it out. M has large values that are projected out by A and B (actually, this was the intent). Dividing M by norm first causes these values to blow up and degrade the precision of the result because the smaller desired result is the leftovers after larger terms cancel.

So this is no longer about numpy or matmul, but it is still relevant to programming so I guess it is good to leave up the question.

PS - to the trailing point of the discussion, this means the logic of my reply to the first comment by @hpaulj breaks down at the point "assuming that a is on the order of elements of A, and so forth". That becomes ambiguous when the large part is projected out by a cofactor.