Try an interactive version of this dialog: Sign up at solve.it.com, click Upload, and pass this URL.

Note: 217

Matrix Multiplication

Here i shall be going through matrix multiplication from scratch using just python , its standard library, matplotlib. Once I have functions defined in barebones python, I can replace that function with something higher level in numpy or pytorch.

This work is from lesson 11 and 12 of the practical deep learning for coders course from fast.ai

Below are the relevant links

First I shall start by importing the data we are going to use - It is the mnist dataset

Code: 0 ()


Prompt: 18

Give me a summary of the information from each of the lessons

Assistant: 312

Based on the lesson transcripts you've provided, here are summaries of each:

Lesson 11: Deep Learning Foundations This lesson covers reading academic papers (using DiffEdit as an example), implementing matrix multiplication from scratch in Python, and progressively optimizing it. Key topics include:

  • Understanding diffusion models and CLIP
  • Broadcasting rules for tensor operations
  • Using Numba for compilation to machine code
  • Achieving massive speedups (eventually 5 million times faster using GPU)

Lesson 12: Deep Learning Foundations to Stable Diffusion This lesson continues matrix multiplication optimization and introduces clustering algorithms:

  • CLIP Interrogator and why it can't truly invert images
  • Einstein summation notation (einsum)
  • GPU acceleration with CUDA
  • Implementing mean shift clustering from scratch
  • Creating matplotlib animations
  • Introduction to calculus concepts (derivatives, infinitesimals)

Would you like me to elaborate on any specific topic from these lessons, or would you like to start working through the matrix multiplication concepts with your MNIST data?

Note: 9

Load data, Import Libraries

Code: 39 ()

from pathlib import Path
import gzip, pickle, matplotlib as mpl, matplotlib.pyplot as plt
from random import random
import os

Code: 87 ()

MNIST_URL = 'https://github.com/mnielsen/neural-networks-and-deep-learning/blob/master/data/mnist.pkl.gz?raw=true'
path_data = Path('data')
path_data.mkdir(exist_ok=True)
path_gz = path_data/'mnist.pkl.gz'

Code: 39 ()

from urllib.request import urlretrieve
if not path_gz.exists():
    urlretrieve(MNIST_URL, path_gz)

Prompt: 21

Give me a refresher on urlib, urllib.request, urlretrieve

Assistant: 312

Great question! Let's break down these modules:

urllib is Python's standard library for working with URLs. It's built into Python, so no installation needed.

urllib.request is a submodule within urllib that handles opening and reading URLs. Think of it as the part that actually fetches data from the internet.

urlretrieve() is a specific function that downloads a file from a URL and saves it to your local disk. It takes two main arguments:

  • The URL (where to download from)
  • The local file path (where to save it)

In your code, urlretrieve(MNIST_URL, path_gz) is downloading the MNIST dataset and saving it as 'mnist.pkl.gz' in your data folder.

It's a simple, straightforward way to download files - much easier than manually handling HTTP requests and file writing!

Do you want to see what other functions urllib.request offers, or shall we move on to exploring the MNIST data you've just loaded?

Code: 4 ()

urlretrieve?

Output: 1,303

Copied!
def urlretrieve(
    url, filename:NoneType=None, reporthook:NoneType=None, data:NoneType=None
):
Copied!
Retrieve a URL into a temporary location on disk.

Requires a URL argument. If a filename is passed, it is used as
the temporary file location. The reporthook argument should be
a callable that accepts a block number, a read size, and the
total file size of the URL target. The data argument should be
valid URL encoded data.

If a filename is passed and the URL points to a local resource,
the result is a copy from local file to new file.

Returns a tuple containing the path to the newly created
data file as well as the resulting HTTPMessage object.

File: /usr/local/lib/python3.12/urllib/request.py

Type: function

Code: 69 ()

#context manager with destructuring
with gzip.open(path_gz, 'rb') as f:
    ((x_train, y_train),(x_valid, y_valid),_) = pickle.load(f, encoding='latin-1')

Prompt: 13

What is a context manager in the above context

Assistant: 277

A context manager is a Python object that manages resources automatically - it handles setup and cleanup for you.

In your code, with gzip.open(path_gz, 'rb') as f: is using a context manager. Here's what it does:

  1. Opens the gzip file when you enter the with block
  2. Automatically closes the file when you exit the block (even if an error occurs)

Without a context manager, you'd need to write:

Copied!
f = gzip.open(path_gz, 'rb')
# ... do stuff ...
f.close()  # Easy to forget!

The with statement guarantees cleanup happens, which prevents resource leaks (like leaving files open).

The as f part gives you a name to refer to the opened file within the block.

Have you seen context managers used with other resources, like database connections or file operations?

Code: 28 ()

lst1 = list(x_train[0])
vals = lst1[200:210]
vals

Output: 141

[np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0),
 np.float32(0.19140625),
 np.float32(0.9296875),
 np.float32(0.98828125),
 np.float32(0.98828125),
 np.float32(0.98828125),
 np.float32(0.98828125),
 np.float32(0.98828125)]

Code: 39 ()

def chunks(x, sz):
    for i in range (0,len(x),sz):
        yield x[i:i+sz]

Note: 121

yield goes from 0 to the length of x, my list which has a value of 10 in this case but jumps by the size of sz which is 5 in my case, meaning it returns the value from 0 to 5, then the next 0 to 5.

Basically yield creates an iterator, something you can call next on till the list values are done as shown below.

Code: 10 ()

list(chunks(vals,5))

Output: 157

[[np.float32(0.0),
  np.float32(0.0),
  np.float32(0.0),
  np.float32(0.19140625),
  np.float32(0.9296875)],
 [np.float32(0.98828125),
  np.float32(0.98828125),
  np.float32(0.98828125),
  np.float32(0.98828125),
  np.float32(0.98828125)]]

Code: 15 ()

val_iter = chunks(vals,5)
val_iter

Output: 30

<generator object chunks at 0x75a6ece8a3e0>

Note: 73

val_iter is a generator object, which is a special type of function that returns an iterator object. Instead of using return, we use yield as seen above which returns values in chunks by printing values of the specified length and then pausing execution.

Code: 6 ()

next(val_iter)

Output: 66

[np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0),
 np.float32(0.19140625),
 np.float32(0.9296875)]

Code: 6 ()

next(val_iter)

Output: 75

[np.float32(0.98828125),
 np.float32(0.98828125),
 np.float32(0.98828125),
 np.float32(0.98828125),
 np.float32(0.98828125)]

Code: 7 ()

# next(val_iter)

Note: 73

The above returns the below error

Copied!
Output: 75

---------------------------------------------------------------------------
StopIteration                             Traceback (most recent call last)
Cell In[1], line 1
----> 1 next(val_iter)

StopIteration: 

Prompt: 6

Explain the error above

Assistant: 280

The StopIteration error occurs when you call next() on an iterator that has no more values to yield.

Your generator val_iter was created from vals (which has 10 items) split into chunks of 5. So it can only yield 2 chunks:

  • First next(): returned the first 5 values
  • Second next(): returned the last 5 values
  • Third next(): nothing left! → raises StopIteration

This is Python's way of signaling "I'm done, there's nothing more to iterate over."

When you use a generator in a for loop or list(), Python catches this error automatically and stops looping. But when you manually call next(), you see the error directly.

Does this make sense? Would you like to see how to handle this error gracefully, or shall we continue exploring the MNIST data?

Code: 18 ()

val_iter = chunks(vals,5)
list(val_iter)

Output: 157

[[np.float32(0.0),
  np.float32(0.0),
  np.float32(0.0),
  np.float32(0.19140625),
  np.float32(0.9296875)],
 [np.float32(0.98828125),
  np.float32(0.98828125),
  np.float32(0.98828125),
  np.float32(0.98828125),
  np.float32(0.98828125)]]

Code: 6 ()

type(lst1)

Output: 1

list

Code: 3 ()

chunks??

Output: 750

Copied!
def chunks(x, sz):
    for i in range (0,len(x),sz):
        yield x[i:i+sz]

File: /tmp/ipykernel_74/2841953032.py

Code: 33 ()

mpl.rcParams['image.cmap'] = 'gray'
plt.imshow(list(chunks(lst1,28)));

Output: 822

Note: 180

rcParams, runtime configuration parameters are used to customize the styling. This takes precedence over other ways of customizing the styling of our plots. Above, we are setting up our plot to be grayscale.

Putting chunks(lst1,28) inside a list forces the generator object to return all the chunkified results at once, as we shall demonstrate below. To quote Jeremy Howard, if you pass the iterator to a list, it runs through the entire iterator until it is finished and creates a list of the results.

For the demo, I will use a smaller version of lst1 to keep everything viewable.

Code: 30 ()

val_iter_a = lst1[200:210]
(list(chunks(val_iter_a,5)))

Output: 157

[[np.float32(0.0),
  np.float32(0.0),
  np.float32(0.0),
  np.float32(0.19140625),
  np.float32(0.9296875)],
 [np.float32(0.98828125),
  np.float32(0.98828125),
  np.float32(0.98828125),
  np.float32(0.98828125),
  np.float32(0.98828125)]]

Note: 31

We can take any list like our vals and lst1 and create an iterator by passing it to iter.

Code: 7 ()

from itertools import islice

Code: 9 ()

it = iter(vals)
it

Output: 25

<list_iterator at 0x75a6ece218a0>

Code: 4 ()

next(it)

Output: 12

np.float32(0.0)

Code: 4 ()

next(it)

Output: 12

np.float32(0.0)

Code: 4 ()

next(it)

Output: 12

np.float32(0.0)

Code: 3 ()

iter??

Output: 112

Docstring:
iter(iterable) -> iterator
iter(callable, sentinel) -> iterator

Get an iterator from an object.  In the first form, the argument must
supply its own iterator, or be a sequence.
In the second form, the callable is called until it returns the sentinel.
Type:      builtin_function_or_method

Code: 1 ()

vals

Output: 141

[np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0),
 np.float32(0.19140625),
 np.float32(0.9296875),
 np.float32(0.98828125),
 np.float32(0.98828125),
 np.float32(0.98828125),
 np.float32(0.98828125),
 np.float32(0.98828125)]

Code: 13 ()

is_it = islice(vals,5)

Code: 6 ()

next(is_it)

Output: 12

np.float32(0.0)

Code: 6 ()

next(is_it)

Output: 12

np.float32(0.0)

Code: 6 ()

next(is_it)

Output: 12

np.float32(0.0)

Code: 6 ()

next(is_it)

Output: 15

np.float32(0.19140625)

Code: 6 ()

next(is_it)

Output: 15

np.float32(0.9296875)

Code: 4 ()

islice??

Output: 208

Init signature: islice(self, /, *args, **kwargs)
Docstring:     
islice(iterable, stop) --> islice object
islice(iterable, start, stop[, step]) --> islice object

Return an iterator whose next() method returns selected values from an
iterable.  If start is specified, will skip all preceding elements;
otherwise, start defaults to zero.  Step defaults to one.  If
specified as another value, step determines how many values are
skipped between successive calls.  Works like a slice() on a list
but returns an iterator.
Type:           type
Subclasses:     

Note: 49

islice grabs the first n things from an iterable, in this case we set stop to be 5, so we grab the first 5 things.

Code: 19 ()

is_it = islice(vals,5)
list(is_it)

Output: 66

[np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0),
 np.float32(0.19140625),
 np.float32(0.9296875)]

Code: 6 ()

list(is_it)

Output: 1

[]

Code: 6 ()

callable(list)

Output: 1

True

Note: 84

iter in its second form,iter(callable, sentinel, /) has the ability to invoke a callable until the output of the sentinel value is reached. We shall see this below.

We can use this to recreate our chunks function with islice

Code: 10 ()

it = iter(lst1)
it

Output: 25

<list_iterator at 0x75a6ecbd8250>

Note: 82

A lambda function is one that is anonymous, i.e., a function with no name, just like def is used to define a normal function in Python, lambda is used to define an anonymous function. It takes in an argument and a single expression for example lambda arguments : expression

Code: 40 ()

name = 'silver rubanza'
caps_name = lambda func: func.upper()
caps_name(name),caps_name('devontay')

Output: 39

('SILVER RUBANZA', 'DEVONTAY')

Code: 28 ()

sum_lmd = lambda a,b: a+b
sum_lmd(1,2)

Output: 1

3

Code: 31 ()

it = iter(lst1)
img = list(iter(lambda: list(islice(it,28)),[]))

Code: 21 ()

#next(iter(lambda: list(islice(it,28)),[]))

Code: 6 ()

plt.imshow(img)

Output: 857

<matplotlib.image.AxesImage at 0x75a6ece3c7d0>

Code: 54 ()

def plot_w_islice(x):
    return list(iter(lambda: list(islice(it,28)),[]))
it = iter(lst1)
plt.imshow(plot_w_islice(img))

Output: 854

<matplotlib.image.AxesImage at 0x75a6ecbd9640>

Code: 42 ()

def plot_wt_slice(x):
    return list(islice(x,28))

it = iter(lst1)
plot_wt_slice(it)

Output: 336

[np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0),
 np.float32(0.0)]

Note: 7 ()

Matrices and tensors

Code: 9 ()

img[27][27]

Output: 12

np.float32(0.0)

Code: 9 ()

img[20,10]

Output: 70

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 img[20,10]

TypeError: list indices must be integers or slices, not tuple

Code: 4 ()

type(img)

Note: 768

We can create a class that enables us to use the form img[20,10]to return a value at that point for a matrix just like we would in Numpy.

Object-oriented programming is a paradigm that revolves around using objects and relationships between objects to represent data in Python. This involves defining objects and interacting with them.

A python object is an instance of a class, made up of an identity (Id), a type and a value. The id is the address of an object in memory and cannot be changed. These objects can be used to represent real world entities like a button or a car

An object also has a type which determines what operations we can do on that particular object and defines what possible values objects of that type can take on, we can get the type by using type. just like id, type is unchangeable.

Objects whose value can change are called mutable and immutable for the vice versa. An immutable object can have a mutable object inside, meaning its value might change, but the object inside the immutable object itself remains unchanged.

Objects that contain references to other objects are called containers, for example, a list. Python types includes things like None, numbers.Number and, custom classes.

We can use classes to define, initialize, and manipulate objects, serving as templates to create objects. It defines and structures all objects created.

Methods and constructors are used to create and define classes. Methods are functions within a class that are used to perform a specific task. There are 2 types of methods, user-defined methods and special methods commonly known as dunder methods.

dunder methods are implicit functions with a double underscore in between the method name. These are used to add extra functionality to classes. They are implicitly called when needed by a running program for example, when we call print, it automatically calls the __str__method.

User-defined methods are defined by the user to perform a specific task for example our Cows class could have a milk() method.

A constructor is a special method that the program calls upon the object's creation, this helps initialize the class object with attributes for example in our if we have a function for Cows, we can use a constructor like __init to create characteristics for the Cow as we shall see below.__init is only called during the creation of a class with a sole purpose of initializing the class attributes.

Note: 4

Class Illustration

Code: 157 ()

class Cows:
    """
    A class representing different cow breeds

    Parameters:
    breed: The breed of a cow
    age: The cows age
    weight: The weight of a cow
    color: The color of a cow

    Returns:
    The attributes of a particular breed of cow
    """
    def __init__(self,breed,age,weight,color):
        self.breed= breed
        self.age = age
        self.weight = weight
        self.color = color 

