improve accuracy of logaddexp.reduce

[juliantaylor]

scipy has the logsumexp function which avoids unnecessary calls to log and thus probably saves a couple ulp of rounding errors.
The underflow avoidance should already be present in numpys current version.

I don't see a reason why we couldn't do this directly in logaddexp.reduce. The only possible issue is that the reduce loop is blocked to a few thousand elements by nditer so the maximum in the kernel is no global maximum. But I don't think that makes a big difference in numeric terms as only the largest terms even contribute to the sum.
But I would appreciate input of an expert in this regard.

A prototype implementation is available here:
https://github.com/juliantaylor/numpy/tree/logaddexp

[njsmith]

I'm not a numerical expert, but it seems like a great idea to me. Certainly
no worse than what we have now.

On Tue, Feb 4, 2014 at 6:23 PM, Julian Taylor notifications@github.comwrote:

scipy has the logsumexp function which avoids unnecessary calls to logand thus probably saves a couple ulp of rounding errors.
The underflow avoidance should also already be present in numpys current
version.

I don't see a reason why we couldn't do this directly in logaddexp.reduce.
The only issue is that the reduce loop is blocked to a few thousand
elements so the max in the kernel is no global max. But I don't think
that makes a big difference in numeric terms as only the largest terms even
contribute to the sum.
But I would appreciate input of an expert in this regard.

A prototype implementation is available here:
https://github.com/juliantaylor/numpy/tree/logaddexp

—
Reply to this email directly or view it on GitHubhttps://github.com//issues/4260
.

Nathaniel J. Smith
Postdoctoral researcher - Informatics - University of Edinburgh
http://vorpus.org

[charris]

I think the blocking is a plus, although probably not noticeable in practice.

[seberg]

I am pretty sure you are aware of it, but just saying that this would only affect reductions along the fast axis. It isn't that I don't like it, but it isn't like we can advertise a better accuracy if we implement it.

[njsmith]

@seberg: By fast axis I assume you mean the contiguous axis? Could we fix
it so custom reduction loops apply even on strided axes? This is useful for
things like pairwise summation as well...

On Tue, Feb 4, 2014 at 11:17 PM, seberg notifications@github.com wrote:

I am pretty sure you are aware of it, but just saying that this would only
affect reductions along the fast axis. It isn't that I like it, but it
isn't like we can advertise a better accuracy if we implement it.

—
Reply to this email directly or view it on GitHubhttps://github.com//issues/4260#issuecomment-34119822
.

Nathaniel J. Smith
Postdoctoral researcher - Informatics - University of Edinburgh
http://vorpus.org

[seberg]

well, it does not have to be contiguous but yes. I don't know the ufunc stuff too well, but NpyIter optimizes the iteration order to some degree, which means that np.sum(arr, axis=1) and np.sum(arr, axis=0) have typically the same speed. We would need some mechanism to signal that we have some special reduction stuff and don't want this optimization. For the simple things like sum, the iteration order is quite a large optimization though.

[njsmith]

Oh, right, the trick where it allocates a single accumulation buffer for
the whole output at once, and then updates it in multiple passes across the
fast axis.

I guess the ideal would be to extend the ufunc reduction loop call
signature so that we can register a reduction function that knows how to
handle a full 2d reduction directly. Something to put on the "someday"
list, I guess...

On Thu, Feb 6, 2014 at 6:28 PM, seberg notifications@github.com wrote:

well, it does not have to be contiguous but yes. I don't know the ufunc
stuff too well, but NpyIter optimizes the iteration order to some degree,
which means that np.sum(arr, axis=1) and np.sum(arr, axis=0) have
typically the same speed. We would need some mechanism to signal that we
have some special reduction stuff and don't want this optimization. For the
simple things like sum, the iteration order is quite a large optimization
though.

—
Reply to this email directly or view it on GitHubhttps://github.com//issues/4260#issuecomment-34353993
.

Nathaniel J. Smith
Postdoctoral researcher - Informatics - University of Edinburgh
http://vorpus.org