## Disabling subnormals in MATLAB

Suppose that we want to compute the following expression, somewhat contrived in complexity for the purpose of example:

$y = \frac{1}{s}\sum_{i=1}^{10^7} \frac{s}{2^{50}}$

The following MATLAB code implements this formula, and measures the time required to evaluate it:

tic;
scale = 1;
y = 0;
for i = 1:10000000
x = scale;
for j = 1:50
x = x * 0.5;
end
y = y + x;
end
y = y / scale;
toc


Now suppose that we execute this same code again, but this time changing the “scale” factor $s$ to a much smaller value: scale = realmin, corresponding to $s=2^{-1022}$, the smallest positive normalized floating-point number. Inspection of the formula above suggests that the value of $y$ should not depend on the changed value of $s$ (as long as it is non-zero); we may spend most of our time working with much smaller numbers, but the end result should be the same.

And indeed, execution of the modified code confirms that we get exactly the same result… but it takes nearly 20 times longer to do so on my laptop than the original version with $s=1$. And there are more complex– and less contrived– calculations where the difference in performance is even greater.

The problem is that these smaller numbers are subnormal, small enough in magnitude to require a slightly different encoding than “normal” numbers which make up most of the range of representable floating-point numbers. Arithmetic operations can be significantly more expensive when required to recognize and manipulate this “special case” encoding.

Depending on your application, there may be several approaches to handling this problem:

1. Rewrite your code to prevent encountering subnormals in the first place. In the above contrived example, this is easy to do: just shift the “scale,” or magnitude, of all values involved in the computation away from the subnormal range (and possibly shifting back only at the end if necessary). This can not only result in faster code, but more accurate results, since subnormal numbers have fewer “significant” mantissa bits in their representation.
2. Disable subnormal numbers altogether, so that for any floating-point operation, input arguments or output results that would otherwise be treated as subnormal are instead “flushed” to zero.

We have seen above how to manage (1). For (2), the following MEX function does nothing but set the appropriate processor flags to disable subnormals. I have only tested this on Windows 7 with an Intel laptop, compiling in MATLAB R2017b with both Microsoft Visual Studio 2015 as well as the MinGW-w64 MATLAB Add-On (edit: a reader has also tried this on Linux Mint 19 with MATLAB R2018b and GCC 7.3.0):

#include <xmmintrin.h>
#include <pmmintrin.h>
#include "mex.h"

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
//_MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);
_mm_setcsr((_mm_getcsr() & ~0x8000) | (0x8000));
//_MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);
_mm_setcsr((_mm_getcsr() & ~0x0040) | (0x0040));
}


After running this MEX function in a MATLAB session, re-running the modified calculation above gets all of the speed back… but now at a different cost: instead of the correct value $y=10^7/2^{50}$, every term in the sum has been flushed to zero, resulting in a final incorrect value of $y=0$.

If there is any moral to this story, it’s that you’re a test pilot. First, this was a very simple test setup; it’s an interesting question whether any MATLAB built-in functions might reset these flags back to the slower subnormal support, and whether it is feasible in your application to reset them back again, possibly repeatedly as needed. And second, even after any algorithm refactoring to minimize the introduction of subnormals, can your application afford the loss of accuracy resulting from flushing them to zero? MathWorks’ Cleve Moler seems to think the answer is always yes. I think the right answer is, it depends.

This entry was posted in Uncategorized. Bookmark the permalink.

This site uses Akismet to reduce spam. Learn how your comment data is processed.