Precision in Mata matrix definitions

[John Mullahy]

Colleagues: When defining Mata matrixes, do you have any favorite tricks you use to avoid the nonzero sum(d) result here:

: d=.25,.1,-.1,-.25

: sum(d)
  -2.77556e-17

: sum(d)==0
  0

: quadsum(d)
  0

: quadsum(d)==0
  1

I suppose a related question is: Is it ever preferable to use sum(.) instead of quadsum(.)?

[Maarten Buis]

You could do this by transforming the matrix in one of integers

: foo = .25, .1, -.1, -.25

: sum(foo)
  -2.77556e-17

: quadsum(foo)
  0

:
: bar = 25, 10, -10, 25

: sum(bar)
  0

However, you can also improve the precision by sorting the matrix such that the smallest absolute numbers come first.

: sum(foo[order(abs(foo)', 1)])==0
  1

Adding the elements in a matrix is based on a running sum, and you loose precision if you add numbers that are different in size. So if you add the smallest numbers (in absolute terms) first, you minimize the difference between the previous sum and the current number, and thus minimize the loss of precision. However, sorting is computationally expensive, so I guess that in most cases quadsum() is more efficient than what I showed here.

[John Mullahy]

Thanks Maarten. I much appreciate these insights.