# Computation with less than 8 bits in gemmlowp

## Introduction

We assume familiarity with gemmlowp's low-precision uint8 computation paradigm,
which is described in [low-precision.md](low-precision.md).

This document is about the possibility of further reducing precision below 8
bits.

That allows to get higher arithmetic throughput on some architectures, at the
cost of decreased accuracy.

## The past, present, and future of less-than-8-bit computation in gemmlowp

A meta note is needed here as to how this fits with the general gemmlowp design.

### The past

Less-than-8-bit computation was initially designed and implemented in gemmlowp
as a drop-in replacement for regular 8bit computation, a plain optimization. The
idea was that to automatically requantize 8bit operands to less-than-8bit during
the O(N^2) packing stage, then take advantage of the lower bit depth during the
O(N^3) compute stage. For large enough matrices, that should be worth it.

### The present

TODO(benoitjacob): update this documentation. This 'present' state just
became the past (February 2017).

At the moment, this less-than-8-bit mode of gemmlowp is not much used in
practice, because the implicit requantization of operands from 8bit to
less-than-8bit turned out to be more expensive than initially expected, both in
terms of speed and accuracy:

1.  Speed: the O(N^2) requantization is only negligible compared to the O(N^3)
    compute kernel when the matrix size N is large enough; in practice, smaller
    matrix sizes turned out to be very important, making the requantization
    approach slower than expected.

2.  Accuracy: As neural networks were optimized for size, their sensitivity to
    numerical accuracy increased. Then the approach of requantizing
    already-quantized data turned out to be more wasteful of accuracy than we
    could afford.

### The future

Less-than-8bit still probably has good prospects; what should be dropped is the
requantization. In other words, in the future, we might have neural networkds
trained right away for some bit depth lower than 8 bits. The resulting values
would probably still be stored as 8 bits (unless the bit depth eventually
becomes very low). Thus, no particular work would be needed in the packing
stage; no overhead or loss of accuracy would be incurred anymore.

In other words: the design of less-than-8-bit kernels is probably useful in the
long run; what is on the way out is requantization and packing/unpacking-stage
aspects.

With that said, the rest of this page retains its old content about the present
approach:

## Public interface

### The BitDepthSetting parameter in the EightBitIntGemm interface

Accessing less-than-8-bit computation via the EightBitIntGemm is very simple:
EightBitIntGemm takes a BitDepthSetting enum which allows to choose among a
fixed set of supported bit-depth combinations.

### The BitDepthParams parameter in the public/gemmlowp.h interface

The public/gemmlowp.h interface exposes more extensive control over
quantization, by means of a BitDepthParams template parameter, which is a type
parameter, carrying information about: 1. The LHS and RHS bit depth, which can
be set arbitrarily and independently; 2. The 'RoundingStrategy', which is the
heuristic used to choose a rounding mode, based on the accumulation size (a.k.a.
the "depth" dimension of the Gemm). Details can be seen in public/bit_depth.h.

### How does BitDepth{Setting,Params} affect input/output uint8 matrix data?

Input/output matrix data is all uint8's, ranging from 0 to 255, regardless of
the BitDepth{Setting,Params}.

So the BitDepth{Setting,Params} is only an internal detail. It only means to
allow gemmlowp to use lower precision internally, but the input/output data
format is unaffected.

As far as the API contract goes, the only thing that the
BitDepth{Setting,Params} does is to relax the accuracy requirement. With
standard 8bit/8bit computation, gemmlowp is required to return the exact result
as specified in [low-precision.md](low-precision.md). With lower bit depths,
gemmlowp is no longer required to return an exact result.

## Implementation

Here we refer to the 3 stages of computation as described in
[design.md](design.md), namely: packing, computation kernel, unpacking.

The general idea is that at the packing stage, we requantize input (Lhs/Rhs)
data to less-than-8-bit depths by scaling them, thus shrinking the range of the
packed matrix entries; for instance, if the Rhs bit depth is to be 5 bits then
packed Rhs matrix entries will be in the range [0 ... 31]. This then allows the
GEMM kernel to use narrower accumulators without risking overflow, thus
achieving higher arithmetic throughput. Finally, at the unpacking stage, it only
remains to scale the result values to compensate for the scalings applied
earlier.

Let us go into more detail for each of those stages:

### Packing stage

The packing stage is where most of the work specific to the BitDepthParams takes
place.

