Matrix multiplication; the holy grail of deep neural networks and fashionable language understanding behemoths. As MLEs or information scientists, our fingers are too fast to sort `tf.matmul`

or `torch.matmul`

and we by no means look again. However don’t inform me you’ve by no means had the millisecond infatuation to know what could be taking place to that matrix when it enters the GPU! For those who did, you’re in the proper place. Be a part of me in a journey by the fascinating intricacies inside a GPU.

I’ll clarify to you ways these compute powerhouses crunch up the numbers. You’ll be taught three little-known spectacular issues GPUs do, after they come face-to-face with matrices. By the tip of this weblog publish, you’ll have understanding of how matrix multiplication works inside GPUs.

GEMM or generalized matrix multiplication is the kernel that’s executed when GPUs carry out matrix multiplication.

`C = a (A.B) + b C`

Right here, `a`

and `b`

are scalars, `A`

is an `MxK`

matrix, `B`

is an `KxN`

matrix, and thus `C`

is an `MxN`

matrix. It’s simple as that! You would possibly surprise why that trailing addition exists. Seems this can be a fairly frequent sample for neural networks (e.g. including bias, making use of ReLU, including residual connections).

For those who’re requested to write down a matrix multiplication algorithm from first ideas, right here’s what you’ll do (except you’re gifted with a GPU in lieu of a mind — wouldn’t that get monetary savings for an MLE!).

`for (int i = 0; i < M; ++i)`

for (int j = 0; j < N; ++j)

for (int okay = 0; okay < Ok; ++okay)

C[i][j] += A[i][k] * B[k][j];

Right here’s an animated visible that reveals you what this does.

However do you know GPUs despise this implementation 🤔? To grasp why that’s the case, it’s good to perceive the GPU reminiscence structure,

For all comparisons and specs, I’ll be utilizing the Nvidia A100 GPU specs.

A GPU has three essential reminiscence ranges,

- World reminiscence or HBM (what you usually discuss with as GPU reminiscence and what you see once you run
`nvidia-smi`

) - Shared reminiscence (an area reminiscence that’s devoted to a single streaming multiprocessor [or SM] and shared between threads operating in that SM)
- Registers (individually allotted to threads to hold out their workload)

That is what it appears like,

The very first thing to notice is that shared reminiscence (known as SRAM to any extent further) is method smaller than the HBM, not to mention registers. So your matrix shouldn’t be going to slot in there (in most events). If we return to our animation, for a single row of `A`

all columns of `B`

must be retrieved, and repeat the method for all rows in `A`

. This implies, the GPU must do many-many reads to compute the output. The HBM (~1.5TB/s) is a number of magnitudes slower than SRAM (~19TB/s).

To place that in numbers, say you need to multiply a `10x20`

and `20x30`

matrix, it’s good to learn columns of `B`

`10x30=300`

occasions. Is there a greater method we are able to do that?

Seems a easy trick can go a good distance right here! Merely flip the order of the loops, in order that `okay`

turns into the outer most loop. And also you’re executed! 😮

`for (int okay = 0; okay < Ok; ++okay) `

for (int i = 0; i < M; ++i)

for (int j = 0; j < N; ++j)

C[i][j] += A[i][k] * B[k][j];

We didn’t contact the precise computation, simply the order of the loops, so we must always get the identical outcome as earlier than. Right here’s what the matrix multiplication appears like now!

You see, we solely carry one *column* of `A`

and one *row* of `B`

at a time and by no means look again. This requires far much less reads than the unique implementation. The one distinction is we have been computing the *interior product* between two vectors earlier than, now we’re computing the *outer product*.

However nonetheless, we want total `C`

in SRAM, which could be too massive to slot in SRAM. What does CUDA do then? That brings us to the second trick.

To not fear! I’m not going to blast you with any complicated arithmetic or Leetcode algorithms. The primary factor to bear in mind is, a matrix is a 2D format of particular person tiles. The next animation does justice to what I’m attempting to clarify.

