doc/quantization.md

# Building a quantization paradigm from first principles

**TLDR:** If you prefer example code over theory, look at
[doc/quantization_example.cc](quantization_example.cc).

## Overview

gemmlowp allows to perform calculations on matrices on uint8 values, but these
matrices are only useful insofar as they somehow approximate matrices of real
numbers. By a _quantization paradigm_ we mean a correspondence between matrices
of quantized 8bit values and matrices of real numbers. The choice of a
quantization paradigm affects the calculations that gemmlowp itself needs to
perform, specifically, it affects how one goes from internal 32bit accumulator
to final 8bit outputs.

The part of gemmlowp transforming internal 32bit accumulator to final
8bit outputs is the "output pipeline" described in [output.md](output.md).

gemmlowp's `GemmWithOutputPipeline` entry point allows specifying an arbitrary
output pipeline, allowing the user to implement their own preferred quantized
arithmetic paradigm.

In the present document, our purpose is to show how, reasoning from first
principles and some domain-specific knowledge of neural networks, we can arrive
naturally at some specific quantization paradigm, and how that can be
implemented using a specific output pipeline.

We also aim to show how that differs from the older, legacy quantization
paradigm implemented by gemmlowp's legacy interfaces and why the change to the
newer quantization paradigm described in this document was useful as far as some
applications of gemmlowp were concerned.

## Quantization as an affine map.

In order for arithmetic on real values to map directly to arithmetic on
quantized uint8 values, the mapping between real and quantized uint8 values must
be affine, which means it must be of the form

```
real_value = A * quantized_value + B             (1)
```

for some constants A, B, or equivalently, of the form

```
real_value = C * (quantized_value + D)           (2)
```

for some constants C, D. Indeed, anything else than such an affine map would
mean that the result of the quantized calculations do no longer readily provide
an approximation to the result of the real-numbers calculation.

## Domain-specific constraint: the real value 0 must be exactly representable.

Here a domain-specific constrain from neural networks appears: for some neural
network layers, it is very useful for optimized implementations that the
real-value 0 be exactly representable.

For instance, in a Convolutional or Pooling layer with padding, it is useful to
be able to implement the padding by zero-padding the input array, so that
optimized loops do not need to become more complex to avoid overrunning the
array bounds.

In order for such zero-padding to be feasible in a quantized implementation of
such layers, it is important that the real value '0' be exactly representable in
quantized form, i.e. that it correspond exactly to some quantized value, which
we call the _zero-point_.

Indeed, if '0' were not exactly representable, then we would have to use some
quantized value for padding, that does not exactly correspond to the real value
'0'. That would typically introduce inaccuracy in the result. In fact, using
always the same such value would be worse: it would introduce _bias_ in the
result.

## The final form of the quantization equation

Now let us phrase what this constraint — that the real value 0 be exactly
representable — means in either quantization equations, (1) and (2).

In equation (1), plugging `real_value = 0` and `quantized_value = zero_point`,
we get:

```
0 = A * zero_point + B
```

equivalently:

```
zero_point = -B / A
```

We are thus left with a rather awkward constraint: the real number `-B / A` must
somehow be guaranteed to be exactly integral, so that the special uint8 value
`zero_point` can be exactly equal to it. Quite awkward!

Now let us look at equation (2). Plugging `real_value = 0` and
`quantized_value = zero_point`, we get:

```
0 = C * (zero_point + D)
```

Conveniently, the constant `C` plays no role anymore, so this equation
simplifies to:

```
0 = zero_point + D
```

In other words, `D = -zero_point`. This suggests rewriting the quantization
equation (2) into the following form (3), which will be the final form that we
will consistently use:

```
real_value = scale * (quantized_value - zero_point)        (3)
```

To go from (2) to (3), we merely renamed `C` to `scale` and `D` to
`-zero_point`.

With this quantization equation (3), the condition that 0 be exactly
representable is vacuously satisfied: `zero_point` is by definition one of the
possible `quantized_value`'s, and equation (3) maps it to a `real_value` of
exactly 0.

Note that the final quantizaton equation (3) depends on two constants, one
integral, the other an arbitrary positive real number:

*   `zero_point` is integral, more specifically is one of the possible quantized
    values (i.e. typically is a uint8 value).
*   `scale` is a positive real number. Thus at this stage we have not yet shown
    how to eliminate all usage of floating-point arithmetic. That will come
    below.

## Quantizing a matrix multiplication

Now that we know — equation (3) — how real numbers are to correspond
to quantized values (typically uint8), we turn to applying this knowledge to
rewriting a multiplication of matrices of real numbers, by the equivalent
multiplication of matrices of quantized values.

Say that we have two matrices of real values `lhs_real_matrix`,
`rhs_real_matrix`. Each entry of their product is the sum (accumulation) of many
products of individual matrix entries, say `lhs_real_value * rhs_real_value`.

Now suppose that we have already quantized these two matrices according to the
above equation (3), with some already-known quantization parameters `lhs_scale`,
`rhs_scale`, `lhs_zero_point`, `rhs_zero_point`, so that their matrix entries
are quantized as

```
lhs_real_value[i] = lhs_scale * (lhs_quantized_value[i] - lhs_zero_point)
rhs_real_value[i] = rhs_scale * (rhs_quantized_value[i] - rhs_zero_point)
```

We then rewrite the matrix product accumulator accordingly:

```
result_real_value
  = Sum_over_i(lhs_real_value[i] * rhs_real_value[i])
  = Sum_over_i(
        lhs_scale * (lhs_quantized_value[i] - lhs_zero_point) *
        rhs_scale * (rhs_quantized_value[i] - rhs_zero_point)
    )
  = lhs_scale * rhs_scale * Sum_over_i(
        (lhs_quantized_value[i] - lhs_zero_point) *
        (rhs_quantized_value[i] - rhs_zero_point)
    )                                                      (4)
```

Now our goal is to represent this result itself as a quantized matrix, i.e.
still according to equation (3), for some pre-established quantization
parameters `result_scale` and `result_zero_point`, as

```
result_real_value = result_scale *
    (result_quantized_value - result_zero_point)
```

Here we need to keep in mind that our goal is to specify what the quantized
matrix multiplication should do, i.e. how to compute `result_quantized_value`.
The last equation above is equivalent to

```
result_quantized_value = result_zero_point +
    result_real_value / result_scale
```

Now we can use equation (4) above to plug into this the expression of
result_real_value in terms of the quantized operands, and we obtain:

```
result_quantized_value = result_zero_point +
    (lhs_scale * rhs_scale / result_scale) *
        Sum_over_i(
            (lhs_quantized_value[i] - lhs_zero_point) *
            (rhs_quantized_value[i] - rhs_zero_point)
        )                                                  (5)
```

Equation (5) is the conclusion of this general discussion of how to specify what
"quantized matrix multiplication" should actually compute, in order to be able
to replace real matrix multiplications.

## Implementation of quantized matrix multiplication

Having obtained the mathematical form (5) of quantized matrix multiplication, we
now turn to its actual implementation.

The inner-most part of (5),

```
int32_accumulator =
    Sum_over_i(
        (lhs_quantized_value[i] - lhs_zero_point) *
        (rhs_quantized_value[i] - rhs_zero_point)
)
```

is the "kernel" accumulation loop. It is where the bulk of the computational
cost goes. Luckily, it only involves integers: the quantized operands matrix
entries, and their `zero_point` quantization parameters. Typically, all of these
values are uint8. Typically, the above differences of uint8 values would be
represented as signed int16; their products as signed int32.