Here, we have to scale input matrix values from their original range of [0 ...
255] to the range specified by the BitDepthParams, which is [0 ... (2^N)-1]
where N is the number of bits for the matrix at hand (Lhs or Rhs). For example,
for a bit depth of 5 bits, we need to scale down to [0 ... 31].

This scaling is what we call "requantization". The pedantic name matches the
fact that this is actually quite nontrivial to do correctly i.e. in such a way
that the result accuracy will be good enough for real-world applications. See
the section below on requantization details.

Concretely, this work happens in PackingRegisterBlock::Pack(), which calls
Requantize(). This is in internal/pack.h. This code can be overridden for
specific architectures, see internal/pack_neon.h.

This requantization work is costly and makes packing slower. This means that, at
least in our approach, less-than-8-bit computation is only interesting for
large-enough, square-enough GEMMs where packing is only a small fraction of the
overall cost. In cases where packing overhead is more prevalent (highly
rectangular cases), less-than-8-bit is probably a waste of time as long as we
treat it as an internal computation detail. What might help there, might be if
we shrunk the input/output data format to lower memory bandwidth usage.

### Computation kernel stage

In principle, the computation kernel stage simply doesn't have to care about the
bit depth at all. In fact, on architectures where we do not have specific
optimized kernels for less-than-8-bit cases, we simply use our standard kernel
there, and that's correct!

However, while the kernel doesn't have to know about the fact that the operands
are on less than 8 bits, it can use that information to make special
optimizations that would be incorrect in the general 8-bit case and become
correct here thanks to the more restricted range of inputs. That's the whole
point of this less-than-8-bit computation idea.

With Lhs entries guaranteed to be smaller than 2^N, and Rhs entries guaranteed
to be smaller than 2^M, each product is thus guaranteed to be smaller than
2^(M+N). Thus, one may accumulate 2^(16-(M+N)) such products and still be
guaranteed that such an accumulator will be smaller than 2^16, and thus can be
stored as a uint16 without risking overflow.

For example, in the L7R5 case, the Lhs enties are on 7 bits (N=7) and the Rhs
entries are on 5 bits (M=5), so each product fits in 12 bits and one can thus
accumulate 16 ( = 2^(16-12)) such products into uint16 accumulators with no risk
of overflow.

This means that a computation kernel may use uint16 accumulators for several
loop iterations (16 in the above example), provided that it is allowed to assume
that inputs are in such restricted range.

After this fixed number of loop iterations, the kernel must accumulate the local
uint16 accumulators back into global uint32 accumulators.

On SIMD architectures with suitable uint16 arithmetic, this in principle allows
to multiply arithmetic throughput by up to 2x, since twice more accumulators now
fit in each SIMD vector register. This is partially offset by the cost of
accumulating back into global uint32 accumulators every several loop iterations,
but our experience on ARM NEON has been that we still get quite close to a 2x
speedup. See internal/kernel_neon.h, specifically
NEON32Kernel12x4Depth2Assuming12BitProducts.

### Unpacking stage

At the unpacking stage, it only remains to scale the result values to compensate
for the scaling of the inputs. This is easier because now we are expanding the
range instead of shrinking it, so we don't need to worry about ways to minimize
a loss of accuracy. We simply need to multiply result values by a constant
fraction, rounding to nearest.

Since the inputs were scaled by factors of (2^lhs_bits - 1)/255 and
(2^rhs_bits - 1)/255 respectively, the scaling of the outputs needs to be by the
following factor:

                 255 * 255
    -----------------------------------
    (2^lhs_bits - 1) * (2^rhs_bits - 1)

This is done by a MultiplyByConstantFraction function, see internal/unpack.h

## Requantization details

Here we go into more detail on the Requantize() function used at the packing
stage to requantize input matrix data. See this function in internal/pack.h.

It depends on the bit depth and on a rounding mode, and requantizes an input
value in [0 ... 255] to the range [0 ... (2^N)-1] specified by the bit depth N.

### Naive, bad rounding, that's plainly biased

Naive and inaccurate ways to achieve this requantization include: 1. By shifting
bits rights by (8-N) bits; 2. By multiplying by ((2^N) - 1) and dividing by 255.

Both of those are biased in some way: 1. has the wrong "derivative", since it
approximates (((2^N) - 1) / 255) by ((2^N) / 256) ; 2. has bias since it
effectively implements rounding towards 0.

