Float16 vs Float32: What GEMM and My Benchmarks Taught Me About GPU Performance

When I first started benchmarking matrix multiplication on my system, I thought I already understood what was going to happen. Float16 uses fewer bits than float32, so it should be faster. That felt obvious. Smaller numbers mean less work, right? But the moment I actually ran the benchmarks across CPU and Apple GPU (MPS), that simple intuition completely broke down, and what I saw forced me to rethink how computation really works at the hardware level.

For a 1024 × 1024 matrix multiplication, my results looked like this: CPU with float32 achieved around 1.66 TFLOPS, which felt reasonable, but CPU with float16 dropped to almost zero effective throughput. On the GPU using MPS, float32 reached around 1.75 TFLOPS, while float16 pushed slightly higher to around 1.94 TFLOPS. The most shocking part was not that GPU float16 was faster, but that CPU float16 was catastrophically slow. At that point, it became clear that this wasn’t just a story about precision. It was a story about how computation interacts with hardware.

It started taking too much time in CPU Fp16 for 2048x2048 matrix

To understand this, I went back to something very basic. If we multiply two tensors using float32 and float16, the operation itself does not change at all.

import torch

a_fp32 = torch.tensor([1.5, 2.0], dtype=torch.float32)
b_fp32 = torch.tensor([3.0, 4.0], dtype=torch.float32)

a_fp16 = torch.tensor([1.5, 2.0], dtype=torch.float16)
b_fp16 = torch.tensor([3.0, 4.0], dtype=torch.float16)

print(a_fp32 * b_fp32)
print(a_fp16 * b_fp16)

In both cases, the same multiplication is happening. The algorithm is identical. What changes is how the numbers are represented in memory. Float32 stores each value using 4 bytes, while float16 uses only 2 bytes. At first glance, that looks like a simple storage optimization, but once you look at matrix multiplication more closely, it becomes something much more fundamental.

Matrix multiplication, or GEMM, is usually written as something very simple like C = A @ B, but internally each element of the result is computed through repeated multiplication and accumulation. If you expand it, each element is formed by taking a row from A and a column from B, multiplying corresponding values, and summing them together. Writing a naive implementation makes this very clear.

def matmul_naive(A, B):
    m = len(A)
    k = len(A[0])
    n = len(B[0])

    C = [[0.0 for _ in range(n)] for _ in range(m)]

    for i in range(m):
        for j in range(n):
            acc = 0.0
            for t in range(k):
                acc += A[i][t] * B[t][j]
            C[i][j] = acc

    return C

Looking at this loop carefully, you notice something important. Every iteration is not just doing a multiplication, it is also loading values from memory. That second part turns out to be the real bottleneck in many cases. Modern hardware is extremely fast at arithmetic, often much faster than it is at moving data. This means that performance is frequently limited not by how fast we can multiply numbers, but by how quickly we can feed those numbers into the computation.

This is exactly where float16 starts to matter. If each value takes half the space, then the amount of data that needs to be moved through the system is also cut roughly in half. You can see this directly by calculating how many bytes GEMM actually moves.

def gemm_bytes(m, n, k, dtype):
    element_size = torch.tensor([], dtype=dtype).element_size()
    return (m * k + k * n + m * n) * element_size

For the same matrix size, float16 reduces the total data movement significantly, while the number of floating point operations remains exactly the same. That means the computation is doing the same amount of work, but with less data flowing through memory. This changes the entire balance of the system.

On the Apple GPU (MPS), this shows up clearly in the benchmarks. When I ran GEMM using float16, the throughput increased slightly compared to float32. The math did not change, but the GPU was able to operate more efficiently because it spent less time waiting on memory. The workload effectively shifted from being partially limited by memory bandwidth toward being more limited by compute, which is exactly where GPUs are strongest.

This idea becomes easier to understand when you think in terms of how much computation you do per byte of data. In GEMM, the number of operations is fixed, but the amount of data depends on the datatype. When you switch from float32 to float16, you reduce the data size while keeping the computation constant. That means you are doing more work per byte of data. In other words, you are using the hardware more efficiently.

The CPU result tells the opposite story. When float16 performance collapses to almost zero, it reveals that the CPU is not executing float16 matrix multiplication efficiently at all. Running a simple benchmark confirms this behavior.

import time
import torch

def benchmark_mm(device, dtype):
    a = torch.randn(1024, 1024, device=device, dtype=dtype)
    b = torch.randn(1024, 1024, device=device, dtype=dtype)

    start = time.perf_counter()
    for _ in range(20):
        torch.mm(a, b)
    end = time.perf_counter()

    return (end - start) / 20

On CPU, float32 is well-optimized and runs through efficient vectorized libraries. Float16, however, often does not have a direct optimized path. In many cases, it is internally converted or handled in a slower way. So instead of benefiting from reduced data movement, the computation itself becomes inefficient, leading to extremely poor performance.

There is another subtle detail that initially confused me. Even when we say we are using float16, the entire computation is not always done in float16. Matrix multiplication involves accumulation, and accumulation is sensitive to precision. If everything were done purely in float16, small rounding errors would quickly build up.

To see this, consider summing many small values:

import torch

values = torch.tensor([0.1] * 1000, dtype=torch.float16)
print(values.sum())

values_fp32 = values.to(torch.float32)
print(values_fp32.sum())

The results differ because float16 cannot represent small increments accurately enough during accumulation. To solve this, modern systems use mixed precision. Multiplications are performed in float16, but accumulation is done in float32. This allows the system to maintain numerical stability while still benefiting from the performance advantages of float16.

Another interesting observation from my benchmarks was that the measured memory bandwidth in GB/s was actually lower for float16 than for float32. At first, this looked like float16 was underperforming, but the opposite is true. Since float16 reduces the total number of bytes being moved, the calculated bandwidth can decrease even while overall performance improves. The system is simply doing the same work with less data.

When I increased the matrix sizes further in my prefill-style benchmarks, performance remained high initially but started to drop slightly at very large sizes. This reflects increasing memory pressure. As the matrices grow, cache reuse becomes less effective, and more data must be fetched from main memory. Even a compute-heavy operation like GEMM eventually becomes influenced by memory bandwidth again.

What all of this made clear to me is that float16 is not inherently faster on its own. It is faster because it changes how the workload interacts with the hardware. It reduces data movement, increases the amount of useful computation per byte, and aligns better with how GPUs are designed to operate.

In the context of large language models, this insight becomes even more important. Almost every major component of a transformer can be expressed as a matrix multiplication. Each layer, each attention head, each projection involves GEMM operations. If each of those operations becomes more efficient, the entire model becomes faster.

x = torch.randn(1, 4096, device="mps", dtype=torch.float16)
W = torch.randn(4096, 4096, device="mps", dtype=torch.float16)

y = x @ W

This is repeated many times across the network. The gains from float16 accumulate across layers, leading to significant improvements in overall inference performance.

Before running these benchmarks, I thought of float16 versus float32 mainly in terms of precision. After working through GEMM and observing real performance, I now see it as a question of alignment. The faster format is not just the smaller one. It is the one that allows the hardware to operate efficiently.

Check out my code:

https://github.com/utsab818/LLMInference