Why is naive matrix multiplication so slow?

The naive triple loop (order i, j, k) is memory-bound, not compute-bound. Its inner loop reads B down a column — each step jumps a full row ahead in memory (a stride of n elements), so almost every access is a fresh cache line and a likely cache miss. Measured on Apple Silicon (macOS) it runs at about 2.64 GFLOP/s at n=1024, while the hardware is capable of hundreds. The multiplies are cheap; waiting on memory is the cost.

Does loop tiling (blocking) always make matrix multiplication faster?

No — it is situational. In our benchmark the single biggest win came from simply reordering the loops (i, k, j) so the inner loop strides contiguously: that alone was 6.3× faster than naive (2.64 → 16.73 GFLOP/s at n=1024), with no blocking at all. Adding cache tiling on top did not beat the reordered loop on this hardware — best tiled result was 14.83 GFLOP/s, about 11% slower than the plain reordered loop. Tiling pays off when the working set genuinely exceeds cache and when combined with register-level blocking; on a big-cache CPU with an auto-vectorizing compiler, a clean reordered loop is already hard to beat.

What is the fastest way to multiply matrices in practice?

Call a tuned BLAS library — numpy's "@" operator, which dispatches to Accelerate, MKL or OpenBLAS. On the same machine that ran at 2.64 GFLOP/s naive, numpy hit about 410 GFLOP/s at n=1024 — roughly 155× the naive loop and 24.5× our best hand-written loop. That gap is not cache tiling; it is SIMD vector instructions, register blocking, multithreading and hand-tuned microkernels. Do not hand-roll matmul for production — reach for the library.

What is the best tile size for blocked matrix multiplication?

It depends on your cache sizes, and you tune it empirically. Sweeping tile sizes at n=1024 on our machine, larger tiles won: performance climbed from about 8.8 GFLOP/s at a 16-wide tile to 14.83 at 256, because a bigger block reuses more data before it is evicted — as long as three tiles still fit in L2 (3 × T² × 8 bytes). Past that point the block overflows cache and speed drops again. There is no universal best value; it is a per-machine knob.

Tiled vs Naive Matrix Multiplication: What Actually Makes It Fast

Tiling — splitting the matrices into cache-sized blocks — is the famous trick for making matrix multiplication fast, and it's the first thing most guides reach for. So we wrote all three versions in C, compiled them the same way (cc -O2), and timed them against the matrix multiply you'd actually use in practice: a BLAS library. The numbers say something the tutorials usually skip — the biggest free win isn't tiling at all. Everything below was measured on Apple Silicon (macOS); the reproduce-this box has the exact loops.

Four ways to multiply the same two matrices

Same 1024×1024 matrices, same machine, four implementations. The scale is logarithmic because the spread is enormous: from the naive loop to the BLAS call is a factor of 155.

n=1024, Apple Silicon (macOS), cc -O2, BLAS via numpy (Accelerate). All four produced the identical result (checksums matched); only the speed differs.

Key finding: the largest free win is loop order, not tiling. Reordering the naive loops to stride contiguously (i, k, j) was 6.3× faster (2.64 → 16.73 GFLOP/s) with no blocking. Our best tiled version (14.83 GFLOP/s) never beat that reordered loop on this hardware. And the real 24.5× that BLAS has over our best loop is SIMD, register blocking and threads — cache tiling is just one rung of that ladder.

The free lunch: loop order

The naive order is for i, for j, for k: C[i][j] += A[i][k]·B[k][j]. The problem is the inner loop over k: A[i][k] walks along a row (contiguous, good), but B[k][j] walks down a column — every step jumps a whole row ahead in memory. At n=1024 that's a 8 KB stride per multiply, so almost every access misses cache. The CPU spends its time waiting on memory, which is why naive sits at just 2.64 GFLOP/s.

Swap the last two loops to for i, for k, for j and the inner loop becomes C[i][j] += a·B[k][j] with a = A[i][k] hoisted out. Now C[i][j] and B[k][j] both march along contiguous rows — unit stride. The hardware prefetcher and the compiler's auto-vectorizer both love that pattern, and the same arithmetic runs 6.3× faster. Not one line of blocking code; just touching memory in the order it's laid out.

