gemmlowp

gemmlowp implements low-precision General Matrix Multiplication (GEMM) using a three-stage computation scheme to balance memory efficiency and arithmetic performance.

Because the library uses uint8 inputs/outputs but accumulates in int32 to maintain precision, it must manage the transition between these types. The process follows these steps:

1. Pack: Reorder input matrix blocks (LHS/RHS) into a layout optimized for cache locality and SIMD loading.
2. Compute: Execute the multiplication using a kernel that operates on the packed blocks and accumulates results into a temporary int32 block.
3. Unpack: Convert the int32 accumulator block back into the uint8 destination matrix.

This approach minimizes the memory footprint of high-precision int32 values by only storing and processing them in small, manageable blocks.

1. Pack lhs/rhs blocks from the input matrices.
2. Compute the product of the packed blocks, using the kernel.
3. Unpack the result block into the output matrix.

To maintain high performance and low memory bandwidth, gemmlowp does not add lhs_offset and rhs_offset to every entry during the GEMM kernel execution. Instead, it uses a mathematical expansion of the matrix product:

$$(lhs + lhs_offset \cdot P) \cdot (rhs + rhs_offset \cdot Q) =$

1. lhs * rhs (The core GEMM kernel only computes this term)
2. lhs_offset * P * rhs (A rank-one update: lhs_offset times the column sums of rhs)
3. lhs * rhs_offset * Q (A rank-one update: rhs_offset times the row sums of lhs)
4. lhs_offset * rhs_offset * depth (A rank-zero update: a constant added to all entries, where depth is the number of columns in lhs)

Implementation Details:

• Packing Stage: gemmlowp computes the sum of each row of the lhs block and the sum of each column of the rhs block. These are stored in a buffer managed by sums_of_each_slice_handle_ in the PackedSideBlock class.
• Compute Kernel Stage: The kernel only computes the first term (lhs * rhs) to maximize speed and register availability.
• Unpacking Stage: The rank-one and rank-zero updates are applied during the unpacking stage (e.g., in UnpackResultImpl) to finalize the result.

Less-than-8-bit computation achieves higher throughput by reducing the range of operands, allowing for narrower accumulators (e.g., uint16) without overflow risk. The process follows three stages:

1. Packing Stage (Requantization)

Input matrices (range [0...255]) are scaled down to the range specified by the bit depth: [0 ... (2^N)-1]. For example, a 5-bit depth scales values to [0...31]. This "requantization" is computationally expensive and can make packing slower, meaning this mode is most effective for large, square matrices where packing overhead is negligible compared to the $O(N^3)$ compute cost.

2. Computation Kernel Stage

Because inputs are restricted to a smaller range, products are also smaller. If LHS entries are $2^N$ and RHS entries are $2^M$, each product is $\le 2^{M+N}$. This allows accumulating $2^{16-(M+N)}$ products into a uint16 accumulator before risking overflow. On SIMD architectures, using uint16 accumulators can potentially double arithmetic throughput compared to uint32 accumulators.

3. Unpacking Stage

To compensate for the scaling applied during packing, the result values are scaled back up. The output is scaled by the following factor:

$$\frac{255 \times 255}{(2^{lhs_bits} - 1) \times (2^{rhs_bits} - 1)}$$

This is typically implemented using a MultiplyByConstantFraction function.

The current recommended way to perform low-precision matrix multiplication in gemmlowp is using an "output pipeline." This design allows you to define a chain of transformations to be applied to internal 32-bit accumulators to produce the final 8-bit outputs. This approach is more flexible and avoids the overflow risks associated with legacy integer multiplication.

To use this paradigm, you provide:

1. A lhs matrix of uint8_t quantized values.
2. A rhs matrix of uint8_t quantized values.
3. An int32 lhs_offset to be added to each entry of the lhs matrix.
4. An int32 rhs_offset to be added to each entry of the rhs matrix.
5. An output pipeline to process the resulting int32 accumulators.

The computation follows these steps:

1. Cast lhs entries to int32 and add lhs_offset.
2. Cast rhs entries to int32 and add rhs_offset.
3. Compute the int32 matrix product.
4. Apply the output pipeline to the int32 accumulators to obtain final outputs.

During the packing stage of less-than-8-bit computation, Requantize() is used to map input matrix data from the standard 8-bit range [0 ... 255] to a smaller bit-depth range [0 ... (2^N)-1].

Choosing the correct rounding method is critical for accuracy:

1. Rounding-to-nearest: Uses the formula dst = (src * maxval + rounding_offset) / 255, where rounding_offset = 127. While theoretically unbiased, it can introduce significant bias in practice when input data is non-uniformly distributed or belongs to a small finite set. This bias causes errors that grow linearly with the GEMM depth.
2. Probabilistic rounding: Uses the same formula but with a random rounding_offset in the range [0 .. 254]. This guarantees zero bias regardless of the input distribution. However, it has a higher error variance (2x higher), meaning the error term grows with the square root of the GEMM depth.

Recommendation: Use probabilistic rounding for large GEMM depths and rounding-to-nearest for smaller GEMM depths.

The gemmlowp architecture is organized into a hierarchical pipeline designed to maximize cache hits and register usage. A typical execution flow follows this nested structure:

1. L2 Block Level: Subdivide matrices into blocks that fit in the L2 cache.
2. Packing: Pack the L2 blocks to ensure efficient traversal.
3. L1/Sub-block Level: Subdivide L2 blocks into smaller sub-blocks that fit in the L1 cache.
4. Kernel Execution: The inner-most loop (the GEMM kernel) performs the actual multiply-accumulate operations using architecture-specific instructions (often written in assembly or SIMD) to maximize register usage.

By organizing data this way, the library minimizes redundant memory accesses, which are the primary bottleneck in matrix multiplication.

allocate(some_lhs_L2_block);
allocate(some_rhs_L2_block);
// new: temp storage for int32 accums
allocate(some_int32_accumulators_block);
for (some_lhs_L2_block) {
  pack(some_lhs_L2_block);
  for (some_rhs_L2_block) {
    pack(some_rhs_L2_block);
    for (some_lhs_sub_block in some_lhs_L2_block) {
      for (some_rhs_sub_block in some_rhs_L2_block) {
        // new: pass int32 accums to kernel
        kernel(&some_int32_accumulators_block,
               some_lhs_sub_block,
               some_rhs_sub_block);
      }
    }
    // new: unpack int32 accums into destination matrix
    unpack(some_int32_accumulators_block);
  }
}