Structure and Interpretation

of Tensor Programs

First Edition

A Whirlwind Tour to Deep Learning and Deep Learning Systems
Runnanochatby buildingteenygradfrom scratch: the bridge from microgradto tinygrad
with the Lean, Python, Rust, and CUDA Rust programming languages!

Made with 🖤🪻 by Jeffrey Zhang, University of Waterloo (BMath)
Made possible by Lambda Labs Research Grant

Dedication

In loving memory of my father, my teacher, and my best friend Dr. Thomas Zhang ND., R.TCMP, R.Ac.
The love I put into this book is but a fraction of the love he gave me.
May you rest in pure land. We’ll meet again dad.

And ye shall know the truth, and the truth shall make you free. John 8:32

We’ll Meet Again by Vera Lynn 1939. Cover by Johnny Cash 2002.

Presenting an early outline of SITP at Toronto School of Foundation Modeling Season 1 (November 2025)

Preface

The Structure and Interpretation of Tensor Programs

This book is aspirationally titled The Structure and Interpretation of Tensor Programs, (henceforth SITP) as it’s goal is to serve a similar role for software 2.0 as The Structure and Interpretation of Computer Programs (henceforth SICP) did for software 1.0. Written by Harold Abelson and Gerald Sussman with Julie Sussman, SICP took learners on a whimsical whirlwind tour throughout the essence of computation starting with the elements of programs with functional programming, higher order functions, data abstraction, streams, and ending with programming their own programming languages with interpreters, compilers, and register machines.

My alma matter was amongst those which took the SICP approach, and as intended, for someone coming into first year college with high school computer science, it blew my mind. After graduating college in 2022, I followed my curiosity for diving deeper into the souls of our machine by going on to developing industrial languages and runtimes.“There is only one project, architecture, operating system and languages, compiler, it’s only one project. It’s all together.” – Boris Babayan. Particularly, I hacked on languages with domain specific cloud compilers and runtimes with cloud provisioners, and cloud garbage collectors. At the end of 2022 though, when ChatGPT was released by OpenAI my mind was blown twice more. As someone programming since high school, I could not believe this at all. After two more years of hacking on cloud languages and runtimes, I started my transition from domain specific cloud compilers from GPS to Terraform to to domain specific tensor compilers from PyTorch to Triton.

The transition started with a tweet showcasing the beginnings of a tensor library evaluating the forward pass of a feed forward network from Andrej Karpathy’s Neural Networks: Zero to Hero course. While it was illuminating to start implementing each individual torch call that the nets from makemore were making, my knowledge felt quite fragmented as I forgot a lot of the foundational mathematics I saw in a single semester, and I wasn’t sure how to bridge myself to industrial deep learning systems like tinygrad, torch, jax, vllm, and sglang. Coloquially speaking, I was a neural network script kiddie.

Shortly after, I decided to take the plunge and started drinking from the firehose all the mathematical foundation I’ve since forgotten. While revisting preliminary foundation like Strang (1988), Nocedal, Wright (1999), Boyd, Lieven, Vandenberghe (2004) and reading deep learning cannon like Russel, Norvig 1995, Sutton, Barto (1992), Hastie Tibshirani (2001), Goodfellow, Bengio, Courtville (2016), Murphy (2022), the one thought I could not get out of my head was where is the SICP for software 2.0? While I found two excellent resources on building your own torch-like autograd by Tianqi Chen at Carnegie Mellon and Sasha Rush at Cornell, I personally would have enjoyed a unified resource that took me from math, to deep learning, to deep learning systems in a single unbroken sequence of thought, and perhaps others would feel similarly. That is the genesis story for this book, whose central research question is the following: What does the SICP for Deep Learning look like?

Frontispiece of Dialogue Concerning the Two Chief World Systems (Galileo Galilei 1632)

WeA modified excerpt from The Structure and Interpretation of Computer Programs Chapter 1: Building Abstractions with Procedures are about to study the idea of a computational process. Computational processes are abstract beings that inhabit computers. As they evolve, processes manipulate other abstract things called data. The evolution of a process is directed by a pattern of rules called a program parameters called a model. People create programs train models to direct processes. In effect, we conjure the spirits of the computer with our spells.

A computational process is indeed much like a sorcerer’s idea of a spirit. It cannot be seen or touched. It is not composed of matter at all. However, it is very real. It can perform intellectual work. It can answer questions. It can affect the world by disbursing money at a bank or by controlling a robot arm in a factory. The programs models we use to conjure processes are like a sorcerer’s spells. They are carefully composed recovered from symbolic numerical expressions in arcane and esoteric parallel programming languages that prescribe the tasks losses we want our processes to perform minimize.

I. Elements of Networks

Not different from learning another mathematical or programming language, they must transition from their finitely discrete structures and deterministic procedures tooling they have grown acustomed to and make the transition to the infintely continuous structures and stochastic procedures. Back then, ancient greek mathematicians were only comfortable with the finiteness of natural numbers like $1$, $2$, and $3$, and had to grapple with the infinite nature of the real numbers such as $\sqrt{2}$, $\pi$, and $e$. Similarly, the programmers of today are being asked to transition from programming algorithms of sets, maps, lists, trees, and graphs to the distributions of scalars, vectors, matrices, tensors, and neural networks.

More coloquially, programmers interested in the deep learning approach to artificial intelligence must make the transition from software 1.0software 1.0 to software 2.0software 2.0 See (Karpathy 2017), a distinction used to differentiate the classical act of programming software line by line, and the newer approach of programming software by specifying a dataset, a neural net architecture with a goal, and searching the space of programs with compute. How to exactly program with this new approach will take the remainder of the book to explain.

While software 2.0 has increased the intelligence and autonomy of our devices throughout the past decade — to name a few, language understanding with Google’s Translate and Apple’s Siri, vision understanding with Tesla Autopilot — at the end of 2022 ChatGPT was released to the world marking the beginning of software 3.0software 3.0 See (Karpathy 2025), enabling the activity of programming with none other than the English language. What may be surprising to realize is that artificial intelligence like ChatGPT is “just” another computer program. However, rather than being implemented in a language like C, Java, or Javascript, it’s implemented in one that goes by the name of PyTorch, a software 2.0 programming language centered around torch.Tensor, a multidimensional array humbly embedded within a Python package.

In this whimsical whirlwind tour dubbed The Structure and Interpretation of Tensor Programs (SITP), we will embark on a quest to build from scratch our own deep neural network like ChatGPT by implementing nanochat and our own deep learning framework like PyTorch by implementing teenygrad using Lean, Python, Rust, and CUDA Rust. After our journey together, we encourage you to modify, extend, and hack on it. Whether you’re an eager high school student, an up coming college student, or a battle-tested industry programmer, SITP has been meticulously designed so that the only prerequisite required is a basic familiarity with the elements of programming, and high school calculus. Any additional experience is helpful, not mandatory.

So with that all said, go on young hacker. Venture forth!

Table of Contents

1. Sequence Learning

$\hookleftarrow$ Table of Contents

In which we transition to the stochastic and infinitely continuous software 2.0 with ngram and linear models using probabiltiy theory, linear algebra, and calculus.

1.1 From Certain to Uncertain Knowledge

$\hookleftarrow$ Table of Contents

In order to build an artificially intelligent conversational machine like ChatGPT, we will need to follow the approach that Philosophical Investigations proposes, rather than the one Tractatus Logico-Philosophicus does. Let us now enter the forest of stochastic and infinitely continuous methods from software 2.0.
You are now ready.

In the previous chapter we saw that the disciplines of computational linguistics and natural language processing model language as a set of strings with the logical certainty of formal language theory.

In contrast, those of statisticsstatistics namely machine learningmachine learning and deep learningdeep learning model language subtly different. That is, as a weighted set of sequences with stochastic uncertainty or randomnessrandomnessThere are few terms in discourse that involve lots of philosophical baggage. This is one of them, namely, whether randomness truly exists or not. Like most treatments of stochastics, we will expand on this when we get to the post-data inverse problem of statistics and differentiate between beliefs and frequencies..

For instance, using an example sentences “Colorless green ideas sleep furiously” and “Furiously sleep ideas green colorless” from the previous chapter, we can evaluate language modellanguage models like nanochat which one has a higher chancechance of occuring?Interestingly enough, Chomsky aruged that the inability to distinguish the two sentences since they were not in the language models training data (for instance, perhaps gpt2) was evidence against such probabilistic language models. As we will see in Part II. Neural Networks, scale solves a lot of problems. We would expect that the former is more likelihoodlikely because it’s only senseless with respect to semantics, whereas the latter is senseless in grammar (which of course implies semantically senseless as well). That is,

$$p(\text{Colorless green ideas sleep furiously})>p(\text{Furiously sleep ideas green colorless})$$

Indeed, evaluating nanochat we can see it assigns the probabilityprobabilities $0.x$ and $0.y$ respectively. While evaluating probabilities for any possible sentence is useful for analytic language understanding, the predominant mode of use is to use nanochat for synthetic language generative AIgeneration. We can evaluate nanochat to produce the chance it assigns to all possible next words $w$ in language $\Omega$ conditional probabilitygiven some text. In other words, to predict the chance of next token $w$ with

$p(w|\text{Colorless green ideas sleep})$

Evaluating nanochat we see

in which nanochat samples the most likely word. So when you ask systems like ChatGPT a questionquestion, it produces an answeranswer word by word by repeating the following inference loop.

1. evaluating the probability of the next word
2. sampling such word
3. appending it to the existing text

Tip

Try completing the generation of the sentence by manually evaluating such loop by 1. clicking the predict button 2. choosing a word and 3. appending such word to the input box.

With all that said, the first task in our journey to implementing an artificially intelligent conversational machine like nanochat is to formalize our intuitions around the notions of uncertainty, chance, and likelihood just introduced. We will do so with the language of probability theoryAs opposed to alternative frameworks such as such as probabilistic logic https://en.wikipedia.org/wiki/Probabilistic_logic or uncertainty quantification.

1.2 Language as a Probability Space of Sequences

$\hookleftarrow$ Table of Contents

Because of how pervasive large language models have become over the past few years, many people are familiar with the ideas that large language models are next token predictors in the same way that digital computers speak 0s and 1s, the internet is an intergalactic highway, and the cloud is some cloud in the computer,

a bigram character-level language model adapted from karpathysyllabus: https://karpathy.ai/zero-to-hero.html, lecture: https://www.youtube.com/watch?v=PaCmpygFfXo, notebook: https://github.com/karpathy/nn-zero-to-hero/blob/master/lectures/makemore/makemore_part1_bigrams.ipynb, makemore: https://github.com/karpathy/makemore/blob/master/makemore.py#L399

dataset = open('./examples/data/names.txt', 'r').read().splitlines()
N = len(dataset)

print("--- TRAINING (counting p(w|h) with python dict ---")
# Histogram (counting frequencies) is the most precise model for training set. it *is* the training set. but it generalizes poorly.
counts_dict = {}
for di in dataset:
  di_normalized = ['<S>'] + list(di) + ['<E>']
  for h,w in zip(di_normalized, di_normalized[1:]): # in the case of bigrams h is a single character, so we can simply zip two strings to get a pair of characters
    # print(h, w)
    counts_dict[(h,w)] = counts_dict.get((h,w), 0) + 1
sorted_counts_dict = sorted(counts_dict.items(), key = lambda x: -x[1])
print("2D (w,h) histogram using python's dict:\n", sorted_counts_dict)

foobarbaz

print("--- TRAINING (counting p(w|h) with numpy ndarray (NxN) ---")
# We will now construct the same 2d histogram, but with numpy's ndarray instead of python's dict
# Because numpy's ndarray uses numerical indices to index into, we need to create a dict[str,int]
# so that when we loop over (w,h) pairs within a word we can update the count at the correct location
import numpy as np
vocab = sorted(list(set(''.join(dataset)))) # construct vocab
c2i = {c:i+1 for i,c in enumerate(vocab)}                      # construct map<char,ord>
c2i['.'] = 0                                                                            # with . as the start token and end token, to remove counting freq of (<E>*) and (*<S>) which are all 0
V = len(c2i)                                                                   # evaluate the vocab len V
C_VV = np.zeros((V,V), dtype=np.int32)            # and use V to construct C_VV

# Now we can proceed
for di in dataset:
  di_normalized = ['.'] + list(di) + ['.']
  for h,w in zip(di_normalized, di_normalized[1:]):
    print(h,w)
    h_index, w_index = c2i[h], c2i[w]                                         # use map<char, ord> to lookup the coordinate index needed for C_VV
    C_VV[h_index, w_index] += 1                                                         # update C_VV
print("2D (xt,xt-1) histogram using numpy dict:\n", C_VV)

# normalize counts C_VV to probs P_VV
C_VVf32 = (C_VV+1).astype(np.float32)             # inductive bias (locally smooth)
s_V1 = C_VVf32.sum(axis=1,keepdims=True)                # reduce along axis=1 because we want p(y|x) not p(x|y)
P_VV = C_VVf32 / s_V1                                                                   # (V, V) / (V, 1) broadcasts

# for P_VV, the elements are the counts of bigrams (h,w) accessed by indexing with ord(h) at axis=0 and ord(w) at axis=1
# now, since numpy's ndarray's are row major order, axis=0 gets printed vertically from up to down while axis=1 gets printed horizontally from left to right
i2c = {i:c for c,i in c2i.items()}  # invert map<char, ord> to map<ord, char> because looping with enumerate provides access to indices
header = '    ' + ' '.join(f'{i2c[y_index]:>4}' for y_index in range(V))
print("2D (ord, ord) histogram using numpy ndarray")
print(header)
for w_index, row in enumerate(C_VV+1):
  h = f'{i2c[w_index]:>4}'
  print(h, ' '.join(f'{count:>4}' for count in row))

foobarbaz

print("\n\n--- INFERENCE (GENERATING a name by 1. evaluating p(W=w|H=h), appending, and repeating ---")
rng = np.random.default_rng(1337)
sample_count = 10

for _ in range(sample_count):
  h, h_index = [], 0
  while True:
    # 1. evaluate p(W=w|h)
    pWcondH_V = P_VV[h_index].squeeze()

    # 2. sampling
    h_index = rng.choice(len(pWcondH_V), size=1, replace=True, p=pWcondH_V)
    sample_char = i2c[h_index.item()]

    # 3. appending the sample to history
    h.append (sample_char)
    if h_index == 0: break
  print(''.join(h))


loglikelihooddataset,n = 0.0, 0
for di in dataset:
  di_normalized = ['.'] + list(di) + ['.']
  for h,w in zip(di_normalized, di_normalized[1:]):
    w_index, h_index = c2i[h], c2i[w] # use map<char, ord> to lookup the coordinate index needed for P_VV
    pycondx = P_VV[w_index, h_index] # maximize likelihood
    logpycondx = np.log(pycondx)     # maximize loglikelihood

    loglikelihooddataset += logpycondx
    n += 1
    # print(f'{x_char}{y_char}: {pycondx:.4f} {logpycondx:.4f}')



nlldataset = -loglikelihooddataset   # minimize -loglikelihood
avgnlldataset = nlldataset / n       # minimize -1/n loglikelihood
print(f'{loglikelihooddataset=}')
print(f'{nlldataset=}')
print(f'{avgnlldataset=}')

Important

we don’t enumerate through the entire sample space nor event space, that would run into similar issues of describing a reality with too many parts to count. we need random variables.

1.3 From Variables to Random Variables

$\hookleftarrow$ Table of Contents

• definitions of rvs using calculus/analysis
• bernouilli an n=1 binomial)
• categorical (an n=1 multinomial)

1.4 Statistics as Inverse Probability

$\hookleftarrow$ Table of Contents

1.5 Iteratively Fitting Logistic Regression with Cross Entropy Loss

$\hookleftarrow$ Table of Contents

# a bigram character-level language model adapted from karpathy's zero to hero (makemore) lectures
# - syllabus: https://karpathy.ai/zero-to-hero.html
# - lecture: https://www.youtube.com/watch?v=PaCmpygFfXo
# - notebook: https://github.com/karpathy/nn-zero-to-hero/blob/master/lectures/makemore/makemore_part1_bigrams.ipynb
# - makemore: https://github.com/karpathy/makemore/blob/master/makemore.py#L399

print("\n\n--- DATA ---")
import torch
import torch.nn.functional as F
dataset = open('./examples/data/names.txt', 'r').read().splitlines()
D = len(dataset)                                                   # D is the length of the dataset. N is reserved for the number of examples, in which each di can comprise many of
vocab = sorted(list(set(''.join(dataset))))                        # construct vocab
c2i = {c:i+1 for i,c in enumerate(vocab)}                          # construct map<char,usize>
c2i['.'] = 0                                                       # with . as the start token and end token, to remove counting freq of (<E>*) and (*<S>) which are all 0
V = len(c2i)                                                       # evaluate the vocab len V

xindicesraw, yindicesraw = [], []
for di in dataset:                                                 # change to dataset[:1] when debugging to limit the number of N examples
  di_normalized = ['.'] + list(di) + ['.']                         # normalize each word di
  for x_char,y_char in zip(di_normalized, di_normalized[1:]):      # loop through the (x,y) "self-supervised" bigrams of each word di with zip
    x_index, y_index = c2i[x_char], c2i[y_char]                    # use map<char, usize> to map representation from discrete characters to discrete integers (ords)
    xindicesraw.append(x_index), yindicesraw.append(y_index)       # append

N = len(xindicesraw)
xindices_N, yindices_N = torch.tensor(xindicesraw), torch.tensor(yindicesraw)   # then, convert python lists into torch tensors
print(f'inputs (usize): {xindices_N}'); print(f'outputs (usize): {yindices_N}')

                                                                   # finally, map representation once again from discrete integers to continuous (but local) one hot vectors