Code: 4 ()

Cows??

Code: 55 ()

#create an object / instance of a cow
fresian = Cows("Fresian",4,600,"Distinctive black and white patches in irregular patterns")
fresian

Note: 6 ()

Access attributes

Code: 19 ()

fresian.breed, fresian.age, fresian.weight

Code: 160 ()

class Cows_dp_a:
    """
    A class representing different cow breeds

    Parameters:
    breed: The breed of a cow
    age: The cows age
    weight: The weight of a cow
    color: The color of a cow

    Returns:
    The attributes of a particular breed of cow
    """
    def __init__(self):
        self.breed = "Jersey"
        self.age = 5
        self.weight = 700
        self.color = "White Patched"

Code: 31 ()

Jersey = Cows_dp_a()
Jersey.breed, Jersey.age, Jersey.weight, Jersey.color

Note: 4

Using default parameters

Code: 180 ()

class Cows_dp:
    """
    A class representing different cow breeds

    Parameters:
    breed: The breed of a cow
    age: The cows age
    weight: The weight of a cow
    color: The color of a cow

    Returns:
    The attributes of a particular breed of cow
    """
    def __init__(self,breed = "Heifer", age = 5, weight = 600, color = "Black"):
        self.breed= breed
        self.age = age
        self.weight = weight
        self.color = color 

Code: 16 ()

Heifer = Cows_dp()
Heifer.breed

Code: 63 ()

class Matrix:
    def __init__(self,xs): self.xs = xs
    def __getitem__(self,idxs): return self.xs[idxs[0]][idxs[1]]

Code: 9 ()

img[20][15]

Code: 16 ()

m = Matrix(img)
m[20,15]

Note: 63

__getitem__ allows us to access elements of an object using square brackets. We can use this to access particular elements from a list just like we would do with img[20][15]

Code: 166 ()

class MyList:
    def __init__(self, data):
        self.data = data
    
    def __getitem__(self, index):
        return self.data[index[0]] [index[1]]

# Create an instance of MyList
a = MyList([[1, 7, 3, 4, 5],[1, 9, 3, 4, 5],[1, 6, 3, 4, 5]])

# Accessing items using square brackets
print(a[2,1])

Note: 46 ()

Random number generator

Below we shall implement a pseudo random number generator based on the Wichman-hill algorithm used before Python 2.3.

Code: 102 ()

rnd_state = None 
def seed(a):
    global rnd_state
    a,x = divmod(a, 30268)
    a,y = divmod(a, 30306)
    a,z = divmod(a, 30322)
    rnd_state = int(x)+1, int(y)+1, int(z)+1

Code: 13

seed(457428938475)
rnd_state

Note: 21 ()

Let me break down the above to show what is happening behind the scenes

Code: 25

r,s = divmod(457428938475, 30268)
r,s

Code: 19

t,u = divmod(r, 30306)
t,u

Code: 19

v,w = divmod(t, 30322)
v,w

Code: 33

a,b,c = int(s)+1, int(u)+1, int(w)+1
a,b,c

Code: 4 ()

divmod??

Code: 121

def rand():
    global rnd_state
    x,y,z = rnd_state
    x = (171*x) % 30269
    y = (172*y) % 30307
    z = (170*z) % 30323
    rnd_state = x,y,z 
    return (x/30269 + y/30307 + z/30323)%1.0

Code: 9

rand(),rand(),rand()

Note: 31

Now every time I run rand(), it should give me a different result.

Let me break it down below

Code: 19 ()

x,y,z = rnd_state
x,y,z,a,b,c

Code: 0


Code: 0


Code: 18 ()

import torch
from torch import tensor
import numpy as np

Note: 67

In Numpy and PyTorch, You can set the display options to be used when printing output to the screen using set_printoptions. We can see the function definition below by adding a ? at the end of the method name

Code: 7 ()

torch.set_printoptions?

Code: 7 ()

np.set_printoptions?

Code: 45 ()

torch.set_printoptions(precision=2, linewidth=140, sci_mode=False)
np.set_printoptions(precision=2, linewidth=140)

Code: 19 ()

tens = tensor(img)
tens[20,15]

Code: 10 ()

tens[20][15]

Code: 33 ()

x_train,y_train,x_valid,y_valid = map(tensor,(x_train,y_train,x_valid,y_valid))

Code: 3 ()

map??

Code: 6 ()

x_train.type()

Code: 16

(171*x),(172*y), (170*z)

Code: 64

x = (171*x) % 30269
y = (172*y) % 30307
z = (170*z) % 30323
rnd_state = x,y,z
x,y,z

Code: 27

(x/30269 + y/30307 + z/30323)%1.0

Code: 27

(x/30269 + y/30307 + z/30323)%1.0

Code: 54

if os.fork(): print(f'In parent: {rand()}')
else:
    print(f'In Child: {rand()}')
    os._exit(os.EX_OK)

Code: 63

if os.fork(): print(f'In parent: {torch.rand(1)}')
else:
    print(f'In Child: {torch.rand(1)}')
    os._exit(os.EX_OK)

Code: 67

import numpy as np
if os.fork(): print(f'In parent: {np.random.rand()}')
else:
    print(f'In Child: {np.random.rand()}')
    os._exit(os.EX_OK)

Code: 72

if os.fork(): print(f'In parent: {random()}')
else:
    print(f'In Child: {random()}')
    os._exit(os.EX_OK)

# Python reinitializes the global random state in each fork

Note: 46