In practice, both of the above requantization functions give results that are
too inaccurate in practice for the neural network that we tried (GoogLeNet).

### Round-to-nearest rounding: unbiased in principle but not in practice

The simplest fix is to avoid the bias in 2. by rounding-to-nearest instead of
rounding towards 0. This can be achieved by doing

dst = (src * maxval + rounding_offset) / 255;

Where maxval = ((2^N) - 1) is the highest requantized value, and the
rounding_offset can be set to

rounding_offset = 127

to achieve rounding-to-nearest (while the above rounding towards 0 corresponded
to rounding_offset = 0).

In principle, rounding-to-nearest is unbiased and optimal in various ways.

In practice though, our input data is not random real numbers, but
already-quantized 8-bit values. That means that even in the best case, there
would be at most 255 different possible input values; in practice, we generally
see the input values distributed non-uniformly in that range, so that a majority
of input values tend to be in a much smaller range. See test/test_data.cc.

Having a large part of the input values in a very small finite set, means that
the corresponding rounding errors are also in a very small finite set, which can
be small enough that the mean of these rounding errors is significantly
different from 0. That rounding-to-nearest is "unbiased" only means that over a
sufficiently large set of input values, the bias would become arbitrarily close
to 0; here, the set of input values is effectively small enough that the
resulting bias is significant.

This leads to biasing the matrix product entries, resulting in an error that
grows linearly with the depth dimension of the GEMM.

### Probabilistic rounding: unbiased even on small finite input distributions

To address that, we can instead use probabilistic rounding. The idea is that for
instance if we have to round the value 3.8 to the nearest integer, we can round
it to 3 with 20% probability and to 4 with probability 80%. If that value 3.8
occurs many times, the mean requantized value will thus tend to 3.8.

This amounts to keeping the above requantization formula,

dst = (src * maxval + rounding_offset) / 255;

but now the rounding_offset is a random value in [0 .. 254].

This guarantees zero bias no matter how small the distribution of input values
is.

On the other hand, the variance of the error term here is higher than with
rounding-to-nearest --- one can check that it is 2x higher.

So the error term coming from the Central Limit Theorem, which grows with the
square root of the accumulator depth i.e. the GEMM depth, will be 2x higher
here.

Still, for large enough GEMM depth, that is better than rounding-to-nearest
which has an error term growing linearly with the GEMM depth.

### Switching between rounding-to-nearest and probabilistic rounding

Thus, for fixed input values and bit depths, we expect that probabilistic
rounding will give more accurate results for large enough GEMM depths, while
rounding-to-nearest will be more accurate for smaller GEMM depths.

That is why use switch between these rounding modes based on GEMM depth, see
ChooseRoundingMode in internal/bit_depth_util.h.

It is based on a constant, kProbabilisticRoundingThreshold, defined in
internal/common.h and empirically determined. See the comment there. It would be
nice to better understand the statistics here and come up with better heuristics
for this switching.

### Choice of pseudorandom number generator

We provide two PRNGs. The first is an 8-bit Xorshift. It is fast, naturally
produces values ranging over an interval of width 255, which is what we need
here (as opposed to an interval of width 256), and turns out, from empirical
tests, to produce better results than a linear congruential generator (LCG).
That's unfortunate, as a 8-bit LCG performs better (we confirmed that on various
ARM devices) but we need as perfect un-biased-ness as we can get.

The second is an "add-mod" sequence generator, which generates non-random values
in the sequence x = (x + 97) % 255. This generates a low-discrepancy sequence
that minimizes the "clumpiness" of the random offsets (Thus, for example,
quantizing a 3x3 matrix will have a maximum additive error of about 200 from the
random offsets). While not random, this sequence performs well empirically for
many quantizations. (For information about why 97 is a good value, see
https://en.wikipedia.org/wiki/Low-discrepancy_sequence#Additive_recurrence and
http://mollwollfumble.blogspot.com/2011/03/subrandom-numbers.html 97/255 = 0.38;
0.382 is the best choice. For discrete numbers, the choice must be relatively
prime to the modulus. 97 is prime, so it is safely relatively prime to 255. 107
is another near-optimal choice.

The low-discrepancy sequence generator is the default.

More details and results are given in a comment on the default PRNG in
internal/pack.h. Interested users can change the PRNG used by setting
DefaultRoundingGenerator in bit_depth_util.h.