“Half Precision” 16-bit Floating Point Arithmetic
The floating point arithmetic format that requires only 16 bits of storage is becoming increasingly popular. Also known as half precision or binary16, the format is useful when memory is a scarce resource.
Contents
Background
The IEEE 754 standard, published in 1985, defines formats for floating point numbers that occupy 32 or 64 bits of storage. These formats are known as binary32 and binary64, or more frequently as single and double precision. For many years MATLAB used only double precision and it remains our default format. Single precision has been added gradually over the last several years and is now also fully supported.
A revision of IEEE 754, published in 2008, defines a floating point format that occupies only 16 bits. Known as binary16, it is primarily intended to reduce storage and memory bandwidth requirements. Since it provides only "half" precision, its use for actual computation is problematic. An interesting discussion of its utility as an image processing format with increased dynamic range is provided by Industrial Light and Magic. Hardware support for half precision is now available on many processors, including the GPU in the Apple iPhone 7. Here is a link to an extensive article about half precision on the NVIDIA GeForce GPU.
Floating point anatomy
The format of a floating point number is characterized by two parameters, p, the number of bits in the fraction and q, the number of bits in the exponent. I will consider four precisions, quarter, half, single, and double. The quarter-precision format is something that I just invented for this blog post; it is not standard and actually not very useful.
The four pairs of characterizing parameters are
p = [4, 10, 23, 52];
q = [3, 5, 8, 11];
With these values of p and q, and with one more bit for the sign, the total number of bits in the word, w, is a power of two.
w = p + q + 1
w =
8 16 32 64
Normalized numbers
Most floating point numbers are normalized, and are expressed as
$$ x = \pm (1+f)2^e $$
The fraction $f$ is in the half open interval
$$ 0 \leq f < 1 $$
The binary representation of $f$ requires at most p bits. In other words $2^p f$ is an integer in the range
$$ 0 \leq 2^p f < 2^p $$
The exponent $e$ is an integer in the range
$$ -b+1 \leq e \leq b $$
The quantity $b$ is both the largest exponent and the bias.
$$ b = 2^{q-1} - 1 $$
b = 2.^(q-1)-1
b =
3 15 127 1023
The fractional part of a normalized number is $1+f$, but only $f$ needs to be stored. That leading $1$ is known as the hidden bit.
Subnormal
There are two values of the exponent $e$ for which the biased exponent, $e+b$, reaches the smallest and largest values possible to represent in q bits. The smallest is
$$ e + b = 0 $$
The corresponding floating point numbers do not have a hidden leading bit. These are the subnormal or denormal numbers.
$$ x = \pm f 2^{-b} $$
Infinity and Not-A-Number
The largest possible biased exponent is
$$ e + b = 2^q-1 $$.
Quantities with this exponent field represent infinities and NaN, or Not-A-Number.
The percentage of floating point numbers that are exceptional because they are subnormal, infinity or NaN increases as the precision decreases. Exceptional exponents are only $2$ values out of $2^q$. For double precision this is $2/2^{11}$, which is less than a tenth of a percent, but for half precision it is $2/2^5$, which is more than 6 percent. And fully one-fourth of all my toy quarter precision floating point numbers are exceptional.
Encode the sign bit with s = 0 for nonnegative and s = 1 for negative. And encode the exponent with an offsetting bias, b. Then a floating point number can be packed in w bits with
x = [s e+b 2^p*f]
Precision and range
epsilon
If a real number cannot be expressed with a binary expansion requiring at most p bits, it must be approximated by a floating point number that does have such a binary representation. This is roundoff error. The important quantity characterizing precision is machine epsilon, or eps. In MATLAB, eps(x) is the distance from x to the next larger (in absolute value) floating point number. With no argument, eps is simply the difference between 1 and the next larger floating point number.
format shortg
eps = 2.^(-p)
eps =
0.0625 0.00097656 1.1921e-07 2.2204e-16
This tells us that double precision is good for about 16 decimal digits of accuracy, single for about 7 decimal digits, half for about 3, and quarter for barely more than one.
realmax
If a real number, or the result of an arithmetic operation, is too large to be represented, it overflows and is replaced infinity. The largest floating point number that does not overflow is
realmax = 2.^b.*(2-eps)
realmax =
15.5 65504 3.4028e+38 1.7977e+308
realmin
Underflow and representation of very small numbers is more complicated. The smallest normalized floating point number is
realmin = 2.^(-b+1)
realmin =
0.25 6.1035e-05 1.1755e-38 2.2251e-308
tiny
But there are numbers smaller than realmin. IEEE 754 introduced the notion of gradual underflow and denormal numbers. In the 2008 revised standard their name was changed to subnormal.
Think of roundoff in numbers near underflow. Before 754 floating point numbers had the disconcerting property that x and y could be unequal, but their difference could underflow so x-y becomes 0. With 754 the gap between 0 and realmin is filled with numbers whose spacing is the same as the spacing between realmin and 2*realmin. I like to call this spacing, and the smallest subnormal number, tiny.
tiny = realmin.*eps
tiny =
0.015625 5.9605e-08 1.4013e-45 4.9407e-324
Floating point integers
flintmax
It is possible to do integer arithmetic with floating point numbers. I like to call such numbers flints. When we write the numbers $3$ and $3.0$, they are different descriptions of the same integer, but we think of one as fixed point and the other as floating point. The largest flint is flintmax.
flintmax = 2./eps
flintmax =
32 2048 1.6777e+07 9.0072e+15
Technically all the floating point numbers larger than flintmax are integers, but the spacing between them is larger than one, so it is not safe to use them for integer arithmetic. Only integer-valued floating point numbers between 0 and flintmax are allowed to be called flints.
Table
Let's collect all these anatomical characteristics together in a new MATLAB table.
T = [w; p; q; b; eps; realmax; realmin; tiny; flintmax]; T = table(T(:,1), T(:,2), T(:,3), T(:,4), ... 'variablenames',{'quarter','half','single','double'}, ... 'rownames',{'w','p','q','b','eps','realmax','realmin', ... 'tiny','flintmax'}); disp(T)
quarter half single double
________ __________ __________ ___________
w 8 16 32 64
p 4 10 23 52
q 3 5 8 11
b 3 15 127 1023
eps 0.0625 0.00097656 1.1921e-07 2.2204e-16
realmax 15.5 65504 3.4028e+38 1.7977e+308
realmin 0.25 6.1035e-05 1.1755e-38 2.2251e-308
tiny 0.015625 5.9605e-08 1.4013e-45 4.9407e-324
flintmax 32 2048 1.6777e+07 9.0072e+15
fp8 and fp16
Version 3.1 of Cleve's Laboratory includes code for objects @fp8 and @fp16 that begin to provide full implementations of quarter-precision and half-precision arithmetic.
The methods currently provided are
methods(fp16)
Methods for class fp16: abs eps isfinite mrdivide rem subsref binary eq le mtimes round svd ctranspose fix lt ne sign tril diag fp16 lu norm single triu disp ge max plus size uminus display gt minus realmax subsasgn double hex mldivide realmin subsindex
These provide only partial implementations because the arithmetic is not done on the compact forms. We cheat. For each individual scalar operation, the operands are unpacked from their short storage into old fashioned doubles. The operation is then carried out by existing double precision code and the results returned to the shorter formats. This simulates the reduced precision and restricted range, but requires relatively little new code.
All of the work is done in the constructors @fp8/fp8.m and @fp16/fp16.m and what we might call the "deconstructors" @fp8/double.m and @fp16/double.m. The constructors convert ordinary floating point numbers to reduced precision representations by packing as many of the 32 or 64 bits as will fit into 8 or 16 bit words. The deconstructors do the reverse by unpacking things.
Once these methods are available, almost everything else is trivial. The code for most operations is like this one for the overloaded addition.
type @fp16/plus.m
function z = plus(x,y) z = fp16(double(x) + double(y)); end
Wikipedia test suite
The Wikipedia page about half-precision includes several 16-bit examples with the sign, exponent, and fraction fields separated. I've added a couple more.
0 01111 0000000000 = 1 0 00101 0000000000 = 2^-10 = eps 0 01111 0000000001 = 1+eps = 1.0009765625 (next smallest float after 1) 1 10000 0000000000 = -2 0 11110 1111111111 = 65504 (max half precision) = 2^15*(2-eps) 0 00001 0000000000 = 2^-14 = r ~ 6.10352e-5 (minimum positive normal) 0 00000 1111111111 = r*(1-eps) ~ 6.09756e-5 (maximum subnormal) 0 00000 0000000001 = r*eps ~ 5.96046e-8 (minimum positive subnormal) 0 00000 0000000000 ~ r*eps/2 (underflow to zero) 0 00000 0000000000 = 0 1 00000 0000000000 = -0 0 11111 0000000000 = infinity 1 11111 0000000000 = -infinity 0 11111 1111111111 = NaN 0 01101 0101010101 = 0.333251953125 ~ 1/3
This provides my test suite for checking fp16 operations on scalars.
clear
zero = fp16(0);
one = fp16(1);
eps = eps(one);
r = realmin(one);
tests = {'1','eps','1+eps','-2','2/r*(2-eps)', ...
'r','r*(1-eps)','r*eps','r*eps/2', ...
'zero','-zero','1/zero','-1/zero','zero/zero','1/3'};
Let's run the tests.
for t = tests(:)' x = eval(t{:}); y = fp16(x); z = binary(y); w = double(y); fprintf(' %18s %04s %19.10g %19.10g %s\n', ... z,hex(y),double(x),w,t{:}) end
0 01111 0000000000 3C00 1 1 1 0 00101 0000000000 1400 0.0009765625 0.0009765625 eps 0 01111 0000000001 3C01 1.000976563 1.000976563 1+eps 1 10000 0000000000 C000 -2 -2 -2 0 11110 1111111111 7BFF 65504 65504 2/r*(2-eps) 0 00001 0000000000 0400 6.103515625e-05 6.103515625e-05 r 0 00000 1111111111 03FF 6.097555161e-05 6.097555161e-05 r*(1-eps) 0 00000 0000000001 0001 5.960464478e-08 5.960464478e-08 r*eps 0 00000 0000000001 0001 5.960464478e-08 5.960464478e-08 r*eps/2 0 00000 0000000000 0000 0 0 zero 0 00000 0000000000 0000 0 0 -zero 0 11111 0000000000 7C00 Inf Inf 1/zero 1 11111 0000000000 FC00 -Inf -Inf -1/zero 1 11111 1111111111 FFFF NaN NaN zero/zero 0 01101 0101010101 3555 0.3333333333 0.3332519531 1/3
Matrix operations
Most of the methods in @fp8 and @fp16 handle matrices. The 4-by-4 magic square from Durer's Melancholia II provides my first example.
clear
format short
M = fp16(magic(4))
M =
16 2 3 13
5 11 10 8
9 7 6 12
4 14 15 1
Let's see how the packed 16-bit elements look in binary.
B = binary(M)
B =
4×4 string array
Columns 1 through 3
"0 10011 0000000000" "0 10000 0000000000" "0 10000 1000000000"
"0 10001 0100000000" "0 10010 0110000000" "0 10010 0100000000"
"0 10010 0010000000" "0 10001 1100000000" "0 10001 1000000000"
"0 10001 0000000000" "0 10010 1100000000" "0 10010 1110000000"
Column 4
"0 10010 1010000000"
"0 10010 0000000000"
"0 10010 1000000000"
"0 01111 0000000000"
Check that the row sums are all equal. This matrix-vector multiply can be done exactly with the flints in the magic square.
e = fp16(ones(4,1)) Me = M*e
e =
1
1
1
1
Me =
34
34
34
34
fp16 backslash
I've overloaded mldivide, so I can solve linear systems and compute inverses. The actual computation is done by lutx, a "textbook" function that I wrote years ago, long before this half-precision project. But now the MATLAB object system insures that every individual arithmetic operation is done on unpacked fp16 numbers.
Let's generate a 5-by-5 matrix with random two-digit integer entries.
A = fp16(randi(100,5,5))
A =
76 71 83 44 49
75 4 70 39 45
40 28 32 77 65
66 5 96 80 71
18 10 4 19 76
I am going to use fp16 backslash to invert A. So I need the identity matrix in half precision.
I = fp16(eye(5))
I =
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
0 0 0 0 1
Now the overloaded backslash calls lutx to compute the inverse.
X = A\I
X =
-0.0044 0.0346 0.0125 -0.0254 -0.0046
0.0140 -0.0116 0.0026 -0.0046 -0.0002
0.0071 -0.0176 -0.0200 0.0238 0.0008
-0.0060 -0.0020 0.0200 0.0006 -0.0125
0.0003 -0.0052 -0.0072 0.0052 0.0174
Compute the residual.
AX = A*X R = I - AX
AX =
1.0000 -0.0011 -0.0002 -0.0007 -0.0001
-0.0001 0.9990 -0.0001 -0.0007 -0.0001
-0.0001 -0.0005 1.0000 -0.0003 -0.0002
-0.0002 -0.0011 -0.0001 0.9995 -0.0002
-0.0000 -0.0001 0.0001 -0.0001 1.0000
R =
0 0.0011 0.0002 0.0007 0.0001
0.0001 0.0010 0.0001 0.0007 0.0001
0.0001 0.0005 0 0.0003 0.0002
0.0002 0.0011 0.0001 0.0005 0.0002
0.0000 0.0001 -0.0001 0.0001 0
Both AX and R are what I expect from arithmetic that is accurate to only about three decimal digits.
Although I get a different random A every time I publish this blog post, I expect that it has a modest condition number.
kappa = cond(A)
kappa = 15.7828
Since A is not badly conditioned, I can invert the computed inverse and expect to get close to the original integer matrix.
Z = X\I
Z = 76.1250 71.0000 83.1250 44.1250 49.1250 75.1250 4.0234 70.1875 39.1250 45.1250 40.0625 28.0000 32.0625 77.0000 65.0625 66.1250 5.0234 96.1875 80.1250 71.1250 18.0156 10.0000 4.0156 19.0156 76.0000
fp16 SVD
I have just nonchalantly computed cond(A). But cond isn't on the list of overload methods for fp16. I was pleasantly surprised to find that matlab\matfun\cond.m quietly worked on this new datatype. Here is the core of that code.
dbtype cond 34:43, dbtype cond 47
34 if p == 2 35 s = svd(A); 36 if any(s == 0) % Handle singular matrix 37 c = Inf(class(A)); 38 else 39 c = max(s)./min(s); 40 if isempty(c) 41 c = zeros(class(A)); 42 end 43 end 47 end
So it is correctly using the singular value decomposition, and I have svd overloaded. The SVD computation is handled by a 433 line M-file, svdtx, that, like lutx, was written before fp16 existed.
Let's compute the SVD again.
[U,S,V] = svd(A)
U =
-0.5210 -0.4841 0.6802 -0.0315 0.1729
-0.4260 -0.2449 -0.3572 -0.4561 -0.6504
-0.4058 0.4683 0.1633 0.6284 -0.4409
-0.5786 0.0268 -0.5620 0.1532 0.5703
-0.2174 0.6968 0.2593 -0.6104 0.1658
S =
267.5000 0 0 0 0
0 71.1875 0 0 0
0 0 55.5000 0 0
0 0 0 37.3750 0
0 0 0 0 16.9531
V =
-0.4858 -0.3108 -0.0175 -0.3306 -0.7471
-0.2063 -0.2128 0.9238 0.2195 0.1039
-0.5332 -0.5205 -0.2920 -0.0591 0.5967
-0.4534 0.2891 -0.2050 0.7993 -0.1742
-0.4812 0.7095 0.1384 -0.4478 0.2126
Reconstruct A from its half precision SVD. It's not too shabby.
USVT = U*S*V'
USVT = 75.9375 71.0000 83.0625 44.0313 49.0000 75.0000 4.0117 70.0625 38.9688 45.0000 40.0313 28.0469 32.0313 77.0625 65.0625 66.0000 4.9688 96.0625 80.0000 71.0000 18.0313 10.0234 4.0156 19.0313 76.0000
Finally, verify that we've been working all this time with fp16 objects.
whos
Name Size Bytes Class Attributes A 5x5 226 fp16 AX 5x5 226 fp16 B 4x4 1576 string I 5x5 226 fp16 M 4x4 208 fp16 Me 4x1 184 fp16 R 5x5 226 fp16 S 5x5 226 fp16 U 5x5 226 fp16 USVT 5x5 226 fp16 V 5x5 226 fp16 X 5x5 226 fp16 Z 5x5 226 fp16 e 4x1 184 fp16 kappa 1x1 8 double
Calculator
I introduced a calculator in my blog post about Roman numerals. Version 3.1 of Cleve's Laboratory also includes a fancier version of the calculator that computes in four different precisions -- quarter, half, single, and double -- and displays the results in four different formats -- decimal, hexadecimal, binary, and Roman.
I like to demonstrate the calculator by clicking on the keys
1 0 0 0 / 8 1 =
because the decimal expansion is a repeating .123456790.
Thanks
Thanks to MathWorkers Ben Tordoff, Steve Eddins, and Kiran Kintali who provided background and pointers to work on half precision.
Comments
To leave a comment, please click here to sign in to your MathWorks Account or create a new one.