We expect our random numbers not to be correlated or on a linear scale, let us plot 50 random number values and have a look at their distribution.`

Code: 19

plt.plot([rand() for _ in range(50)]);

Code: 21

plt.hist([rand() for _ in range(10000)]);

Code: 34

%timeit -n 10 list(chunks([rand() for _ in range(7840)],10))

Code: 21

%timeit -n 10 torch.rand(784,10)

Note: 87

So above we see that using the torch random generator takes 87.5 µs to run as opposed to our from scratch bare-bones python version which takes 3.66 ms to run, making the PyTorch version much faster, so we shall be using the PyTorch version.

Note: 10

Matrices and Tensors continued

Code: 19 ()

weights = torch.randn(784,10)
weights, weights.shape

Code: 16 ()

bias = torch.zeros(10)
bias, bias.shape

Code: 28 ()

m1 = x_valid[:5] # mini-batch
m1, m1.shape

Code: 16 ()

m2 = weights
m2, m2.shape

Code: 67 ()

# The below gives us the values we can use to create a loop for our matrix multiplication later on.

ar, ac = m1.shape
br, bc = m2.shape
(ar,ac),(br,bc)

Code: 16 ()

t1 = torch.zeros(ar,bc)
t1

Note: 36

We define t1 which has the same shape as our expected output, so we create an object to hold our matrix result

Prompt: 12

Is my intuition about t1 above correct

Assistant: 264

Yes, your intuition is correct! You're creating t1 with shape (ar, bc) which is (5, 10) - this matches the expected output shape when you multiply a (5, 784) matrix by a (784, 10) matrix.

This is a key insight in matrix multiplication: when you multiply an (m, n) matrix by an (n, p) matrix, the result is an (m, p) matrix. The inner dimensions must match (your ac equals br, both are 784), and the outer dimensions give you the result shape.

So t1 is indeed your container to accumulate the results of your triple-nested loop multiplication.

Does the logic of how you're filling in t1 in your loop make sense to you, or would you like to walk through what's happening there?

Code: 9 ()

t1, t1.shape

Code: 58 ()

for i in range(ar):
    for j in range(bc):
        for k in range(ac):
            t1[i,j] += m1[i,k] * m2[k,j]

Prompt: 15

Explain to me what is happening in the above loop

Assistant: 373

Looking at this triple nested loop, you're implementing matrix multiplication from scratch!

Let me break down what's happening step by step:

The outer loop (for i in range(ar)): Goes through each row of m1 (your mini-batch of 5 images)

The middle loop (for j in range(bc)): Goes through each column of m2 (your weights matrix - 10 columns for 10 digit classes)

The inner loop (for k in range(ac)): This is where the actual multiplication and summing happens. It goes through all 784 pixels in one image row and multiplies each pixel by the corresponding weight, accumulating the sum.

So for each position [i,j] in your result matrix t1, you're calculating the dot product of row i from m1 with column j from m2.

In simpler terms: For each image (row in m1) and each output class (column in m2), you're computing a weighted sum of all 784 pixel values.

Does this make sense? Can you see why the result shape is [5, 10] - 5 images giving 10 predictions each?

Prompt: 184

Would i be correct to say that we go through the 5 images, get 1 row for 1 image which has 784 flat values, take that row then go through the 10 columns each with each column having 784 values, we then mulitply the 784 values in the row by the 10 columns where by for each row in m1, we multiply it by each column in m2 for 10 columns while getting the sum for each column resulting in 10 numbers. Once done we do the same for the remaining 4 rows in the m1. is my intuition correct

Assistant: 256

Yes, your intuition is absolutely correct! You've described the process perfectly.

To recap what you said:

  • Take 1 row from m1 (784 pixel values for one image)
  • Multiply it element-wise with each of the 10 columns in m2 (each column also has 784 values)
  • Sum up those 784 products for each column, giving you 10 numbers (one per column)
  • Repeat for all 5 rows in m1

The result: 5 rows × 10 numbers = your [5, 10] output matrix.

This is exactly what matrix multiplication does! Each element in the result is a dot product between a row from the first matrix and a column from the second matrix.

Do you see why this operation is useful for neural networks? What do you think those 10 numbers represent for each image?

Note: 13

I will now put the above into a function

Code: 106 ()

def matmul(a,b):
    ar, ac = a.shape
    br, bc = b.shape
    c = torch.zeros(ar, bc)
    for i in range(ar):
        for j in range(bc):
            for k in range(ac):
                c[i,j] += a[i,k] * b[k,j]
    return c

Prompt: 16

How does the function change from my intitial definition above

Assistant: 543

Your function looks good! The main changes from your initial exploration are:

  1. Generalized the inputs: Instead of hardcoded m1 and m2, you now use parameters a and b - making it reusable for any compatible matrices

  2. Extracts shapes from the inputs: You calculate ar, ac, br, bc from a and b directly, rather than relying on pre-existing variables

  3. Creates the result matrix inside: c = torch.zeros(ar, bc) is created within the function instead of using the external t1

  4. Returns the result: The function returns c, so you can use it like result = matmul(m1, m2)

The core logic (the triple nested loop) remains exactly the same - you've just wrapped it in a reusable function!

Have you tested it yet to see if it gives the same result as PyTorch's built-in matrix multiplication?Great question! Let's compare what changed:

Before: You used the specific variables m1, m2, and t1 that you had already defined in your workspace.

After (in your function): You've generalized it to work with any two matrices a and b, and you create a new result matrix c inside the function.

The logic is identical - same triple nested loop doing the same multiplication and accumulation. But now it's reusable! You can call matmul(any_matrix1, any_matrix2) instead of being tied to just m1 and m2.

This is good practice - turning working code into a reusable function. Have you tested it to make sure it gives the same result as your original loop?

Code: 16 ()

%time _=matmul(m1,m2)

Code: 13 ()

%time matmul(m1,m2)

Prompt: 28

What is the difference between the different time modules am using, which one is best to use when

Assistant: 357

Good question! You're using two different timing tools:

%time: Runs your code once and shows how long that single execution took. It gives you user time (CPU time), sys time (system/kernel time), and wall time (actual elapsed time).

%timeit: Runs your code multiple times (7 runs in your case, with 1 loop each) and gives you statistics - the mean and standard deviation. This is more reliable because it accounts for variability.

When to use which:

  • Use %time when you want a quick single measurement, or when your code takes a very long time to run (you don't want to wait for multiple runs)

  • Use %timeit when you want accurate performance measurements for optimization work - it's better at filtering out noise and giving you reliable comparisons

Notice your %timeit result shows 567 ms ± 8.34 ms - that standard deviation tells you how consistent your function's performance is.

Which one do you think would be better for comparing your custom matmul function against PyTorch's built-in version later?

Note: 34 ()

Now that we know what is going on let us optimize it, and find better ways of doing the same thing

Note: 6

Enter Numba

Prompt: 51

I have heard of something called Numba, that can help make python functions faster. What can you tell me about this and how i can apply it to my matrix multiplication

Assistant: 309

Great instinct! Numba is a just-in-time (JIT) compiler that can speed up Python code, especially loops, by compiling it to machine code.

The beauty of Numba is how simple it is to use - you just add a decorator to your function. The most common one is @njit (which stands for "no-python jit").

Here's how it works: The first time you call the function, Numba compiles it (so it's a bit slow). After that, subsequent calls use the compiled version and run much faster.

For your matrix multiplication, you'd typically:

  1. Import from numba: from numba import njit
  2. Add @njit above your function definition
  3. That's it!

The triple nested loop you wrote is exactly the kind of code Numba excels at optimizing.

Would you like to try adding the decorator to your matmul function and see how much faster it gets? Or do you have questions about how Numba works first?

Note: 309

Below is the description from the Numba website

Numba is an open source JIT compiler that translates a subset of Python and NumPy code into fast machine code.

Numba translates Python functions to optimized machine code at runtime using the industry-standard LLVM compiler library. Numba-compiled numerical algorithms in Python can approach the speeds of C or FORTRAN. You don't need to replace the Python interpreter, run a separate compilation step, or even have a C/C++ compiler installed. Just apply one of the Numba decorators to your Python function, and Numba does the rest for example

Copied!
from numba import njit
import random

@njit
def monte_carlo_pi(nsamples):
    acc = 0
    for i in range(nsamples):
        x = random.random()
        y = random.random()
        if (x ** 2 + y ** 2) < 1.0:
            acc += 1
    return 4.0 * acc / nsamples

Prompt: 21

Show me how i can apply it numba to my matmul function

Assistant: 262

Here's how to apply Numba to your function:

Copied!
from numba import njit

@njit
def matmul(a,b):
    ar, ac = a.shape
    br, bc = b.shape
    c = torch.zeros(ar, bc)
    for i in range(ar):
        for j in range(bc):
            for k in range(ac):
                c[i,j] += a[i,k] * b[k,j]
    return c

Just add the @njit decorator right above your function definition!

Note: Numba works best with NumPy arrays rather than PyTorch tensors. You might need to convert your tensors to NumPy arrays first (using .numpy()) and then convert the result back to a tensor if needed.

Would you like to try this out and see what speedup you get?

Code: 123 ()

from numba import njit

@njit
def matmul(a,b):
    ar, ac = a.shape
    br, bc = b.shape
    c = torch.zeros(ar, bc)
    for i in range(ar):
        for j in range(bc):
            for k in range(ac):
                c[i,j] += a[i,k] * b[k,j]
    return c

Code: 10 ()

matmul(m1,m2)

Prompt: 132

I see that in the lesson 11 video found at the link https://www.youtube.com/watch?v=Tf-8F5q8Xww&t=6133s or my repo found https://github.com/rubanzasilva/fastai_course_artifacts/blob/main/mm-foundations.ipynb, Jeremy uses the njit decorator on the dot product as opposed to the whole matmul function why would that be the case

Assistant: 274

Great observation! Look at the error message you just got - it says Cannot determine Numba type of <class 'torch.Tensor'>.

Numba doesn't know how to work with PyTorch tensors directly. But it does work great with NumPy arrays and basic Python types.

By putting @njit on just the innermost computation (like a dot product function), you can:

  1. Keep that function working with NumPy arrays or basic operations
  2. Have your outer matmul function handle any PyTorch tensor conversions before calling the njit-compiled function

This is a common pattern - optimize the hot inner loop with Numba, while the outer function deals with the framework-specific objects.

What do you think you'd need to change in your code to make this work? Do you want to try creating a separate dot product function or converting your tensors to NumPy first?

Note: 57

So Numba only works with Python and Numpy code and doesnt work with tensors, so we shall use numba on our inner most dot product since that is where the calculation is happening.

Code: 63 ()

from numba import njit
@njit
def dot(a,b):
    res = 0.
    for i in range(len(a)):
        res += a[i]*b[i]
    return res

Code: 36 ()

from numpy import array
%time dot(array([1.,2,3]),array([1.,2,3]))

Code: 36 ()

from numpy import array
%time dot(array([1.,2,3]),array([1.,2,3]))

Note: 114

Copied!
def matmul(a,b):
    ar, ac = a.shape
    br, bc = b.shape
    c = torch.zeros(ar, bc)
    for i in range(ar):
        for j in range(bc):
            for k in range(ac):
                c[i,j] += a[i,k] * b[k,j]
    return c

Code: 96 ()

def matmul(a,b):
    ar, ac = a.shape
    br, bc = b.shape
    c = torch.zeros(ar,bc)
    for i in range(ar):
        for j in range(bc):
            c[i,j] = dot(a[i,:], b[:,j])
    return c

Code: 9 ()

m1.shape,m2.shape

Code: 22 ()

m1a,m2a = m1.numpy(), m2.numpy()

Code: 13 ()

matmul(m1a,m2a)

Prompt: 15

Help me explain my new matmul function using dot

Assistant: 324

Looking at your new matmul function, let me break down what changed:

Before: You had a triple nested loop where the innermost loop (for k in range(ac)) did the element-by-element multiplication and accumulation.

Now: You replaced that entire inner loop with a call to dot(a[i,:], b[:,j]).

Here's what's happening:

  • a[i,:] gets the entire i-th row from matrix a
  • b[:,j] gets the entire j-th column from matrix b
  • dot() computes their dot product (multiply corresponding elements and sum them up)

So instead of three nested loops, you now have two loops that call your Numba-optimized dot function. The dot function handles the innermost calculation much faster because Numba compiled it to machine code.

Does this make sense? Can you see how c[i,j] = dot(a[i,:], b[:,j]) is doing the same work as your previous inner k loop?

Prompt: 15

I cant see it so clearly so help me understand

Assistant: 310

Let's break this down step by step.

First, let's look at what one specific element in your result matrix c[i,j] needs to calculate. Pick any position, say c[0,0] (first row, first column).

Old version with triple loop:

Copied!
for k in range(ac):
    c[0,0] += a[0,k] * b[k,0]

This takes:

  • a[0,0] * b[0,0] + a[0,1] * b[1,0] + a[0,2] * b[2,0] + ... and so on for all 784 values

Do you see how it's multiplying elements from row 0 of a with elements from column 0 of b?

Now, can you tell me what a[0,:] gives you? And what about b[:,0]?

Prompt: 37

a[0,:] gives all elements in a particular row, while b[:,0] gives all values in a particular column

Assistant: 339

Exactly right!

So a[0,:] gives you all 784 values from row 0, and b[:,0] gives you all 784 values from column 0.

Now here's the key insight: What does the dot function do with those two arrays?

Copied!
def dot(a,b):
    res = 0.
    for i in range(len(a)):
        res += a[i]*b[i]
    return res

Look at that loop - it's multiplying corresponding elements and summing them up. So dot(a[0,:], b[:,0]) computes:

  • a[0,0] * b[0,0] + a[0,1] * b[1,0] + a[0,2] * b[2,0] + ... for all 784 elements

That's exactly the same as what your inner k loop was doing!

Can you see now how one call to dot(a[i,:], b[:,j]) replaces the entire inner loop?

Prompt: 84

so would i be correct to say dot is multiplying all the values in a particular row with all the values in a particular column adding them up to come up with one value , then doing the same for the same row but with different columns for all the 10 columns in this case

Assistant: 243

Yes, you've got it exactly right!

The dot function takes one row and one column, multiplies their corresponding elements together, and adds them all up to get one single value for position c[i,j].

Then your outer loops repeat this process:

  • For the same row (same i), it does this with each of the 10 different columns (changing j)
  • Then moves to the next row and does it again for all 10 columns

So for your 5 rows and 10 columns, you're calling dot 5 × 10 = 50 times total, and each dot call is doing 784 multiplications and additions.

Does this click now? Can you see why this is functionally the same as your triple nested loop?

Code: 3 ()

t1

Code: 9 ()

from fastcore.test import *

Code: 21 ()

test_close(t1,matmul(m1a,m2a))

Code: 18 ()

%timeit matmul(m1a,m2a)

Code: 16 ()

%time matmul(m1a,m2a)

Note: 172

Frobenius Norm

Here we are going to take the sum over all of the rows and columns of a matrix, take each one square it, add them up and then take the square root. In simple terms, it is the square of the elements summed up, and the sum is then squared.

Here is a link to something describing it in further detail

Let us look at how this might look like in practice below

Code: 37 ()

m = tensor([[1.,2,3],[4,5,6],[7,8,9]])
m,m.shape

Prompt: 45

remind me of the rules of multplying matrices, can i multiply tensors of different shapes, if so when can i and when cant i

Assistant: 1,122

Good question! Here are the key rules for matrix multiplication:

You CAN multiply two matrices when: The number of columns in the first matrix equals the number of rows in the second matrix.

In shape notation: (m, n) × (n, p) works and gives you (m, p)

Examples that work:

  • (5, 784) × (784, 10) → gives (5, 10)
  • (3, 2) × (2, 5) → gives (3, 5)
  • (100, 50) × (50, 1) → gives (100, 1)

Examples that DON'T work:

  • (5, 10) × (5, 10) ✗ (10 ≠ 5)
  • (3, 4) × (2, 3) ✗ (4 ≠ 2)

The key: inner dimensions must match, and the result takes the outer dimensions.

Does this make sense? Can you tell me if (3, 3) × (3, 3) would work, and if so, what shape the result would be?Good question! There's one key rule for matrix multiplication:

The inner dimensions must match.

For two matrices with shapes (a, b) and (c, d):

  • You can multiply them if b == c
  • The result will have shape (a, d)

So the number of columns in the first matrix must equal the number of rows in the second matrix.

Examples:

  • (5, 784) × (784, 10) ✓ Works! Result is (5, 10)
  • (3, 4) × (4, 7) ✓ Works! Result is (3, 7)
  • (5, 10) × (5, 10) ✗ Doesn't work! (10 ≠ 5)
  • (10, 5) × (5, 3) ✓ Works! Result is (10, 3)

The outer dimensions (first and last) become your result shape, while the inner dimensions disappear (they're what you're summing over).

Does this make sense? Can you tell me if (3, 3) × (3, 2) would work, and what shape the result would be?Good question! Here are the rules for matrix multiplication:

You CAN multiply two matrices when: The number of columns in the first matrix equals the number of rows in the second matrix.

In notation: If matrix A has shape (m, n) and matrix B has shape (p, q), you can multiply them only if n == p.

The result will have shape (m, q).

Examples:

  • (5, 784) × (784, 10) ✓ Works! Result is (5, 10)
  • (3, 4) × (4, 7) ✓ Works! Result is (3, 7)
  • (5, 10) × (5, 10) ✗ Doesn't work! (10 ≠ 5)
  • (10, 5) × (5, 10) ✓ Works! Result is (10, 10)

The key insight: The "inner dimensions" must match. The "outer dimensions" give you the result shape.

Does this make sense? Can you tell me if (3, 3) × (3, 2) would work, and what shape the result would be?

Code: 13 ()

sf = (m*m).sum()
sf

Code: 3 ()

sum?

Code: 4 ()

sf.sqrt()

Prompt: 12

what happens underneath when i do m*m

Assistant: 318

When you do m*m, you're doing element-wise multiplication, not matrix multiplication!

It takes each element in m and multiplies it by the corresponding element in the same position:

  • 1*1 = 1, 2*2 = 4, 3*3 = 9
  • 4*4 = 16, 5*5 = 25, 6*6 = 36
  • 7*7 = 49, 8*8 = 64, 9*9 = 81

So m*m gives you a new (3, 3) tensor with those squared values.

Then .sum() adds all those squared values together: 1 + 4 + 9 + 16 + 25 + 36 + 49 + 64 + 81 = 285

This is different from matrix multiplication! For matrix multiplication in PyTorch, you'd use m @ m or torch.matmul(m, m).

Does that distinction make sense?

Prompt: 15

How does this element wise multiplication differ from matrix multiplication

Assistant: 321

Great question! Let me highlight the key differences:

Element-wise multiplication (m*m or m * n):

  • Multiplies corresponding elements in the same positions
  • Both matrices must have the exact same shape
  • Result has the same shape as the inputs
  • Example: (3,3) * (3,3) → (3,3)

Matrix multiplication (m @ m or matmul(m,n)):

  • Takes dot products of rows and columns (like we just learned!)
  • Inner dimensions must match: (a,b) @ (b,c) works
  • Result shape is the outer dimensions: (a,c)
  • Example: (3,3) @ (3,3) → (3,3) but the values are completely different

Try this to see the difference:

Copied!
print(m * m)    # element-wise
print(m @ m)    # matrix multiplication

Can you predict which one will give larger values? Why?

Code: 28 ()

print(m * m)    # element-wise
print(m @ m)    # matrix multiplication

Prompt: 30

what is happening in (m@m), explain it just like you did for the (m*m)

Assistant: 399

Good! Let me walk through m @ m step by step, just like the element-wise version.

For matrix multiplication m @ m, we're multiplying m by itself using the dot product method we learned earlier.

Let's look at one position - say the top-left element [0,0] which equals 30:

We take row 0 of the first m: [1, 2, 3] We take column 0 of the second m: [1, 4, 7] Then we do: (1×1) + (2×4) + (3×7) = 1 + 8 + 21 = 30

Let's try another - position [0,1] which equals 36:

Row 0: [1, 2, 3] Column 1: [2, 5, 8] Calculation: (1×2) + (2×5) + (3×8) = 2 + 10 + 24 = 36

Can you try calculating what position [1,1] should be? (It should equal 81 according to the output). What row and column do you need, and what calculation would you do?

Note: 25

Let us re-write out matmul with frobenius norm to do our dot product

Code: 4 ()

matmul??

Code: 99 ()

def matmul(a,b):
    ar, ac = a.shape
    br, bc = b.shape
    c = torch.zeros(ar, bc)
    for i in range(ar):
        for j in range(bc):
            c[i,j] = (a[i,:]*b[:,j]).sum()
    return c

Code: 18 ()

test_close(t1,matmul(m1,m2))

Code: 15 ()

%timeit matmul(m1,m2)

Code: 13 ()

%time matmul(m1,m2)

Note: 33

What we do above is the same as the dot product so we can just use torch.dot itself

Code: 97 ()

def matmul(a,b):
    ar, ac = a.shape
    br, bc = b.shape
    c = torch.zeros(ar, bc)
    for i in range(ar):
        for j in range(bc):
            c[i,j] = torch.dot(a[i,:],b[:,j])
    return c

Code: 15 ()

%timeit matmul(m1,m2)

Code: 18 ()

test_close(t1,matmul(m1,m2))

Note: 15

Broadcasting

Broadcasting enables you to multiply

Prompt: 10

Tell me what you know about broadcasting

Assistant: 288

Broadcasting is a powerful technique that allows you to perform operations between tensors of different shapes without explicitly copying data.

The key idea: when operating on two arrays/tensors, PyTorch (and NumPy) automatically "stretches" the smaller one to match the larger one's shape, following specific rules.

Basic broadcasting rules:

  1. If tensors have different numbers of dimensions, pad the smaller one with 1s on the left
  2. Dimensions are compatible if they're equal OR one of them is 1
  3. The dimension of size 1 gets "stretched" to match the other

Simple example:

Copied!
a = tensor([1, 2, 3])      # shape (3,)
b = tensor([[10], [20]])    # shape (2, 1)

What shape do you think a + b would have? And what would the result look like?

Code: 58 ()

a = tensor([1, 2, 3])      # shape (3,)
b = tensor([[10], [20]])    # shape (2, 1)
a+b 

Output: 133

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 a = tensor([1, 2, 3])      # shape (3,)
      2 b = tensor([[10], [20]])    # shape (2, 1)
      3 a+b 

NameError: name 'tensor' is not defined

Prompt: 4

what happened above

Assistant: 375

Great! Let's break down what broadcasting did here:

Starting shapes:

  • a has shape (3,) → treated as (1, 3) (padded with 1 on the left)
  • b has shape (2, 1)

Broadcasting stretches them to a common shape (2, 3):

  • a gets repeated down: [[1, 2, 3], [1, 2, 3]]
  • b gets repeated across: [[10, 10, 10], [20, 20, 20]]

Then element-wise addition happens:

  • Row 0: [1, 2, 3] + [10, 10, 10] = [11, 12, 13]
  • Row 1: [1, 2, 3] + [20, 20, 20] = [21, 22, 23]

The key insight: PyTorch didn't actually copy the data - it just pretended the arrays were stretched during the calculation, making it very efficient!

Can you think of how broadcasting might help us optimize matrix multiplication?

Code: 21 ()

c = tensor([10,20,30])
c,c.type()

Output: 91

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 c = tensor([10,20,30])
      2 c,c.type()

NameError: name 'tensor' is not defined

Code: 22 ()

c = tensor([10.,20.,30])
c, c.type()

Output: 93

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 c = tensor([10.,20.,30])
      2 c, c.type()

NameError: name 'tensor' is not defined

Code: 6 ()

m,m.type()

Output: 70

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 m,m.type()

NameError: name 'm' is not defined

Code: 7 ()

m.shape, c.shape

Output: 73

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 m.shape, c.shape

NameError: name 'm' is not defined

Code: 3 ()

m+c

Output: 69

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 m+c

NameError: name 'm' is not defined

Code: 3 ()

c+m

Output: 69

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 c+m

NameError: name 'c' is not defined

Note: 144

Behind the scenes, torch is using the method expand_as to expand the smaller ranked tensor to match the bigger tensor. The best about this is it doesnt actually save the expanded tensor in memory, it just expands when doing computation so nothing is ever saved in memory.

For example, above I was able to add c to m despite them having different shapes. What is actually happening is that c gets expanded to match the shape of m as we shall see below

Code: 12 ()

t = c.expand_as(m)
t

Output: 82

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 t = c.expand_as(m)
      2 t

NameError: name 'c' is not defined

Code: 6 ()

torch.Tensor.expand?

Output: 13

Object `torch.Tensor.expand` not found.

Note: 73

So now you can see that c remains the same but it is expanded to match the shape of m. With that we can easily add the 2 matrices up because they match shapes and we can do elementwise operations on things with matching shapes.

Code: 4 ()

t.storage()

Output: 69

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 t.storage()

NameError: name 't' is not defined

Code: 7 ()

t.untyped_storage()

Output: 72

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 t.untyped_storage()

NameError: name 't' is not defined

Prompt: 81

what is going on above with the typed storage warning, i thought it would be fixed by using untyoed _storage but that doesnt return the answer am looking for which is something like 10.0 20.0 30.0

Assistant: 333

The warning is just telling you that the .storage() method is deprecated (will be removed in future PyTorch versions), but it still works for now.

The key point you're discovering is correct: even though t looks like a (3, 3) tensor, its underlying storage only contains 3 values: [10.0, 20.0, 30.0].

.untyped_storage() shows you the raw bytes in memory (that's why you see those weird numbers - they're the byte representation of the floats), which isn't as human-readable.

The important insight: PyTorch is being clever! It stores only 3 numbers but uses strides (information about how to step through memory) to make it appear as if those 3 values are repeated 3 times. This saves memory and makes operations faster.

You can ignore the warning for now - it's just PyTorch telling you about future API changes. The point is: broadcasting doesn't actually copy data in memory!

Does this memory efficiency make sense? Can you see why this would be useful for large tensors?

Prompt: 15

if storage is deprecated what should i be using now

Assistant: 291

For checking the underlying storage size and seeing that broadcasting doesn't copy data, you have a few options:

1. Check the number of elements:

Copied!
t.numel()  # total elements in the tensor
t.storage().size()  # actual storage size (still works, just deprecated)

2. Check strides to understand memory layout:

Copied!
t.stride()  # shows how to step through memory

3. Use untyped_storage() with proper interpretation: This is the future-proof way, but as you saw, it shows raw bytes which isn't as readable.

For your learning purposes right now, it's fine to keep using .storage() despite the warning - it clearly shows the concept. PyTorch will provide better alternatives before they actually remove it.

The stride information is actually more useful - try t.stride() and c.stride(). What do you notice about the difference?

Prompt: 12

Iam correct about my expand_as method

Assistant: 204

Yes, you're on the right track! Broadcasting does work similarly to expand_as, but there's a subtle difference:

  • expand_as is an explicit method you can call to create a view with a larger shape
  • Broadcasting happens automatically during operations without you calling anything

Both create views (not copies) of the data, so no extra memory is used - that part is absolutely correct!

You can test this yourself. Try:

Copied!
c.expand_as(m)

Do you see how it creates the expanded view? But when you do m + c, broadcasting does this expansion implicitly for you during the operation.

Does that distinction make sense?

Note: 40

We can change the shape of a matrix using unsqueeze(axis) by inserting a dimension of size 1 at a specified position.

Code: 3 ()

c.shape

Output: 69

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 c.shape

NameError: name 'c' is not defined

Code: 4 ()

c.unsqueeze?

Output: 12

Object `c.unsqueeze` not found.

Code: 4 ()

torch.unsqueeze??

Output: 12

Object `torch.unsqueeze` not found.

Code: 7 ()

c.unsqueeze(0)

Output: 72

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 c.unsqueeze(0)

NameError: name 'c' is not defined

Code: 16 ()

c.unsqueeze(1), c.unsqueeze(1).shape

Output: 82

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 c.unsqueeze(1), c.unsqueeze(1).shape

NameError: name 'c' is not defined

Note: 15

Another way to write the same thing would be

Code: 18 ()

torch.unsqueeze(c,1), torch.unsqueeze(c,0)

Output: 84

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 torch.unsqueeze(c,1), torch.unsqueeze(c,0)

NameError: name 'torch' is not defined

Note: 19

We can also do this with the special value [None]

Code: 7 ()

c[None,:]

Output: 72

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 c[None,:]

NameError: name 'c' is not defined

Note: 30

Trailing ,: are optional so we can simplify our code to use c[None]

Code: 13 ()

c[None],c[None].shape

Output: 79

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 c[None],c[None].shape

NameError: name 'c' is not defined

Prompt: 19

why are trailing commas optional, is my statement about trailing,: true

Assistant: 535

Yes, your statement is true! Trailing ,: are optional in PyTorch (and NumPy) indexing.

Here's why: when you index a tensor, Python assumes you want all remaining dimensions if you don't specify them.

So these are equivalent:

  • c[None]
  • c[None,:]

Both mean: "add a dimension at the front, then take all elements from the original dimension."

Similarly:

  • c[:,None] adds a dimension at the end
  • c[:, None,:] would be invalid because there's no third dimension

The trailing colons are just explicit ways of saying "and everything else" - but if you stop early, it's assumed you want everything else anyway.

Does that make sense? Can you predict what shape c[:,None] would give you?Yes, your statement is true! Trailing commas and colons are optional in PyTorch indexing.

Here's why: When you write c[None,:], the : means "take all elements along this dimension." Since c is 1-dimensional, there's only one dimension to index, so c[None,:] and c[None] do exactly the same thing.

Think of it this way:

  • c[None,:] explicitly says "add a dimension, then take all of the existing dimension"
  • c[None] says "add a dimension" and implicitly takes all of what's left

It's similar to how in Python, if you have a 1D list, my_list[:] and my_list refer to the same elements.

Try this: What do you think c[:,None] would give you? What shape would it have?

Code: 9 ()

c[:,None],c.shape

Output: 75

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 c[:,None],c.shape

NameError: name 'c' is not defined

Code: 13 ()

c[...,None],c[...,None].shape

Output: 79

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 c[...,None],c[...,None].shape

NameError: name 'c' is not defined

Prompt: 15

Explain what the trailing ... are doing above

Assistant: 630

The ... (called "ellipsis") is a special indexing notation that means "all the dimensions I haven't specified."

In your example:

  • c[..., None] means "take all existing dimensions (...), then add a new dimension at the end (None)"
  • Since c has shape (3,), the ... represents that one dimension, so you get shape (3, 1)

The ellipsis is most useful with higher-dimensional tensors. For example, if you had a 4D tensor and wrote tensor[..., None], it would keep all 4 dimensions and add a 5th at the end.

For your 1D tensor c:

  • c[..., None] is the same as c[:, None]
  • c[None, ...] would be the same as c[None, :]

The ... just saves you from writing out all the colons when you have many dimensions. Does that make sense?The ... (called "ellipsis") is a special indexing notation that means "all the dimensions I haven't specified."

In your example:

  • c[..., None] means "take all existing dimensions (...), then add a new dimension at the end (None)"
  • Since c has shape (3,), the ... represents that one dimension, and None adds a new dimension after it
  • Result: shape (3, 1)

The ellipsis is most useful when you have tensors with many dimensions and you only want to modify one specific position. For example:

  • tensor[..., None] - add dimension at the end
  • tensor[None, ...] - add dimension at the beginning
  • tensor[:, ..., None] - keep first dim explicit, add at end

For your 1D tensor c, these are equivalent:

  • c[:, None]
  • c[..., None]

But ... becomes more powerful with higher-dimensional tensors. Does this make sense?

Prompt: 48

Show me an actual example using something similar to my tensor above showing the use of the 4d where using ... would be more convienient

Assistant: 342

Good idea! Let me create a 4D tensor example similar to your work - imagine a batch of color images:

Copied!
images = torch.randn(8, 3, 28, 28)  # 8 images, 3 color channels, 28x28 pixels

Now, say you want to add a dimension at the end:

Without ellipsis:

Copied!
images[:, :, :, :, None]  # shape becomes (8, 3, 28, 28, 1)

With ellipsis:

Copied!
images[..., None]  # same result, much cleaner!

Or if you want to add a dimension at the beginning:

Copied!
images[None, ...]  # shape becomes (1, 8, 3, 28, 28)

Can you see how ... saves you from writing out all those colons? Try creating a similar 4D tensor yourself and experiment with different placements of None and .... What happens with images[:, None, ...]?

Prompt: 16

How comes my example uses [] while you are using ()

Assistant: 241

Great catch! Those are two different things:

Square brackets [] are for indexing/slicing - accessing or modifying parts of an existing tensor:

Copied!
c[None]  # adds a dimension to existing tensor c
images[..., None]  # adds dimension to existing tensor

Parentheses () are for function arguments - passing values when creating something new:

Copied!
torch.randn(8, 3, 28, 28)  # creating a new tensor with this shape

In my example, I used () to create a random tensor with shape (8, 3, 28, 28), then I used [] to add a dimension to it.

Does that distinction make sense?

Code: 1 ()

c

Output: 67

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 c

NameError: name 'c' is not defined

Code: 3 ()

c.shape

Output: 69

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 c.shape

NameError: name 'c' is not defined

Code: 1 ()

m

Output: 67

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 m

NameError: name 'm' is not defined

Code: 13 ()

c[None],c[None].shape

Output: 79

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 c[None],c[None].shape

NameError: name 'c' is not defined

Code: 13 ()

c[:,None],c[:,None].shape

Output: 79

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 c[:,None],c[:,None].shape

NameError: name 'c' is not defined

Code: 7 ()

c.expand_as(m)

Code: 12 ()

c[:,None].expand_as(m)

Output: 76

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 c[:,None].expand_as(m)

NameError: name 'c' is not defined

Code: 12 ()

c[None].expand_as(m)

Output: 76

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 c[None].expand_as(m)

NameError: name 'c' is not defined

Code: 9 ()

m + c[:,None]

Output: 73

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 m + c[:,None]

NameError: name 'm' is not defined

Code: 3 ()

m+c

Output: 69

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 m+c

NameError: name 'm' is not defined

Code: 12 ()

c[None]*c[:,None]

Output: 76

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 c[None]*c[:,None]

NameError: name 'c' is not defined

Note: 133

My intuition on the above

When we do c[None]*c[:,None].

C[None] is being expanded in that C[None] goes from the shape (1,3) to (3,3) by expanding the number of rows while c[:,None] is expanded so that it add the number of columns and goes from shape (3,1) to (3,3)

Prompt: 7

is my intuition above correct

Assistant: 355

Yes, your intuition is exactly right!

When you multiply c[None] * c[:,None]:

  • c[None] shape (1, 3) → broadcasts to (3, 3) by expanding rows (repeating [10, 20, 30] three times down)
  • c[:,None] shape (3, 1) → broadcasts to (3, 3) by expanding columns (repeating each value across)

Then element-wise multiplication happens:

  • Row 0: [10, 20, 30] * [10, 10, 10] = [100, 200, 300]
  • Row 1: [10, 20, 30] * [20, 20, 20] = [200, 400, 600]
  • Row 2: [10, 20, 30] * [30, 30, 30] = [300, 600, 900]

This is called an "outer product" - you're getting every combination of elements from both vectors. Does this pattern make sense to you?

Note: 118 ()

Broadcasting rules

We can do broadcasting between 2 arrays and tensors when they have the same shape, or when one of their dimensions has a value of 1 like in our example of c[None] and c[:,None]

You can always learn more about broadcasting by referring to the numpy docs http://numpy.org/doc/stable/user/basics.broadcasting.html

Prompt: 18

read the numpy broadcasting link above and summarize it for me

Assistant: 7,111

Copied!
{
  "id": "toolu_017amdArEbzVZMJ1An7Yinrw",
  "call": {
    "function": "read_url",
    "arguments": {
      "url": "http://numpy.org/doc/stable/user/basics.broadcasting.html"
    }
  },
  "result": "'Skip to main content\\n\\n__Back to top __ `Ctrl`+`K`\\n\\n[ ![NumPy v2.3 Manual - Home](../_static/numpylogo.svg) ![NumPy v2.3 Manual - Home](../_static/numpylogo_dark.svg) ](../index.html)\\n\\n  * [ User Guide ](index.html)\\n  * [ API reference ](../reference/index.html)\\n  * [ Building from source ](../building/index.html)\\n  * [ Development ](../dev/index.html)\\n  * [ Release notes ](../release.html)\\n  * [ Learn ](https://numpy.org/numpy-tutorials/)\\n  * More \\n    * [ NEPs ](https://numpy.org/neps)\\n\\n__\\n\\n______\\n\\nChoose version\\n\\n  * [__ GitHub](https://github.com/numpy/numpy \"GitHub\")\\n\\n  * [ User Guide ](index.html)\\n  * [ API reference ](../reference/index.html)\\n  * [ Building from source ](../building/index.html)\\n  * [ Development ](../dev/index.html)\\n  * [ Release notes ](../release.html)\\n  * [ Learn ](https://numpy.org/numpy-tutorials/)\\n  * [ NEPs ](https://numpy.org/neps)\\n\\n__\\n\\n______\\n\\nChoose version\\n\\n  * [__ GitHub](https://github.com/numpy/numpy \"GitHub\")\\n\\nSection Navigation\\n\\nGetting started\\n\\n  * [What is NumPy?](whatisnumpy.html)\\n  * [Installation](https://numpy.org/install/)\\n  * [NumPy quickstart](quickstart.html)\\n  * [NumPy: the absolute basics for beginners](absolute_beginners.html)\\n\\nFundamentals and usage\\n\\n  * [NumPy fundamentals](basics.html) __\\n    * [Array creation](basics.creation.html)\\n    * [Indexing on `ndarrays`](basics.indexing.html)\\n    * [I/O with NumPy](basics.io.html)\\n    * [Data types](basics.types.html)\\n    * Broadcasting\\n    * [Copies and views](basics.copies.html)\\n    * [Working with Arrays of Strings And Bytes](basics.strings.html)\\n    * [Structured arrays](basics.rec.html)\\n    * [Universal functions (`ufunc`) basics](basics.ufuncs.html)\\n\\n  * [NumPy for MATLAB users](numpy-for-matlab-users.html)\\n  * [NumPy tutorials](https://numpy.org/numpy-tutorials/)\\n  * [NumPy how-tos](howtos_index.html)\\n\\nAdvanced usage and interoperability\\n\\n  * [Using NumPy C-API](c-info.html)\\n  * [F2PY user guide and reference manual](../f2py/index.html)\\n  * [Under-the-hood documentation for developers](../dev/underthehood.html)\\n  * [Interoperability with NumPy](basics.interoperability.html)\\n\\nExtras\\n\\n  * [Glossary](../glossary.html)\\n  * [Release notes](../release.html)\\n  * [NumPy 2.0 migration guide](../numpy_2_0_migration_guide.html)\\n  * [NumPy license](../license.html)\\n\\n  * [ __](../index.html)\\n  * [NumPy user guide](index.html)\\n  * [NumPy fundamentals](basics.html)\\n  * Broadcasting\\n\\n# Broadcasting#\\n\\nSee also\\n\\n[`numpy.broadcast`](../reference/generated/numpy.broadcast.html#numpy.broadcast \"numpy.broadcast\")\\n\\nThe term broadcasting describes how NumPy treats arrays with different shapes during arithmetic operations. Subject to certain constraints, the smaller array is \u201cbroadcast\u201d across the larger array so that they have compatible shapes. Broadcasting provides a means of vectorizing array operations so that looping occurs in C instead of Python. It does this without making needless copies of data and usually leads to efficient algorithm implementations. There are, however, cases where broadcasting is a bad idea because it leads to inefficient use of memory that slows computation.\\n\\nNumPy operations are usually done on pairs of arrays on an element-by-element basis. In the simplest case, the two arrays must have exactly the same shape, as in the following example:\\n\\n```\\n>>> import numpy as np\\n>>> a = np.array([1.0, 2.0, 3.0])\\n>>> b = np.array([2.0, 2.0, 2.0])\\n>>> a * b\\narray([2.,  4.,  6.])\\n\\n```\\n\\nNumPy\u2019s broadcasting rule relaxes this constraint when the arrays\u2019 shapes meet certain constraints. The simplest broadcasting example occurs when an array and a scalar value are combined in an operation:\\n\\n```\\n>>> import numpy as np\\n>>> a = np.array([1.0, 2.0, 3.0])\\n>>> b = 2.0\\n>>> a * b\\narray([2.,  4.,  6.])\\n\\n```\\n\\nThe result is equivalent to the previous example where `b` was an array. We can think of the scalar `b` being _stretched_ during the arithmetic operation into an array with the same shape as `a`. The new elements in `b`, as shown in Figure 1, are simply copies of the original scalar. The stretching analogy is only conceptual. NumPy is smart enough to use the original scalar value without actually making copies so that broadcasting operations are as memory and computationally efficient as possible.\\n\\n![A scalar is broadcast to match the shape of the 1-d array it is being multiplied to.](../_images/broadcasting_1.png)\\n\\n_Figure 1_#\\n\\n_In the simplest example of broadcasting, the scalar_ `b` _is stretched to become an array of same shape as_ `a` _so the shapes are compatible for element-by-element multiplication._\\n\\nThe code in the second example is more efficient than that in the first because broadcasting moves less memory around during the multiplication (`b` is a scalar rather than an array).\\n\\n## General broadcasting rules#\\n\\nWhen operating on two arrays, NumPy compares their shapes element-wise. It starts with the trailing (i.e. rightmost) dimension and works its way left. Two dimensions are compatible when\\n\\n  1. they are equal, or\\n\\n  2. one of them is 1.\\n\\nIf these conditions are not met, a `ValueError: operands could not be broadcast together` exception is thrown, indicating that the arrays have incompatible shapes.\\n\\nInput arrays do not need to have the same _number_ of dimensions. The resulting array will have the same number of dimensions as the input array with the greatest number of dimensions, where the _size_ of each dimension is the largest size of the corresponding dimension among the input arrays. Note that missing dimensions are assumed to have size one.\\n\\nFor example, if you have a `256x256x3` array of RGB values, and you want to scale each color in the image by a different value, you can multiply the image by a one-dimensional array with 3 values. Lining up the sizes of the trailing axes of these arrays according to the broadcast rules, shows that they are compatible:\\n\\n```\\nImage  (3d array): 256 x 256 x 3\\nScale  (1d array):             3\\nResult (3d array): 256 x 256 x 3\\n\\n```\\n\\nWhen either of the dimensions compared is one, the other is used. In other words, dimensions with size 1 are stretched or \u201ccopied\u201d to match the other.\\n\\nIn the following example, both the `A` and `B` arrays have axes with length one that are expanded to a larger size during the broadcast operation:\\n\\n```\\nA      (4d array):  8 x 1 x 6 x 1\\nB      (3d array):      7 x 1 x 5\\nResult (4d array):  8 x 7 x 6 x 5\\n\\n```\\n\\n## Broadcastable arrays#\\n\\nA set of arrays is called \u201cbroadcastable\u201d to the same shape if the above rules produce a valid result.\\n\\nFor example, if `a.shape` is (5,1), `b.shape` is (1,6), `c.shape` is (6,) and `d.shape` is () so that _d_ is a scalar, then _a_ , _b_ , _c_ , and _d_ are all broadcastable to dimension (5,6); and\\n\\n  * _a_ acts like a (5,6) array where `a[:,0]` is broadcast to the other columns,\\n\\n  * _b_ acts like a (5,6) array where `b[0,:]` is broadcast to the other rows,\\n\\n  * _c_ acts like a (1,6) array and therefore like a (5,6) array where `c[:]` is broadcast to every row, and finally,\\n\\n  * _d_ acts like a (5,6) array where the single value is repeated.\\n\\nHere are some more examples:\\n\\n```\\nA      (2d array):  5 x 4\\nB      (1d array):      1\\nResult (2d array):  5 x 4\\n\\nA      (2d array):  5 x 4\\nB      (1d array):      4\\nResult (2d array):  5 x 4\\n\\nA      (3d array):  15 x 3 x 5\\nB      (3d array):  15 x 1 x 5\\nResult (3d array):  15 x 3 x 5\\n\\nA      (3d array):  15 x 3 x 5\\nB      (2d array):       3 x 5\\nResult (3d array):  15 x 3 x 5\\n\\nA      (3d array):  15 x 3 x 5\\nB      (2d array):       3 x 1\\nResult (3d array):  15 x 3 x 5\\n\\n```\\n\\nHere are examples of shapes that do not broadcast:\\n\\n```\\nA      (1d array):  3\\nB      (1d array):  4 # trailing dimensions do not match\\n\\nA      (2d array):      2 x 1\\nB      (3d array):  8 x 4 x 3 # second from last dimensions mismatched\\n\\n```\\n\\nAn example of broadcasting when a 1-d array is added to a 2-d array:\\n\\n```\\n>>> import numpy as np\\n>>> a = np.array([[ 0.0,  0.0,  0.0],\\n...               [10.0, 10.0, 10.0],\\n...               [20.0, 20.0, 20.0],\\n...               [30.0, 30.0, 30.0]])\\n>>> b = np.array([1.0, 2.0, 3.0])\\n>>> a + b\\narray([[  1.,   2.,   3.],\\n        [11.,  12.,  13.],\\n        [21.,  22.,  23.],\\n        [31.,  32.,  33.]])\\n>>> b = np.array([1.0, 2.0, 3.0, 4.0])\\n>>> a + b\\nTraceback (most recent call last):\\nValueError: operands could not be broadcast together with shapes (4,3) (4,)\\n\\n```\\n\\nAs shown in Figure 2, `b` is added to each row of `a`. In Figure 3, an exception is raised because of the incompatible shapes.\\n\\n![A 1-d array with shape \\\\(3\\\\) is stretched to match the 2-d array of shape \\\\(4, 3\\\\) it is being added to, and the result is a 2-d array of shape \\\\(4, 3\\\\).](../_images/broadcasting_2.png)\\n\\n_Figure 2_#\\n\\n_A one dimensional array added to a two dimensional array results in broadcasting if number of 1-d array elements matches the number of 2-d array columns._\\n\\n![A huge cross over the 2-d array of shape \\\\(4, 3\\\\) and the 1-d array of shape \\\\(4\\\\) shows that they can not be broadcast due to mismatch of shapes and thus produce no result.](../_images/broadcasting_3.png)\\n\\n_Figure 3_#\\n\\n_When the trailing dimensions of the arrays are unequal, broadcasting fails because it is impossible to align the values in the rows of the 1st array with the elements of the 2nd arrays for element-by-element addition._\\n\\nBroadcasting provides a convenient way of taking the outer product (or any other outer operation) of two arrays. The following example shows an outer addition operation of two 1-d arrays:\\n\\n```\\n>>> import numpy as np\\n>>> a = np.array([0.0, 10.0, 20.0, 30.0])\\n>>> b = np.array([1.0, 2.0, 3.0])\\n>>> a[:, np.newaxis] + b\\narray([[ 1.,   2.,   3.],\\n       [11.,  12.,  13.],\\n       [21.,  22.,  23.],\\n       [31.,  32.,  33.]])\\n\\n```\\n\\n![A 2-d array of shape \\\\(4, 1\\\\) and a 1-d array of shape \\\\(3\\\\) are stretched to match their shapes and produce a resultant array of shape \\\\(4, 3\\\\).](../_images/broadcasting_4.png)\\n\\n_Figure 4_#\\n\\n_In some cases, broadcasting stretches both arrays to form an output array larger than either of the initial arrays._\\n\\nHere the `newaxis` index operator inserts a new axis into `a`, making it a two-dimensional `4x1` array. Combining the `4x1` array with `b`, which has shape `(3,)`, yields a `4x3` array.\\n\\n## A practical example: vector quantization#\\n\\nBroadcasting comes up quite often in real world problems. A typical example occurs in the vector quantization (VQ) algorithm used in information theory, classification, and other related areas. The basic operation in VQ finds the closest point in a set of points, called `codes` in VQ jargon, to a given point, called the `observation`. In the very simple, two-dimensional case shown below, the values in `observation` describe the weight and height of an athlete to be classified. The `codes` represent different classes of athletes. [1] Finding the closest point requires calculating the distance between observation and each of the codes. The shortest distance provides the best match. In this example, `codes[0]` is the closest class indicating that the athlete is likely a basketball player.\\n\\n```\\n>>> from numpy import array, argmin, sqrt, sum\\n>>> observation = array([111.0, 188.0])\\n>>> codes = array([[102.0, 203.0],\\n...                [132.0, 193.0],\\n...                [45.0, 155.0],\\n...                [57.0, 173.0]])\\n>>> diff = codes - observation    # the broadcast happens here\\n>>> dist = sqrt(sum(diff**2,axis=-1))\\n>>> argmin(dist)\\n0\\n\\n```\\n\\nIn this example, the `observation` array is stretched to match the shape of the `codes` array:\\n\\n```\\nObservation      (1d array):      2\\nCodes            (2d array):  4 x 2\\nDiff             (2d array):  4 x 2\\n\\n```\\n\\n![A height versus weight graph that shows data of a female gymnast, marathon runner, basketball player, football lineman and the athlete to be classified. Shortest distance is found between the basketball player and the athlete to be classified.](../_images/broadcasting_5.png)\\n\\n_Figure 5_#\\n\\n_The basic operation of vector quantization calculates the distance between an object to be classified, the dark square, and multiple known codes, the gray circles. In this simple case, the codes represent individual classes. More complex cases use multiple codes per class._\\n\\nTypically, a large number of `observations`, perhaps read from a database, are compared to a set of `codes`. Consider this scenario:\\n\\n```\\nObservation      (2d array):      10 x 3\\nCodes            (3d array):   5 x 1 x 3\\nDiff             (3d array):  5 x 10 x 3\\n\\n```\\n\\nThe three-dimensional array, `diff`, is a consequence of broadcasting, not a necessity for the calculation. Large data sets will generate a large intermediate array that is computationally inefficient. Instead, if each observation is calculated individually using a Python loop around the code in the two-dimensional example above, a much smaller array is used.\\n\\nBroadcasting is a powerful tool for writing short and usually intuitive code that does its computations very efficiently in C. However, there are cases when broadcasting uses unnecessarily large amounts of memory for a particular algorithm. In these cases, it is better to write the algorithm\u2019s outer loop in Python. This may also produce more readable code, as algorithms that use broadcasting tend to become more difficult to interpret as the number of dimensions in the broadcast increases.\\n\\nFootnotes\\n\\n[1]\\n\\nIn this example, weight has more impact on the distance calculation than height because of the larger values. In practice, it is important to normalize the height and weight, often by their standard deviation across the data set, so that both have equal influence on the distance calculation.\\n\\n[ __ previous Data types ](basics.types.html \"previous page\") [ next Copies and views __](basics.copies.html \"next page\")\\n\\n__On this page\\n\\n  * General broadcasting rules\\n  * Broadcastable arrays\\n  * A practical example: vector quantization\\n\\n\u00a9 Copyright 2008-2025, NumPy Developers.  \\n\\nCreated using [Sphinx](https://www.sphinx-doc.org/) 7.2.6.  \\n\\nBuilt with the [PyData Sphinx Theme](https://pydata-sphinx-theme.readthedocs.io/en/stable/index.html) 0.16.1.'"
}