x1hots_NV, y1hots_NV = F.one_hot(xindices_N,num_classes=V).float(), F.one_hot(yindices_N,num_classes=V).float()
print(f'inputs (1hot): {x1hots_NV.shape} {x1hots_NV}'); print(f'outputs (1hot): {x1hots_NV.shape} {x1hots_NV}')





print("\n\n--- ARCH ---")
g = torch.Generator().manual_seed(1337+1)                          # avgnll: 3.5 -> 4.9
W_VV = torch.randn((V, V), generator=g, requires_grad=True)        # V neurons

print("\n\n--- TRAINING ---")
K = 100
lr = 50
print(f'{D=}, {N=}, {K=}')
for k in range(K): # gradient descent
  # TRAINING FORWARD {k} --- (including softmax) with N examples from dataset D')
  logits_NV = x1hots_NV @ W_VV                                     # batch matmul print(f'logits_NV: {logits_NV.shape}, {logits_NV}');
  counts_NV = logits_NV.exp()                                      # equivalent to C_VV print(f'counts_NV: {counts_NV.shape}, {counts_NV}')
  probsycondx_NV = counts_NV / counts_NV.sum(dim=1, keepdims=True) # normalize, completing the evaluation of softmax
  # print(f'(forward) probsycondx_NV: {probsycondx_NV.shape}\n {probsycondx_NV}')

                                                                   # vectorized evaluation of loss in TEST section below
  indices_N = torch.arange(N)                                      # in order to evaluate loss of first N examples, use torch.arange(N) NOT xindices_N
  loss = -probsycondx_NV[indices_N, yindices_N].log().mean()       # however we do use yindices_N to find the corresponding likelihood predictions of the *targets*
  print(f'{k=}, {loss=}')                                          # pluck the likelihoods, then eval the log, average it, and take the inverse
                                                                   # ...to get negative log likelihood

  # TRAINING BACKWARD {k} --- (including softmax) with N examples from dataset D')
  W_VV.grad = None                                                 # same as zeroing gradients
  loss.backward()                                                  # eval f'(x)
  # print(f'(backward): {W_VV.grad.shape}, {W_VV.grad}')

  # TRAINING STEP {k} ---\n')
  W_VV.data += -lr*W_VV.grad







print("\n\n--- INFERENCE ---")





print("\n\n--- TEST ---")
i2c = {i:c for c,i in c2i.items()}                               # invert map<char, ord> to map<ord, char> for decoding

# logpD,n = 0.0, 0
# nlls = torch.zeros(N)                                            # replace sum and count tracking for average with .mean() (we won't be pushing logpycondx though only -logpycondx)
# for i in range(N):                                               # loop through the 5 input-outputs (x^(i), y^(i)) of the first word di in dataset
#   xord, yord = xindices_N[i].item(), yindices_N[i].item()        # map (x^(i), y^(i))'s index to ordinal
#   xchar, ychar = i2c[xord], i2c[yord]                            # map (x^(i), y^(i))'s ordinal to char

#   pYcondx_V = probsycondx_NV[i]
#   pycondx_ = pYcondx_V[yord]
#   logpycondx_ = torch.log(pycondx_)
#   nll = -logpycondx_
#   print(f'(x(^{i}),y(^{i})): ({xchar},{ychar}) ---> py(^{i})condx(^{i})hat: {pycondx_:.4f}, logpy(^{i})condx(^{i}): {logpycondx_:.4f}, nll: {nll:.4f}')

#   nlls[i] = nll
  # logpD += logpycondx
  # n += 1

# nllD = nlls.sum()
# avgnllD = nlls.mean()
# print(f'{nllD=}, {avgnllD=}')

1.6 Directly Fitting Linear Regression with Least Squared Error

$\hookleftarrow$ Table of Contents

1.7 Summary

$\hookleftarrow$ Table of Contents

1.8 Bibliographic Notes

$\hookleftarrow$ Table of Contents