Where tiling helps — and where it doesn't

Tiling computes the answer one small block at a time so a T×T patch of each matrix stays hot in cache while it's reused. That's a real effect, and you can watch it work: sweeping the tile size, throughput climbs steadily as the block grows — bigger tiles reuse more before eviction, right up until three tiles (3·T²·8 bytes) stop fitting in L2. But here's the honest part — on this machine the whole tiled family stayed below the plain reordered loop.

Tiling clearly gets better with bigger blocks — but every bar sits under the dashed line, the same loops reordered without any tiling at all.

Why the famous trick underperformed here: this CPU has a large L2, and the compiler already auto-vectorizes the clean reordered inner loop into wide SIMD stores. A hand-written blocked loop adds index arithmetic and loop overhead the reordered version doesn't pay, and the cache pressure it relieves wasn't the bottleneck yet. Tiling earns its keep on smaller-cache CPUs, at much larger matrices, or — crucially — as one layer of a multi-level blocking scheme (registers → L1 → L2), which is exactly what real libraries do. As a bolt-on to a loop that already streams cache-friendly, it's a wash. The lesson isn't "tiling is useless" — it's "measure on your hardware; the famous optimization isn't automatically the winning one."

The 24.5× that's left over

Our best hand loop reached 16.73 GFLOP/s. numpy reached about 410 — still 24.5× more, from the same silicon. None of that remaining gap is cache tiling. It's SIMD (one instruction multiplying a whole vector of numbers at once), register blocking (keeping a tiny sub-block of the result in registers so it never touches memory mid-accumulation), multithreading across cores, and microkernels hand-tuned per CPU. Cache tiling is one rung on that ladder, not the ladder. The practical takeaway is the oldest one in performance work: for real matmul, call the library — the matrix multiply tool here uses exactly that path — and spend your own effort only where no library exists.

Reproduce this

The three loops are below — compile with cc -O2, time each, and divide 2·n³ by the seconds for GFLOP/s. Unlike our deterministic studies, these are wall-clock timings: your absolute numbers will differ by machine, but the shape — naive slow, reorder ~6.3×, tiling situational, BLAS far ahead — reproduces.

// naive — inner loop reads B down a column (stride n): cache-hostile
void ijk(const double*A,const double*B,double*C,int n){
  for(int i=0;i<n;i++)for(int j=0;j<n;j++){
    double s=0;
    for(int k=0;k<n;k++) s+=A[i*n+k]*B[k*n+j];
    C[i*n+j]=s;
  }
}

// reordered — inner loop over j is unit-stride: the 6.3x win
void ikj(const double*A,const double*B,double*C,int n){
  for(int i=0;i<n*n;i++)C[i]=0;
  for(int i=0;i<n;i++)for(int k=0;k<n;k++){
    double a=A[i*n+k];
    for(int j=0;j<n;j++) C[i*n+j]+=a*B[k*n+j];
  }
}

// tiled — same math, blocked into T-sized patches for cache reuse
void tiled(const double*A,const double*B,double*C,int n,int T){
  for(int i=0;i<n*n;i++)C[i]=0;
  for(int ii=0;ii<n;ii+=T)for(int kk=0;kk<n;kk+=T)for(int jj=0;jj<n;jj+=T)
    for(int i=ii;i<ii+T&&i<n;i++)for(int k=kk;k<kk+T&&k<n;k++){
      double a=A[i*n+k];
      for(int j=jj;j<jj+T&&j<n;j++) C[i*n+j]+=a*B[k*n+j];
    }
}

The BLAS baseline is one line of Python: C = A @ B with numpy, timed the same way. Chart geometry and derived ratios are baked by scripts/gen-matmul-study.mjs from the measured rates.

Do the matrix math

Matrix Multiply — the product of two matrices, computed in your browser via the fast path.
Matrix Inverse and Solve Ax = b — Gaussian elimination on your own numbers.
Transpose and Rank — the rest of the math tools.

Share this X Facebook Reddit