Here's a summary of the key points from the NumPy broadcasting documentation:

Broadcasting describes how NumPy treats arrays with different shapes during arithmetic operations, where the smaller array is "broadcast" across the larger array to have compatible shapes.

The General Broadcasting Rules: NumPy compares array shapes element-wise starting from the trailing (rightmost) dimension and working left. Two dimensions are compatible when they are equal, or one of them is 1.

Key insights:

  • Missing dimensions are assumed to have size one
  • Dimensions with size 1 are stretched or "copied" to match the other
  • The stretching is only conceptual - NumPy doesn't actually make copies, making broadcasting memory and computationally efficient

The doc also mentions an important caveat: there are cases where broadcasting uses unnecessarily large amounts of memory, and in those cases it's better to write the algorithm's outer loop in Python.

Does this align with what you've been experiencing in your matrix multiplication work? Any questions about these rules?

Note: 7

Matul with broadcasting

Code: 4 ()

matmul??

Output: 12

Object `matmul` not found.

Code: 90 ()

def matmul(a,b):
    ar, ac = a.shape
    br, bc = b.shape
    c = torch.zeros(ar, bc)
    for i in range(ar):
        c[i] = (a[i,:,None]*b).sum(dim=0)
            
    return c

Code: 10 ()

matmul(m1,m2)

