This work was funded by the A*CRC and carried out in collaboration with Dr. John Gustafson and Dr. Marek Michalewicz. It was presented at the SuperComputing Frontiers 2016 in Singapore by Prof. Edelman.
The unum is a floating point format proposed by John Gustafson, proposed as an alternative to the now ubiquitious IEEE 754 formats. The proposal and justification are explained in his somewhat ambitiously-titled book The end of error, but the following slides give a good overview:
The two defining features of the unum format are:
For performing computation with the format, Gustafson proposes using interval arithmetic with a pair of unums, what he calls an ubound, providing the guarantee that the resulting interval contains the exact solution (he also suggests a more complicated ubox method, but I won’t discuss this further).
The Unums.jl package provides an implementation of the unum format and ubound arithmetic in Julia.
New instances can be constructed from existing numbers
julia> x = Unum22(1)
julia> y = Unum22(2.1)
As 2.1 (or, more correctly, the Float64 nearest 2.1) cannot be represented as an exact Unum, the resulting Unum represents an open interval.
The low-level bit representation of an unum is similar to that of an IEEE754 floating point number, but with extra places for the ubit and the size of the exponent and fraction:
|sign|exponent|fraction|ubit|exponent size|fraction size|
This is available via the print_bits function:
Standard arithmetic operations are defined, and Unums will automatically be promoted to Ubounds:
One of the Julia’s most powerful features are its generic methods. In particular, there are several linear algebra routines, such as LU and QR factorizations, which can be applied to an array of any numeric type, such as Quaternions, simply by defining the necessary elementary operations (+, -, *, / and sqrt) and a few simple methods such as abs and copysign.
We can make use of them here with Unums:
julia> X = map(Unum34,randn(10,6))
(0.76873779296875,0.7687454223632812) … (-1.520599365234375,-1.5205841064453125)
(2.5906982421875,2.590728759765625) … (-0.7313995361328125,-0.7313919067382812)
julia> Q,R = qr(X)
(-0.203277587890625,-0.2031707763671875) … (0.5075531005859375,0.8198089599609375)
(-0.6849288940429688,-0.684814453125) … (-0.4050254821777344,-0.11310958862304688)
(0.264190673828125,0.2642326354980469) (0.362701416015625,0.576568603515625) ,
(-3.782867431640625,-3.78271484375) … (0.4373779296875,0.4379119873046875)
0.0 … (-2.2342529296875,-2.19671630859375) )
We can check that the matrix reconstructed from the factors contains the original matrix X:
julia> Y = Q*R
(0.7685317993164062,0.7689743041992188) … (-1.89501953125,-1.15325927734375)
(2.5904541015625,2.59100341796875) … (-1.071014404296875,-0.38515472412109375)
Unfortunately, as with any interval arithmetic approach, the problem is that the intervals grow as the number of operations increase. In this case, as each column of the orthogonal matrix Q depends on the previous ones, the widths increase with column number:
julia> W = map(Unums.width,Q)
0.000106812 0.000665665 0.00137138 0.00675392 0.0597467 0.312256
6.48499e-5 0.000686646 0.00142288 0.0072937 0.0620613 0.31245
9.05991e-6 0.000112534 0.000808716 0.00331306 0.0285912 0.138903
1.52588e-5 0.000366211 0.000694275 0.00553513 0.0377731 0.184459
1.62125e-5 0.000196457 0.00105286 0.00414467 0.0401134 0.196717
0.000114441 0.000652313 0.00124741 0.00605774 0.0500031 0.291916
2.86102e-5 0.000366211 0.000835419 0.00431061 0.0451851 0.215823
1.81198e-5 0.000141144 0.000701904 0.00354385 0.0270338 0.133575
7.24792e-5 0.000396729 0.000806808 0.00379539 0.0399399 0.194012
4.19617e-5 0.000411987 0.000822067 0.00435543 0.0425873 0.213867
julia> using Gadfly
julia> spy(W, Scale.color_log10)
Internally, Unums.jl uses Julia’s BigFloat type (which wraps the MPFR library). This provides convenient per-operation rounding control and inexact flags, which are necessary
for computing the correct intervals and determining the state of the u-bit. As a result,
the performance is long way behind the corresponding machine floating point operations:
julia> X = ones(Ubound34,100,100);
julia> t = @elapsed X*X
So that’s around 105 “unops” (unum operations per second), compared with
i.e. 6×1010 flops (floating point operations per second). Although it is
not going to be possible to approach such speeds for unums without dedicated hardware
support, it should be possible to see some gain via a lower-level implementation based on integer operations and bit-twiddling, as this would avoid many of the bottlenecks involved with the use of BigFloats, such as memory allocation and garbage collection.
I think the key lesson here is that although interval arithmetic is an extremely powerful
tool, intervals cannot be used as simple substitutes for numbers. John Gustafson proposes
using fused operations for polynomial evaluation and dot products, and it should be possible
to test these ideas out in Julia, for example by overloading the dot function,
or via macros to detect reused bindings in expressions such as x*(2.0+x).
Fully leveraging the power of interval arithmetic does require a slightly different mindset, however,
using specialised algorithms such as interval Newton’s method.
Unfortunately, such approaches need to be designed from scratch, and as a result we can no
longer fully leverage Julia’s generic functions.
I didn’t really discuss the utility of variable-width encoding. As it stands, the Unum
format is far from the most efficient storage format, with a lot of redundant encodings.
This problem that has been recognised by John Gustafson in his recent presentation on
“Unums version 2.0”.
Given the constraints of existing computing architectures, it may be that a byte-aligned
format (similar to UTF-8) might be more amenable to optimization.
Need help with Julia?
We also provide training and consulting services
and build open source or proprietary packages
for our customers on a consulting basis. Email us:
Julia Computing was founded by all the creators
of the language to provide commercial support
to Julia users. We are based in Boston, New York,
San Francisco, London and Bangalore with
customers across the world.
© 2016 - 2020 Julia Computing, Inc. All Rights Reserved.