CODEWITHSUNDEEP
pythonarraysnumpy
Nico Schlömer
Nico Schlömer

Reputation: 58831

math.fsum for arrays of multiple dimensions

I have a numpy array of dimension (i, j) in which I would like to add up the first dimension to receive a array of shape (j,). Normally, I'd use NumPy's own sum

import numpy

a = numpy.random.rand(100, 77)
numpy.sum(a, axis=0)

but in my case it doesn't cut it: Some of the sums are very ill-conditioned, so the computed results only have a few correct digits.

math.fsum is fantastic at keeping the errors at bay, but it only applies to iterables of one dimension. numpy.vectorize doesn't do the job either.

How to efficiently apply math.fsum to an array of multiply dimensions?

Upvotes: 2

Views: 291

Answers (2)

Nico Schlömer
Nico Schlömer

Reputation: 58831

Check out the signature keyword to vectorize. 

_math_fsum_vec = numpy.vectorize(math.fsum, signature='(m)->()')

Unfortunately, it's slower than the for solution:

enter image description here

Code to reproduce the plot:

import math
import numpy
import perfplot


_math_fsum_vec = numpy.vectorize(math.fsum, signature='(m)->()')


def fsum_vectorize(a):
    return _math_fsum_vec(a.T).T


def fsum_for(a):
    return numpy.array([math.fsum(row) for row in a.T])


perfplot.save(
    'fsum.png',
    setup=lambda n: numpy.random.rand(n, 100),
    kernels=[fsum_vectorize, fsum_for],
    n_range=[2**k for k in range(12)],
    logx=True,
    logy=True,
    )

Upvotes: 0

jsmolka
jsmolka

Reputation: 800

This one works fast enough for me.

import numpy
import math

a = numpy.random.rand(100, 77)
a = numpy.swapaxes(a, 0, 1)
a = numpy.array([math.fsum(row) for row in a])

Hopefully it's the axis you are looking for (returns 77 sums).

Upvotes: 2

Related Questions