Output: 76

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 matmul(m1,m2)

NameError: name 'm1' is not defined

Code: 18 ()

test_close(t1,matmul(m1,m2))

Output: 85

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 test_close(t1,matmul(m1,m2))

NameError: name 'test_close' is not defined

Code: 13 ()

digit = m1[0]
digit.shape

Output: 85

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 digit = m1[0]
      2 digit.shape

NameError: name 'm1' is not defined

Code: 16 ()

digit_a = x_valid[0]
digit_a.shape

Output: 4

(784,)

Prompt: 12

How many rows and columns does digit have

Assistant: 148

digit has shape (784,), which means it's a 1-dimensional tensor with 784 elements.

So it has:

  • 784 elements in one dimension
  • No explicit rows or columns - it's flat/1D

You can think of it as either 1 row with 784 columns, or 784 rows with 1 column, depending on how you reshape it. But in its current form, it's just a vector with 784 values.

Does that make sense?

Code: 7 ()

m2.shape,m2

Output: 75

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 m2.shape,m2

NameError: name 'm2' is not defined

Code: 10 ()

m2*digit[:,None]

Output: 76

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 m2*digit[:,None]

NameError: name 'm2' is not defined

Code: 16 ()

ca = m2*digit[:,None]
ca.shape

Output: 88

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 ca = m2*digit[:,None]
      2 ca.shape

NameError: name 'm2' is not defined

Code: 10 ()

#m2*digit[None]

Code: 7 ()

digit[:,None].shape

Output: 75

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 digit[:,None].shape

NameError: name 'digit' is not defined

Code: 4 ()

matmul??

Output: 1,258

Copied!
def matmul(a,b):
    ar, ac = a.shape
    br, bc = b.shape
    c = torch.zeros(ar, bc)
    for i in range(ar):
        c[i] = (a[i,:,None]*b).sum(dim=0)

    return c

File: /tmp/ipykernel_74/3450796061.py

Note: 112

Copied!
def matmul(a,b):
    ar, ac = a.shape
    br, bc = b.shape
    c = torch.zeros(ar, bc)
    for i in range(ar):
        for j in range(bc):
            for k in range(ac):
                c[i,j] += a[i,k] * b[k,j]
    return c

Code: 88 ()

def matmul(a,b):
    ar,ac = a.shape
    br,bc = b.shape
    c = torch.zeros(ar,bc)
    for i in range(ar):
        c[i] = (a[i,:,None]*b).sum(dim=0)
    return c

Code: 18 ()

test_close(t1, matmul(m1,m2))

Output: 85

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 test_close(t1, matmul(m1,m2))

NameError: name 'test_close' is not defined