It is out of scope of the present doc to discuss how to avoid the overhead of
having to subtract these `zero_point` constants in this inner loop; refer to
[this section of
low-precision.md](low-precision.md#efficient-handling-of-offsets) for that. The
gist of it is that a mathematical trick allows us to take the handling of these
`zero_point` constants out of this accumulation loop, so that it simplifies to

```
int32_accumulator =
    Sum_over_i(
      lhs_quantized_value[i] *
      rhs_quantized_value[i]
    )                                                      (6)
```

Anyway, the result is a `int32_accumulator` that we now plug back into the rest
of (5):

```
result_quantized_value = result_zero_point +
    (lhs_scale * rhs_scale / result_scale) * int32_accumulator       (7)
```

The difficulty here is of course that `(lhs_scale * rhs_scale / result_scale)`
is a positive real number, not an integer in general. It is a constant, though.
So what we have to implement here is the (approximate) scaling of a int32 value
by some arbitrary positive constant multiplier.

Moreover, it is safe to assume that this positive constant multiplier is smaller
than one — each of the `scale` values here is typically smaller than one,
as we are typically mapping the `[0..255]` quantized uint8 value range to an
interval of real values that is much narrower than that, typically within
`[-10,10]` in most neural networks. For example, a neural network using Relu6
activation functions will typically have real activation values in the interval
[0,6].

So how do we implement the multiplication of a int32 value by a positive real
constant that is smaller than one? Typically, by multiplying by a fixed-point
constant multiplier in the normalized interval `[1/2,1)`, and right-shifting
the result to achieve the correct multiplier.

At this point we have obtained the int32 value of the product

```
(lhs_scale * rhs_scale / result_scale) * int32_accumulator
```

Looking at (7), it only remains to add to it the integral value
`result_zero_point`, and we are done.

## How this is implemented in gemmlowp

The different parts of gemmlowp implementing aspects of the above discussion
are:

*   The packing stage (see [packing.md](packing.md)) implements the special
    mathematical trick to handle `lhs_offset`, `rhs_offset` that we alluded to
    above, see [this section of
    low-precision.md](low-precision.md#efficient-handling-of-offsets) for
    details. Thanks to is, the rest of the calculation can proceed as if
    `lhs_offset`, `rhs_offset` were 0.

*   The compute/kernel stage (see [kernel.md](kernel.md)) performs the core
    accumulation loop producing the `int32_accumulator`, see equation (6) above.

*   The unpacking stage feeds into the output pipeline (see
    [output.md](output.md)), which implements the rest of the evaluation of the
    above equation (5), that we discussed in the previous section.

Now, the point of gemmlowp's flexible output-pipelines mechanism (see
[output.md](output.md)) is to support different quantization paradigms, so we
now have to specify which particular flavor of output pipeline corresponds to
the particular quantization paradigm that we detailed above in this document.

The specific output pipeline stage implementing the present quantization
paradigm, i.e. implementing the precise computation detailed in the previous
section (equation (5)), is
`OutputStageQuantizeDownInt32ByFixedPoint`.

Please refer to the comment explaining it in
[public/output_stages.h](../public/output_stages.h).

## How this differs from the older legacy gemmlowp quantization paradigm

The difference between the older legacy quantization paradigm described in
[low-precision.md](low-precision.md) and the newer one described in this
document boils down to the difference between the legacy output stage
implementing it, `OutputStageQuantizeDownInt32ToUint8Scale`, and the new output
stage implementing the new paradigm,
`OutputStageQuantizeDownInt32ByFixedPoint`.

Please refer to the comments in
[public/output_stages.h](../public/output_stages.h) for details about these two
output stages and how they differ.

Issues with the old output stage `OutputStageQuantizeDownInt32ToUint8Scale` are:

1.  The int32 accumulators (inputs to the output stage) undergo a plain int32
    multiplication with a int32 multiplier, which may overflow. By contrast, in
    the newer `OutputStageQuantizeDownInt32ByFixedPoint`, this
    integer multiplication becomes a fixed-point multiplication and cannot
    overflow.

    *   In practice, to limit the risk of overflow, this pushes users to choose
        smaller values for this integer multiplier, which means limited
        multiplicative accuracy, which may cause multiplicative bias depending
        on how it is used.

2.  Note how the order of multiplying by the multipler and adding the
    `result_offset` are swapped. This reflects a quantizatin equation of the
    form (1) above, as opposed to the form (2)/(3) that the new quantization
    paradigm uses. As a result, it is essentially impossible to guarantee that 0
    is an exactly-representable value, which as discussed above is an issue at
    least in some convolutional neural network applications.

## Example code illustrating the new quantization paradigm

Example code showing how to perfom a quantized matrix multiplication in the
quantization paradigm discussed here is in
[doc/quantization_example.cc](quantization_example.cc).