The primary resources consulted in the writing of this chapter have been (Deisenroth, Faisal, Ong 2020), (Hastie, Tibshirani, and Friedman 2009), and (Murphy 2022); For probability theory, (Bertsekas and Tsitsiklis 2008), (Chan 2021), (Wasserman 2004), and (Durrett (2019); For linear algebra, (Strang 2026), (Strang 2019), and (Axler 2026); For optimization, (Boyd and Vandenberghe 2004) and (Kochenderfer and Wheeler 2019).

1.9 Problems

$\hookleftarrow$ Table of Contents

Intermezzo One: The Language of Logic, Sets and Functions

Dependent type theoryHarper’s forward to (Friedman and Christiansen 2018)’s https://thelittletyper.com/ (…) is a wonderfully beguiling, and astonishingly effective, unification of mathematics and programming. In type theory when you prove a theorem you are writing a program to meet a specification-and you can even run it when you are done! A proof of the fundamental theorem of arithmetic amounts to a program for factoring numbers. And it works the other way as well: every program is a proof that its specification is sensible enough to be implementable. Typе theory is a hacker’s paradise.

We begin with the latter, namely introducing the mathematical foundations that underlie the stochastic and infintely continuous needed for programming software 2.0 with the Lean programming language and proof assistant. Once we’ve introduced and formalized the basic mathematical preliminaries with Lean, we will only formally construct concepts in Intermezzo chapters after intuitively using or implementing such concept with the teenygrad deep learning framework.

If you’d like to start with the former, you are encouraged to consult Appendix A and Appendix B, which historically retrace the transitions from symbolic software 1.0 to stochastic software 2.0 and from classical to constructive mathematics respectively.

With that said, down the rabbit hole we go.

print(type(True))
print(type(42))
print(type(4.2))
print(type("42"))
print(type([42]))
print(type({"what's the meaning of life?", 42}))
print(type(dict(one="01", two="10", three="11")))

Tip

Run code snippets by clicking the play button located at the top right.

/-- Doubles a natural number. -/
def double (n : Nat) : Nat := n + n

theorem double_eq (n : Nat) : double n = 2 * n := by
  simp [double, Nat.two_mul]

#check And
#check Or
#check Or

#eval double 21

variable {α : Type*}
variable (s t u : Set α)
open Set

example (h : s ⊆ t) : s ∩ u ⊆ t ∩ u := by
  rw [subset_def, inter_def, inter_def]
  rw [subset_def] at h
  simp only [mem_setOf]
  rintro x ⟨xs, xu⟩
  exact ⟨h _ xs, xu⟩

example (h : s ⊆ t) : s ∩ u ⊆ t ∩ u := by
  simp only [subset_def, mem_inter_iff] at *
  rintro x ⟨xs, xu⟩
  exact ⟨h _ xs, xu⟩

Note

Because we are dealing with the domain of language, we will denote alphabets of symbols and strings of symbols with $w\in V$ and $\mathbf{s}=(w1,w2,\cdots ,w_t),w_i\in V$ respectively rather than $a\in \Sigma$ and $\mathbf{a}=(a1,a2,\cdots ,a_t),a_i\in \Sigma$ to denote the fact that our alphabets and strings are modeling vocabularies and sentences of words within the domain of language. The alphabet and string formalism of formal language theory can be applied to other domains that admit sequences of tokens i.e biology with protein folding.

Chapter 0 deals with formal language theory and set theory as if it were a natural language.

It introduces the “basic words” of the language, suggests how to compose “words” into “sentences,” and appeals to your knowledge of algebra for an intuitive understanding of these “sentences.” While this kind of introduction works to some extent, truly effective communication requires some formal study.

Intermezzo Two: The Language of Probability and Matrix Calculus

Fixed-Size DataA modified excerpt from How to Design Programs (Felleisen et al., 2014) Intermezzo 1: Beginning Student Language. Chapter 0 deals with BSL formal language theory and set theory as if it were a natural language. It introduces the “basic words” of the language (which in turn, models natural language), suggests how to compose “words” into “sentences,” and appeals to your knowledge of algebra sets as a collection of objects for an intuitive understanding of these “sentences.” While this kind of introduction works to some extent, truly effective communication requires some formal study.

In the previous subchapter, we mentioned that the statistical disciplines of machine learning and deep learning model language (thus the term language model) as a weighted set of sequences with stochastic uncertainty. This is done with probability theory, which models uncertaintyuncertainty in an experimentexperiment generically with two sets and a function.

The first set is the set of all possible outcomes $\omega$outcomes of an experiment and is called the sample space $\Omega$sample space, denoted by $\omega$ and $\Omega$ respectively, corresponding to the blue dots in the diagram above. In the case of modeling the phenomena of language as a random experiment, we have to select the level of granularity we care about. As in, should the set of outcomes contain elements of letters, words, sentences, paragraphs, or documents?

For language models like nanogpt, individual outcomes $\omega$ are sentences, and so $\Omega$ is the set of all sentences of length $m$ stringed into a sequence from a vocabulary of words $V$. That is, $\omega =[v1,v2,\cdots ,v_m]$.

For instance, t aking the vocabulary to be $V:=\{\text{colorless},\text{furiously},\text{green},\text{ideas},\text{sleep}\}$ and $m:=5$, enumerating the sample space $\Omega$ which is the set of all possible length 5 sentences stringed together from $V$ looks like:

$$\begin{array}{rl}\Omega :=\{& \\ {\omega }_1& =[v1=\text{colorless},v2=\text{green},v3=\text{ideas},v4=\text{sleep},v5=\text{furiously}]\\ {\omega }_2& =[v1=\text{colorless},v2=\text{green},v3=\text{ideas},v4=\text{sleep},v5=\text{colorless}]\\ {\omega }_3& =[v1=\text{colorless},v2=\text{green},v3=\text{ideas},v4=\text{sleep},v5=\text{green}]\\ {\omega }_4& =[v1=\text{colorless},v2=\text{green},v3=\text{ideas},v4=\text{sleep},v5=\text{ideas}]\\ {\omega }_5& =[v1=\text{colorless},v2=\text{green},v3=\text{ideas},v4=\text{sleep},v5=\text{sleep}]\\ & \cdots \\ {\omega }_{21}& =[v1=\text{furiously},v2=\text{sleep},v3=\text{ideas},v4=\text{green},v5=\text{colorless}]\\ {\omega }_{22}& =[v1=\text{furiously},v2=\text{sleep},v3=\text{ideas},v4=\text{green},v5=\text{green}]\\ {\omega }_{23}& =[v1=\text{furiously},v2=\text{sleep},v3=\text{ideas},v4=\text{green},v5=\text{furiously}]\\ {\omega }_{24}& =[v1=\text{furiously},v2=\text{sleep},v3=\text{ideas},v4=\text{green},v5=\text{ideas}]\\ {\omega }_{25}& =[v1=\text{furiously},v2=\text{sleep},v3=\text{ideas},v4=\text{green},v5=\text{sleep}]\\ \}& \end{array}$$

where each ${\omega }_i$ corresponds to a dot within Figure 1.x’s blue sample space.

Caution

If we wanted a probability model so that we could ask the chance or likelihood of a single word occuring, we could shorten the length of our sentences to a single word with $m:=1$ while maintaining the same vocabulary $V:=\{\text{colorless},\text{furiously},\text{green},\text{ideas},\text{sleep}\}$, but then the outcomes $\omega \in \Omega$ are reduced to the vocabulary $v\in V$ and if we wanted a model that is expressive enough to span the entirety of a natural language like English, we would need to

“how is v_m an event?” is answered by: it isn’t; v_m is a value; the notation smuggles the value into the event {W_m = v_m} through the random variable

the formal language theoretic vocabulary is now seen as the sample spaceIn the case of language modeling, it’s somewhat straightforward how to construct the sample space. Namely, each element is a word from the vocabulary. However in other experiments, constructing spaces with the correct amount of detail, fidelity or granularity is somewhat of an art. For instance, you wouldn’t want to construct a sample space $\Omega :=\{\text{win},\text{lose}\}$ to model the chance of winning the lottery.:

The second set is the set of all possible subsets of the sample space $\Omega$, corresponding to all possible subsets drawn with colored rectangles. This is the set of all possible events of the experiment called the event space $\mathcal{F}$event space, denoted by $A$ and $\mathcal{F}$ respectively. That is, $A\subseteq \Omega$. For instance, two possible events $A,B\in \mathcal{F}$ are “all words that start with the letter ‘a’” and “all words that start with the letter ‘b’”:

$$\begin{array}{rl}A& :=\{\text{aardvark}\}\\ B& :=\{\text{banana, bongo}\}\end{array}$$

where the first event $A$ has a single outcome $\omega \in \Omega$ whereas the second event $B$ has two. Finally, the function is the probability law $\mathbb{P}:\mathcal{F}\to [0,1]$probability law with type $\mathbb{P}:\mathcal{F}\to \mathbb{R}$, mapping any event in the event space to a number on the real number line.

$$\begin{array}{rl}\mathbb{P}(A)& :=0.00001\\ \mathbb{P}(B)& :=0.00001\\ & \cdots \end{array}$$

and so on, for all events $E\in \mathcal{F}$. Note that we are using $:=$ to denote that the mapping above is one possible assignment $\mathbb{P}$. For instance, two different language models can have two different assignments. That is, probability is the weight, mass, or more generally the measure of a set relative to the sample space.

Note

Curious why $\mathbb{P}$ is defined on subsets $E\subset \Omega$ rather than elements $\omega \in \Omega$? To learn more, consult Intermezzo Two: The Language of Probability and Matrix Calculus.

Together, this triplet of the sample space, event space and probability law comprise a probability space $(\Omega ,\mathcal{F},\mathbb{P})$probability space $(\Omega ,\mathcal{F},\mathbb{P})$, in which $\mathbb{P}$ must satisfy three axioms:

1. non-negativity: $\mathbb{P}(\mathbb{E})\ge 0$
2. normalization: $\mathbb{P}(\Omega )=1$
3. additivity: $A\cap B=\varnothing ⟹\mathbb{P}(A\cup B)=\mathbb{P}(A)+\mathbb{P}(B)$

The first two axioms are quite easy to accept. Namely, non-negativity ensures all probabilities are positive, and normalization ensures that the probability of something happening is 1. This value can be any value in which everyone agrees that it represents near certainty such as 10, or 100. In practice it is always 1. The third axiom of additivity ensures we can evaluate the probability of disjoint events by summing the measure of individual events. For instance, because the event $B$ which contains all words that start with the letter “b” is composed of disjoint events (corresponding to the left figure below) the axiom of additivity applies:

$$\begin{array}{rlrl}\mathbb{P}(\{\text{banana},\text{bongo}\})& =\mathbb{P}(\{\text{banana}\}\cup \{\text{bongo}\})& & \\ & =\mathbb{P}(\{\text{banana}\})+\mathbb{P}(\{\text{bongo}\})& & \text{[by additivity]}\\ & =0.1+0.2& & \\ & =0.3& & \end{array}$$

However the event $C$ which contains all foobarbaz is composed of overlapping events (corresponding to the right figure above), so the axiom of addivitiy does not apply. We can evaluate the union of non-disjoint (overlapping) events with the following corollary, known as the sum rulesum rule:

$$\mathbb{P}(A\cup B)=\mathbb{P}(A)+\mathbb{P}(B)-\mathbb{P}(A\cap B)$$

a conditional probability $\mathbb{P}(A\mid B)$conditional probability is definedFor the mathematically inclined, conditional probability must be defined rather than deduced with Kolmogorov’s axiomatization of probability as sets becuase it concentrates on the notion of additive measure. For a separate axiomatic approach in treating conditional probability as a first-class citizen, see Conditional Measures and Applications (Rao 1993). as

$$\mathbb{P}(A\mid B)\stackrel{△}{=}\frac{\mathbb{P}(A\cap B)}{\mathbb{P}(B)}$$

since probability is the weight, mass, or more generally the measure of a set relative to the sample space. not defined P(B)=0 (observed near uncertainty) false positive/negative?

the product rule $\mathbb{P}(A\cap B)$product rule is

$$\mathbb{P}(A\cap B)=\mathbb{P}(B)\mathbb{P}(A\mid B)$$

two events $A$ and $B$ are independenceindependent when

$$\mathbb{P}(A\cap B)=\mathbb{P}(A)\mathbb{P}(B)$$

the bayes’ rulebayes’ rule when Finally, in other scenarios where $P(A\mid B)$ is given but $P(B\mid A)$ is not, this is amenable by rewriting the numerator and denominator of conditional probability using the product rule and sum rule:

$$\begin{array}{rlrl}P(A\mid B)& =\frac{P(A\cap B)}{P(B)}& & \text{[by definition of conditional]}\\ & =\frac{P(A)P(B\mid A)}{P(B)}& & \text{[by product rule on joint]}\end{array}$$

$$\underset{\text{posterior}}{\underbrace{p(x\mid y)}}=\frac{\stackrel{\text{prior}}{\overbrace{p(x)}}\stackrel{\text{likelihood}}{\overbrace{p(y\mid x)}}}{\underset{\text{evidence}}{\underbrace{\int p(y\mid x)p(x)dx}}}.$$

and is referred to as bayes’ rule. The scenario where $P(A\mid B)$ is accessible but $P(B\mid A)$ is not comes up many times practice, where $A$ is some latent event and $B$ is an observable event. This happens in science (why bayes rule is also referred to as the logic of science) where $P(A\mid B)$ is used to model hypothesis and evidence where the posterior $P(H\mid E)$ (latent given observed) is updated by evaluating the product of the prior $P(H)$ and the likelihood (observed given latent) $P(E\mid H)$ normalized by the evidence $P(E)$.

two law of total probabilitylaw of total probability when

sam’s broken thermometer (page 293)

2. From IPL’s Array to APL’s Multidimensional Array

$\hookleftarrow$ Table of Contents

2.1 From Virtual to Physical Machines (and Shapes)

$\hookleftarrow$ Table of Contents

• justify native for eager performance
• pyo3 https://github.com/j4orz/ateenysitp/blob/master/ARCHITECTURE.md#level-1-teenygrads-build-configuration-and-development-environment
• native components using cpython as encapsulation boundary
• freethreaded python eliminates multithreading problem

/// SGEMM with the classic BLAS signature (row-major, no transposes):
/// C = alpha * A * B + beta * C
fn sgemm(
  m: usize, n: usize, k: usize,
  alpha: f32, a: &[f32], lda: usize,
  b: &[f32], ldb: usize,
  beta: f32, c: &mut [f32], ldc: usize) {
  assert!(m > 0 && n > 0 && k > 0, "mat dims must be non-zero");
  assert!(lda >= k && a.len() >= m * lda);
  assert!(ldb >= n && b.len() >= k * ldb);
  assert!(ldc >= n && c.len() >= m * ldc);

  for i in 0..m {
    for j in 0..n {
      let mut acc = 0.0f32;
      for p in 0..k { acc += a[i * lda + p] * b[p * ldb + j]; }
      let idx = i * ldc + j;
      c[idx] = alpha * acc + beta * c[idx];
    }
  }
}

fn main() {
  use std::time::Instant;

  for &n in &[16usize, 32, 64, 128, 256] {
    let (m, k) = (n, n);
    let (a, b, mut c) = (vec![1.0f32; m * k], vec![1.0f32; k * n], vec![0.0f32; m * n]);

    let t0 = Instant::now();
    sgemm(m, n, k, 1.0, &a, k, &b, n, 0.0, &mut c, n);
    let secs = t0.elapsed().as_secs_f64().max(std::f64::MIN_POSITIVE);
    let gflop = 2.0 * (m as f64) * (n as f64) * (k as f64) / 1e9;
    let gflops = gflop / secs;

    println!("m=n=k={n:4} | {:7.3} ms | {:6.2} GFLOP/s", secs * 1e3, gflops);
  }
}

Thus far in our journey we’ve successfully made the transition from programming software 1.0 which involves specifying of discrete data structures — like sets, associations, lists, trees, graphs, — with the determinism of types to programming software 2.0 which involves reovering continous data structures — like vectors, matrices, and tensors — with the stochasticity of probabilities. If you showed the mathematical equations for your models, loss functions, and optimizers to our mathematical cousinsIn particular, the physicists whom realized describing reality as minimizing free-energy was fruitful, such as Helmholtz with energy, Gibbs with enthalpy, and Boltzmann with entropy. who also made the same transition, they would understand the mathematical equations and even algorithmsas until recently, algorithms were understood to be a sequence of instructions to be carried out by a human. i.e Egyptian Multiplication, Euclid’s Greatest Common Divisor. But they wouldn’t understand how the code we’ve programmed thus far with Python and Numpy is being automatically calculated by the mechanical machine we call computers. For that we need to transition from users of the numpy framework to implementors our own framework, which we’ll call teenygrad. This requires us moving from the beautiful abstraction of mathematics to the assembly heart and soul of our machines. This book uses RustTo backtrack to familiarize yourself with Rust, take a look at the experimental version of The Rust Programming Language, by from Will Crichton and Shriram Krishnamurthi, originally written by Steve Klabnik, Carol Nichols, and Chris Krycho., but feel free to follow along with C/C++In that case, take a look at “The C Programming Language by Brian Kernighan and Dennis Ritchie. — what’s fundamental to the purpose of accelerating said basic linear algebra subroutines on the CPU is is choosing a language that 1. forgoes the dynamic memory management overhead of garbage collection and 2. has a compiled implementation which allows us to analyze and tune the actual instructions which an underlying processor executes. Using Rust will allow us to start uncovering what is really happening under the hood of accelerated Tensors To sidetrack to familiarize yourself with systems programming in the context of software 1.0 , read the classic of Computer Systems A Programmer’s Perspectivewith instruction set architectures, microarchitecture, memory hierarchies

• instruction set architecture
• computation: control path, data path
• communication: memory, input/output

In the last chapter we developed teenygrad.Tensor by virtualizing physical 1D storage buffers into arbitrarily-shaped logical ND arrays by specifying an iteration space with strides, and started our journey in accelerating the basic linear algebra subroutines defined on the Tensor with Rust, which speeds up the throughput performance of SGEMM from 10MFLOP/S to 1GFLOP/S. Simply executing a processor’s native code — referred to as assembly code — rather than Python’s bytecode results in an increase of two orders of magnitude.

Now that the code that is being executed is native assembly, we can dive deeper into the architecture of the machine in order to reason and improve the performance of teenygrad’s BLAS.

2.2 Accelerating the Communication of Hierarchies

Loop Reordering, Register and Cache Blocking

$\hookleftarrow$ Table of Contents

2.3 Accelerating the Computation of Pipelines

$\hookleftarrow$ Table of Contents

Instruction Level Parallelism via Loop Unrolling

2.4 From Abstract to Numerical Linear Algebra

2.6 Summary

One quick way to summarize the milestones in high performance computing, compilers and architecture is to list the Turing Award winners: Alan Perlis (1966) for his influence on advanced programming techniques and compiler construction, including his role in designing ALGOL and establishing the discipline of programming languages as a field; John Backus (1977) for designing FORTRAN — the first high-level language to achieve widespread practical adoption — and formalizing language syntax through Backus-Naur Form; Tony Hoare (1980) for axiomatic semantics, giving programmers a formal logical framework for reasoning about program correctness; Niklaus Wirth (1984) for designing a sequence of clean, teachable languages — EULER, ALGOL-W, Pascal, and Modula — that shaped how programming languages are structured and implemented; John Cocke (1987) for pioneering optimizing compilers and the Reduced Instruction Set Computer (RISC) architecture, showing that simpler instruction sets allow faster hardware; William Kahan (1989) for fundamental contributions to numerical analysis, most consequentially the IEEE 754 floating-point standard that made reliable numerical computation reproducible across hardware; Frederick Brooks (1999) for landmark contributions to computer architecture — most notably the IBM System/360 — and for articulating the enduring lessons of large-scale software engineering; Frances Allen (2006) for foundational contributions to the theory and practice of optimizing compilers, including dataflow analysis and the program dependence graph; John Hennessy and David Patterson (2017) for a systematic, quantitative approach to designing and evaluating computer architectures, whose RISC principles underpin billions of processors and the open RISC-V standard; and Alfred Aho and Jeffrey Ullman (2020) for foundational contributions to programming language theory and compiler construction, most durably codified in the Dragon Book; and finally,Jack Dongarra (2021) for pioneering the numerical libraries — BLAS, LAPACK, and MPI — that became the substrate of high-performance scientific computing and modern deep learning accelerators;

2.7 Bibliographic Notes

2.8 Problems

rough

After SITP’s first two chapters, you now understand how to write software 2.0 programs which learn with numpy, and it’s core multidimensional np.ndarray data structure. For data, instead of using locally discrete representations, software 2.0 programs use dimensionally stochastic representations. — this is why in Chapter 1 you learned the technique of 1. Representing Data with High Dimensional Stochasticity in numpy And for functions, instead of constructing maps with explicit instructions, software 2.0 programs recover such maps with differential optimization — this is why in Chapter 2 you learned the technique of 2. Learning Functions from Data with Parameter Estimation in numpy.

Throughout the first two chapters, a question you may have been asking yourself is exactly how numpy as a library is implemented. While there are many aspects to numpy (sidenote?), the primary components we care about is the core np.ndarray data structure and it’s acceleration of basic linear algebra subroutines through the usage of low-level software known as the BLAS. In this chapter, we will finally understand how numpy works under the hood by implementing our very own np.ndarray with teenygrad.Tensor, and it’s subset of the BLAS with teenygrad.eagkers. In Part II. Neural Networks’s Chapter 5. Accelerating Sequence Models on GPU in teenygrad, you will evolve teenygrad’s implementation by adding neural network primitives, gradient-based optimization, and gpu acceleration with cuda cores in order to support the various inductive biases of deep neural networks explored in the “era of research”, which roughly spans from 2012-2020.

However, before we begin with teenygrad’s implementation, there is another transition we must make moving from programming machine learning algorithms to programming machine learning systems themselves. Namely, the understanding that comes with the shift from abstract linear algebra to numerical linear algebra, which is that the “linear” in computational linear algebra is moreso aspirational. The first place to understand that is in the transition from real number arithmetic to floating point arithmetic.

Although the discipline of machine learning heavily relies on the language of linear algebra, as you now know from the previous two chapters, the primary vector spaces that are used are d-dimensional coordinate spaces ${\mathbb{R}}^d$ in which we can construct random vectors and their high dimensional joint distributions to represent data. Recall the set-theoretic definitions of ${\mathbb{R}}^2$, ${\mathbb{R}}^3$, and so on, up until ${\mathbb{R}}^d$. Because these definition are constructive, they translate quite easily to code:

from typing import Self

class Rd():
  def __init__(self: Self, data: list[T]) -> None:
    self.data = data

  def __add__(self: Self, other: Self) -> Self:
    return Rd([xi + yi for xi, yi in zip(self.data, other.data)])

But then that begs the question, what is the definition of the set $\mathbb{R}$?

class Real():
  def __init__(self: Self) -> None: raise NotImplementedError("todo")
  def __add__(self: Self, other: Self) -> Self: raise NotImplementedError("todo")

While we now have a correct implementation of numpy’s multidimensional np.ndarray with teenygrad.Tensor, there is still exists a large difference between the two, primarily boiling down to performance. Let’s take the example of solving a linear system $A\mathbf{x}=b$ where $A\in {\mathbb{R}}^{n\times m}$, $x\in {\mathbb{R}}^n$, and $b\in {\mathbb{R}}^m$, with the the number of equations (the size of the dataset) $n$ and the number of unknowns (the number of dimensions) $m$ being large. (todo, or linear least squares?):

import numpy as np

import teenygrad as tg

The running time with np.ndarray is X, whereas with tg.Tensor it’s Y. This is because although you are programming numpy with Python, underneath the hood it’s core basic linear algebra subroutines are accelerated using native software, from a library called (for lack of a better term), the basic linear algebra subroutines* (abbreviated here on in as the BLAS). The BLAS has a longstanding history dating back

That is, programs whose instructions are executed by the physical machine, rather than a virtual machine. In order to accelerate tg.Tensor, we must dive deeper in gaining a mechanical empathy for the soul of the machine. This is where programmers differ from mathematical cousins — not only do we concern ourselves about correctness, but also that of performance — making us

You have to be kind to the car computer. You feel the poor thing groaning underneath you. If you’re going to push a piece of machinery to the limit, and expect it to hold together, you have to have some sense of where that limit is. Out there, is the perfect lap speed of light. No mistakes: every gear change load, every corner store. Perfect. Do you see it? Most people don’t.

To build SITP’s capstone framework project of teenygrad, you will need to accelerate the training and inference of the capstone model project nanochat with many core parallel processors known as GPUs in Part II’s Chapter 5. Accelerating Sequence Models on GPU in teenygrad with CUDA Rust and PTX. To prepare you for such massively parallel programming, in the remainder of Chaoter 3 you will accelerate the basic linear algebra subroutines of teenygrad’s tg.Tensor in order deeply understand the following question:

Q: How do high performance software libraries accelerate their speeds?
A: By accelerating the computation of pipelines and the communication of hierarchies.

Historically speaking, the way in which biological humans first started communicating with digital computers was by writing programs directly in a machine language.

machine programming language: assemblers high level programming languages: fortran, cobol, algol systems programming language: C (only survivor) virtual machine and jits: euluer/pascal, java, javascript, python

The core idea behind the borrow checker is that variables have three kinds of permissions on their data:

• Read (R): data can be copied to another location.
• Write (W): data can be mutated in-place.
• Own (O): data can be moved or dropped.

shapes

from typing import Self
import array, math
import teenygrad

class InterpretedTensor:
  @classmethod
  def arange(cls, end: int, requires_grad: bool=False) -> Self: return InterpretedTensor((end,), list(range(end)), requires_grad=requires_grad)
  @classmethod
  def zeros(cls, shape: tuple[int, ...]) -> Self:
    numel = math.prod(shape)
    tensor = InterpretedTensor((numel,), [0.0]*numel).reshape(shape)
    return tensor
  @classmethod
  def ones(cls, shape: tuple[int, ...]) -> Self:
    numel = math.prod(shape)
    tensor = InterpretedTensor((numel,), [1.0]*numel).reshape(shape)
    return tensor
  
  def __init__(self, shape: tuple[int, ...], storage: list[float], inputs: tuple[Self, ...]=(), requires_grad: bool=False) -> None:
    self.shape: tuple[int, ...] = shape
    self.stride: tuple[int, ...] = [math.prod(shape[i+1:]) for i in range(len(shape))] # row major, and math.prod([]) produces 1
    self.storage: list[float] = storage

    self.inputs: tuple[Self, ...] = inputs
    self._backward = lambda: None # callers override after init with self captured in closure
    self.grad: InterpretedTensor = InterpretedTensor.zeros(shape) if requires_grad else None # python can recursively type (no need for Box<_>) bc everything is a heap-allocated reference
  @property
  def numel(self): return math.prod(self.shape) # np (and thus jax) call this .size
  @property
  def ndim(self): return len(self.shape)
  @property
  def T(self) -> Self:
    assert self.ndim == 2
    m, n = self.shape
    t = InterpretedTensor((n, m), self.storage)
    t.stride = [self.stride[1], self.stride[0]]
    return t

  def reshape(self, shape: tuple[int, ...]) -> Self:
    self.shape = shape
    self.stride = [math.prod(shape[i+1:]) for i in range(len(shape))] # math.prod([]) produces 1
    return self
  
  def __repr__(self) -> str:
    return f"InterpretedTensor({self.chunk(self.storage, self.shape)})"
  @staticmethod
  def chunk(flat, shape):
    if len(shape) == 1: return flat[:shape[0]]
    size = len(flat) // shape[0]
    return [InterpretedTensor.chunk(flat[i*size:(i+1)*size], shape[1:]) for i in range(shape[0])]
  
  # backward f'(x)
  @staticmethod
  def topo(node: InterpretedTensor, seen: set[InterpretedTensor], output: list[InterpretedTensor]) -> None:
    if node in seen: return
    seen.add(node)
    for input in node.inputs: InterpretedTensor.topo(input, seen, output)
    output.append(node)

  def backward(self) -> None:
    seen, topologically_sorted_expression_graph = set(), []
    InterpretedTensor.topo(self, seen, topologically_sorted_expression_graph)

    self.grad = InterpretedTensor.ones(self.shape) # base case
    for tensor in reversed(topologically_sorted_expression_graph): tensor._backward()
  
  # forwards f(x)
  def __radd__(self, other: Self) -> Self: return self.__add__(other)
  def __add__(self, other: Self) -> Self:
    n, alpha = self.numel, 1
    x, y, z = array.array('f', self.storage), array.array('f', other.storage), array.array('f', [0.0]*(n))
    teenygrad.eagkers.cpu.saxpy(n, alpha, x, y) # y=axpy
    requires_grad = self.grad is not None or other.grad is not None
    output_tensor = InterpretedTensor(self.shape, list(y), (self, other), requires_grad=requires_grad)
    def _backward():
      self.grad += output_tensor.grad
      other.grad += output_tensor.grad
    output_tensor._backward = _backward
    return output_tensor

  def __rmul__(self, other: Self) -> Self: return  self.__mul__(other)
  def __mul__(self, other: Self) -> Self:
    n = self.numel
    x, y, z = array.array('f', self.storage), array.array('f', other.storage), array.array('f', [0.0]*n)
    teenygrad.eagkers.cpu.smul(n, x, y, z)
    requires_grad = self.grad is not None or other.grad is not None
    output_tensor = InterpretedTensor(self.shape, list(z), (self, other), requires_grad=requires_grad)
    def _backward():
      self.grad += output_tensor.grad * other
      other.grad += output_tensor.grad * self
    output_tensor._backward = _backward
    return output_tensor

  def __neg__(self) -> Self:
    n = self.numel
    x, y = array.array('f', self.storage), array.array('f', [0.0]*n)
    teenygrad.eagkers.cpu.saxpy(n, -1, x, y)
    requires_grad = self.grad is not None
    output_tensor = InterpretedTensor(self.shape, list(y), (self,), requires_grad=requires_grad)
    def _backward():
      self.grad += -output_tensor.grad
    output_tensor._backward = _backward
    return output_tensor

  def __sub__(self, other: Self) -> Self:
    n = self.numel
    x, y = array.array('f', other.storage), array.array('f', self.storage)
    teenygrad.eagkers.cpu.saxpy(n, -1, x, y)
    requires_grad = self.grad is not None or other.grad is not None
    output_tensor = InterpretedTensor(self.shape, list(y), (self, other), requires_grad=requires_grad)
    def _backward():
      self.grad += output_tensor.grad
      other.grad += -output_tensor.grad
    output_tensor._backward = _backward
    return output_tensor

  def tanh(self) -> Self:
    n = self.numel
    x, y = array.array('f', self.storage), array.array('f', [0.0]*n)
    teenygrad.eagkers.cpu.stanh(n, x, y)
    requires_grad = self.grad is not None
    output_tensor = InterpretedTensor(self.shape, list(y), (self,), requires_grad=requires_grad)
    def _backward():
      self.grad += output_tensor.grad * (InterpretedTensor.ones(self.shape) - output_tensor * output_tensor) # f(x) = tanh(x) ==> f'(x) = 1 - tanh(x)^2
    output_tensor._backward = _backward
    return output_tensor

  def __rmatmul__(self, other: Self) -> Self: return other.__matmul__(self) # GEMM does not commute: AB != BA
  def __matmul__(self, other: Self) -> Self:

Intermezzo Three: The Language of Numerical Analysis

The Creation of Adam, Michelangelo 1508-1512.

II. Neural Networks

In part one of The Structure and Intepretation of Tensor Programs you have developed a solid foundation in the mathematical preliminaries and statistical models used throughout the machine learning approach to artificial intelligence. It’s amazing how close you are to of deep learning without you even knowing. By the end of Part II. Neural Networks, we will be one step closer in achieving our quest of building our own ChatGPT by reproducing GPT2Presented in Language Models are Unsupervised Multitask Learners (Raford et al. 2019), following Andrej Karpathy’s nanogpt. But before we get there, there is some more work for us to do.

So in Chapter 4. Learning Sequences via Deep Neural Networks with teenygrad, you will increase the expressivity of the linear models implemented in Part I by non-linearities to get a class of models known as deep neural networks. In order to build ourselves up to nanogpt, we will hold the training goal of learning sequences constant, and incrementally implement more expressive neural networks architectures following the nets in Andrej Karpathy’s makemore, starting from feedforward neural networks (FNNs), to convolutional neural networks (CNNs), to recurrent neural networks (RNNs), and finally, transformer neural networks (GPTs). These various neural network architectures implement different inductive biases, which were all explored during the 2012-2019 time period of what is coloqially known as the age of research.

Then, in Chapter 5. Accelerating Sequence Models on GPU in teenygrad, you will evolve teenygrad from a numerical linear algebra library implemented in Part I to a full blown batteries-included deep learning framework like PyTorch. This means implementing the optimizers for neural networks whose evaluations are accelerated with massively parallel processors and whose gradients are automatically evaluated with an automatic differentiation engine.

After completing part two, you will be ready for Part III. Scaling Networks of the book where we finally achieve our quest of building our own ChatGPT. Part three follows the 2020-2025 time period of what is colloquially as the age of scaling where researchers focused on scaling up the generality of generative pretrained transformers by adding assistant-like behavior in a midtraining phase with reinforcement learning with human feedbackOriginally presented in Introducing ChatGPT (OpenAI 2022), and reproduced by open source in Llama 2: Open Foundation and Fine-Tuned Chat Models (Touvron et al. 2023) and by adding reasoning-like behavior in a posttraining phase with foobarbazOriginally presented in Introducing OpenAI o1 (OpenAI 2024), and reproduced by open source in DeepSeek-R1: Incentivizing Reasoning Capability in LLMs via Reinforcement Learning

II. Neural Networks

• 3. Automated Sequence Learning via Deep Neural Networks with teenygrad
  • 3.1 From Supervised to Self-Supervised Learning with Sequences
  • 3.2 Learning Sequences with Linear Models
  • 3.3 From Linear to Non-Linear Learning with Deep Neural Networks
  • 3.4 Learning Sequences with Feedforward Neural Networks (FNNs)
  • 3.5 Learning Sequences with Convolutional Neural Networks (CNNs)
  • 3.6 Learning Sequences with Recurrent Neural Networks (RNNs)
  • 3.7 Learning Sequences with Transformer Neural Networks (GPTs)
• 4. Accelerating Deep Neural Networks on GPUs in teenygrad
  • 4.1 From Numerical Linear Algebra Libraries to Deep Learning Frameworks
  • 4.2 Network Primitives with teenygrad.nn
  • 4.3 Automatic Differentiation with Tensor.forward() and Tensor.backward()
  • 4.4 Gradient-Based Optimization with teenygrad.optim
  • 4.4 From SIMD to SIMT: From Multi Core to Many Core Processors
  • 4.5 Accelerating GEMV on GPU with CUDA Rust and PTX via Rooflines
  • 4.6 Accelerating GEMM on GPU with Data Reuse
  • 4.7 Accelerating GEMM on GPU with Scheduling
• 6. Speedrunning GPT2 with Advanced Linear Algebra
  • 6.1 From Abstract to Numerical Linear Algebra

4. Learning Sequences from Data with Deep Neural Networks in torch

$\hookleftarrow$ Table of Contents

(todo, some explanation) on learning sequences

4.1 From Supervised to Self-Supervised Learning with Sequences

4.2 Learning Sequences with Linear Models

4.3 From Linear to Non-Linear Learning with Deep Neural Networks

$\hookleftarrow$ Table of Contents

Recall that the task of house price prediction and sentiment classification which can be modelled by functions of the form $f:{\mathbb{R}}^d\to \mathbb{R}$ and $f:{\mathbb{R}}^d\to [0,1]$ respectively. The simplest inductive bias was made, in a which a linear relationship was assumed to hold between the input and output spaces, and where the output was subsequently modeled as an inner product $<,>:{\mathbb{R}}^d\times {\mathbb{R}}^d\to \mathbb{R}$ between an input vector and a weight vector. For the case of regression, we have $y=f(\mathbf{x};\theta ):={\theta }^{\top }\mathbf{x}$, and for the case of classification, we have $y=f(\mathbf{x};\theta ):=\sigma ({\theta }^{\top }\mathbf{x})$, todo:glm,exp. The key entry point into the function class of deep neural networks is that of logistic regression, because the log odds produced by the inner product (which are indeed affine) required a mapping into a valid probability via sigmoid function $\sigma :\mathbb{R}\to [0,1]$, where $\sigma (x):=\frac{1}{1+e(-x)}$ which is in fact not linear nor affine.

The next natural question to ask then is whether the logistic regression model is considered a deep neural network? The answer is that technically yes, it can be considered a degenerative deep neural network with $0$ hidden layersIn the same way that a list can be considered a degenerative binary tree or graph. These so-called hidden layers automate the construction of the representation through learning, so that the model not only discovers the mapping from representation to output, but also the representation itselfIn the same way that for certain computations the positional representation of arabic numerals are more suitable compared to that of roman numerals, and polar coordinates over cartesian coordinates. See [A Representational Analysis of Numeration Systems (Zhang, Norman 1995)](A representational analysis of numeration systems Author links open overlay panel). The functions of deep neural networks will take the form of $f(\mathbf{x}):=g\circ \varphi (\mathbf{x})$ with $g$ being a linear classifier on feature extractor $\varphi :{\mathbb{R}}^{d_1}\to {\mathbb{R}}^{d_n}$, $\varphi (\mathbf{x}):=h^{[l]}\circ h^{[l-1]}\circ \cdots h^{[1]}\circ h^{[0]}(\mathbf{x})$ where $l$ is the number of compositional layers, and each intermediate function $h^{[i]}$ has the form of $h^{[i]}:{\mathbb{R}}^{d_i}\to {\mathbb{R}}^{d_{i+1}}$. Each hidden layer successively and graduallyIn the same way of Grothendieck’s preferred style of mathematics described in Récoltes et Semailles: I can illustrate the second approach with the same image of a nut to be opened. The first analogy that came to my mind is of immersing the nut in some softening liquid, and why not simply water? From time to time you rub so the liquid penetrates better, and otherwise you let time pass. The shell becomes more flexible through weeks and months—when the time is ripe, a touch of the hand is enough, and the shell opens like a perfectly ripened avocado! A different image came to me a few weeks ago. The unknown thing to be known appeared to me as some stretch of earth or hard marl, resisting penetration. One can go at it with pickaxes or crowbars or even jackhammers: this is the first approach, that of the “chisel” (with or without a hammer). The other is the sea. The sea advances insensibly and in silence, nothing seems to happen, nothing moves, the water is so far off you hardly hear it… yet it finally surrounds the resistant substance. lifts the complexity and abstraction of the data’s representationChris Olah, cofounder of Anthropic and the lead of it’s interpretability research wrote an excellent article on how software 2.0’s representation learning loosely correspond to software 1.0’s types in Neural Networks, Types, and Functional Programming.

Together, these two aspects of learning non-linear, representations form the essence of deep learning.

Let’s now turn out attention to the function bodies of these $h^{[i]}$’s with a deep neural netork of $1$ hidden layer, carrying out the task of price regression and sentiment classification so that $f$ has the form $f:{\mathbb{R}}^d\to \mathbb{R}$. With the statistical learning foundation from part one, we will simply present the forward pass $f(\mathbf{x})$, the loss function $\mathcal{L}(d,\theta )$, and the backward pass ${\ell }^{\prime }(y,\hat{y},\theta )$.

Forward Pass

The functions of deep neural networks will take the form of $f(\mathbf{x}):=g\circ \varphi (\mathbf{x})$ with $g$ being a linear classifier on feature extractor

$$\begin{array}{rl}\varphi & :{\mathbb{R}}^{d_1}\to {\mathbb{R}}^{d_n}\\ \varphi (\mathbf{x})& :=h^{[l]}\circ h^{[l-1]}\circ \cdots h^{[1]}\circ h^{[0]}(\mathbf{x})\end{array}$$ where $l$ is the number of compositional layers, and each intermediate function $h^{[i]}$ has the form of $h^{[i]}:{\mathbb{R}}^{d_i}\to {\mathbb{R}}^{d_{i+1}}$.

Loss Function

Backward Pass

TODO

• figure/diagram/lecun circuits ->algebraic/symbolic equations->torch code
  • first train net for price regression and classification
    • intuition of automating feature engineering
  • XOR: playground.tensorflow
  • change code below to XOR
  • mention that 2.1 treats backward pass as a black box, which is the magic of the abstraction
  • mention to readers that if they want to, they can read 2.1.3, 2.2, and then back to 2.1.4

TODO: generate prose/exposition by repaging lecture/text in ram

class MLP(nn.Module):
  """
  takes the previous block_size tokens, encodes them with a lookup table,
  concatenates the vectors and predicts the next token with an MLP.

  Reference:
  Bengio et al. 2003 https://www.jmlr.org/papers/volume3/bengio03a/bengio03a.pdf
  """

  def __init__(self, config):
    super().__init__()
    self.block_size = config.block_size
    self.vocab_size = config.vocab_size
    self.wte = nn.Embedding(config.vocab_size + 1, config.n_embd) # token embeddings table
    # +1 in the line above for a special <BLANK> token that gets inserted if encoding a token
    # before the beginning of the input sequence
    self.mlp = nn.Sequential(
      nn.Linear(self.block_size * config.n_embd, config.n_embd2),
      nn.Tanh(),
      nn.Linear(config.n_embd2, self.vocab_size)
    )

  def get_block_size(self):
    return self.block_size

  def forward(self, idx, targets=None):
    # gather the word embeddings of the previous 3 words
    embs = []
    for k in range(self.block_size):
      tok_emb = self.wte(idx) # token embeddings of shape (b, t, n_embd)
      idx = torch.roll(idx, 1, 1)
      idx[:, 0] = self.vocab_size # special <BLANK> token
      embs.append(tok_emb)

    # concat all of the embeddings together and pass through an MLP
    x = torch.cat(embs, -1) # (b, t, n_embd * block_size)
    logits = self.mlp(x)

    # if we are given some desired targets also calculate the loss
    loss = None
    if targets is not None:
      loss = F.cross_entropy(logits.view(-1, logits.size(-1)), targets.view(-1), ignore_index=-1)

    return logits, loss

4.4 Learning Sequences with Feedforward Neural Networks

$\hookleftarrow$ Table of Contents

class MLP(nn.Module):
  """
  takes the previous block_size tokens, encodes them with a lookup table,
  concatenates the vectors and predicts the next token with an MLP.

  Reference:
  Bengio et al. 2003 https://www.jmlr.org/papers/volume3/bengio03a/bengio03a.pdf
  """

  def __init__(self, config):
    super().__init__()
    self.block_size = config.block_size
    self.vocab_size = config.vocab_size
    self.wte = nn.Embedding(config.vocab_size + 1, config.n_embd) # token embeddings table
    # +1 in the line above for a special <BLANK> token that gets inserted if encoding a token
    # before the beginning of the input sequence
    self.mlp = nn.Sequential(
      nn.Linear(self.block_size * config.n_embd, config.n_embd2),
      nn.Tanh(),
      nn.Linear(config.n_embd2, self.vocab_size)
    )

  def get_block_size(self):
    return self.block_size

  def forward(self, idx, targets=None):
    # gather the word embeddings of the previous 3 words
    embs = []
    for k in range(self.block_size):
      tok_emb = self.wte(idx) # token embeddings of shape (b, t, n_embd)
      idx = torch.roll(idx, 1, 1)
      idx[:, 0] = self.vocab_size # special <BLANK> token
      embs.append(tok_emb)

    # concat all of the embeddings together and pass through an MLP
    x = torch.cat(embs, -1) # (b, t, n_embd * block_size)
    logits = self.mlp(x)

    # if we are given some desired targets also calculate the loss
    loss = None
    if targets is not None:
      loss = F.cross_entropy(logits.view(-1, logits.size(-1)), targets.view(-1), ignore_index=-1)

    return logits, loss

4.5 Learning Sequences with Convolutional Neural Networks

$\hookleftarrow$ Table of Contents

# Near copy paste of the layers we have developed in Part 3

# -----------------------------------------------------------------------------------------------
class Linear:
  
  def __init__(self, fan_in, fan_out, bias=True):
    self.weight = torch.randn((fan_in, fan_out)) / fan_in**0.5 # note: kaiming init
    self.bias = torch.zeros(fan_out) if bias else None
  
  def __call__(self, x):
    self.out = x @ self.weight
    if self.bias is not None:
      self.out += self.bias
    return self.out
  
  def parameters(self):
    return [self.weight] + ([] if self.bias is None else [self.bias])

# -----------------------------------------------------------------------------------------------
class BatchNorm1d:
  
  def __init__(self, dim, eps=1e-5, momentum=0.1):
    self.eps = eps
    self.momentum = momentum
    self.training = True
    # parameters (trained with backprop)
    self.gamma = torch.ones(dim)
    self.beta = torch.zeros(dim)
    # buffers (trained with a running 'momentum update')
    self.running_mean = torch.zeros(dim)
    self.running_var = torch.ones(dim)
  
  def __call__(self, x):
    # calculate the forward pass
    if self.training:
      if x.ndim == 2:
        dim = 0
      elif x.ndim == 3:
        dim = (0,1)
      xmean = x.mean(dim, keepdim=True) # batch mean
      xvar = x.var(dim, keepdim=True) # batch variance
    else:
      xmean = self.running_mean
      xvar = self.running_var
    xhat = (x - xmean) / torch.sqrt(xvar + self.eps) # normalize to unit variance
    self.out = self.gamma * xhat + self.beta
    # update the buffers
    if self.training:
      with torch.no_grad():
        self.running_mean = (1 - self.momentum) * self.running_mean + self.momentum * xmean
        self.running_var = (1 - self.momentum) * self.running_var + self.momentum * xvar
    return self.out
  
  def parameters(self):
    return [self.gamma, self.beta]

# -----------------------------------------------------------------------------------------------
class Tanh:
  def __call__(self, x):
    self.out = torch.tanh(x)
    return self.out
  def parameters(self):
    return []

# -----------------------------------------------------------------------------------------------
class Embedding:
  
  def __init__(self, num_embeddings, embedding_dim):
    self.weight = torch.randn((num_embeddings, embedding_dim))
    
  def __call__(self, IX):
    self.out = self.weight[IX]
    return self.out
  
  def parameters(self):
    return [self.weight]

# -----------------------------------------------------------------------------------------------
class FlattenConsecutive:
  
  def __init__(self, n):
    self.n = n
    
  def __call__(self, x):
    B, T, C = x.shape
    x = x.view(B, T//self.n, C*self.n)
    if x.shape[1] == 1:
      x = x.squeeze(1)
    self.out = x
    return self.out
  
  def parameters(self):
    return []

# -----------------------------------------------------------------------------------------------
class Sequential:
  
  def __init__(self, layers):
    self.layers = layers
  
  def __call__(self, x):
    for layer in self.layers:
      x = layer(x)
    self.out = x
    return self.out
  
  def parameters(self):
    # get parameters of all layers and stretch them out into one list
    return [p for layer in self.layers for p in layer.parameters()]

# original network
# n_embd = 10 # the dimensionality of the character embedding vectors
# n_hidden = 300 # the number of neurons in the hidden layer of the MLP
# model = Sequential([
#   Embedding(vocab_size, n_embd),
#   FlattenConsecutive(8), Linear(n_embd * 8, n_hidden, bias=False), BatchNorm1d(n_hidden), Tanh(),
#   Linear(n_hidden, vocab_size),
# ])

# hierarchical network
n_embd = 24 # the dimensionality of the character embedding vectors
n_hidden = 128 # the number of neurons in the hidden layer of the MLP
model = Sequential([
  Embedding(vocab_size, n_embd),
  FlattenConsecutive(2), Linear(n_embd * 2, n_hidden, bias=False), BatchNorm1d(n_hidden), Tanh(),
  FlattenConsecutive(2), Linear(n_hidden*2, n_hidden, bias=False), BatchNorm1d(n_hidden), Tanh(),
  FlattenConsecutive(2), Linear(n_hidden*2, n_hidden, bias=False), BatchNorm1d(n_hidden), Tanh(),
  Linear(n_hidden, vocab_size),
])

# parameter init
with torch.no_grad():
  model.layers[-1].weight *= 0.1 # last layer make less confident

parameters = model.parameters()
print(sum(p.nelement() for p in parameters)) # number of parameters in total
for p in parameters:
  p.requires_grad = True

4.6 Learning Sequences with Recurrent Neural Networks

$\hookleftarrow$ Table of Contents

class RNNCell(nn.Module):
  """
  the job of a 'Cell' is to:
  take input at current time step x_{t} and the hidden state at the
  previous time step h_{t-1} and return the resulting hidden state
  h_{t} at the current timestep
  """
  def __init__(self, config):
    super().__init__()
    self.xh_to_h = nn.Linear(config.n_embd + config.n_embd2, config.n_embd2)

  def forward(self, xt, hprev):
    xh = torch.cat([xt, hprev], dim=1)
    ht = F.tanh(self.xh_to_h(xh))
    return ht

class RNN(nn.Module):
  def __init__(self, config, cell_type):
    super().__init__()
    self.block_size = config.block_size
    self.vocab_size = config.vocab_size
    self.start = nn.Parameter(torch.zeros(1, config.n_embd2)) # the starting hidden state
    self.wte = nn.Embedding(config.vocab_size, config.n_embd) # token embeddings table
    if cell_type == 'rnn':
        self.cell = RNNCell(config)
    elif cell_type == 'gru':
        self.cell = GRUCell(config)
    self.lm_head = nn.Linear(config.n_embd2, self.vocab_size)

  def get_block_size(self):
    return self.block_size

  def forward(self, idx, targets=None):
    device = idx.device
    b, t = idx.size()

    # embed all the integers up front and all at once for efficiency
    emb = self.wte(idx) # (b, t, n_embd)

    # sequentially iterate over the inputs and update the RNN state each tick
    hprev = self.start.expand((b, -1)) # expand out the batch dimension
    hiddens = []
    for i in range(t):
      xt = emb[:, i, :] # (b, n_embd)
      ht = self.cell(xt, hprev) # (b, n_embd2)
      hprev = ht
      hiddens.append(ht)

    # decode the outputs
    hidden = torch.stack(hiddens, 1) # (b, t, n_embd2)
    logits = self.lm_head(hidden)

    # if we are given some desired targets also calculate the loss
    loss = None
    if targets is not None:
      loss = F.cross_entropy(logits.view(-1, logits.size(-1)), targets.view(-1), ignore_index=-1)

    return logits, loss

4.7 Learning Sequences with Generative Pretrained Transformers

$\hookleftarrow$ Table of Contents

class NewGELU(nn.Module):
  """
  Implementation of the GELU activation function currently in Google BERT repo (identical to OpenAI GPT).
  Reference: Gaussian Error Linear Units (GELU) paper: https://arxiv.org/abs/1606.08415
  """
  def forward(self, x):
    return 0.5 * x * (1.0 + torch.tanh(math.sqrt(2.0 / math.pi) * (x + 0.044715 * torch.pow(x, 3.0))))

class CausalSelfAttention(nn.Module):
  """
  A vanilla multi-head masked self-attention layer with a projection at the end.
  It is possible to use torch.nn.MultiheadAttention here but I am including an
  explicit implementation here to show that there is nothing too scary here.
  """

  def __init__(self, config):
    super().__init__()
    assert config.n_embd % config.n_head == 0
    # key, query, value projections for all heads, but in a batch
    self.c_attn = nn.Linear(config.n_embd, 3 * config.n_embd)
    # output projection
    self.c_proj = nn.Linear(config.n_embd, config.n_embd)
    # causal mask to ensure that attention is only applied to the left in the input sequence
    self.register_buffer("bias", torch.tril(torch.ones(config.block_size, config.block_size))
                                  .view(1, 1, config.block_size, config.block_size))
    self.n_head = config.n_head
    self.n_embd = config.n_embd

  def forward(self, x):
    B, T, C = x.size() # batch size, sequence length, embedding dimensionality (n_embd)

    # calculate query, key, values for all heads in batch and move head forward to be the batch dim
    q, k ,v  = self.c_attn(x).split(self.n_embd, dim=2)
    k = k.view(B, T, self.n_head, C // self.n_head).transpose(1, 2) # (B, nh, T, hs)
    q = q.view(B, T, self.n_head, C // self.n_head).transpose(1, 2) # (B, nh, T, hs)
    v = v.view(B, T, self.n_head, C // self.n_head).transpose(1, 2) # (B, nh, T, hs)

    # causal self-attention; Self-attend: (B, nh, T, hs) x (B, nh, hs, T) -> (B, nh, T, T)
    att = (q @ k.transpose(-2, -1)) * (1.0 / math.sqrt(k.size(-1)))
    att = att.masked_fill(self.bias[:,:,:T,:T] == 0, float('-inf'))
    att = F.softmax(att, dim=-1)
    y = att @ v # (B, nh, T, T) x (B, nh, T, hs) -> (B, nh, T, hs)
    y = y.transpose(1, 2).contiguous().view(B, T, C) # re-assemble all head outputs side by side

    # output projection
    y = self.c_proj(y)
    return y

class Block(nn.Module):
  """ an unassuming Transformer block """

  def __init__(self, config):
    super().__init__()
    self.ln_1 = nn.LayerNorm(config.n_embd)
    self.attn = CausalSelfAttention(config)
    self.ln_2 = nn.LayerNorm(config.n_embd)
    self.mlp = nn.ModuleDict(dict(
        c_fc    = nn.Linear(config.n_embd, 4 * config.n_embd),
        c_proj  = nn.Linear(4 * config.n_embd, config.n_embd),
        act     = NewGELU(),
    ))
    m = self.mlp
    self.mlpf = lambda x: m.c_proj(m.act(m.c_fc(x))) # MLP forward

  def forward(self, x):
    x = x + self.attn(self.ln_1(x))
    x = x + self.mlpf(self.ln_2(x))
    return x

class Transformer(nn.Module):
  """ Transformer Language Model, exactly as seen in GPT-2 """

  def __init__(self, config):
    super().__init__()
    self.block_size = config.block_size

    self.transformer = nn.ModuleDict(dict(
        wte = nn.Embedding(config.vocab_size, config.n_embd),
        wpe = nn.Embedding(config.block_size, config.n_embd),
        h = nn.ModuleList([Block(config) for _ in range(config.n_layer)]),
        ln_f = nn.LayerNorm(config.n_embd),
    ))
    self.lm_head = nn.Linear(config.n_embd, config.vocab_size, bias=False)

    # report number of parameters (note we don't count the decoder parameters in lm_head)
    n_params = sum(p.numel() for p in self.transformer.parameters())
    print("number of parameters: %.2fM" % (n_params/1e6,))

  def get_block_size(self):
    return self.block_size

  def forward(self, idx, targets=None):
    device = idx.device
    b, t = idx.size()
    assert t <= self.block_size, f"Cannot forward sequence of length {t}, block size is only {self.block_size}"
    pos = torch.arange(0, t, dtype=torch.long, device=device).unsqueeze(0) # shape (1, t)

    # forward the GPT model itself
    tok_emb = self.transformer.wte(idx) # token embeddings of shape (b, t, n_embd)
    pos_emb = self.transformer.wpe(pos) # position embeddings of shape (1, t, n_embd)
    x = tok_emb + pos_emb
    for block in self.transformer.h:
        x = block(x)
    x = self.transformer.ln_f(x)
    logits = self.lm_head(x)

    # if we are given some desired targets also calculate the loss
    loss = None
    if targets is not None:
        loss = F.cross_entropy(logits.view(-1, logits.size(-1)), targets.view(-1), ignore_index=-1)

    return logits, loss

5. Accelerating Sequence Models on GPU in teenygrad with CUDA Rust

$\hookleftarrow$ Table of Contents

5.1 From Numerical Linear Algebra to Deep Learning Frameworks

$\hookleftarrow$ Table of Contents

5.2 Network Primitives and Optimizers with teenygrad.nn and teenygrad.optim

$\hookleftarrow$ Table of Contents

Consider the function $f:\mathbb{R},\mathbb{R}\to \mathbb{R}$ where $f:=\exp (x_1)\ast {\sin }^2{x}_2$, and translate it to it’s computational counterpart in python with one-dimensional Tensors:

import picograd as pg

def f(x1: pg.Tensor, x2: pg.Tensor) -> pg.Tensor:
  a = pg.exp(x1)
  b = pg.sin(x2)
  c = b**2
  d = a*c
  return d

Figure 1. Python source for the function $f:\mathbb{R},\mathbb{R}\to \mathbb{R}$ where $f:=\exp (x_1)\ast {\sin }^2{x}_2$

Here we’ve broken up the function to render the subexpressions more clearly. But this isn’t necessary — automatic differentiation will work if the function was expressed in one line. In part one, the development of picograd followed that of numpy — an array programming language similar to Matlab but embedded in the host language of Python, that could evaluate functions of the form $f(\mathbf{x})$ where Tensor objects stored their values $f(\mathbf{x})$ with the value: field and the function types that produced their values with Op. For instance, evaluating the specified function f from above with 9 and 10

if __name__ == "__main__":
  print(f(9, 10))

populates the Tensor.value fields. In part one of the book we verified this with a REPL-interface, but we can also represent the entire expression being evaluated with a graph of vertices and edges where the vertices are Tensors (along with their Ops and values) and the edges are their data dependencies:

Here you can see that even if the function was specified in one line, the graph of the expression always parses into Tensor vertices, and data dependency edges. You may have noticed the Tensor.grad fields, which supposedly store the values of derivatives $f^{\prime }(\mathbf{x})$. The question now remains in how to populate these fields.

Taking a step back to differential calculus, deriving the derivative of $f:=\exp (x_1)\ast {\sin }^2{x}_2$ involves the application of the chain rule where $f(x):=g(h(x))⟹f^{\prime }(x)=g^{\prime }(h(x))h^{\prime }(x)$. Evaluating the derivative of the function with respect to its inputs $f^{\prime }(x_1)$ and $f^{\prime }(x_2)$ results in

$$f^{\prime }(x_1)=\frac{\partial }{\partial x_1}[\exp (x_1)\ast {\sin }^2{x}_2]=\frac{\partial }{\partial x_1}[\exp (x_1)]\ast {\sin }^2{x}_2=\exp (x_1)(1)\ast ({\sin }^2{x}_2)$$ $$f^{\prime }(x_2)=\frac{\partial }{\partial x_2}[\exp (x_1)\ast {\sin }^2{x}_2]=\exp (x_1)\ast \frac{\partial }{\partial x_2}[{\sin }^2{x}_2]=\exp (x_1)\ast 2\sin x_2\cos x_2$$

symbolic and numeric differentiattion symbolic differentiation has performance issues since a large unrolled expression must be constructed in order to differentiate[^0], whereas numerical differentiation has correctness issues since evaluating finite differences requires evaluating functions to a precision point resulting in numerical instability. (trace through EXAMPLE for both. talking nets widrow)

To populate the Tensor.grad fields, the simplest idea would be to literally translate the manual derivation of the derivative into code. The translation from math to code involves a design decision: should we evaluate from outputs to inputs (symbolically outside-in, graphically right-to-left) or from inputs to outputs (symbolically inside-out, graphically left-to-right)? Although the former order seems more natural with symbolic expressions, there’s nothing illegal about the latter.

import picograd as pg

def f(x1: pg.Tensor, x2: pg.Tensor) -> pg.Tensor:
  a = pg.exp(x1)
  b = pg.sin(x2)
  c = b**2
  d = a*c
  return d

# dict[f(x), f'(x)] of local derivatives (adjoints)
dd_da, dd_dc = [c, a] # d(a,c):=a*c ==> d'(a)=c, d'(c)=a
da_dx1 = pg.exp(x1) # a(x1):=exp(x1) ==> a'(x1)=exp(x1)
dc_db = 2*b # c(b):=b^2 ==> c'(b)=2b
db_dx2 = pg.cos(x2) # b(x2):=sin(x2) ==> b'(x2)=cos(x2)

# outputs to inputs: outside-in symbolically, right-to-left graphically
dd_dd = pg.Tensor(1) # base case
dd_da, dd_dc = [dd_dd*dd_da, dd_dd*dd_dc]
dd_dx1 = dd_da*da_dx1 # DONE for the x1->d path

dd_db = dd_dc*dc_db
dd_dx1 = dd_db*db_dx2 # DONE for x2->path

# inputs to outputs: inside-out symbolically, left-to-right graphically
dx1_dx1, dx2_dx2 = [pg.Tensor(1), pg.Tensor(1)] # base case
da_dx1 = da_dx1*dx1_dx1
dd_dx1 = dd_da*da_dx1 # DONE for the x1->d path

db_dx2 = db_dx2*dx2_dx2
dc_dx2 = dc_dc*db_dx2
dd_dx2 = dd_dc*dc_dx_2 # DONE for the x2->d path

Do you notice any difference in the number of evaluations between the two orders?

The outputs-to-input ordering takes 6 arithmetic operations (including the destructuring), whereas the input-to-output ordering take 7 arithmetic operations. This is because the former can reuse dd_dd as a dynamic programming solution to a subproblem for the two inputs, whereas the latter cannot. And taking a step back, we only want to reuse the output because the shape of the function is of $f:{\mathbb{R}}^2\to \mathbb{R}$. Alternatively, if $f$ had type $f:\mathbb{R}\to {\mathbb{R}}^2$, then the input-to-output ordering would be able to reuse results. This distinction is referred to as “forward-mode” vs “reverse-mode”, and reflects the fact that for some function $f:{\mathbb{R}}^n\to {\mathbb{R}}^m$ the time complexity of forward-mode differentiation is proportional to $n$, whereas that of forward-mode differentiation is proportional to $m$. If the expression graph fans-in so that $n>m$, reverse-mode is preferred. If the expression graph fans-out so that $m>n$, forward-mode is preferred. However, if we take a step with a graph-theory lens, we can see that the derivative is the sum of paths, where each path is a product of local derivatives from the input source to the output sink. From a combinatorics perspective, we are calculating all the possible (ors) ways (ands) on how the inputs perturb the output. That is:

$$\frac{\partial f}{\partial x_i}=\sum \prod \frac{\partial v_j}{\partial v_k}$$

and as long as the operations along this path are associative — $(AB)C=A(BC)$ — then we can choose the order in how we perform these path products to minimize the number of operations. Finding the optimal ordering is an NP-hard problem because ____. For instance, if the expression graph is diamond-shaped, evaluating the derivative with forward-mode for the left-half and reverse-mode for the right-half would be more performant. In practice, we use reverse-mode as a heuristic, since most of the functions that are differentiated (so they can be optimized) in the field of machine learning are neural networks of the form $f:{\mathbb{R}}^n\to \mathbb{R}$

How can we generalize this into an algorithm?
All we need are 1. mappings from $f(\mathbf{x})\to f^{\prime }(\mathbf{x})$ and 2. a topological sort

For the derivative rules, the same way that optimizing compilers implement an optimization “manually” once which then gets reused many times, the authors of deep learning frameworks also implement derivatives manually which then become reused many times through automatic differentiation. In theory, we can differentiate any expression with f’(x) with only a few derivative rules for addition and multiplication, but in practice most frameworks provide sugar for complex derivatives.

For topological sort, we can simply reversed the ordering produced by a depth-first-search:

def toposort(self):
  order: list[Op] = []
  visited: set[Op] = set()

  def dfs(node: Op) -> None:
    if node in visited: return
    visited.add(node)
    for src in node.src: dfs(src)
    order.append(node)

  dfs(self)
  return order

class Tensor():
  def backward():
    for t in reversed(topo):
      t.backward()

We will now use this idea to modify the interpretation of our deep learning framework to not only evaluate $f(\mathbf{x})$, but $f^{\prime }(\mathbf{x})$ as well. This is done by dynamically overloading the operators at runtime[^0] to trace the expression graph

chain_rules = PatternMatcher([
  (Pattern(OpCode.MATMUL, name="input"), lambda output_grad, input: (_____,)),
  (Pattern(OpCode.MATVEC, name="input"), lambda output_grad, input: (_____,)),
  (Pattern(OpCode.RECIPROCAL, name="input"), lambda output_grad, input: (-output_grad * input * input,)),
  (Pattern(OpCode.SIN, name="input"), lambda output_grad, input: ((math.pi/2 - input.src[0]).sin() * output_grad,)),
  (Pattern(OpCode.LOG2, name="input"), lambda output_grad, input: (output_grad / (input.src[0] * math.log(2)),)),
  (Pattern(OpCode.EXP2, name="input"), lambda output_grad, input: (input * output_grad * math.log(2),)),
  (Pattern(OpCode.SQRT, name="input"), lambda output_grad, input: (output_grad / (input*2),)),
  (Pattern(OpCode.ADD), lambda output_grad: (1.0*output_grad, 1.0*output_grad)),
  (Pattern(OpCode.MUL, name="input"), lambda output_grad, input: (input.src[1]*output_grad, input.src[0]*output_grad)),
])

class Tensor:
  def _forward(self, f:Callable, *other:Tensor) -> Tensor: #extra_args=(), **kwargs)
    out_tensor = evaluator.eval_uop([self, other], out_uop)

  def backward(self, grad:Tensor|None=None) -> Tensor:
    """
    backward performs by collecting tensors, computing gradients with automatic differentiation, and updating said tensors.
    """
    # 1. collect all tensors that requires grad by topologically sorting the graph of uops and filter
    all_uops = self.uop.toposort()
    tensors_require_grad: list[Tensor] = [t for tref in all_tensors if (t:=tref()) is not None and t.uop in all_uops and t.requires_grad]
    uops_require_grad = [t.uop for t in tensors_require_grad]
    assert grad is not None or self.shape == tuple(), "when no gradient is provided, backward must be called on a scalar tensor"
    if not (self.is_floating_point() and all(t.is_floating_point() for t in tensors_require_grad)): raise RuntimeError("only float Tensors have gradient")
    
    # 2. compute the gradient with a map of tensors to partials
    if grad is None: grad = Tensor(1.0, dtype=self.dtype, device=self.device, requires_grad=False) # base case is 1.0
    tens2grads = Tensor._automatically_differentiate(self.uop, grad.uop, set(uops_require_grad)) # skipping materializing zerod grads for now
    grads = [Tensor(g, device=t.device) for t,g in zip(tens2grads.keys, tens2grads.values)] # initialize tensor grads on device
    
    # 3. update the tensors that require grad with the gradient's partials
    for t,g in zip(tensors_require_grad, grads):
      assert g.shape == t.shape, f"grad shape must match tensor shape, {g.shape!r} != {t.shape!r}"
      t.grad = g if t.grad is None else (t.grad + g) # accumulate if t.grad exists
    return self

  @staticmethod
  def _automatically_differentiate(root:Op, root_grad:Op, targets:set[Op]) -> dict[Op, Op]:
    """
    _differentiate backpropagates partials on a topologically sorted expression graph with the chain rule
    and produces the gradient in the form of a map of ops to their partials (which, in turn, are ops)
    """
    tens2grads = {root: root_grad}

    # 1. topological sort
    in_target_path: dict[Op, bool] = {}
    for u in root.toposort(): in_target_path[u] = any(x in targets or in_target_path[x] for x in u.src)
    dfs = list(root.toposort()) # lambda node: node.op not in {OpCode.DETACH, OpCode.ASSIGN} and in_target_path[node])) # don't flow through DETACH/ASSIGN or anything not in target path

    # 2. backpropagation with the chain rule
    for tensor in reversed(dfs):
      if tensor not in tens2grads: continue

      local_grads: tuple[Op|None, ...]|None = cast(tuple[Op, ...]|None, chain_rules.rewrite(tensor, ctx=tens2grads[tensor]))
      if local_grads is None: raise RuntimeError(f"failed to compute gradient for {tensor.op}\n\nin {str(tensor)[0:1000]}...")
      assert len(local_grads) == len(tensor.src), f"got {len(local_grads)} gradient, expected {len(tensor.src)}"

      for tensor,local_grad in zip(tensor.src, local_grads): # <--------------------- MOOOSE: why are we accumulating inside ad()? don't we do it in backward()??
        if local_grad is None: continue
        if tensor in tens2grads: tens2grads[tensor] = tens2grads[tensor] + local_grad # accumulate if tensor exists
        else: tens2grads[tensor] = local_grad # o/w initialize

To implement automatic differentiation with Tensor.backward(), there is a design decision to be made — the choice of implementing it dynamically or just-in-time[^3], similar to the decision of how to implement types for general programming languages[^4]. This stands in contrast to the alternative of performing a just-in-time, source-to-source transformation.

Let’s now move onto automatically differentiating the functions of neural networks, specifically the FFN language model from earlier. (johnson/ryan adams ordering) n^2 vs n^3

5.3 Automatic Differentiation with Tensor.forward() and Tensor.backward()

$\hookleftarrow$ Table of Contents

5.4 From SIMD of Multi-Core Latency-Oriented Processors to SIMT of Many-Core Throughput-Oriented Processors with PTX

$\hookleftarrow$ Table of Contents

5.5 Accelerating GEMV on GPU with CUDA Rust via Rooflines

$\hookleftarrow$ Table of Contents

#![allow(unused)]
fn main() {
// gpu_host.rs
use cudarc::{driver::{self, PushKernelArg}, nvrtc};
use src_device::T; // shared type with device code
static PTX: &str = include_str!(concat!(env!("OUT_DIR"), "/gpu_device.ptx")); // Embed the PTX code as a static string.

pub fn cudars_helloworld() -> Result<(), Box<dyn std::error::Error>> {
  // initialize device context and stream via driver api
  let process = driver::CudaContext::new(0)?; // device 0
  let queue = process.default_stream();
  
  // load ptx via nvrtc
  let dylib = process.load_module(nvrtc::Ptx::from_src(PTX))?;
  let add_kernel = dylib.load_function("add")?;

  // allocate on device
  let (a, b): ([T; _], [T; _]) = ([1.0, 2.0, 3.0, 4.0], [2.0, 3.0, 4.0, 5.0]);
  let (a_gpu, b_gpu, mut c_gpu) = (queue.clone_htod(&a)?, queue.clone_htod(&b)?, queue.alloc_zeros::<T>(a.len())?);
  let (a_len, b_len) = (a_gpu.len(), b_gpu.len());

  let cfg = driver::LaunchConfig { grid_dim: (1, 1, 1), block_dim: (4, 1, 1), shared_mem_bytes: 0, };
  unsafe {
    queue
    .launch_builder(&add_kernel).arg(&a_gpu).arg(&a_len).arg(&b_gpu).arg(&b_len).arg(&mut c_gpu)
    .launch(cfg)?;
  }
  queue.synchronize()?;

  let c = queue.clone_dtoh(&c_gpu)?;
  println!("c from cuda is = {:?}", c);
  Ok(())
}
}

#![allow(unused)]
fn main() {
// gpu_device.rs
use cuda_std::kernel;
use crate::T;

#[allow(improper_ctypes_definitions)]
#[kernel] pub unsafe fn add(a: &[T], b: &[T], c: *mut T) {
  let i = cuda_std::thread::index_1d() as usize;
  if i < a.len() {
    let elem = unsafe { &mut *c.add(i) };
    *elem = a[i] + b[i];
  }
}

#[allow(improper_ctypes_definitions)]
#[kernel] pub unsafe fn saxpy(a: &[T], b: &[T], c: *mut T) {
  let i = cuda_std::thread::index_1d() as usize;
  todo!()
}

#[allow(improper_ctypes_definitions)]
#[kernel] pub unsafe fn smul(a: &[T], b: &[T], c: *mut T) {
  let i = cuda_std::thread::index_1d() as usize;
  todo!()
}

#[allow(improper_ctypes_definitions)]
#[kernel] pub unsafe fn stanh(a: &[T], b: &[T], c: *mut T) {
  let i = cuda_std::thread::index_1d() as usize;
  todo!()
}
}

5.6 Accelerating _ of GEMM on GPU with Data Reuse

$\hookleftarrow$ Table of Contents

5.7 Accelerating _ of GEMM on GPU with Scheduling

$\hookleftarrow$ Table of Contents

III. Scaling Networks

Part 3 covers the age of scaling by building a distributed tensor compiler using tinygrad IR

Table of Contents

• 6. Scaling LLMs
• 7. Accelerating LLMs
• 8. Biology of LLMs
• 9. Physics of LLMs

Afterword

From nanochat to llama and deepseek

From teenygrad to torch and jax

Appendix

Table of Contents

A. From Symbolic Software 1.0 to Stochastic Software 2.0

$\hookleftarrow$ Table of Contents

In which we historically retrace the development and failure of the discretely symbolic approach to build artificially intelligent machines with common sense and motivate the need to transition from logical and finitely discrete software 1.0 to stochastic and infintely continuous software 2.0.

A.1 From Psychology to Artificial Intelligence

The transition between the two representations happened relatively recently at a summer worshop at Dartmouth in 1956. There, a group of researchers unsatisfied with the theories that the discipline of psychology were using to explain the phenomena of the mind and it’s intelligence came together to discuss a different approach, namely, one where the computer is the instrument for conducting scientific experiment. Although seemingly trivial from the modern perspective where most if not all sciences use the computer, they were were arguably the first with motivation arisen from the epistemological: using the computer as basis for the science of mind (and all sciences in general) strengthened it’s explanations from the observationally simple like stimulus-response to the constructively complex such as neural networksPractical applications are often a result of inquiry that is philosophical and gradiose with no immediately obvious economic value. Namely, computers with Hilbert wanting to automate mathematics as beers, tables and mugs; language models with McCarthy, Minsky, Newell and Simon wanting to mechanize and naturalize the mind.. That is, constructive because explaining via computer means simulating the phenomena by programming processes with procedures. And, complex because computers allow for the simulating of many things at onceParaphrasing Minsky, “Under certain conditions mathematical analysis can describe complex phenomena where the parts of the system can be treated as individual and independently random (i.e statistical thermodynamics), but there is no reason to suspect that intelligence is the result of averaging out many events.”. The proposal for the workshop states:

We propose that a (…) study of artificial intelligence be carried out (…) the study is to proceed on the basis of the conjecture that every aspect of learning or any other feature of intelligence can in principle be so precisely described that a machine can be made to simulate it.

Besides the intellectual pursuit of finding better explanations for a clearer picture of reality, using the computer also means something quite practically profound. If the explanation it comes up with are accurate, we will have artificial systems that exhibit behavior which we would attribute intelligence to. This is what the Turing Test posited and predicted in 195X (todo, read computing machinery and intelligence). artificial intelligenceartificial intelligence natural language processingnatural language processing computational linguisticscomputational linguistics. And this is why we have ChatGPT.

In this book we embark on a quest to build from scratch our own deep neural network like ChatGPT by implementing nanochat and our own deep learning framework like PyTorch by implementing teenygrad capable of running nanochat itself. These systems by nature are stochastic and infintelyTurns out not quite infinte, as we will see in chapter 3. continuous software 2.0 rather than the logical and finitely discrete software 1.0 and are implemented not by programming algorithms and their procedures line by line with sets, maps, lists, trees, and graphs, but rather, by searching the space of programs by providing a goal to calculus, which then optimizes said goal — in the case of ChatGPT, producing a probability distribution over tokens — with the linear algebra of tensors. However, there was a time where the dominant approach involved using software 1.0 and in chapter 0 we will build various systems using such techniques to display their shortcomings, understanding the underlying philosophical principle, and ultimately motivating the need for software 2.0The art of programming software 1.0 is necessary however on your quest to learn software 2.0! PyTorch is embedded and implemented within Python afterall. For instance, to those who spent countless nights learning esotoric spells such as that of dynamic programming to enter the kingdoms of our feudal lords only to create web page buttons should not fret as it turns out that dynamic programming over a graph is in fact the beating heart of all deep learning frameworks.

If you’d like to revisit the fundamentals of programming, we recommend the Data Centric Introduction to Computing, which begins with the teaching language Pyret and graduates to Python. You can then take a look at the documentation of the Python Tutorial, Python Language Reference and Python Standard Library.

briefly mention intelligence must be told knowledge before learning it and such knowledge should be represented a symbolic logic

• mccarthy’s excerpts from “programs with common sense”
• minsky’s “descriptive languages and problem solving”
• newell and simon’s symbolic hypothesis
• newell knowledge level

Are you ready to begin?

A.2 Weizenbaum Cheats Turing’s Test with the Pattern Matching of ELIZA

Humans, it seems, know things: and what they know helps them do things. The early approach to artificial intelligence using logical and finitely discrete techniques from software 1.0 focused on building systems that reasoningreasoned over an internal representationrepresentation of knowledge. Iteratively deepening software 1.0’s symbolic perspective of such terms will be the focus of this first chapterBy the end of the book you will have come to understand the software 2.0 perspective of such terms..

Although there was various flavors of the symbolic approach to AI — game playing, puzzle solving, problem solving to name a few — our focus is on building conversational machines within the realm of natural language processing and computational linguistics, to pass something like the aforementioned Turing Test. With that said, what is the simplest way to build a conversational machine with the logical and finitely discrete techniques from software 1.0?

Warning

Pause and think!

What if we represented words with strs and, produce answers as output with if statements conditioned on questions as input?

That’s effectively what the ELIZASee A Computer Program For the Study of Natural Language Communication Between Man and Machine (Weizenbaum 1966) system doesDo you feel dissapointed after learning ELIZA’s trick? The trick with all explanations is that after the explanation, no trick remains. See Matter, Mind, and Models (Minsky 1965), and The Nature of Explanation (Craik 1952). You might feel the same way by the end of the book even after learning how nanochat and teenygrad work under the hood. Don’t say we didn’t warn you! to imitate a Rogerian psychotherapist.

For instance, a question that tends to get asked is the meaning of life. Answering such question seems quite grandiose for now, so let’s have our system produce some random string for now. Let’s document and interatively test our exampleThe implementation of all functions implemented in the book will start with examples, following the principled design of programs with the How to Design Programs’ Design Recipe. We still find it useful in the era of agentic coding. with Python’s convinent standard library module docttest:

def eliza(input: str) -> str:
  """
  >>> eliza("What's the meaning of life?")
  '42'
  """
  if input == "What's the meaning of life?": return "42"
  else: raise NotImplementedError("")

  if __name__ == "__main__":
    import doctest
    doctest.testmod()

Let’s start with other questions that are perhaps less grandiose but as equally important in which a patient might ask a psychotherapist. For instance, a patient reporting to a therapist that they are unhappy or upset. Given that ELIZA imitates a Rogerian psychotherapist which follows the principle of person-centered therapy — that is, no immediate rejection (todo: read wiki) — we might expect on an a priori basisWe can also empirically confirm the following question-answer pairs by consulting the appendix in (Weizenbaum 1963). that our system responds like so:

def eliza(input: str) -> str:
  """
  >>> eliza("What's the meaning of life?")
  '42'

  >>> eliza("I am unhappy")
  'Why do you say you are unhappy?'
  """
  if input == "What's the meaning of life?": return "42"
  else:                                       raise NotImplementedError("")

if __name__ == "__main__":
  import doctest
  doctest.testmod()

Evaluating the tests fails as expected. How should we implement the function body for eliza() so that they pass?

Warning

Pause and think!

The most naive way to make them pass is to add an if-then rule for each example, following the question-answer pair for the meaning of life:

def eliza(input: str) -> str:
  """
  >>> eliza("What's the meaning of life?")
  '42'

  >>> eliza("I am unhappy")
  'Why do you say you are unhappy?'
  """
  if input == "What's the meaning of life?": return "42"
  elif input == "I am unhappy":              return "Why do you say you are unhappy?"
  else:                                       raise NotImplementedError("")

if __name__ == "__main__":
  import doctest
  doctest.testmod()

Clearly eliza lacks any true understanding of word meaning found in natural language, for it’s simply reflecting the prompt back to the user. In the paper:

The ELIZA program itself is merely a translating processor in the technical programming sense. Gorn [2] in a paper on language systems says: ‘Given a language which already possesses semanticssemantic content, then a translating processor, even if it operates only syntaxsyntactically, generates corresponding expressions of another language to which we can attribute as “meanings” (possibly multiple — the translator may not be one to one) the “semantic intents” of the generating source expressions; whether we find the result consistent or useful or both is, of course, another problem.’

The classic linguistics example to distinguish syntactic form and semantic meaning comes from Syntactic Structures (Chomsky 1957):

• Furiously sleep ideas green colorless
• Colorless green ideas sleep furiously

where the first sentence is gramatically incorrect whereas the second, while gramatically correct, is semantically meaningless(todo).. With ELIZA however, it doesn’t outright produce sentences that are as meaningless as the second sentence, but as the interaction with the chatbot progresses, the mirage of such semantic understanding unveils itself, and most people start to understand the gist of ELIZA’s gimmicks.

But even if only operating syntactically speaking, another issue is that the implementation clearly does not scale, for in the case where our patient prompts eliza with cases not handled — that they are upset for instance — our implementation immediately fails:

def eliza(input: str) -> str:
  """
  >>> eliza("What's the meaning of life?")
  '42'

  >>> eliza("I am unhappy")
  'Why do you say you are unhappy?'

  >>> eliza("I am upset")
  'Why do you say you are upset?' # <-- FAIL

  >>> eliza("Why does Alice hate me?")
  'Why do you say Alice hates you?' # <-- FAIL
  """
  if input == "What's the meaning of life?": return "42"
  elif input == "I am unhappy":              return "Why do you say you are unhappy?"
  else:                                       raise NotImplementedError("")

if __name__ == "__main__":
  import doctest
  doctest.testmod()

We can add another conditional statement to handle such case, but then the the patient can come in reporting yet another sentiment after that. Although enumerating through the entire space of possible questions a patient could ask is indeed intractable, perhaps we could collapse said space with a few conditional statements that provided reuse within each branch. For instance, in the case where a patient’s prompt takes the syntactical form “I am BLAH”, eliza can respond with “Why do you say you are BLAH?” independent of BLAH’s semantic meaning. In another case with “BLAH hates me”, eliza can respond with “Why do you say BLAH hates you?”.

def eliza(input: str) -> str:
  """
  >>> eliza("What's the meaning of life?")
  '42'

  >>> eliza("I am unhappy")
  'Why do you say you are unhappy?'

  >>> eliza("I am upset")
  'Why do you say you are upset?' # <-- FAIL

  >>> eliza("Why does Alice hate me?")
  'Why do you say Alice hates you?' # <-- FAIL
  """
  if input == "What's the meaning of life?": return "42"
  elif input == "I am BLAH":                 return "Why do you say you are BLAH?"
  elif input == "BLAH hates me":             return "Why do you say BLAH hates me?"
  else:                                       raise NotImplementedError("")

if __name__ == "__main__":
  import doctest
  doctest.testmod()

In order to implement the code sketch above, some formal language theory is needed. The theory models languagelanguage as a set of strings, where each string is a sequence of elements from some finite alphabetalphabet. Even if such set is infinite, the set itself can be characterized with a finite set of rules. The core interest of such theory are the syntactical aspects of languages, namely the membership problemmembership problem. That is, to determine based off structural form whether a given string is in a language or not. Given that ELIZA is merely operating with the syntactic structure of the “Rogerian psychotherapist language”, that is not a problem.

In our case, we’d like to somehow define the language of all strings that take the form “I am BLAH”, and then match all strings that are inside that set. How do we characterize such a set? We can do so with a regular expressionregular expression, which defines the regular languageregular language of said strings, “I am BLAH”. A regular expression is one which can include the following elements

• a literal character drawn from some alphabet $\Sigma$
• the empty string $ϵ$
• the Kleene star $R^{\ast }$, where $R$ is a regular expression
• concatenation $RS$, where $R$ and $S$ are regular expressions
• alternation $R|S$, where $R$ and $S$ are regular expressions
• and parentheses $(R)$, where $R$ is a regular expression

With Python particularly, regular expressions are available via standard library’s re module with a two step process.

1. The first step is passing a regular expression to re.compile() to produce a re.Pattern object.
2. Then, the second step is to match against said pattern with an input string via Pattern.search(string), Pattern.match(string) or Pattern.fullmatch(string) which returns a corresponding re.Match object or None.

However if the pattern is only going to be matched against a single time without any reuse, you can evaluate the re.Pattern and re.Match objects with a single function invocation rather than two. For example:

import re
pattern = re.compile(r"I am unhappy")
result1 = pattern.fullmatch(r"I am unhappy")    # match
result2 = pattern.fullmatch(r"foobar")          # no match
result3 = pattern.fullmatch(r"You are unhappy") # close, but still no match
print(f'{result1=}')
print(f'{result2=}')
print(f'{result3=}')

result1_singlestep = re.fullmatch(r"I am unhappy", "I am unhappy") # match, with a single call

In our case where we’d like to characterize all strings that take the form “I am BLAH”, we need to use the Kleene star which effectively acts as a wildcard, and subsequently capture the BLAH with Math.groups()

import re
result = re.fullmatch("I am (.*)", input)
print("captured: {0}?".format(*result.groups()))

Using our new machinery with regular expression in eliza()’s implementation, we now have:

import re

def eliza(input: str) -> str:
  """
  >>> eliza("What's the meaning of life?")
  '42'

  >>> eliza("I am unhappy")
  'Why do you say you are unhappy?'

  >>> eliza("I am upset")
  'Why do you say you are upset?'
  
  >>> eliza("Why does Alice hate me?")
  'Why do you say Alice hates you?'
  """
  if match := re.fullmatch("What's the meaning of life?", input): return "42"
  elif match := re.fullmatch("I am (.*)", input):                 return "Why do you say you are {0}?".format(*m.groups())
  elif match := re.fullmatch("(.*) hate (.*)", input):            return "Why do you think {0} hates {1}?".format(*m.groups())
  else:                                                           raise NotImplementedError("")

print(eliza("I am very unhappy these days")) # Why do you say you are unhappy?

Tests pass! Let’s add some more rules, including a catchall rule in the else branch where any phrase uttered in which eliza() does not recognize will be responded with "Please go on.". It’s effectively an escape hatch for eliza()’s utter lack of semantic understanding with word meaning.

todo you -> me me -> you

import re

def eliza(input: str) -> str:
  """
  >>> eliza("What's the meaning of life?")
  '42'

  >>> eliza("I am unhappy")
  'Why do you say you are unhappy?'

  >>> eliza("I am upset")
  'Why do you say you are upset?'
  """
  if match := re.fullmatch("What's the meaning of life?", input): return "42"
  elif match := re.fullmatch("I am (.*)", input):                 return "Why do you say you are {0}?".format(*m.groups())
  elif match := re.fullmatch("It seems that (.*)", input):        return "What makes you think {0}?".format(*m.groups())
  elif match := re.fullmatch("(.*) hate (.*)", input):            return "Why do you think {0} hates {1}?".format(*m.groups())
  else:                                                           return "Please go on." # <-- the magic trick

print(eliza("I am very unhappy these days")) # Why do you say you are unhappy?
print(eliza("I am very unhappy these days")) # How long have you been very unhappy these days?
print(eliza("It seems that you hate me"))    # What makes you think you hate me?

With the last transformation rule, you can see how brittle ELIZA’s so-called “understanding”, or semantics, truly is. The primary reason a simple pattern matcher over strings can be endowed with human understanding (in other words, why the magic works) is because of the psychiatric context — especially the Rogerian one with person-centered therapy — where users are effectively talking with oneselvesIn a 1978 interview, “Well, I would deny that that there’s any important sense, non-negligible sense in which the program understands. It certainly creates the illusion of understanding. there’s no question about that. But we have to understand that that illusion is an attribution that the person conversing with the program contributes to the conversation. It’s not a function of the program itself.”. The paper goes on to report that:

This mode of conversation was chosen because the psychiatric interview is one of the few examples of categorized dyadic natural language communication in which one of the participating pair is free to assume the pose of knowing almost nothing of the real world. If, for example, one were to tell a psychiatrist “I went for a long boat ride” and he responded “Tell me about boats”, one would not assume that he knew nothing about boats, but that he had some purpose in so directing the subsequent conversation. It is important to note that this assumption is one made by the speaker. Whether it is realistic or not is an altogether separate question. In any case, it has a crucial psychological utility in that it serves the speaker to maintain his sense of being heard and understood.

Important

Our eliza() will not qualitatively improve it’s breadth of common sense nor it’s depth of understanding by simply adding another if-then rule. This is because it’s trying to describe a reality with too many parts to count.

Describing a reality with too many parts is the philosophical principle and problem that logical and finitely discrete techniques from software 1.0 ultimately run into. But perhaps it’s too soon to jump the software 1.0 ship to the stochastic and infintely continuous methods of software 2.0? Afterall, ELIZA’s representation are only strs, and it performs no reasoning whatsoever. Maybe all we need to build a natural language processing system that can match the capability of nanochat is with a stronger syntactic and semantic analysis?

A.3 Wood’s Winograd Challenge with the Translation of LUNAR

(INSERT LUNAR EXAMPLE)

Computers are being used today to take over many of our jobs. They can perform millions of calculations in a second, handle mountains of data, and perform routine office work much more efficiently and accurately than humans. But when it comes to telling them what to do, they are tyrants. They insist on being spoken to in special computer languages, and act as though they can’t even understand a simple English sentence.

Let us envision a new way of using computers so they can take instructions in a way suited to their jobs. We will talk to them just as we talk to a research asisstant, librarian, or secretary, and they will carry out our commands and provide us with the information we ask for. If our instructions aren’t clear enough, they will ask for more information before they do what we want, and this dialog will all be in English.

Procedures as a Representation for Data in a Computer Program for Understanding Natural Language (Winograd’s Dissertation 1971)

change segue to be more historical (minsky’s semantic informationation processing, the intros to wood and winograd’s dissertations)

Starting where the previous chapter left off, how can we build a natural language processing system with a stronger analysis of syntactic (grammar) and semantic (meaning) of the English language?

Warning

Pause and think! When questions become more difficult to think about from first principles perhaps such as this one, one heuristic is to use a combination of both history and theories from other disciplines as a guide.

That’s effectively what the LUNAR and SHRDLU systemsPresented across a series of papers. See (Woods 1969), (Woods 1970), (Woods 1971), (Woods 1972), and (Winograd 1971), respectively. do. The former, implements a natural language processing system to interface with a database containing chemical analysis data on lunar rock and soil from the from the Apollo 11 moon missions so that non-expert geologists could express their questions in the natural language of English rather than a formal one such as SQL. The latter, implements the same but rather than interface with a database, it does so with a 3D graphically simulated world with blocksCalled Blocks World, a toy environment which was home to many of MIT’s symbolic AI projects.. The LUNAR and SHRDLU systems were amongst the first natural language proecssing systems to implement a deeper level of word understanding via natural language compilers, compared to systems like ELIZA. Like many of the AI systems at the time, although seemingly applied to the narrow domains of toy environments, these researchers were after general principles.

(quote something from woods, winograd, or minsky’s semantic information processing)

A.3.1 From Linguistics to Logic

Although we might not know exactly how that’s done, we might be somewhat familiar with the distinction of “lower level” vs “higher level” languages where the former are closer to the machine (whatever that may mean) such as C++ and the latter are English-like (whatever that may mean) like Python. Whatever is going underneath there, it seems like these languages understand at least something about our intent in order to preserve it through the translation from a higher level language to a lower level one. Perhaps we can use the same techniques that formal programming languages use but extend them one level “higher” by applying them to the natural language of English itself?

With the approach settled, the problem reduces down into answering the question of how do we implement a compiler? Depending on your appetite, you can try to tackle this question from first principles as well. We however don’t have the ability nor time to perform such a feat, and so we will consult the existing literatureWe recommend Programming Languages: Application and Interpretation (Krishmaurthi 2022).

The first idea to understand is that a compiler is a function which takes a string in and produces a string out (programs are written in text editors afterall) where the input string is referred to as the source program and the output string as the target program. While a traditional compiler might take in C as the source language and produce x86 as the target language, LUNAR takes in a question expressed in English and produces a query expressed in SQL. (todo: remove SQL) Just like a human translator, a compiler has a two step process:

1. First, a compiler needs to understand the message expressed in the source language
2. Then, said compiler can translate the message expressed in the target language

A.3.1 Analysis: Understanding the Source

Starting with the first step, that awfully sounds like formal language theory’s membership problem which we encountered in the previous chapter. That is, determining whether a given string is in a language (subsequently modeled as a set of strings) or not. Perhaps we can reuse regular expressions? Why not try the machinery we already have?

Clearly, regular expressions have no problem describing languages

memory to model long range dependencies i.e a^nb^n (push down, recursion (via stack)) a lexical analysislexical analysis are defined by context-free grammar loremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsum a lexical grammarlexical grammar are defined by context-free grammar loremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsum

a syntactic analysissyntactic analysis are defined by context-free grammar loremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsum a syntactic grammarsyntactic grammar are defined by context-free grammar loremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsum a syntax treesyntax tree are defined by context-free grammar loremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsum a derivationderivation loremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsum a parserparsing is the problem of finding a derivation for a string in a grammar (recognizer)

a context-free languagecontext-free language are defined by context-free grammar loremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsum a context-free grammarcontext-free grammar loremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsumloremipsum

It first converts the source program from a concrete syntaxconcrete str to a more abstract syntaxabstract data structure that has more “understanding” of what’s being spoken to it compared to a simple silly str. Because this abstract representation of the program is not the final representation (namely x86), it’s referred to as the <span

A.3.2 Synthesis: Translating to the Target

While there do exist compilers that perform this translation in a single step, most compilers (and interpreters for that matter) allocate an intermediate representationintermediate representation of the program in order to better analyze, understand, and perhaps optimize said programSaid compilers that translate in a single step are single-pass compilers, primarily from the past when memory bottlenecks prevented the allocation of data structures that represented the entire source program. This is why you can’t invoke a function above it’s definition in C, unless you explicitly provide a forward declaration.. The translation step from source to intermediate representation and from intermediate representation to target are called parsingparsing and generationgeneration respectivelyWith all industrial language implementations, there is almost always a third step in between the two, which optimizes the program.

def compile(source: str) -> str:
  """translates a source program expressed in C to a target program expressed in x86

  """
  ir = parse(source)
  target = generate(ir)
  return target

from dataclasses import dataclass
@dataclass
class IR():
  # ???
def parse(source: str) -> IR: raise NotImplementedError("")
def generate(ir: IR) -> str: raise NotImplementedError("")

def lunar(english_source: str) -> str:
  """translates a source question expressed in english to a target question expressed in SQL
  (todo, examples)
  """
  ir = parse(source)
  target = generate(ir)
  return target

from dataclasses import dataclass
@dataclass
class IR():
  # ???

def parse(source: str) -> IR: raise NotImplementedError("")
def generate(ir: IR) -> str: raise NotImplementedError("")

Evaluating the code above fails as expected. Before implementing the respective function bodies of parse() and generate(), we need to design the IR data structure.

A compiler however is only a program translatortranslator, namely from the C language to the x86 language. In order to produce a final answer you will need a program evaluatorevaluator for x86, which is any Intel processor that understands x86. Similarly with LUNAR, after translating the query from English to SQL, a database that understands SQL is neededWe will learn about the internals of the former evaluator (a hardware processor) in Chapter 3 when accelerating linear algebra in order to implement teenygrad. The later evaluator (a database processor) however is unfortunately out of scope — if you’d like to learn more about database internals, we recommend Andy Pavlo’s CMU 15-445/645.. Following the Design Recipe, let’s start with examples for lunar().

def lunar(english_source: str) -> str:
  """translates a source question expressed in english to a target question expressed in SQL

  """
  ir = parse(source)
  target = generate(ir)
  return target

from dataclasses import dataclass
@dataclass
class IR():
  # ???
def parse(source: str) -> IR: raise NotImplementedError("")
def generate(ir: IR) -> str: raise NotImplementedError("")

carnap, montague, frege (philosophical principles)

(Woods 1972 BBN Report)

Although the goal of accepting an input request in any phrasing which a user might ask is one which will require additional grammar development and semantic work, the system has already achieved considerable progress towards this goal, and the components and organization which we have used in building the system permit conditinoal gradual evolution towards its achievement.

In Progress in Natural Language Understanding — An Application to Lunar Geology (Woods 1973)

The advent of computer networks such as the ARPA net has significantly increased the opportunity for access by a single researcher to a variety of different computer facilities and data bases, thus raising expectations of a day when it will be a common occurrence rather than an exception that a scientist will casually undertake to use a computer facility located 3000 miles away and whose languages, formats, and convetions are unknown to him. In this foreseeable future, learning and remembering the number of different languages and convetions that such a scientist would have to know will require significant effort — much greater than that now required to learn the conventions of his local computing center.

• end the chapter with lighthill
• end chapter with winograd’s two phd students sergey brin and larry page duck semantics with information retrieval and search engine.
• motivate expert systems (feigenbaym)
• the other reaction is lenat with CYC (let’s increase the knowledge base)

Important

Although lunar()’s depth of understanding the natural language of English is deeper than eliza()s, it came at the cost of a shallower breadth in it’s common sense. That is, simply adding more syntactic grammar or semantic interpretation rules will not enable lunar() to start helping patients in the psychiatric setting because it will always be limited to the chemical analysis data in it’s knowledge base from the Apollo 11 moon missions. Although in different ways, LUNAR too, like ELIZA, is trying to describe a reality with too many parts to count.

A.4 Lenat’s Advice Taker with the Frames of CYC

• minsky’s programs with common sense
• mccarthy’s ontology of frames

A.5 From A Logical to Distributional Semantics

feigenbaum’s concept learning lenat’s bitter lesson: from the tractatus to the investigations

lighthill report (1973)

then expert systems in ’80s (feigenbaum and raj reddy), expert systems being abandoned in ’90s, creating the second winter. Parallel Distributed Processing (Rumelhart and McClelland, 1986)

although obvious posthoc that neural networks, this was all predicated with foresight by wittgenstein.

vector semantics by (Osgood et al. 1957) distributional semantics (Harris 1954)

• from feigenbaum/reddy to pdp
• from the organon (knowledge representation and reasoning with upper ontologies and deductive inference) to norvum organon (occam’s razor)
• from the tractatus to the investigations is effectively the transition from software 1.0 to sofware 2.0
• to understand the claude, we must return to claude
• data science begins where computer science begins

question answering systems eventually incorporated the web as it’s knowledge base, and the field of information retrieval emerged. https://start.csail.mit.edu/index.php

A.6 Summary

One quick way to summarize the software 1.0 approach to AI is to list the first six Turing Award winners for AI: Marvin Minsky (1969) and John McCarthy (1971) for defining the foundations of the field based on representation and reasoning; Allen Newell and Herbert Simon (1975) for symbolic models of problem solving and human cognition; Ed Feigenbaum and Raj Reddy (1994) for developing expert systems that encode human knowledge to solve real-world problems. Although not comprehensive, we explored the flavor of logically and discretely finite methods that the software 1.0 approach to AI employed focused on natural language processing and computational linguistics by implementing pattern matching with ELIZA, compilation with LUNAR, and inference with CYC. There were many other earlier approaches to embedding machines with intelligence such as game playing, solving math problems, and ___, which can be found in the secondary resources listed in the bibliographic notes.

The remainder of the book is spent focused on the software 2.0 approach to AI: Judea Pearl (2011) for developing probabilistic reasoning techniques that deal with uncertainty in a principled manner; Yoshua Bengio, Geoffrey Hinton, and Yann LeCun (2019) for making “deep learning” (multilayer neural networks) a critical part of modern computing; and finally, Richard Sutton, Andrew Barto (2024) for pioneering reinforcement learning in which agents learn by maximizing reward via trial and error in which we will implement FFNs, CNN, RNNs, and GPTs.

A.7 Bibiliographic Notes

• The primary resources used throughout this chapter, as listed on the sidenotes, were (Weizenbaum 1966), a short paper on ELIZA; (Minsky 1968)’s anthology of dissertations from his doctoral students; (Winograd 1971) and (Woods 1972), doctoral dissertations for SHRDLU and LUNAR (the latter is a BBN report, based off his PhD thesis)
• For a historical perspective, refer to (Gardner 1958) and (Nilsson 2010)’s Part II (Chapter 6, 7), Part III (Chapter 11, 13), Part IV (Chapter 18), Part VII (Chapter 26);
• For classical natural language processing techniques including syntactic, semantic, and logical analysis, refer to (Eisenstein 2018)’s Chapter 9, 10, 11, 12, 13 and (Jurafsky and Martin 2026)’s Appendix E, F, G, H and (Russel and Norvig 2022)’s Part III (Chapter 7, 8, 9, 10) and Part VI (Chapter 24).

0.8 Problems

Intermezzo One: The Language of Sets, Functions, Logic

Fixed-Size DataA modified excerpt from How to Design Programs (Felleisen et al., 2014) Intermezzo 1: Beginning Student Language. Chapter 0 deals with BSL formal language theory and set theory as if it were a natural language. It introduces the “basic words” of the language (which in turn, models natural language), suggests how to compose “words” into “sentences,” and appeals to your knowledge of algebra sets as a collection of objects for an intuitive understanding of these “sentences.” While this kind of introduction works to some extent, truly effective communication requires some formal study.

Although these deterministic and finitely discrete methods of software 1.0 run into the Bitter Lesson — that is, describing a reality with too many parts to count — we will need the machinery that underlies such techniques, namely that of sets, functions, logic because the same machinery is the common unifying foundation for all of mathematics, including the set of stochastic and infinitely continuous mathematics we need for our journey up ahead with software 2.0. While we’ve introduced such languages of set theory, functions, and logic throughout chapter 0 by appealing to your intuition as a programmer, as per the opening exercept, truly effective communication and understanding requires the formalformal study of mathematics using the axiomaxiomatic method.

Historically speaking, (elements of euclid…)

The good news as a programmer is that you have an advantage to learning mathematics, because the essence of both activities are in fact one and the same with the Univalent FoundationsSee https://en.wikipedia.org/wiki/Univalent_foundations, and https://ncatlab.org/nlab/show/univalent+foundations+for+mathematics.

Let us begin.
You will understand in due time.

/-- Doubles a natural number. -/
def double (n : Nat) : Nat := n + n

theorem double_eq (n : Nat) : double n = 2 * n := by
  simp [double, Nat.two_mul]

#check And
#check Or
#check Or

#eval double 21

variable {α : Type*}
variable (s t u : Set α)
open Set

example (h : s ⊆ t) : s ∩ u ⊆ t ∩ u := by
  rw [subset_def, inter_def, inter_def]
  rw [subset_def] at h
  simp only [mem_setOf]
  rintro x ⟨xs, xu⟩
  exact ⟨h _ xs, xu⟩

example (h : s ⊆ t) : s ∩ u ⊆ t ∩ u := by
  simp only [subset_def, mem_inter_iff] at *
  rintro x ⟨xs, xu⟩
  exact ⟨h _ xs, xu⟩

Note

Because we are dealing with the domain of language, we will denote alphabets of symbols and strings of symbols with $w\in V$ and $\mathbf{s}=(w1,w2,\cdots ,w_t),w_i\in V$ respectively rather than $a\in \Sigma$ and $\mathbf{a}=(a1,a2,\cdots ,a_t),a_i\in \Sigma$ to denote the fact that our alphabets and strings are modeling vocabularies and sentences of words within the domain of language. The alphabet and string formalism of formal language theory can be applied to other domains that admit sequences of tokens i.e biology with protein folding.

Chapter 0 deals with formal language theory and set theory as if it were a natural language.

It introduces the “basic words” of the language, suggests how to compose “words” into “sentences,” and appeals to your knowledge of algebra for an intuitive understanding of these “sentences.” While this kind of introduction works to some extent, truly effective communication requires some formal study.

B. From Classical to Constructive Mathematics

In which we historically retrace the development of the foundations of mathematics, from the geometry of ancient greek mathematics, to the logicism and formalism of Göttingen, and finally, to the intuitionism and constructivism of the Valley.

B.1 Ancient Mathematics

B.2 Classical Mathematics

B.3 Constructive Mathematics

C. From Sequential to Parallel Processors

In which we historically retrace the development of processors from sequential to parallel processors

Bibliography

Abelson, H., & Gerald Jay Sussman. (1996). Structure and Interpretation of Computer Programs. MIT Press. https://mitp-content-server.mit.edu/books/content/sectbyfn/books_pres_0/6515/sicp.zip/index.html

Amanzhol Salykov. (2025, January 12). Advanced Matrix Multiplication Optimization on NVIDIA GPUs. Salykova. https://salykova.github.io/sgemm-gpu

Ansel, J., Yang, E., He, H., Gimelshein, N., Jain, A., Voznesensky, M., Bao, B., Bell, P., Berard, D., Evgeni Burovski, Chauhan, G., Anjali Chourdia, Constable, W., Alban Desmaison, DeVito, Z., Ellison, E., Feng, W., Gong, J., Gschwind, M., & Hirsh, B. (2024). PyTorch 2: Faster Machine Learning Through Dynamic Python Bytecode Transformation and Graph Compilation. ACM, ASPLOS’2024. https://doi.org/10.1145/3620665.3640366

Bach, F. (2024). Learning Theory from First Principles. MIT Press. https://www.di.ens.fr/~fbach/ltfp_book.pdf

Bakhvalov, D. (2024). Performance Analysis and Tuning on Modern CPUs. https://github.com/dendibakh/perf-book

Bertsekas, D. P., & Tsitsiklis, J. N. (2008). Introduction to Probability. Athena Scientific.

Bryant, R. E., & O’hallaron, D. R. (2016). Computer Systems : a Programmer’s Perspective. Pearson.

Boehm, S. (2022, December 31). How to Optimize a CUDA Matmul Kernel for cuBLAS-like Performance: a Worklog. Siboehm.com. https://siboehm.com/articles/22/CUDA-MMM

Bright, P., Edelman, A., & Johnson, S. G. (2025). Matrix Calculus (for Machine Learning and Beyond). ArXiv.org. https://arxiv.org/abs/2501.14787

Chan, S. H. (2021). Introduction to Probability for Data Science. Michigan Publishing. https://probability4datascience.com/

Chen, T., Moreau, T., Jiang, Z., Zheng, L., Yan, E., Cowan, M., Shen, H., Wang, L., Hu, Y., Ceze, L., Guestrin, C., & Krishnamurthy, A. (2018). TVM: An Automated End-to-End Optimizing Compiler for Deep Learning. https://doi.org/10.48550/arxiv.1802.04799

Dao, T. (2023, July 17). FlashAttention-2: Faster Attention with Better Parallelism and Work Partitioning. ArXiv.org. https://arxiv.org/abs/2307.08691

Dao, T., Fu, D. Y., Ermon, S., Rudra, A., & Ré, C. (2022, June 23). FlashAttention: Fast and Memory-Efficient Exact Attention with IO-Awareness. ArXiv.org. https://doi.org/10.48550/arXiv.2205.14135

Darve, E., & Wootters, M. (2021). Numerical Linear Algebra with Julia. SIAM. https://ericdarve.github.io/NLA

Demmel, J. W. (1997). Applied numerical linear algebra. Society For Industrial And Applied Mathematics.

Dongarra, J., J. Du Croz, Sven Hammarling, & Duff, I. S. (1990). A Set of Level 3 Basic Linear Algebra Subprograms. ACM Transactions on Mathematical Software, 16(1), 1–17. https://doi.org/10.1145/77626.79170

Dongarra, J., J. Du Croz, Sven Hammarling, & Hanson, R. J. (1988a). Algorithm 656: An Extended Set of Basic Linear Algebra Subprograms: Model Implementation and Test Programs. ACM Transactions on Mathematical Software, 14(1), 18–32. https://doi.org/10.1145/42288.42292

Dongarra, J., J. Du Croz, Sven Hammarling, & Hanson, R. W. (1988b). An Extended Set of FORTRAN Basic Linear Algebra Subprograms. ACM Transactions on Mathematical Software, 14(1), 1–17. https://doi.org/10.1145/42288.42291

Fisler, K., Krishnamurthi, S., Lerner, B. S., & Politz, J. G. (2025). A Data-Centric Introduction to Computing. Dcic-World.org. https://dcic-world.org/

Fog, A. (n.d.). Software Optimization: C++ and Assembly. Agner.org. https://agner.org/optimize/

Goodfellow, I., Bengio, Y., & Courville, A. (2016). Deep Learning. The MIT Press. https://www.deeplearningbook.org/

Gordić, A. (2025, October 29). Inside NVIDIA GPUs: Anatomy of high performance matmul kernels - Aleksa Gordić. Aleksagordic.com. https://www.aleksagordic.com/blog/matmul

Goto, K., & Geijn, R. A. van de. (2008). Anatomy of High-Performance Matrix Multiplication. ACM Transactions on Mathematical Software, 34(3), 1–25. https://doi.org/10.1145/1356052.1356053

Goto, K., & Van De Geijn, R. (2008). High-Performance Implementation of the Level-3 BLAS. ACM Transactions on Mathematical Software, 35(1), 1–14. https://doi.org/10.1145/1377603.1377607

Güneş Baydin, A., Pearlmutter, B., Siskind, J., Baydin, G., Radul, A., & Mark, J. (2018). Automatic Differentiation in Machine Learning: a Survey. Journal of Machine Learning Research, 18, 1–43. https://www.jmlr.org/papers/volume18/17-468/17-468.pdf

Hastie, T., Tibshirani, R., & Friedman, J. (2009). The Elements of Statistical learning, Second Edition: Data mining, inference, and Prediction (2nd ed.). Springer. https://hastie.su.domains/ElemStatLearn/

Harris, C. R., Millman, K. J., van der Walt, S. J., Gommers, R., Virtanen, P., Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N. J., Kern, R., Picus, M., Hoyer, S., van Kerkwijk, M. H., Brett, M., Haldane, A., del Río, J. F., Wiebe, M., Peterson, P., & Gérard-Marchant, P. (2020). Array Programming with numpy. Nature, 585(7825), 357–362. https://doi.org/10.1038/s41586-020-2649-2

Hennessy, J. L., Patterson, D. A., & Christos Kozyrakis. (2025). Computer Architecture. Morgan Kaufmann.

Hwu, W.-M. W., Kirk, D. B., & Hajj, I. E. (2022). Programming Massively Parallel Processors: A Hands-on Approach. Morgan Kaufmann.

James, G., Witten, D., Hastie, T., Tibshirani, R., & Taylor, J. (2023). An Introduction to Statistical Learning. Springer. https://www.statlearning.com/

Jurafsky, D., & H. Martin, J. (2026). Speech and Language Processing. Stanford.edu. https://web.stanford.edu/~jurafsky/slp3/

Klein, P. N. (2013). Coding the Matrix : Linear Algebra Through Applications to Computer Science. Newtonian Press. https://codingthematrix.com/

Kwon, W., Li, Z., Zhuang, S., Sheng, Y., Zheng, L., Cody Hao Yu, Gonzalez, J. E., Zhang, H., & Stoica, I. (2023). Efficient Memory Management for Large Language Model Serving with PagedAttention. https://doi.org/10.1145/3600006.3613165

Lambert, N. (2026). RLHF Book. Rlhfbook.com. https://rlhfbook.com/

Lawson, C. L., Hanson, R. J., Kincaid, D. R., & Krogh, F. T. (1979). Basic Linear Algebra Subprograms for Fortran Usage. ACM Transactions on Mathematical Software, 5(3), 308–323. https://doi.org/10.1145/355841.355847

Marc Peter Deisenroth, A Aldo Faisal, & Cheng Soon Ong. (2020). Mathematics for Machine Learning. Cambridge University Press. https://mml-book.github.io/book/mml-book.pdf

Minsky, M., Papert, S., & Léon Bottou. (2017). Perceptrons : An Introduction to Computational Geometry. The MIT Press.

Matthias Felleisen, Robert Bruce Findler, Flatt, M., & Shriram Krishnamurthi. (2018). How to Design Programs: An Introduction to Programming and Computing. The MIT Press. https://htdp.org/

Mitchell, T. M. (1997). Machine learning. Mcgraw-Hill. https://www.cs.cmu.edu/~tom/files/MachineLearningTomMitchell.pdf

Murphy, K. P. (2022). Probabilistic Machine Learning: An Introduction. MIT Press.

Myers, M., Van De Geijn, P., & Van De Geijn, R. (2021). Linear Algebra: Foundations to Frontiers. https://www.cs.utexas.edu/~flame/laff/laff/LAFF-2.00M.pdf

Nakatsukasa, Y. (n.d.). Numerical Linear Algebra. Retrieved March 17, 2026, from https://courses.maths.ox.ac.uk/pluginfile.php/105965/mod_resource/content/35/NLA_lecture_notes.pdf

Ng, A., & Ma, T. (2023). CS229 Lecture Notes. https://cs229.stanford.edu/main_notes.pdf

Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., & Bai, J. (2019). PyTorch: An Imperative Style, High-Performance Deep Learning Library. Neural Information Processing Systems. https://papers.nips.cc/paper_files/paper/2019/hash/bdbca288fee7f92f2bfa9f7012727740-Abstract.html

Ragan-Kelley, J., Barnes, C., Adams, A., Paris, S., Durand, F., & Amarasinghe, S. (2013). Halide. Proceedings of the 34th ACM SIGPLAN Conference on Programming Language Design and Implementation. https://doi.org/10.1145/2491956.2462176

Raschka, S. (2024). Build a Large Language Model (From Scratch). Manning.

Raschka, S. (2026). Build a Reasoning Model (From Scratch). Simon and Manning.

Roberts, D. A., Yaida, S., & Hanin, B. (2022). The Principles of Deep Learning Theory: An Effective Theory Approach to Understanding Neural Networks. Cambridge University Press. https://deeplearningtheory.com/

Russell, S., & Norvig, P. (2021). Artificial Intelligence: A Modern Approach (4th ed.). Prentice Hall. https://aima.cs.berkeley.edu/

Slotin, S. (n.d.). Algorithms for Modern Hardware - Algorithmica. En.algorithmica.org. https://en.algorithmica.org/hpc/

seb-v. (2025, January 20). Optimizing Matrix Multiplication on RDNA3: 50 TFlops and 60% Faster Than rocBLAS. Seb-v. https://seb-v.github.io/optimization/update/2025/01/20/Fast-GPU-Matrix-multiplication.html

Shah, J., Bikshandi, G., Zhang, Y., Thakkar, V., Ramani, P., & Dao, T. (2024, July 24). FlashAttention-3: Fast and Accurate Attention with Asynchrony and Low-precision. ArXiv.org. https://arxiv.org/abs/2407.08608

Shalizi, C. R. (n.d.). Advanced Data Analysis from an Elementary Point of View. https://www.stat.cmu.edu/~cshalizi/ADAfaEPoV/ADAfaEPoV.pdf

Shalizi, C. R. (2015). Modern Regression Lecture Notes. Cmu.edu. https://www.stat.cmu.edu/~cshalizi/mreg/15/

Shankhdhar, P. (2024, November 29). Outperforming cuBLAS on H100: a Worklog. Substack.com. https://cudaforfun.substack.com/p/outperforming-cublas-on-h100-a-worklog

Smith, T. M., Geijn, R. van de, Smelyanskiy, M., Hammond, J. R., & Zee, F. G. V. (2014, May 1). Anatomy of High-Performance Many-Threaded Matrix Multiplication. IEEE Xplore. https://doi.org/10.1109/IPDPS.2014.110

Spector, B., Singhal, A., Arora, S., & Re, C. (2024, May 12). GPUs Go Brrr. Stanford.edu. https://hazyresearch.stanford.edu/blog/2024-05-12-tk

Spector, B. F., Arora, S., Singhal, A., Fu, D. Y., & Ré, C. (2024, October 27). ThunderKittens: Simple, Fast, and Adorable AI Kernels. ArXiv.org. https://arxiv.org/abs/2410.20399

Sul, S., & Ré, C. (2026). ThunderKittens 2.0: Even Faster Kernels for Your GPUs. Stanford.edu. https://hazyresearch.stanford.edu/blog/2026-02-19-tk-2

Spector, B., Juravsky, J., Sul, S., Dugan, O., Lim, D., Fu, D., Arora, S., & Ré, C. (2025). Look Ma, No Bubbles! Designing a Low-Latency Megakernel for Llama-1B. Stanford.edu. https://hazyresearch.stanford.edu/blog/2025-05-27-no-bubbles

Spector, B., Juravsky, J., Sul, S., Lim, D., Dugan, O., Arora, S., & Ré, C. (2025). We Bought the Whole GPU, So We’re Damn Well Going to Use the Whole GPU. Stanford.edu. https://hazyresearch.stanford.edu/blog/2025-09-28-tp-llama-main

Stillwell, J. (2010). Mathematics and its history. Springer New York.

Strang, G. (2023). Introduction to Linear Algebra. Wellesley-Cambridge Press. https://math.mit.edu/~gs/linearalgebra/

Sutton, R. S., & Barto, A. (2018). Reinforcement learning: An introduction (2nd ed.). The MIT Press. http://incompleteideas.net/book/the-book-2nd.html

Tillet, P., Kung, H.-T., & Cox, D. G. (2019). Triton: An Intermediate Language and Compiler for Tiled Neural Network Computations. https://doi.org/10.1145/3315508.3329973

Trefethen, L. N., & Bau, D. (1997). Numerical Linear Algebra. Society For Industrial And Applied Mathematics.

Valiant, L. (2014). Probably Approximately Correct: Nature’s Algorithms for Learning and Prospering in a Complex World. Basic Books, A Member Of The Perseus Books Group.

Wasserman, L. (2010). All of Statistics. Springer.

Zheng, L., Yin, L., Xie, Z., Sun, C., Huang, J., Yu, C. H., Cao, S., Kozyrakis, C., Stoica, I., Gonzalez, J. E., Barrett, C., & Sheng, Y. (2023). SGLang: Efficient Execution of Structured Language Model Programs. ArXiv.org. https://arxiv.org/abs/2312.07104

Zadouri, T., Hoehnerbach, M., Shah, J., Liu, T., Thakkar, V., & Dao, T. (2026). FlashAttention-4: Algorithm and Kernel Pipelining Co-Design for Asymmetric Hardware Scaling. ArXiv.org. https://arxiv.org/abs/2603.05451