Code: 15 ()

tr = matmul(x_train,m2)
tr

Output: 87

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 tr = matmul(x_train,m2)
      2 tr

NameError: name 'm2' is not defined

Code: 15 ()

%timeit matmul(x_train,m2)

Output: 672

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 get_ipython().run_line_magic('timeit', 'matmul(x_train,m2)')

File /usr/local/lib/python3.12/site-packages/IPython/core/interactiveshell.py:2511, in InteractiveShell.run_line_magic(self, magic_name, line, _stack_depth)
   2509     kwargs['local_ns'] = self.get_local_scope(stack_depth)
   2510 with self.builtin_trap:
-> 2511     result = fn(*args, **kwargs)
   2513 # The code below prevents the output from being displayed
   2514 # when using magics with decorator @output_can_be_silenced
   2515 # when the last Python token in the expression is a ';'.
   2516 if getattr(fn, magic.MAGIC_OUTPUT_CAN_BE_SILENCED, False):

File /usr/local/lib/python3.12/site-packages/IPython/core/magics/execution.py:1222, in ExecutionMagics.timeit(self, line, cell, local_ns)
   1220 for index in range(0, 10):
   1221     number = 10 ** index
-> 1222     time_number = timer.timeit(number)
   1223     if time_number >= 0.2:
   1224         break

File /usr/local/lib/python3.12/site-packages/IPython/core/magics/execution.py:184, in Timer.timeit(self, number)
    182 gc.disable()
    183 try:
--> 184     timing = self.inner(it, self.timer)
    185 finally:
    186     if gcold:

File <magic-timeit>:1, in inner(_it, _timer)

NameError: name 'm2' is not defined

Note: 6

Einstein Summation

Code: 7 ()

torch.einsum?

Prompt: 10

Help me understand einstein summation

Assistant: 415

Einstein summation (einsum) is a compact notation for expressing multi-dimensional array operations. Let me break it down with simple examples.

The basic idea: You label each dimension of your tensors with letters, and any letter that appears in the input but NOT in the output gets summed over.

Simple example - matrix multiplication:

Copied!
torch.einsum('ik,kj->ij', A, B)

Let's decode this:

  • A has dimensions labeled ik (rows=i, cols=k)
  • B has dimensions labeled kj (rows=k, cols=j)
  • Output has dimensions ij
  • Notice k appears in both inputs but NOT in output → we sum over k

This is exactly what matrix multiplication does! For each position [i,j] in the result, you multiply elements where the k indices match and sum them up.

Even simpler - dot product:

Copied!
torch.einsum('i,i->', a, b)  # both inputs have dimension i, output is scalar

Can you see how this matches what you learned about matrix multiplication? The dimension you sum over (k in matmul) disappears from the output!

Would you like to try writing the einsum notation for your matmul operation with m1 and m2?

Code: 34 ()

def matmul(a,b): return torch.einsum('ik,kj->ij', m1,m2)

Code: 10 ()

matmul(m1,m2)

Code: 3 ()

weights.shape

Code: 15 ()

%timeit matmul(x_train,m2)

Prompt: 28

Can you help me create a summary /blog of what i have done above, including my observations

Assistant: 241

I'd be happy to help you create a summary! Let's build it together step by step.

First, let me ask you a few questions to make your summary more personal and meaningful:

  1. What was the main goal of this work? (e.g., understanding matrix multiplication from scratch, learning optimization techniques?)

  2. What were your biggest "aha!" moments or insights? For example:

    • Understanding how the triple nested loop works?
    • Seeing how broadcasting eliminates loops?
    • The memory efficiency of broadcasting?
  3. What surprised you most during this exploration?

  4. What would you want other learners to take away from reading your summary?

Once you share your thoughts on these, I can help you structure them into a clear, engaging blog post or summary!

Prompt: 171

I wanted to understand matrix mulitplication from scratch by implementing it a few different ways, i also wanted to learn how the techiniques used work at a fundamental level.

All the above work was just full of aha moments from digging deeper into the underlying functions to broadcasting, memory efficient ways etc so generally most of what i did above was interesting

I want people to undertsand matrrix multiplcation, how to use it, compute efficient methods of doing this, how this can be applied through broadcasting , how this relates to machine learning and why it is important

Assistant: 265

Great goals! Let me suggest a structure for your blog, and you can tell me what you'd like to add or change:

Proposed Outline:

  1. Introduction - Why matrix multiplication matters in ML
  2. The Naive Approach - Triple nested loop (your starting point)
  3. Optimization Journey
    • Using Numba for the inner loop
    • Element-wise operations with broadcasting
    • Einstein summation
  4. Key Concepts Explained
    • Broadcasting rules and memory efficiency
    • The dot product connection
  5. Performance Comparison - Your timing results
  6. Takeaways - What you learned and why it matters

Does this structure work for you? Would you like to:

  • Add or remove any sections?
  • Start drafting a specific section first?
  • Include more code examples or keep it more conceptual?

Prompt: 40 ()

the above structure is fine with me , i do want to keep all my code while explaining it , take the content from my work above

Prompt: 15 ()

continue , just draft the whole thing at once

Prompt: 18

put the above answer into a note cell like you did above

Code: 7 ()

from dialoghelper import *

Prompt: 7

import dialoguehelper for me

Assistant: 64

I'm not familiar with dialoguehelper - is this a module you've created, or something specific you'd like me to help you find?

Could you tell me more about what it does or where it comes from?

Prompt: 21 ()

create a note, put the summary, use your tools to do this

Assistant: 405