The results of the inexperienced block 💚 is the sunshine blue strip of A 💙 and the sunshine yellow strip of B 💛. Taking this a step additional, to compute the output, you possibly can carry one block of that strip of A and one block from B’s strip at a time, compute the output and accumulate the outcome within the inexperienced field.

This offers us a versatile framework the place we are able to load an arbitrary dimension block (or tile) of A and B and nonetheless compute the ultimate reply. We don’t must cease there, we are able to hold recursively dividing the issue to even smaller issues. i.e. the matrix is damaged into tiles, tiles are damaged into fragments, and fragments to particular person values.

And this lends itself properly to the method execution structure in a GPU. There are three layers to a kernel execution in a GPU. For simplicity, we’ll say a SM runs a single thread block (though in apply they execute them concurrently, to cut back one thing often known as the tail effect).

- Threads
- Warps (a set of 32 threads)
- Thread blocks (a set of a number of warps)

The precise variety of threads in a thread block depends upon a selected structure. For instance, an A100 has the following specifications.

- Most of 2048 threads per SM
- Most of 1024 threads per block
- Most of 32 thread blocks per SM

## Sidebar #2: Magic of the ability of two

Going again to the tiling, It has been discovered that (heuristically) a matrix tile of dimension `256x128`

per thread block provides affordable effectivity for many issues. Subsequently it’s a standard tile dimension utilized by CUDA.

You might need heard a few finest apply of holding batch dimension, hidden dimension dimension as powers of two. That is the place this comes from! When your matrix dimensions are of powers of two, it will likely be absolutely divisible to a set of tiles with no the rest. If not, it makes your code less efficient.

GPU computations are extra environment friendly when your matrix dimensions are within the energy of two

What occurs when it’s not an influence of two?

## Sidebar #2: Tile quantization

What occurs is an impact often known as *tile quantization*. In different phrases, when you’ve got a tile row dimension of 128 however your matrix has 257 parts in a row, you’ll needn’t two, however three tiles in a row (i.e. 256+1). That is illustrated under.

Drawback with that is that, the thread block does the identical quantity of computation whatever the helpful information residing in it. So, you’re taking the chance to do helpful computation away out of your GPU, resulting in inefficiencies.

An identical impact is named wave quantization, the place the matrix is over-sized and the SMs collectively can not match it without delay. Then the GPU must do the computation in 2 “waves”. Nonetheless that is much less of a priority for contemporary GPUs as they leverage concurrency to cut back wave quantization.

Tile quantization occurs when a thread block has to spill information partially, wave quantization occurs when SMs must spill information.

The ultimate trick is kernel fusion. Most of the time, it’s quicker to do all of the computations in a single kernel than having two kernels known as one after the opposite. Why? As a result of one kernel wants to write down the info to HBM and different must learn that again. We already talked about how gradual that is. A greater method is simply mix the 2 operations into one. Some examples are,

In order it’s seen right here (I’m positive Pytorch has the same glossary), there are a lot of fused kernels provided by TensorFlow that mixes commodious operations in to a single kernel. In code, it means one thing like this,

`for (int m = 0; m < M; m += Mtile) `

for (int n = 0; n < N; n += Ntile)

for (int okay = 0; okay < Ok; ++okay)

float tmp = 0

for (int i = 0; i < Mtile; ++i)

for (int j = 0; j < Ntile; ++j)

int row = m + i;

int col = n + j;

tmp += A[row][k] * B[k][col];

// Do different issues

C[row][col] = tmp

In different phrases, we maintain dearly to our `tmp`

variable till after we’ve completed all our computations. Then solely we’ll write the outcome again to `C`

.

That’s it of us. I hope this was an pleasurable tour by the weeds of a GPU. For those who’re within the audio-visual model right here’s the hyperlink to my YouTube video.