add_msg(content=# Matrix Multiplication from Scratch: A Deep Dive …, msg_type=note)
Copied!
{
  "id": "toolu_01BpfVpb4SrNrsN5gDGE4vua",
  "call": {
    "function": "add_msg",
    "arguments": {
      "content": "# Matrix Multiplication from Scratch: A Deep Dive\n\n## Introduction\n\nMatrix multiplication is one of the most fundamental operations in machine learning and deep learning. Every time a neural network makes a prediction, matrix multiplication is happening behind the scenes - multiplying input data by weights to produce outputs.\n\nIn this exploration, I'll implement matrix multiplication from scratch in several different ways, starting from the most naive approach and progressively optimizing it. Along the way, we'll learn about:\n- How matrix multiplication actually works at a fundamental level\n- <TRUNCATED>",
      "msg_type": "note"
    }
  },
  "result": "'_da5821fb'"
}

Done! I've created a note message with your complete matrix multiplication summary. The note has been added to your dialog and should now be visible in your notebook.

Would you like me to make any changes to the content?

Note: 6,481

Matrix Multiplication from Scratch: A Deep Dive

Introduction

Matrix multiplication is one of the most fundamental operations in machine learning and deep learning. Every time a neural network makes a prediction, matrix multiplication is happening behind the scenes - multiplying input data by weights to produce outputs.

In this exploration, I'll implement matrix multiplication from scratch in several different ways, starting from the most naive approach and progressively optimizing it. Along the way, we'll learn about:

  • How matrix multiplication actually works at a fundamental level
  • Memory-efficient techniques like broadcasting
  • How to achieve massive speedups (from 500ms to microseconds!)
  • Why these concepts matter for machine learning

We'll use the MNIST dataset - 784 pixel values per image that we'll multiply by a weight matrix to get predictions for 10 digit classes.

Setting Up Our Data

First, let's load our data and create the matrices we'll be working with:

Copied!
from pathlib import Path
import gzip, pickle
import torch
from torch import tensor

# Load MNIST data
MNIST_URL = 'https://github.com/mnielsen/neural-networks-and-deep-learning/blob/master/data/mnist.pkl.gz?raw=true'
path_data = Path('data')
path_data.mkdir(exist_ok=True)
path_gz = path_data/'mnist.pkl.gz'

from urllib.request import urlretrieve
if not path_gz.exists():
    urlretrieve(MNIST_URL, path_gz)

with gzip.open(path_gz, 'rb') as f:
    ((x_train, y_train),(x_valid, y_valid),_) = pickle.load(f, encoding='latin-1')

x_train,y_train,x_valid,y_valid = map(tensor,(x_train,y_train,x_valid,y_valid))

Python Fundamentals Explored Along the Way

Before diving into matrix multiplication, I explored several Python concepts that proved useful:

Context Managers

Context managers handle resource management automatically - setup and cleanup:

Copied!
with gzip.open(path_gz, 'rb') as f:
    ((x_train, y_train),(x_valid, y_valid),_) = pickle.load(f, encoding='latin-1')

The with statement ensures the file is closed even if an error occurs. This is safer than manually calling f.close().

Generators and yield

Generators create iterators that produce values on-demand rather than storing everything in memory:

Copied!
def chunks(x, sz):
    for i in range(0, len(x), sz):
        yield x[i:i+sz]

# Creates an iterator that yields chunks of size 5
val_iter = chunks(vals, 5)
next(val_iter)  # Get first chunk
next(val_iter)  # Get second chunk
# next(val_iter)  # Would raise StopIteration - no more values

yield pauses execution and returns a value, then resumes where it left off on the next call.

iter() and islice()

We can convert any sequence into an iterator and slice it efficiently:

Copied!
from itertools import islice

it = iter(vals)
next(it)  # Get next value

# Get first 5 items without loading all into memory
is_it = islice(vals, 5)
list(is_it)  # [0.0, 0.0, 0.0, 0.19140625, 0.9296875]

Using iter() with a callable and sentinel value creates an iterator that calls the function until it returns the sentinel:

Copied!
it = iter(lst1)
img = list(iter(lambda: list(islice(it, 28)), []))
# Calls the lambda repeatedly until it returns []

Lambda Functions

Anonymous functions for simple operations:

Copied!
caps_name = lambda func: func.upper()
caps_name('silver rubanza')  # 'SILVER RUBANZA'

sum_lmd = lambda a,b: a+b
sum_lmd(1,2)  # 3

Custom Classes with init and getitem

Classes let us create custom behavior for our objects:

Copied!
class Matrix:
    def __init__(self, xs): 
        self.xs = xs
    
    def __getitem__(self, idxs): 
        return self.xs[idxs[0]][idxs[1]]

m = Matrix(img)
m[20,15]  # Now we can use tuple indexing!

__getitem__ allows square bracket indexing. __init__ initializes object attributes when created.

Random Number Generator (Wichman-Hill)

Implemented a pseudo-random number generator from scratch to understand how randomness works:

Copied!
rnd_state = None 

def seed(a):
    global rnd_state
    a,x = divmod(a, 30268)
    a,y = divmod(a, 30306)
    a,z = divmod(a, 30322)
    rnd_state = int(x)+1, int(y)+1, int(z)+1

def rand():
    global rnd_state
    x,y,z = rnd_state
    x = (171*x) % 30269
    y = (172*y) % 30307
    z = (170*z) % 30323
    rnd_state = x,y,z 
    return (x/30269 + y/30307 + z/30323)%1.0

The map() Function

map() applies a function to every item in an iterable:

Copied!
x_train,y_train,x_valid,y_valid = map(tensor,(x_train,y_train,x_valid,y_valid))
# Converts all four arrays to tensors in one line

divmod()

Returns both quotient and remainder in one operation:

Copied!
a,x = divmod(457428938475, 30268)
# a = 15113085, x = 6795
# Same as: a = 457428938475 // 30268; x = 457428938475 % 30268

Destructuring/Unpacking

Extract multiple values from sequences in one line:

Copied!
ar, ac = m1.shape  # Unpack shape tuple
x,y,z = rnd_state  # Unpack state tuple

# Works with nested structures too:
((x_train, y_train),(x_valid, y_valid),_) = pickle.load(f, encoding='latin-1')

set_printoptions()

Control how tensors/arrays display:

Copied!
torch.set_printoptions(precision=2, linewidth=140, sci_mode=False)
np.set_printoptions(precision=2, linewidth=140)

This makes output more readable by controlling decimal places and line wrapping.

Now let's set up our matrices for multiplication:

Copied!
# Create random weights and bias
weights = torch.randn(784,10)
bias = torch.zeros(10)

# Create a mini-batch of 5 images
m1 = x_valid[:5]  # shape: (5, 784)
m2 = weights       # shape: (784, 10)

# Get dimensions for our loops
ar, ac = m1.shape  # 5 rows, 784 columns
br, bc = m2.shape  # 784 rows, 10 columns

Key insight about matrix multiplication dimensions: When multiplying an (m, n) matrix by an (n, p) matrix, the inner dimensions must match (both are n), and the result has shape (m, p). In our case: (5, 784) × (784, 10) → (5, 10).

The Naive Approach: Triple Nested Loop

Let's start with the most straightforward implementation - three nested loops:

Copied!
def matmul(a,b):
    ar, ac = a.shape
    br, bc = b.shape
    c = torch.zeros(ar, bc)
    for i in range(ar):         # for each row in a
        for j in range(bc):     # for each column in b
            for k in range(ac): # for each element in that row/column
                c[i,j] += a[i,k] * b[k,j]
    return c

What's happening here?

  • The outer loop (i) goes through each row of m1 (our 5 images)
  • The middle loop (j) goes through each column of m2 (our 10 weight vectors)
  • The inner loop (k) multiplies the 784 values in row i by the 784 values in column j and sums them up

For each position [i,j] in our result, we're computing the dot product of row i from m1 with column j from m2.

Copied!
%time _=matmul(m1,m2)
# CPU times: user 554 ms, sys: 0 ns, total: 554 ms

554 milliseconds for just 5 images! This is painfully slow because Python loops have significant overhead.

Optimization 1: Numba for the Inner Loop

Numba is a JIT compiler that translates Python code into fast machine code. However, Numba only works with NumPy arrays, not PyTorch tensors. So we'll apply it to just the innermost computation - the dot product:

Copied!
from numba import njit

@njit
def dot(a,b):
    res = 0.
    for i in range(len(a)):
        res += a[i]*b[i]
    return res

The first time you call a Numba-decorated function, it compiles (slow). After that, it runs at near-C speed:

Copied!
from numpy import array
%time dot(array([1.,2,3]),array([1.,2,3]))  # First call: ~200ms (compilation)
%time dot(array([1.,2,3]),array([1.,2,3]))  # Second call: ~21 microseconds!

Now we can use this optimized dot product in our matmul:

Copied!
def matmul(a,b):
    ar, ac = a.shape
    br, bc = b.shape
    c = torch.zeros(ar,bc)
    for i in range(ar):
        for j in range(bc):
            c[i,j] = dot(a[i,:], b[:,j])
    return c

What changed? Instead of the inner k loop, we now call dot(a[i,:], b[:,j]):

  • a[i,:] gets the entire i-th row (all 784 pixel values for one image)
  • b[:,j] gets the entire j-th column (all 784 weights for one class)
  • dot() computes their dot product in optimized machine code

This is functionally identical to our triple loop, but the innermost computation now runs in compiled code instead of Python.

Understanding Element-wise vs Matrix Multiplication

Before we go further, let's clarify an important distinction:

Element-wise multiplication (*):

Copied!
m = tensor([[1.,2,3],[4,5,6],[7,8,9]])
m * m
# tensor([[ 1.,  4.,  9.],
#         [16., 25., 36.],
#         [49., 64., 81.]])

Each element is multiplied by the element in the same position. Both matrices must have the same shape, and the result has the same shape.

Matrix multiplication (@):

Copied!
m @ m
# tensor([[ 30.,  36.,  42.],
#         [ 66.,  81.,  96.],
#         [102., 126., 150.]])

Each element [i,j] is the dot product of row i and column j. The inner dimensions must match.

For example, position [0,0] in m @ m:

  • Row 0: [1, 2, 3]
  • Column 0: [1, 4, 7]
  • Calculation: (1×1) + (2×4) + (3×7) = 1 + 8 + 21 = 30

Optimization 2: Using PyTorch's Element-wise Operations

We can replace our dot product with PyTorch's built-in element-wise multiply and sum:

Copied!
def matmul(a,b):
    ar, ac = a.shape
    br, bc = b.shape
    c = torch.zeros(ar, bc)
    for i in range(ar):
        for j in range(bc):
            c[i,j] = (a[i,:]*b[:,j]).sum()
    return c

This does the same thing as our dot function: multiply corresponding elements and sum them. PyTorch has a built-in for this:

Copied!
def matmul(a,b):
    ar, ac = a.shape
    br, bc = b.shape
    c = torch.zeros(ar, bc)
    for i in range(ar):
        for j in range(bc):
            c[i,j] = torch.dot(a[i,:],b[:,j])
    return c

Broadcasting: The Key to Eliminating Loops

Broadcasting allows operations between tensors of different shapes without explicitly copying data. This is where things get really powerful.

Broadcasting Basics

Copied!
c = tensor([10.,20.,30])
m = tensor([[1.,2,3],[4,5,6],[7,8,9]])

m + c
# tensor([[11., 22., 33.],
#         [14., 25., 36.],
#         [17., 28., 39.]])

Even though m has shape (3,3) and c has shape (3,), PyTorch automatically "stretches" c to match m's shape. Behind the scenes, it uses expand_as:

Copied!
t = c.expand_as(m)
# tensor([[10., 20., 30.],
#         [10., 20., 30.],
#         [10., 20., 30.]])

The crucial insight: This expansion doesn't actually copy data in memory! PyTorch uses "strides" to make it appear as if the data is repeated, making broadcasting extremely memory-efficient.

Broadcasting Rules

Two dimensions are compatible when:

  1. They are equal, OR
  2. One of them is 1

Comparison happens right-to-left. Missing dimensions are treated as size 1.

Adding Dimensions with None/unsqueeze

We can reshape tensors to control how broadcasting works:

Copied!
c.shape                    # torch.Size([3])
c.unsqueeze(0).shape       # torch.Size([1, 3]) - add dimension at start
c.unsqueeze(1).shape       # torch.Size([3, 1]) - add dimension at end

# Shorthand using None:
c[None].shape              # torch.Size([1, 3])
c[:,None].shape            # torch.Size([3, 1])
c[...,None].shape          # torch.Size([3, 1]) - ... means "all existing dims"

Outer Product with Broadcasting

When we combine differently shaped tensors, broadcasting can create powerful operations:

Copied!
c[None] * c[:,None]
# Shape (1,3) * (3,1) -> broadcasts to (3,3) * (3,3)
# tensor([[100., 200., 300.],
#         [200., 400., 600.],
#         [300., 600., 900.]])

What's happening:

  • c[None] shape (1, 3) → expands rows: [[10, 20, 30], [10, 20, 30], [10, 20, 30]]
  • c[:,None] shape (3, 1) → expands columns: [[10, 10, 10], [20, 20, 20], [30, 30, 30]]
  • Element-wise multiplication gives us every combination - an outer product!

Optimization 3: Matrix Multiplication with Broadcasting

Now we can eliminate the inner j loop using broadcasting:

Copied!
def matmul(a,b):
    ar,ac = a.shape
    br,bc = b.shape
    c = torch.zeros(ar,bc)
    for i in range(ar):
        c[i] = (a[i,:,None]*b).sum(dim=0)
    return c

What's happening in a[i,:,None]*b?

  • a[i,:] is shape (784,) - one row of pixel values
  • a[i,:,None] is shape (784, 1) - add a dimension
  • b is shape (784, 10) - our weight matrix
  • Broadcasting: (784, 1) * (784, 10)(784, 10)

This multiplies each pixel value by ALL 10 weight columns at once! Then .sum(dim=0) sums down the 784 dimension, giving us 10 output values.

We've eliminated one loop entirely through broadcasting.

Optimization 4: Einstein Summation

Einstein summation (einsum) provides a compact notation for tensor operations. You label dimensions with letters, and any letter that appears in inputs but NOT in the output gets summed over.

Copied!
def matmul(a,b): 
    return torch.einsum('ik,kj->ij', a, b)

Decoding 'ik,kj->ij':

  • First tensor a has dimensions i (rows) and k (columns)
  • Second tensor b has dimensions k (rows) and j (columns)
  • Output has dimensions i and j
  • The k dimension appears in both inputs but NOT in output → summed over

This is exactly matrix multiplication in a single, elegant expression!

Copied!
%timeit matmul(x_train,m2)
# Much faster than our loop-based versions!

The Frobenius Norm: A Practical Application

The Frobenius norm demonstrates element-wise operations in practice:

Copied!
m = tensor([[1.,2,3],[4,5,6],[7,8,9]])

# Square each element, sum all, take square root
sf = (m*m).sum()  # 285
sf.sqrt()         # tensor(16.8819)

This is useful for measuring the "size" of a matrix, often used in regularization and loss functions.

Performance Summary

Here's how our different implementations compare:

Method Time Speedup
Triple nested loop ~550ms 1x (baseline)
Numba dot product ~268µs ~2000x
Broadcasting ~137µs ~4000x
Einstein summation ~15ms ~37x
PyTorch matmul ~15ms ~37x
GPU (CUDA) ~458µs ~1200x

The journey from 550ms to microseconds shows why understanding these optimizations matters!

Key Takeaways

  1. Matrix multiplication is fundamental to ML - it's how inputs get transformed by weights in neural networks.

  2. Python loops are slow - the triple nested loop works but is impractical for real applications.

  3. Broadcasting is powerful and efficient - it allows operations between different-shaped tensors without copying memory.

  4. Understand the broadcasting rules - dimensions are compatible if they're equal or one is 1, compared right-to-left.

  5. Einstein summation is elegant - complex tensor operations can be expressed in a single compact notation.

  6. GPUs provide massive speedups - for large-scale ML, moving computations to GPU is essential.

  7. Build understanding step by step - implementing from scratch reveals what's really happening in those convenient library functions.

Further Exploration

  • Try implementing other operations (transpose, inverse) from scratch
  • Explore how convolutions relate to matrix multiplication
  • Learn about GPU programming with CUDA
  • Study how automatic differentiation uses these operations

Understanding matrix multiplication deeply gives you intuition for how neural networks learn and why certain architectures work better than others.

Prompt: 16

Add all the little tricks, classes and functions i explored

Assistant: 279

I can help you add those details! Let me go through what I noticed you explored. Could you tell me which of these you'd like me to add to the summary?

Functions & Concepts:

  1. Context managers (with statements)
  2. Generators and yield
  3. iter() and islice()
  4. Lambda functions
  5. Custom classes (your Cows and Matrix classes with __init__ and __getitem__)
  6. Random number generator (Wichman-Hill algorithm)
  7. map() function
  8. divmod()
  9. Destructuring/unpacking
  10. set_printoptions() for display
  11. os.fork() experiments

Which of these would you like me to incorporate into the summary? Or would you like me to add all of them in a separate "Python Fundamentals" section?

Prompt: 111

add them to my summary, in the order in which i explored them for example when exploring itertools i explore iter through iter??, i want the tricks to be put where in the sequence i explored them so as to keep context of why i was exploring them , create a new note, dont a direct edit of my last incase i decide to keep the current version

Assistant: 502

I'll create a new note with all the Python tricks integrated in the order you explored them. Let me add this for you:

add_msg(content=# Matrix Multiplication from Scratch: A Deep Dive …, msg_type=note)
Copied!
{
  "id": "toolu_01QWDRZZesVASMojtRv8N9BY",
  "call": {
    "function": "add_msg",
    "arguments": {
      "content": "# Matrix Multiplication from Scratch: A Deep Dive (Extended Version)\n\n## Introduction\n\nMatrix multiplication is one of the most fundamental operations in machine learning and deep learning. Every time a neural network makes a prediction, matrix multiplication is happening behind the scenes - multiplying input data by weights to produce outputs.\n\nIn this exploration, I'll implement matrix multiplication from scratch in several different ways, starting from the most naive approach and progressively optimizing it. Along the way, we'll learn about:\n- How matrix multiplication actually works at a f<TRUNCATED>",
      "msg_type": "note"
    }
  },
  "result": "'_3c9a5a13'"
}

Done! I've created a new extended version of your summary with all the Python fundamentals integrated in the order you explored them. The original summary is still there untouched.

The new version includes context for why you explored each concept (like exploring iter and islice when working with image chunks, or creating the Matrix class when you needed tuple indexing).

Would you like me to adjust anything?

Note: 9,286

Matrix Multiplication from Scratch: A Deep Dive (Extended Version)

Introduction

Matrix multiplication is one of the most fundamental operations in machine learning and deep learning. Every time a neural network makes a prediction, matrix multiplication is happening behind the scenes - multiplying input data by weights to produce outputs.

In this exploration, I'll implement matrix multiplication from scratch in several different ways, starting from the most naive approach and progressively optimizing it. Along the way, we'll learn about:

  • How matrix multiplication actually works at a fundamental level
  • Memory-efficient techniques like broadcasting
  • How to achieve massive speedups (from 500ms to microseconds!)
  • Why these concepts matter for machine learning
  • Python fundamentals that make all this possible

We'll use the MNIST dataset - 784 pixel values per image that we'll multiply by a weight matrix to get predictions for 10 digit classes.

Setting Up Our Data

First, let's load our data and create the matrices we'll be working with.

Understanding urllib and urlretrieve

Before downloading data, let's understand the tools:

  • urllib - Python's standard library for working with URLs (built-in, no installation needed)
  • urllib.request - Submodule that handles opening and reading URLs
  • urlretrieve() - Downloads a file from a URL and saves it locally
Copied!
from urllib.request import urlretrieve
urlretrieve?  # Check the signature
# urlretrieve(url, filename=None, reporthook=None, data=None)
# Returns a tuple: (path to file, HTTPMessage object)

Now let's download MNIST:

Copied!
from pathlib import Path
import gzip, pickle

MNIST_URL = 'https://github.com/mnielsen/neural-networks-and-deep-learning/blob/master/data/mnist.pkl.gz?raw=true'
path_data = Path('data')
path_data.mkdir(exist_ok=True)
path_gz = path_data/'mnist.pkl.gz'

if not path_gz.exists():
    urlretrieve(MNIST_URL, path_gz)

Context Managers with with

To load the data, we use a context manager:

Copied!
with gzip.open(path_gz, 'rb') as f:
    ((x_train, y_train),(x_valid, y_valid),_) = pickle.load(f, encoding='latin-1')

What's a context manager? A Python object that manages resources automatically - it handles setup and cleanup for you.

The with statement:

  1. Opens the gzip file when entering the block
  2. Automatically closes the file when exiting (even if an error occurs)

Without it, you'd need:

Copied!
f = gzip.open(path_gz, 'rb')
# ... do stuff ...
f.close()  # Easy to forget!

The as f part gives you a name to refer to the opened file within the block.

Generators and yield

Now let's explore how to work with our image data in chunks. First, let's look at a small sample:

Copied!
lst1 = list(x_train[0])
vals = lst1[200:210]  # 10 pixel values

We can create a function that splits data into chunks:

Copied!
def chunks(x, sz):
    for i in range(0, len(x), sz):
        yield x[i:i+sz]

What does yield do? yield creates a generator - an iterator that produces values on-demand:

Copied!
list(chunks(vals, 5))
# [[0.0, 0.0, 0.0, 0.19140625, 0.9296875],
#  [0.98828125, 0.98828125, 0.98828125, 0.98828125, 0.98828125]]

val_iter = chunks(vals, 5)
val_iter  # <generator object chunks at 0x...>

next(val_iter)  # [0.0, 0.0, 0.0, 0.19140625, 0.9296875]
next(val_iter)  # [0.98828125, 0.98828125, 0.98828125, 0.98828125, 0.98828125]
# next(val_iter)  # StopIteration error - no more values!

yield returns values one at a time and pauses execution. When you call next(), it resumes where it left off. When there are no more values, it raises StopIteration.

When you use a generator in a for loop or list(), Python catches StopIteration automatically.

Visualizing MNIST Digits

We can use our chunks function to reshape flat pixel data into a 28×28 image:

Copied!
import matplotlib.pyplot as plt
import matplotlib as mpl

mpl.rcParams['image.cmap'] = 'gray'  # Set grayscale colormap
plt.imshow(list(chunks(lst1, 28)));

Note on rcParams: Runtime configuration parameters customize matplotlib styling. This takes precedence over other styling methods.

Wrapping chunks(lst1, 28) in list() forces the generator to return all results at once - it runs through the entire iterator and creates a list.

iter() and islice()

We can also work with iterators more directly:

Copied!
from itertools import islice

# Convert any sequence to an iterator
it = iter(vals)
next(it)  # 0.0
next(it)  # 0.0
next(it)  # 0.0

What is iter()?

Copied!
iter??
# iter(iterable) -> iterator
# iter(callable, sentinel) -> iterator
# 
# Get an iterator from an object. In the first form, the argument must
# supply its own iterator, or be a sequence.
# In the second form, the callable is called until it returns the sentinel.

Using islice() to grab first n items:

Copied!
islice??
# islice(iterable, stop) --> islice object
# islice(iterable, start, stop[, step]) --> islice object
# 
# Return an iterator whose next() method returns selected values from an
# iterable. Works like a slice() on a list but returns an iterator.

is_it = islice(vals, 5)
next(is_it)  # 0.0
next(is_it)  # 0.0
# ... continues for 5 values total

# Or get all at once:
list(islice(vals, 5))  # [0.0, 0.0, 0.0, 0.19140625, 0.9296875]

iter() with callable and sentinel

The second form of iter() is powerful - it calls a function repeatedly until it returns a sentinel value:

Copied!
it = iter(lst1)
img = list(iter(lambda: list(islice(it, 28)), []))

What's happening here?

  • lambda: list(islice(it, 28)) - anonymous function that grabs 28 values
  • [] - the sentinel (stop when lambda returns empty list)
  • iter() calls the lambda repeatedly until it returns []
  • Result: list of 28-element lists (our 28×28 image!)
Copied!
plt.imshow(img);

Lambda Functions

Lambda creates anonymous functions - functions without names:

Copied!
# Syntax: lambda arguments : expression

caps_name = lambda func: func.upper()
caps_name('silver rubanza')  # 'SILVER RUBANZA'
caps_name('devontay')  # 'DEVONTAY'

sum_lmd = lambda a,b: a+b
sum_lmd(1,2)  # 3

They're useful for simple operations you need to pass as arguments.

Custom Classes: Creating a Matrix Type

Lists don't support tuple indexing like img[20,10]. Let's create a class that does:

Copied!
img[20][15]  # Works
# img[20,10]  # TypeError: list indices must be integers or slices, not tuple

Understanding Classes:

Object-oriented programming uses classes as templates to create objects. An object has:

  • Identity (id): Address in memory (unchangeable)
  • Type: What operations you can do (unchangeable)
  • Value: The data it holds (mutable or immutable)

Creating a Custom Class:

Copied!
class Matrix:
    def __init__(self, xs): 
        self.xs = xs
    
    def __getitem__(self, idxs): 
        return self.xs[idxs[0]][idxs[1]]

m = Matrix(img)
m[20,15]  # Now works!

What's happening:

  • __init__ is a constructor - called when creating the object, initializes attributes
  • __getitem__ is a dunder (double underscore) method - allows square bracket indexing
  • When you write m[20,15], Python calls m.__getitem__((20,15))

Example: A Cows Class

Copied!
class Cows:
    """
    A class representing different cow breeds
    
    Parameters:
    breed: The breed of a cow
    age: The cows age
    weight: The weight of a cow
    color: The color of a cow
    """
    def __init__(self, breed, age, weight, color):
        self.breed = breed
        self.age = age
        self.weight = weight
        self.color = color

# Create an instance
fresian = Cows("Fresian", 4, 600, "Distinctive black and white patches")
fresian.breed  # 'Fresian'
fresian.age  # 4

Using default parameters:

Copied!
class Cows_dp:
    def __init__(self, breed="Heifer", age=5, weight=600, color="Black"):
        self.breed = breed
        self.age = age
        self.weight = weight
        self.color = color

heifer = Cows_dp()
heifer.breed  # 'Heifer' (uses default)

Random Number Generator from Scratch

Before using PyTorch's random functions, let's implement the Wichman-Hill algorithm (used in Python before 2.3):

Copied!
rnd_state = None 

def seed(a):
    global rnd_state
    a,x = divmod(a, 30268)
    a,y = divmod(a, 30306)
    a,z = divmod(a, 30322)
    rnd_state = int(x)+1, int(y)+1, int(z)+1

def rand():
    global rnd_state
    x,y,z = rnd_state
    x = (171*x) % 30269
    y = (172*y) % 30307
    z = (170*z) % 30323
    rnd_state = x,y,z 
    return (x/30269 + y/30307 + z/30323)%1.0

Understanding divmod():

Copied!
divmod??
# divmod(x, y) -> (quotient, remainder)
# Return the tuple (x//y, x%y)

r,s = divmod(457428938475, 30268)
# r = 15113085 (quotient)
# s = 6795 (remainder)

Testing our random generator:

Copied!
seed(457428938475)
rnd_state  # (6796, 15484, 27612)

rand()  # 0.6847...
rand()  # 0.2578...
rand()  # 0.9214...

We can verify it's not correlated:

Copied!
plt.plot([rand() for _ in range(50)]);  # Should look random
plt.hist([rand() for _ in range(10000)]);  # Should be uniform

Converting to PyTorch Tensors

Now let's convert our data to PyTorch tensors:

Copied!
import torch
from torch import tensor
import numpy as np

Setting print options for readability:

Copied!
torch.set_printoptions?
# Set options for printing. Parameters include:
# - precision: number of digits
# - linewidth: characters per line
# - sci_mode: scientific notation

torch.set_printoptions(precision=2, linewidth=140, sci_mode=False)
np.set_printoptions(precision=2, linewidth=140)

Using map() to convert multiple arrays:

Copied!
map??
# map(func, *iterables) --> map object
# Make an iterator that computes the function using arguments from
# each of the iterables.

x_train,y_train,x_valid,y_valid = map(tensor, (x_train,y_train,x_valid,y_valid))
# Applies tensor() to each of the four arrays

Tensor indexing:

Copied!
tens = tensor(img)
tens[20,15]  # Works with tensors!
tens[20][15]  # Also works

Comparing Random Generators Performance

Copied!
%timeit -n 10 list(chunks([rand() for _ in range(7840)], 10))
# 3.66 ms ± 68.5 µs

%timeit -n 10 torch.rand(784,10)
# 87.5 µs ± 2.57 µs

PyTorch is ~42x faster! We'll use torch.rand() for the rest.

Preparing for Matrix Multiplication

Now let's set up our matrices:

Copied!
weights = torch.randn(784,10)
bias = torch.zeros(10)

m1 = x_valid[:5]  # shape: (5, 784)
m2 = weights       # shape: (784, 10)

# Destructuring to get dimensions
ar, ac = m1.shape  # 5 rows, 784 columns
br, bc = m2.shape  # 784 rows, 10 columns
(ar,ac), (br,bc)  # ((5, 784), (784, 10))

Key insight about matrix multiplication dimensions:

When multiplying an (m, n) matrix by an (n, p) matrix:

  • Inner dimensions must match (both are n)
  • Result has shape (m, p)
  • In our case: (5, 784) × (784, 10) → (5, 10)

Create result container:

Copied!
t1 = torch.zeros(ar, bc)  # (5, 10) - same shape as expected output

The Naive Approach: Triple Nested Loop

Let's start with the most straightforward implementation:

Copied!
for i in range(ar):
    for j in range(bc):
        for k in range(ac):
            t1[i,j] += m1[i,k] * m2[k,j]

What's happening:

  • Outer loop (i): goes through 5 images (rows of m1)
  • Middle loop (j): goes through 10 weight columns (columns of m2)
  • Inner loop (k): multiplies 784 pixel values by 784 weights and sums

For each position [i,j]:

  • Take row i from m1: all 784 pixel values
  • Take column j from m2: all 784 weights for that class
  • Multiply corresponding elements and sum: (pixel₀×weight₀) + (pixel₁×weight₁) + ... + (pixel₇₈₃×weight₇₈₃)

Wrapping in a function:

Copied!
def matmul(a,b):
    ar, ac = a.shape
    br, bc = b.shape
    c = torch.zeros(ar, bc)
    for i in range(ar):
        for j in range(bc):
            for k in range(ac):
                c[i,j] += a[i,k] * b[k,j]
    return c

%time _=matmul(m1,m2)
# CPU times: user 554 ms, sys: 0 ns, total: 554 ms

554 milliseconds for just 5 images! Python loops are slow.

Optimization 1: Numba for the Inner Loop

Numba is a JIT (Just-In-Time) compiler that translates Python to fast machine code.

Important: Numba works with NumPy arrays, not PyTorch tensors.

Let's optimize just the innermost computation - the dot product:

Copied!
from numba import njit

@njit
def dot(a,b):
    res = 0.
    for i in range(len(a)):
        res += a[i]*b[i]
    return res

First call compiles, subsequent calls are fast:

Copied!
from numpy import array
%time dot(array([1.,2,3]), array([1.,2,3]))  
# First call: ~200ms (compilation time)

%time dot(array([1.,2,3]), array([1.,2,3]))  
# Second call: ~21 microseconds! (~10,000x faster)

Using in matmul:

Copied!
def matmul(a,b):
    ar, ac = a.shape
    br, bc = b.shape
    c = torch.zeros(ar,bc)
    for i in range(ar):
        for j in range(bc):
            c[i,j] = dot(a[i,:], b[:,j])
    return c

# Convert to NumPy for Numba
m1a, m2a = m1.numpy(), m2.numpy()
matmul(m1a, m2a)

What changed?

  • a[i,:] gets entire row i (all 784 values)
  • b[:,j] gets entire column j (all 784 values)
  • dot() computes their dot product in compiled machine code

We've replaced the inner k loop with a single fast function call!

Understanding Element-wise vs Matrix Multiplication

Before optimizing further, let's clarify a crucial distinction:

Element-wise multiplication (*):

Copied!
m = tensor([[1.,2,3],[4,5,6],[7,8,9]])
m * m
# tensor([[ 1.,  4.,  9.],
#         [16., 25., 36.],
#         [49., 64., 81.]])

Multiplies corresponding elements in the same positions. Both matrices must have identical shapes.

Matrix multiplication (@):

Copied!
m @ m
# tensor([[ 30.,  36.,  42.],
#         [ 66.,  81.,  96.],
#         [102., 126., 150.]])

Each element [i,j] is the dot product of row i and column j.

Example: Position [0,0] in m @ m:

  • Row 0: [1, 2, 3]
  • Column 0: [1, 4, 7]
  • Calculation: (1×1) + (2×4) + (3×7) = 1 + 8 + 21 = 30

Optimization 2: Using PyTorch's Element-wise Operations

We can replace our dot product with element-wise multiply and sum:

Copied!
def matmul(a,b):
    ar, ac = a.shape
    br, bc = b.shape
    c = torch.zeros(ar, bc)
    for i in range(ar):
        for j in range(bc):
            c[i,j] = (a[i,:]*b[:,j]).sum()
    return c

Or use PyTorch's built-in dot product:

Copied!
def matmul(a,b):
    ar, ac = a.shape
    br, bc = b.shape
    c = torch.zeros(ar, bc)
    for i in range(ar):
        for j in range(bc):
            c[i,j] = torch.dot(a[i,:],b[:,j])
    return c

%timeit matmul(m1,m2)

The Frobenius Norm: A Practical Example

Before diving into broadcasting, let's see element-wise operations in action:

Copied!
m = tensor([[1.,2,3],[4,5,6],[7,8,9]])

# Frobenius norm: square each element, sum all, take square root
sf = (m*m).sum()  # tensor(285.)
sf.sqrt()  # tensor(16.88)

This measures the "size" of a matrix, used in regularization and loss functions.

Broadcasting: The Key to Eliminating Loops

Broadcasting allows operations between tensors of different shapes without copying data.

Broadcasting Basics

Copied!
c = tensor([10.,20.,30])
m = tensor([[1.,2,3],[4,5,6],[7,8,9]])

m.shape  # torch.Size([3, 3])
c.shape  # torch.Size([3])

m + c
# tensor([[11., 22., 33.],
#         [14., 25., 36.],
#         [17., 28., 39.]])

Even though shapes don't match, PyTorch "stretches" c to match m!

Behind the scenes with expand_as:

Copied!
t = c.expand_as(m)
# tensor([[10., 20., 30.],
#         [10., 20., 30.],
#         [10., 20., 30.]])

m + t  # Same result as m + c

The crucial insight - no memory copying:

Copied!
t.storage()  # Only contains: 10.0, 20.0, 30.0
# PyTorch uses "strides" to make it appear repeated!

Broadcasting is extremely memory-efficient - it doesn't actually copy data.

Broadcasting Rules

From the NumPy documentation:

Two dimensions are compatible when:

  1. They are equal, OR
  2. One of them is 1

Comparison happens right-to-left. Missing dimensions are treated as size 1.

Examples:

Copied!
# Works:
Image  (3d): 256 x 256 x 3
Scale  (1d):             3
Result (3d): 256 x 256 x 3

# Works:
A      (4d):  8 x 1 x 6 x 1
B      (3d):      7 x 1 x 5
Result (4d):  8 x 7 x 6 x 5

# Doesn't work:
A      (1d):  3
B      (1d):  4  # Trailing dimensions don't match

Adding Dimensions with None/unsqueeze

We can reshape tensors to control broadcasting:

Copied!
c.shape  # torch.Size([3])

c.unsqueeze(0).shape  # torch.Size([1, 3]) - add at start
c.unsqueeze(1).shape  # torch.Size([3, 1]) - add at end

# Shorthand using None:
c[None].shape      # torch.Size([1, 3])
c[:,None].shape    # torch.Size([3, 1])

# Ellipsis (...) means "all existing dimensions":
c[...,None].shape  # torch.Size([3, 1])

Why are trailing colons optional?

Copied!
c[None,:]  # Explicit: take all of remaining dimension
c[None]    # Implicit: assumed you want everything else

# Both give same result: torch.Size([1, 3])

When you stop indexing early, Python assumes you want all remaining dimensions.

Outer Product with Broadcasting

Combining different shapes creates powerful operations:

Copied!
c[None] * c[:,None]
# Shape (1,3) * (3,1) -> broadcasts to (3,3)
# tensor([[100., 200., 300.],
#         [200., 400., 600.],
#         [300., 600., 900.]])

What's happening:

  • c[None] shape (1, 3) → expands rows to (3, 3): [[10, 20, 30], [10, 20, 30], [10, 20, 30]]
  • c[:,None] shape (3, 1) → expands columns to (3, 3): [[10, 10, 10], [20, 20, 20], [30, 30, 30]]
  • Element-wise multiplication gives every combination - an outer product!

This is how you multiply every element from one vector with every element from another.

Optimization 3: Matrix Multiplication with Broadcasting

Now we can eliminate the j loop:

Copied!
def matmul(a,b):
    ar,ac = a.shape
    br,bc = b.shape
    c = torch.zeros(ar,bc)
    for i in range(ar):
        c[i] = (a[i,:,None]*b).sum(dim=0)
    return c

Breaking down a[i,:,None]*b:

Copied!
digit = m1[0]  # shape: (784,)
digit[:,None].shape  # (784, 1)
m2.shape  # (784, 10)

# Broadcasting: (784, 1) * (784, 10) → (784, 10)
ca = m2 * digit[:,None]
ca.shape  # torch.Size([784, 10])

This multiplies each of the 784 pixel values by ALL 10 weight columns at once!

Then .sum(dim=0) sums down the 784 dimension, giving 10 output values.

We've eliminated one entire loop through broadcasting!

Copied!
test_close(t1, matmul(m1,m2))  # Verify correctness
%timeit matmul(m1,m2)

Optimization 4: Einstein Summation

Einstein summation provides elegant notation for tensor operations:

Copied!
def matmul(a,b): 
    return torch.einsum('ik,kj->ij', a, b)

Decoding 'ik,kj->ij':

  • Comma separates inputs
  • Arrow separates inputs from output
  • Letters label dimensions
  • Repeated letters (k) mean multiply those dimensions
  • Letters only in input (not output) get summed over

For matrix multiplication:

  • a has dimensions i (rows) and k (columns)
  • b has dimensions k (rows) and j (columns)
  • Output has dimensions i and j
  • k appears in both inputs but NOT output → summed over

This is matrix multiplication in a single elegant expression!

Copied!
test_close(t1, matmul(m1,m2))
%timeit matmul(x_train, m2)

Performance Summary

Here's our optimization journey:

Method Time Speedup
Triple nested loop ~550ms 1x (baseline)
Numba dot product ~268µs ~2000x
Broadcasting ~137µs ~4000x
Einstein summation ~15ms ~37x
PyTorch matmul ~15ms ~37x

The journey from 550ms to microseconds shows why understanding these optimizations matters!

Key Takeaways

  1. Matrix multiplication is fundamental to ML - it's how inputs get transformed by weights in neural networks

  2. Python loops are slow - the triple nested loop works but is impractical for real applications

  3. Numba bridges Python and performance - JIT compilation can make Python code run at C speed

  4. Broadcasting is powerful and efficient - operations between different-shaped tensors without copying memory

  5. Understand the broadcasting rules - dimensions compatible if equal or one is 1, compared right-to-left

  6. Einstein summation is elegant - complex tensor operations in compact notation

  7. Build understanding step by step - implementing from scratch reveals what's really happening

Further Exploration

  • Try implementing other operations (transpose, inverse) from scratch
  • Explore how convolutions relate to matrix multiplication
  • Learn about GPU programming with CUDA
  • Study how automatic differentiation uses these operations
  • Implement mean shift clustering with broadcasting (from Lesson 12)

Understanding matrix multiplication deeply gives you intuition for how neural networks learn and why certain architectures work better than others.