Structure and Interpretation

of Tensor Programs

First Edition

A Whirlwind Tour to Deep Learning Alchemy and Deep Learning Systems
Run 🤖nanochatby building borschtfrom scratch: a hackable subset of 🔥torchWhat’s in a name? Borscht is a sour soup with meat from Eastern Europe, most often associated with the variant originating in Ukraine with red beetroots . It’s a permutation-like pun on torch, given how the deep learning framework borscht aims to be a rough subset of torch. You want a torch.Tensor? You get the borscht.Tensor! Also, my mother in law’s borscht is to die for: Мне очень нравится ваш борщ, миссис Ким)

A new breadth of mathematical computer science for the world of artificial intelligence
Made with 🖤🪻 by Jeffrey Zhang, University of Waterloo (BMath)
Made possible by Lambda Labs Research Grant

Citation:

@book{zhang2026sitp,
  author = "Jeffrey David Zhang",
  title = "Structure and Interpretation of Tensor Programs",
  year = 2027,
  url = "https://sitp.ai"
}

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 and course is but a fraction of the love he gave me.
May you rest in pure land. I am happy that it was in our karma to meet and learn during this lifetime.

《定风波·莫听穿林打叶声》（元丰五年）　苏东坡　　大山及友人 英译
“Calming the Waves” (1082) by Su Dongpo, tr. Dashan & Friends

Course Thesis

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

The Structure and Interpretation of The AI Curriculum

My alma matter was amongst those which took the SICP approach, and as intended, for someone coming into first year college have taken some elementary computer science courses from high school, it blew my mind (Lambda Calculus!). 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 for some time, 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 domain specific tensor compilers from TORCH to TRITON.It turns out there is effectively zero overlap except that they are both compilers in that they take some source program and produce some target program.

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 personally forgot lots of foundation from collegeI had the definition of a gradient and basis in the basement of my subconscious at best. Thank you Professor.West for sharing your juggling, and Professor. Wolczuk for sharing your Ukrainian floor dancing, and I wasn’t sure how to bridge myself to industrial deep learning systems like tinygradHere’s a PR I was able to land in tinygrad making logsumexp numerically stable. It felt like I was contributing to LLVM without having taken a course on compiler construction to learn the basics of parsing, optimization, and lowering., torch, jax, vllm, and sglang.

Shortly after, I decided to take the plunge and started drinking from the firehose of deep learning canon: Hastie et al., Murphy, Goodfellow et al., you name it. 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’s 10-414 Carnegie Mellon and Alexander Rush’s CS5781 at Cornell, I personally would have enjoyed a more 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?

Pur shortly, SITP, is a function that maps learners from the product type of an introductory computer science book such as DCIC (Fisler et al. 2025) and an introductory systems programming book such as TRPL (Klabnik et al. 2026) with experimental modifications (Crichton et al. 2023) to a frontier model report such as the frontier technical report of The Llama 3 Herd of Models (Meta 2024) and frontier systems book such as Hugging Face’s Ultra Scale Book (Tazi et al. 2025).

$$\begin{array}{rl}\text{SITP}& :\text{DCIC}\times \text{TRPL}\to \text{LLAMA3}\times \text{USPB}\\ \text{SITP}& :=\text{HTDP}(\text{PTDS},\text{LADR},\text{NA},\text{NO},\\ & \text{SLP},\text{DL},\text{RL},\\ & \text{TAOHPC},\text{TAODCP},\text{PMPP})\end{array}$$

where the function $\text{HTDP}(\cdot )$ is the systematic and principled design of a curriculum implemented by chimera-like professors and lecturers whom are usually interested and capable in the intersection of both programming language theory (PLT) and computer science education (CSE). In the same way that we use non-intersective adjectives in the distinctions of fake gun vs gun, non-tight language model vs language model, there is also the teaching professor vs professor“I think of the position as having “twice” the teaching and “half” the research as an assistant professor (of course individual approaches can and do differ!).” – @JoePolitz, September 5, 2024. That is, SITP applies the form of HTDP developed by such teaching professors to the substance of deep learning. What then, is such a form?

The How to Design Programs HTDP text This was the actual text used at my alma matter. Thank you to Professor Ragde for bringing it over, and to Professor Vasiga for bringing it to life. was also seminal in that it was the first computer science text to treat curriculum design as scientific research and an engineering problem (Felleisen et al. 2004), building off the shoulders of the giant that SICP (Abelson and Sussman 1996) is.

First, the book discusses explicitly how programs should be constructed. Second, to tame the complexity of programming, it defines a series of teaching languages based on Scheme that represent five distinct knowledge levels through which students pass during their first course. The levels correspond to the complexity of data definitions that the program design guidelines use. Third, the book uses exercises to reinforce the explicit guidelines on program design; few, if any, exercises are designed for the sake of domain knowledge. Finally, the book uses more accessible forms of domain knowledge than SICP. Because of this shift in emphasis, we gave our book the title How to Design Programs (HTDP).

the concreteness fading (Fyfe et al. 2014) of declarative concepts in mathematics (Hestenes et al. 1992) by presenting notions just in time by defining them by their context, and iteratively refining (progressive disclosure) such concepts from the informal to formal ending in HtDP-style IntermezzosWhich present the formal syntax and semantics of a language..

that transfer (Bransford and Schwartz 1999) to procedural skills of computation (Price et al. 2021) in the context of deep learning. Such transfer is possible because the computational discipline of deep learning has lots of structural similarity with mathematics, and is why we have books such as Deisenroth et al. (2020), and subsequent courses follwing such book.

When it comes to the discipline of deep learning, there are two large jumps of transfer a learner must make. The first one is that a learner must learn the mathematics for machine learning, while although elementary, has a non-trivial breadth of probability theory, linear algebra, differential calculus, and if the learner has the time, to deepen such concepts with information theory, measure theory, and analysis. One of the primary texts that seeks to provide a unified treatment in such mathematics for machine learning is (Deisenroth et al. 2020)“For instance, Princeton’s COS302 offered by Ryan Adams, which turns out to be the lab that created created HIPS/autograd (led by Matthew Johnson, Dougal Maclaurin, and David Duvenaud), which inspired PyTorch! See https://soumith.ch/blog/2023-12-17-pytorch-design-origins.md.html#/origins., which provides an excellent two part structure with the foundational mathematics, followed by the applied. The core motivation is that any gaps in mathematics can be filled, preparing a learner for a text like such (Goodfellow et al. 2016) to deepen their understanding in neural network architecture.

More recently, with the advent of scaling large language models (Kaplan et al. 2020) and the slowing of Moore’s Law (Moore 1965), another jump of transfer the learner must make is how to efficiently map their architecture and optimizers to massively parallel hardware at the level of a single graphics processing unit (GPU) and multi GPU. The two canonical texts for this are (Hwu et al. 2022) for systems, and more recently (Tazi et al. 2025)“Alternatively, JAX Scaling Book (link.).

As someone opening up PMPP, it was extremely frustrating to read the appendix that briefly covered the process of deep learning without a proper foundation in the modeling of architecture and optimizers.

When I realized the two large jumps I was making from re-learning foundational mathematics to frontier architecture, and from modeling frontier architecture to accelerating frontier architecture, my “aha” moment was realizing the vast bridges a learner was being asked to cross, and realizing Shriram’s taste of finding some intermediary point to smoothen such a progression. This was the motivation for the explicit Design Recipe in HtDP to go from blank page to a well-designed program, and more recently, property-based testing to go from unit tests to model checking and SAT solvers. To return to the short description, SITP that maps learners from the product type of an introductory computer science book and an introductory systems programming book to a frontier model report such as the frontier technical report and frontier systems book. That is,

$$\text{SITP}:\text{DCIC}\times \text{TRPL}\to \text{LLAMA3}\times \text{USPB}$$

unity/intertwined mathematics and programming

• Evaluation (Benchmarks): Stochastics of probability theory $(\Omega ,\mathcal{F},\mathbb{P})$
• Specification (Architecture): Dimensionality of linear algebra $(V,⟨\cdot ,\cdot ⟩)$
• Implementation (Optimizers): Approximation of differential calculus $\frac{d\mathcal{L}}{dw}$, $\nabla \mathcal{L}$, $\mathbb{J}$, $\mathbb{H}$
• Computation (Systems): Rooflines of processors $\text{max}(compute,communication)$

For instance, the traditional ordering one might approach to the discipline of deep learning is to synthetically deduce definitions starting from what is considered to be elementary (i.e your choice of foundations which you pay lip service to as a working mathematician: ZFC or DTT). The benefits of this approach is that the transfer (Bransford and Schwartz 1999) of concepts is maximized — afterall, studying mathematics is simply generic programming. While this is no doubt the standard for most mathematical texts, this is what was special about the textbook How to Design Programs, an explicitly curriculum designed and engineered for transfer (Felleisen et al. 2004) expositing programming based off the principles of set theory (whatever that may mean).

all in the context of autoregressive language modeling, culminating in the transformers architecture

Jeffrey Zhang
Waterloo, Ontario
August 2026

Course Design Thesis

“Civilization advances by extending the number of important operations which we can perform without thinking of them.” Thanks, Alfred North Whitehead (1861-1947)

> 7:22AM Let me first look at what's in the working directory to scope this.
> Thought for 432s

> 7:30AM The directory is empty and lacks numpy, torch, and pytest,
  so I'll build "nanotorch" — a NumPy-based PyTorch clone with a
  familiar API (Tensor, autograd, nn, optim, data loading, saving/loading), 
  testing against real PyTorch where possible. Starting with environment setup.
> Bash Check CPU count, memory, and Python package tooling
> Machine: 8-core M3 with 8GB RAM, so a workflow can run about 6 agents at once.
  Next I'll create a virtualenv with numpy and pytest, plus real PyTorch to use as the reference in tests.

> Bash Create virtualenv and install numpy and pytest
> Bash Install real PyTorch as a test oracle
> Thought for 2s

> Real PyTorch 2.14.1 is installed as the oracle.
  Before writing the design spec, I'll check a few numpy behaviours the core will depend on.
> Bash Verify numpy reshape/promotion features and torch promotion
> Thought for 607s

> Environment's set with numpy 2.5.4, pytest 9.1.1, and PyTorch 2.14.1 as the oracle.
  Now I'll verify the PyTorch behaviours against the real library before drafting the design contract.
> Bash Probe PyTorch semantics for the design spec
> These results match the spec I had planned. Next I'll write the design contract (docs/DESIGN.md)
  and pyproject.toml, and make the package importable from the venv.
> Write pyproject.toml
> Bash Make repo importable in venv and confirm torch sources exist

> 7:40AM Now the design contract. It's long on purpose: it pins down the core semantics
  and file ownership so parallel agents produce one coherent library.
> Write DESIGN.md
> Thought for 233s

> With the design contract done, I'll now build the core foundation
  first—Tensor, dtypes/type promotion, autograd, views/in-place ops, indexing,
  creation, random, test helpers, and the package skeleton—then subject
  it to five red-team reviews against real PyTorch with fixer rounds until tests pass,
  since everything else depends on this layer.
> Thought for 141s

> 7:50AM Workflow 1 is running in the background: one core implementer,
  then 5 red-team reviewers in parallel, then up to 3 fixer rounds.
  The core is the largest single piece, so this will take a while.
  I'll be notified when it finishes, then review the core myself
  before launching the 9 parallel op/nn families of Workflow 2.

“Everything is vague to a degree you do not realize till you have tried to make it precise.” Thanks, Bertrand Russel (1872-1970)

DESIGN.md

# nanotorch — design and contributor contract

nanotorch is a from-scratch reimplementation of the PyTorch Python API on top of NumPy.
Code written for PyTorch should run unchanged with `import nanotorch as torch` for the
supported subset and produce the same results (up to floating-point rounding and RNG streams).

* **We implement:** tensor semantics (dtypes and type promotion, broadcasting, views,
  in-place ops, indexing), the autograd engine (dynamic graph, higher-order gradients,
  hooks, version counters), every derivative formula, `nn`, `optim`, `utils.data`,
  serialization, printing.
* **NumPy provides:** strided n-d array storage and vectorised kernels (ufuncs, BLAS matmul).

This document is the contract every contributor (human or agent) follows. Section 5 (op
authoring guide) is finalised by the core implementation and must stay in sync with the code.

## 0. Ground rules

1. **Fidelity target: PyTorch 2.x.** The oracle is the real `torch` 2.14.1 installed in
   `.venv`. When in doubt about a behaviour, run torch and match it: names, signatures
   (positional order, keyword names, defaults), output shapes, output dtypes, exception
   *types* (message text should be close but need not be identical).
   torch's own Python sources (`.venv/lib/python3.14/site-packages/torch/...`, e.g.
   `nn/modules/module.py`, `optim/adam.py`, `_tensor_str.py`, `nn/init.py`) are a legitimate
   reference for semantics. Write clean, minimal implementations; when you substantially
   adapt PyTorch code, add a `# Adapted from PyTorch (BSD-3-Clause)` comment.
2. **nanotorch never imports torch.** Only tests import torch (through `pytest.importorskip`).
3. **NumPy (>= 2.1) is the only runtime dependency.** No scipy. Pure-Python fallbacks are
   fine for special functions (implement accurate vectorised approximations).
4. Python >= 3.10 syntax (we run 3.14). Type hints welcome, not required.
5. **Non-goals** (do not implement, do not test): CUDA/MPS execution (CPU only; `device`
   objects for other types can be constructed but allocating on them raises), complex
   dtypes, bfloat16 (placeholders that raise a clear error on use), sparse / quantized /
   named / nested tensors, `torch.compile` / `jit` / `fx` (no-op shims only), distributed,
   multi-process DataLoader workers (run in-process), ONNX, `vmap`, FFT, exact RNG-stream
   parity with PyTorch (same *distributions*, different numbers).
6. Repr strings mirror PyTorch exactly where reasonable, including the `torch.` prefix:
   `torch.float32`, `torch.Size([2, 3])`, `device(type='cpu')`,
   `tensor([1., 2.], requires_grad=True)`, `<MulBackward0 object at 0x...>`.

(i.e Why another deep learning book?)

Learning Outcomes (Tasks) learning probes, observable, signals, channels for feedback interactive learning -> no longer producing code (hieroglyphics)

concept inventory (Hestenes et al. 1992)

Software 2.0, following (Krishnamurthi and Fisler 2020)

• Lectures
• Readings (with concreteness fading following (Fyfe et al. 2014)) (with interactive quizzes following (Crichton and Krishnamurthi 2024)) with SITP
  • Declarative concepts
    • property based testing
    • tradeoff between x and y
    • system design/system invariants
• Labs (with Agentic Coding following (Fisler et al. 2026))
  • Procedural skills
    • write a good design doc/test suite
    • code review
• centered around the Teaching Array Programming Language borscht.Tensor with language levels following the form of (Felleisen et al. 2004) but with the substance of NETSOR and NETSOR$\top$ from (Yang 2019) and (Yang 2020) as well as PTX, CUDA, CuTeDSL, Cutile.

Future Directions

• integrate a trimmed version of such curricula with highschool mathematics, following (Felleisen and Krishnamurthi 2009) and Guzdial
• collect profile data to suggest pedagogic interventions to imporove learner outcomes
• actually encapsulate language levels like racket-like #lang with julia metaprogramming
• map out the semantic space of language features with a Standard Model of APLs following (Krishnamurthi 2008)

Syllabus

Any sufficiently advanced technology is indistinguishable from magic. Thanks, Arthur C. Clarke (1917-2008)

All text (every single word) has been written by a human, for a human, as this book is not a reference book, but rather a cognitive helmet of understanding for a human. You are highly encouraged if not expected to read with an LLM, but be aware (that is, to beware) of how much you cognitively offload and outsource your understanding.

Week 1. Hello Tensor 👋

• Watch ▶️ Lecture 1
• Read 📖 Chapter 1.1 From Certain to Uncertain Knowledge,
• Read 📖 Chapter 1.2 Next Token Prediction is Classification and Compression
• Write Lab 1. Elements of Tensors

Optional Reading

• Read An Example of Statistical Investigation of the Text Eugene Onegin Concerning the Connection of Samples in Chains (Markov 1913)
• Read A Mathematical Theory of Communication (Shannon 1948)
• Read A Proposal for the Dartmouth Summer Research Project on Artificial Intelligence (McCarthy et al. 1955)

Week 2. I Perceive You 👀

• Watch ▶️ Lecture 2
• Read 📖 Chapter 1.3 Parameterizing Classification with Discriminant Functions
• Read 📖 Chapter 1.4 Parameterizing Classification with Probabilistic Models
• Write Lab 2.

Optional Reading

• Read The Perceptron: A Probabilistic Model for Information Storage and Organization in the Brain (Rosenblatt 1958)
• Read Statistical Modeling: The Two Cultures (Breiman 2001)
• Read The Epic Story of Maximum Likelihood (Stigler 2008)

Week 3. MetaCircularEvaLuAtor! ↩️

• Watch ▶️ Lecture 3
• Read 📖 Chapter 2.1 Virtualizing Shapes with Strides
• Read 📖 Chapter 2.2 Accelerating Basic Linear Algebra on CPUs
• Read 📖 Chapter 2.3 From Performance Primitives to Linear Solvers
• Write ✍️ Lab 3. Tensor++ and Matrix Decompositions

Optional Reading

• Read A Set of Level 3 Basic Linear Algebra Subprograms (Dongarra et al. 1990)
• Read Anatomy of High-Performance Matrix Multiplication (Goto and van de Geijn 2008a)
• Read The History of Fortran I, II, and III (Backus 1978)
• Read A History of MATLAB (Moler and Little 2020)
• Read Array Programming with NumPy (Harris et al. 2020)

Week 4. GPUs Go Brrr 🏎️

• Watch ▶️ Lecture 4
• Read 📖 Chapter 2.4 Accelerating Matrix Multiplication on GPUs
• Write ✍️ Lab 4. Accelerating Matmul with Tiles and Layouts

Week 5. King - Man + Queen = 👸

• Watch ▶️ Lecture 5
• Read 📖 Chapter 3.1 Learning Static Representations with FFNs
• Read 📖 Chapter 3.2 Iterative Optimization
• Write Lab 5. Embedding Projector

Week 6. Plato’s Cave 𐦂𖨆𐀪𖠋

• Watch ▶️ Lecture 6
• Read 📖 Chapter 3.3 Learning Representations with CNNs
• Read 📖 Chapter 3.4 Learning Representations with RNNs
• Read 📖 Chapter 3.5 Learning Representations with GPTs
• Write ✍️ Lab 6. Vision Language Action Model

Week 7. AutoNabla (≧∇≦)

• Watch ▶️ Lecture 7
• Read 📖 Chapter 4.1 Automatic Differentiation
• Read 📖 Chapter 4.2 Gradient Based Optimization
• Write ✍️ Lab 7. Custom Autograd Backwards for VLA

Week 8. Don’t Think 🤔

Week 9. Just Do 🏃‍♂️

Week 10. Compilers, They’re InTense 😰

Frontispiece to Galileo’s Dialogue Concerning the Two Chief World Systems (Stefano della Bella 1632)

I. Elements of Networks

No different from learning another mathematical or programming language, they must transition from their finitely discrete structures and deterministic procedures 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 on serial machines to optimizing the distributions of scalars, vectors, matrices, tensors, and neural networks on massively parallel machines. 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 feeding a dataset to a neural net architecture that has a goal, and searching the space of programs with compute. How to exactly program with this approach will take the remainder of this 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 Translate and Apple 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.

WeA modified excerpt from The Structure and Interpretation of Computer Programs §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.

So 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 torch by implementing borscht. Whether you’re an eager high school student, a rising college student, or a battle-tested industrial programmer, SITP has been meticulously designed so that the only prerequisite required is a basic familiarity with the elements of programming, and high school algebra. Any additional experience is helpful, not mandatory. If you’d like to revisit your programming fundamentals, we recommend A Data-Centric Introduction to Computing (Fisler et al. 2025) (DCIC) which is a great introduction to computing with a focus on data scienceSee (Krishnamurthi and Fisler 2020) that starts with the Pyret teaching language and the table data structure (the torch.Tensor in disguise!), and smoothly graduates to Python, the starting point of our adventure herein.

I. Table of Contents

La perspective curieuse (Jean-François Niceron 1638)

1. Self-Supervised Sequence Learning with Single-Layer Networks

In which we transition from the deterministic and discretely finite software 1.0 to the stochastic and infinitely continuous software 2.0 by rediscovering ngrams by Markov and Shannon as well as perceptrons by Rosenblatt using the languages of probability theory and linear algebra afforded by the borscht.Tensor. At the end, we will have rediscovered the early period where researchers in artificial intelligence were uncertain between methods of connectionism and logical symbolism.

1.1 From Certain to Uncertain Knowledge

$\uparrow$ Table of Contents §1.2 $\to$

In which we have our first encounter with the language modeling of GPT-2 via borscht; the primary notions from probability and linear algebra thus used in modeling next-token prediction, namely random variables, their distributions, and their computational equivalent of the multidimensional array.

1.1.1 From Software 1.0 to Software 2.0

$\hookleftarrow$ Table of Contents (§1.1 From Certain to Uncertain Knowledge)

Although ChatGPT, Claude and friends are relatively new to our cosmos, the idea of generating sentences with next token prediction is surprisingly not new and dates back to the work of Russian mathematician Andrey Andreyevich Markov in 1913, and shortly after Claude Shannon in 1948. So why is humanity’s so-called tech tree late to such technology? Predominantly for two reasons.

Philosophically, because of the eternal tension between the discretediscrete and continuouscontinuous methods in describing our reality. Programmers were reluctant to use stochastic and infinitely continuous techniques, favoring those that were logical and finitely discrete. However, … like a physicist position of every particle in a vacuum, even with a set of initial equations for position and momentum with equations for change, describing reality with too many parts to countbitter lesson. This created the need for statistical mechanics, describing particles with probability distributionsprobability distributions.

Practically, it’s because of the fossil-fuel like subsidy of datadata provided by the aforementioned intergalactic computer network we call the web, the computecompute provided by massively parallel processors originally designed for video games we call graphic processing units (henceforth GPUs), which can be used efficiently by a neural network architecturearchitecture called transformer. This is why ChatGPT and Claude are referred to as large language modelslarge language model.

Large Language Models explained briefly (3Blue1Brown 2024)

Large language models are trained using methods from the discipline of deep learningdeep learning, which in turn, are based in the statistical machine learningmachine learning approach to artificial intelligence. This means, in order to produce such a probability distribution over possible next words, ChatGPT, Claude and others use mathematical machinery from the areas of probability theory (clearly), linear algebra, and calculus.

That is, the testingtesting resulting from the evaluation of next-token prediction is stated using the stochasticitystochasticity of probability theory; the specificationspecification of an architecture’s parameterized inductive bias is stated using the dimensionalitydimensionality of linear algebra rather than the discreteness of combinatorics; the implementationimplementation of an architecture’s realized parameters is iteratively optimizedoptimized using the approximation of differential calculus rather than the control flow of structured programming; and the computationcomputation that carries out such an optimization is speed of light acceleration using the rooflines of processor pipelines and memory hierarchies rather than the more abstract cost model of worst case complexity. Put shortly, software 2.0 is built by programming continuous mathematics on highly parallelized and distributed machines — how to exactly program with this new approach will take the remainder of the book to iteratively deepen our understanding of evalsevals, architecturearchitecture, optimizersoptimizers, and systemssystems throughout SITP.

Starting with Part I. Elements of Networks, in §1. Self-Supervised Sequence Learning with Single-Layer Networks we introduce the goal or evaluation of next-token prediction and train our first architectures with linear models which correspond to a single layer of a neural net. This is followed by §2. The Massively Parallel Multidimensional Tensor where we start the implementation of our deep learning system borscht.

Then, in Part II. Neural Networks, particularly §3. Self-Supervised Sequence Learning with Deep Networks, we deepen our understanding of architecture by transitioning from linear models corresponding to single layer networks to non-linear models of multi layer networks, followed by §4. The Engine of Automatic Differentiation where we evolve the implementation of borscht to support automatic differentiation on gpu-accelerated neural network primitives. Finally, in Part III. Scaling Networks (expand).

We will introduce such mathematical primitives by keeping our language models simple at first in Part I. Elements of Networks — namely what is called the ngram model and linear models in which the aforementioned Markov and Shannon were working on around a century ago — before diving into the design of neural network architectures (including transformers with the attention operator) in Part II. Deep Neural NetworksPhysicists learn the standard model before general relativity. Similarly, we must learn the foundations of machine learning before deep learning..

1.1.2 From Variables to Random Variables

$\hookleftarrow$ Table of Contents (§1.1 From Certain to Uncertain Knowledge)

But rather than produce a distribution with two outcomes like a coin, six outcomes like a die, or fifty two outcomes like cards, it produces one for $|V|$ outcomes, where $V$ is the set of words in some vocabularyvocabulary. The size of GPT-2’s vocabulary is 50257, and the list of probabilities you see in Demo 1.1 are the top 10 most likely. As a first approximation, language models can be viewed as as an urn containing a ball labeled with each word in the vocabulary, where some balls are weighted heavier than others.

To be more precise, because we are passing an input sentence, such a distribution is a conditional distributionconditional distribution. That is, GPT-2 is producing the distribution of the next word conditioned on the sequence of words passed in as input, and is denoted by

$p(y\underset{\text{given}}{\underbrace{|}}\text{input sentence})$

and with Demo 1.1, evaluating GPT-2 is equivalent to $p(y|\text{Hello world, nice to})$. That is, language modelslanguage models are next-token predictors modeled by conditional distributions. More accurately, each number in that list of probabilities is the chance that GPT-2 assigns a random variablerandom variable taking on a specific outcome. A random variable is like a deterministic variable in that it can take on values, but rather than being assigned one specific value, it can possibly take on many with a certain degree of chance. A random variable’s distribution then corresponds to an array or list of values which we call a distribution because it distributes the truth or state of the random variable across many outcomes.

In the case of GPT-2, the conditional distribution $p(y|x^{(1)},\dots ,x^{(t)})$’s random variable is denoted as $Y|X^{(1)},\dots ,X^{(t)}$. To make random variables more explicit, some people will denote distributions with their random variables included, which in our case of a conditional distribution over next words, is $p(Y=y|X^{(1)}=x^{(1)},\dots ,X^{(t)}=x^{(t)})$. However, this brings us to our next point which is that with conditional distributions, the random variable being conditioned on is actually no longer random, because it is assumed that it has already taken on a value, in this case where $x^{(0)}=\text{Hello}$, $x^{(1)}=\text{world}$, $x^{(2)}=\text{nice},x^{(3)}=\text{to}$. So, with $p_{Y|X}(y|\mathbf{x})$ more generally, there is no randomness associated with the history $x^{(1)},\dots ,x^{(t)}$. Let us print the distribution over the next word that GPT-2 produces given the history $\text{Hello world, nice to}$ from Demo 1.1:

1.1_gpt2output_torch.ipynb

import torch

# The same GPT-2 output distribution p(w | "Hello world, nice to") from above, top 10 of 50257 words
# (renormalized so this top-10 slice sums to 1)
vocab = ["see", "meet", "hear", "have", "be", "know", "you", "talk", "say", "get"]
probs = torch.tensor([0.387361, 0.154993, 0.152304, 0.127857, 0.057572, 0.043027, 0.039115, 0.015646, 0.013935, 0.008190])

print("Asking GPT-2 what is p(y|Hello world, nice to):")
for i in torch.argsort(probs, descending=True):
    print(f"{vocab[i]} : {probs[i]:.4f}")
print(f"Total sum of p(y|Hello world, nice to): {probs.sum():.4f}")

Asking GPT-2 what is p(w|Hello world, nice to):
see : 0.3874
meet : 0.1550
hear : 0.1523
have : 0.1279
be : 0.0576
know : 0.0430
you : 0.0391
talk : 0.0156
say : 0.0139
get : 0.0082
Total sum of p(w|Hello world, nice to): 1.0000

1.1_gpt2output_torch.ipynb

import matplotlib.pyplot as plt

order = torch.argsort(probs, descending=True)
plt.figure(figsize=(6, 3))
plt.bar([vocab[i] for i in order], probs[order], color="#0072B2")
plt.ylabel("p(y|Hello world, nice to)")
plt.title("GPT-2 next-token distribution (top 10)")
plt.tight_layout()
plt.show()

Where the state of $Y|X$ is being spread across 50257 possible outcomes of GPT-2’s vocabulary. Each index i of the probs array corresponds to a value in the token outcomes tokens[i], with each probs[i] corresponding to the chance, or possibility, of the random variable taking on the value tokens[i]. That is, what connects the weighted array tokens and the state space probs is the index i. It is more accurate then, to describe a random variablerandom variable Stanford CS109 | Random Variables and Expectation | Lecture 6 as states represented by indices $i\in \mathbb{N}$, and the type of a distribution $p$ ends up being some function $p:\mathbb{N}\to [0,1]$ that sends states of indices to their probabilities, such that for all $n\in \mathbb{N},p(n)\ge 0$ and $\sum _{n\in \mathbb{N}}=1$. In general, any list of values between 0-1 and normalizes to 1 is a valid distribution. These properties are referred to as non-negativity and normalization respectively. We can validate that probs is a valid distribution by verifying that each probability is between 0-1 and the distribution normalizes to 1However, since we are only including the top 10 most likely words out of 50256, we’ve re-normalized probs so that the total sum of probs is 1. The original sum is 0.8181, with the other 1-0.8181=0.1819 spread across the other 50256-10=50246 unseen words.. Moreover, since the space of possible outcomes has size of 50257, the random variable $Y|X$ is said to be a discrete random variablediscrete random variable with a probability mass functionprobability mass function as opposed to a continuous random variable with a probability density functionprobability density function.

experimentexperiment trialtrial frequencyfrequency

g = torch.Generator()

n = 10_000_000
target_e = "see"

# p must normalize to 1, so we renormalize our top-10 slice of the distribution
p = probs / probs.sum()
# draw all n samples at once, according to p(y|"Hello world, nice to")
# (tensors hold numbers, not strings, so we sample vocab indices rather than words)
outcomes = torch.multinomial(p, n, replacement=True, generator=g)

i_e = vocab.index(target_e)
prob_e = p[i_e].item()
freq_e = (outcomes == i_e).float().mean().item()
print(f"the expected probability of \"{target_e}\" is {prob_e}")
print(f"the actual frequency of \"{target_e}\" in {n} trials is {freq_e}")

the expected probability of "see" is 0.38736099004745483
the actual frequency of "see" in 10000000 trials is 0.38721099495887756

1.1.3 Probability as Array Programming

$\hookleftarrow$ Table of Contents (§1.1 From Certain to Uncertain Knowledge)

1.1_gpt2output_torch.ipynb

print("gpt-2's output distribution is stored with a torch.Tensor rather than a vanilla python list")
print(f"probs.shape {type(probs.shape)}: {probs.shape}")
print(f"probs.dtype {type(probs.dtype)}: {probs.dtype}")

print("\ngpt-2's individual probabilities are accessed by indexing into")
print(f'0 ({vocab[0]}) -> {probs[0]}')
print(f'1 ({vocab[1]}) -> {probs[1]}')
print(f'2 ({vocab[2]}) -> {probs[2]}')

gpt-2's output distribution is stored with a torch.Tensor rather than a vanilla python list
probs.shape <class 'torch.Size'>: torch.Size([10])
probs.dtype <class 'torch.dtype'>: torch.float32
gpt-2's individual probabilities are accessed by indexing into
0 (see) -> 0.38736099004745483
1 (meet) -> 0.1549929976463318
2 (hear) -> 0.15230399370193481

The two primary properties of a borscht.Tensorborscht.Tensor is it’s .shapeborscht.Tensor.shape and .dtypeborscht.Tensor.dtype. evaluating to (10,) and dtype64 respectively means it’s a vector with support in ${\mathbb{R}}^{10}$ whose real values are being approximately represented with double precision floating point numbers.

Since distributions can be viewed as vectors which in turn can be viewed as lists, probability can be viewed as none other than array programming.

More generally, a borscht.Tensor with .shape of (V), corresponds to a vector with a support $\mathbf{x}\in {\mathbb{R}}^{|V|}$, where $V$ is the vocabulary. (todo: context-length is dimensionality of $\mathbf{h}$). For convenience sake, we analyzed the top 10 probabilities with a borscht.Tensor whose .shape was (10,), which corresponded to a vector $\mathbf{x}\in {\mathbb{R}}^{10}$.

The visual lecture from 3Blue1Brown’s mentions that vector can be viewed as a fancy word for a list — in the context of deep learning we use the term tensor rather than vector.

Vectors | Chapter 1, Essence of linear algebra (3Blue1Brown 2016)

Viewing tensors as lists is taking the data-oriented perspective — for instance probs being some vector $v\in {\mathbb{R}}^{10}$. However, if you squint your eyes and view the probs array from a more functional perspective, it is also a function $p:\mathbb{N}\to \mathbb{R}$statistical learningAlso, if you squint your eyes at the hashmap data structure (in Python a dict) with the functional perspective, it’s a function $map:\text{str}\to T$. In fact, there was a named tensor proposal, Tensor Considered Harmful (Rush 2019), for treating and viewing tensors more like dicts rather than tuples.. The notation ${\mathbb{R}}^{10}$ can be viewed as a function $f:\{1,\dots ,10\}\to \mathbb{R}$, or more generally for ${\mathbb{R}}^n$, the function $f:\mathbb{N}\to \mathbb{R}$. The two notations are equivalent!

missing notebook: /opt/build/repo/sitp/notebooks/1.1_gpt2output.ipynb

Important

Probability can be seen as array programming! Every valid distribution which is some function $p:\mathbb{N}\to [0,1]$ whose values are between 0 and 1 and normalizes to 1 has a corresponding borscht.Tensor with .shape of (todo).

outcomes -> events

• todo eventevent
• todo probability of ORprobability of OR
• todo sum rulesum rule

import numpy as np

# The same <span class="smallcaps">GPT-2</span> output distribution p(w | "Hello world, nice to") from above, top 10 of 50257 words
tokens = ["see", "meet", "hear", "have", "be", "know", "you", "talk", "say", "get"]
probs = np.array([0.3169, 0.1268, 0.1246, 0.1046, 0.0471, 0.0352, 0.0320, 0.0128, 0.0114, 0.0067])

tokens_start_with_h_or_s_and_contains_a_probs = []
for token, prob in zip(tokens, probs):
    stripped = token.strip()
    starts_with_h = stripped.startswith("h")
    starts_with_s = stripped.startswith("s")
    contains_a = "a" in stripped
    if (starts_with_h or starts_with_s) and contains_a:
        tokens_start_with_h_or_s_and_contains_a_probs.append(prob)

prob_starts_with_h_or_s_and_contains_a = sum(tokens_start_with_h_or_s_and_contains_a_probs)

print(f"probability that word starts with letter h: {prob_starts_with_h}")
print(f"probability that word starts with letter s: {prob_starts_with_s}")
print(f"probability that word starts with letter h or s: {prob_starts_with_h + prob_starts_with_s}")
print(f"probability that word starts with letter h or s, and contains letter a: {prob_starts_with_h_or_s_and_contains_a}")

• todo: numpy vectorization

So when you ask an LLM a question, it generates a full answer (with many sentences) word by word by repeating the following loop:

1. evaluating the probability of the next word conditioned on the input
2. selecting (or sampling) a word. halt if the word is the special END word.
3. appending it to the existing text, and evaluating the probability again with the modified input

Now that you understand the basics of language modeling, the trillion dollar question is how to produce such a conditional distribution $p(w|\mathbf{h})$? In some sense that’s all there is to large language models.

1.1.4 Are You Pythonic?

$\hookleftarrow$ Table of Contents (§1.1 From Certain to Uncertain Knowledge)

1.2 Next Token Prediction is Classification and Compression

$\leftarrow$ §1.1 $\uparrow$ Table of Contents §1.3 $\to$

In which The spelled-out intro to language modeling: building makemore (Karpathy) we introduce the terminology of software 2.0 in the context of language modeling; the joint, marginal, conditional with multidimensional arrays, reductions, and broadcasting; implement our first language model by maximizing the likelihood of biased ngrams; evaluate said model’s limits with information theoretical and probabilistic quantities such as perplexity, entropy, bias variance tradeoff, the curse of dimensionality; and finally, understand intelligence more generally as compression.

1.2.1 The Estimation of Software 2.0

$\hookleftarrow$ Table of Contents (§1.2 Next Token Prediction is Classification and Compression)

As mentioned in the previous chapter, large language models are trained using methods from the discipline of deep learningdeep learning, which in turn are based in machine learningmachine learning, which in turn are based in statistical learningstatistical learningActually, the implementation of learning is not just limited to the stochastic and infinitely continuous methods of software 2.0, more classically known as connectionism. It’s also possible to implement with logical and finitely discrete methods of software 1.0 (symbolism), albeit with much less success. See inductive logic programming also known as reverse deduction https://en.wikipedia.org/wiki/Inductive_logic_programming., which in turn are based in statistical inferencestatistical inference, which in turn employs the use of statistical estimationstatistical estimation. What are the meanings of all these terms? Roughly speaking, people in the business of estimation like AI researchers at frontier labs are given data (in statistical jargon, observe samples) such as the internet and would like to estimate distributions (of the population) from said data, namely a $p_{Y|X}(y|x)$ like GPT-2. Because distributions are functions, statistical estimation is also known as function approximationfunction approximation.

Estimation of distributions is one possible activity of making inferences from data — as opposed to simply describing data See https://en.wikipedia.org/wiki/Descriptive_statistics. — another useful form of inference is hypothesis testingSee https://en.wikipedia.org/wiki/Statistical_hypothesis_test., the process of deciding whether observed data provide sufficient evidence to reject a null hypothesis. Estimation however is the primary statistical tool used as opposed to hypothesis testing, and once a distribution like $p_{Y|X}(y|x)$ is estimated from data by research labs, users can use such a distribution in order to generate sentences — hence this most recent wave of statistical methods from software 2.0 being coloquially known as genaigenai, other previous monikers being big databig data, data miningdata mining, data sciencedata science, and pattern recognitionpattern recognition — with an inference loop, namely predicting a distribution over the next token, sampling and appending a token to the history and repeating.

Generally speaking, the difference in direction of estimating a distribution from data versus generating data from a distribution is the primary distinction between probability and statistics, although the line as we will soon see gets quite blurryIn the same way traditional data and function blurs with software 1.0, since functions are also data.. This distinction between the two is more broadly known as the difference between analysisanalysis and synthesissynthesisSee https://plato.stanford.edu/entries/analytic-synthetic/. Researchers analyze the observed data of the internet and recover an estimated distribution whereas users synthesize data and generate sentences. Sometimes they are referred to with direction in time: synthesis is the forward direction whereas analysis is the backward or inverse. Throughout the book, we will incrementally explore increasingly expressive languge models that can do a better job at recovering $p_{Y|X}(y|x)$ from data.

The availibility of said data determines the regime of statistical learning. That is, if sample data is available in the form of input-output pairs $d:=\{({\mathbf{x}}^{(i)},{\mathbf{y}}^{(i)}){\}}_{i=1}^n$, then the regime is known to be supervised learningsupervised learning or self-supervised learningself-supervised learning, for the machine that is attempting to learn a distribution via induction is being supervised with the expected outputs, whereas the regime where such data is unavailable is known as unsupervised learningunsupervised learning. Where in unsupervised learning the goal of inference is that of understanding like clustering or dimensionality reduction, in the regime of supervised learning the goal is the aforementioned statistical estimation of a function that predicts the expected output given the input. For instance, with the following self-supervised dataset

missing notebook: /opt/build/repo/sitp/notebooks/1.2_bigram.ipynb

Moreover, the regime of supervised learning is further refined depending on the kind of output the estimated function $f:\mathcal{X}\to \mathcal{Y}$is predicting. In the case of a categorical output, the supervised learning is known as classificationclassification for the machine’s goal (the estimated function’s output) is trying to classify the input’s category. Modeling languagelanguage models as next-token prediction is in fact an instance of classification since it is classifying the next most likely token $y\in \mathcal{Y}=V$ given an input sentence $(x^{(1)},\dots ,x^{(t)})\in \mathcal{X}=V^{\ast }$ This stands in contrast to the case when the output that the estimated function is predicting is quantitative, which is known as regressionregression.

Somewhat more confusingly, depending on the domain or context of statistical modeling, the terminology for inputs $x$ and outputs $y$ might be further refined. That is, inputs and outputs are also respectively referred to as independent variableindependent variable and dependent variablesdependent variable, as well as predictorspredictors and responsesresponses in the context of experimental design. They can also be referred to as featuresfeatures and targetstargets. The difficulty of naming is not limited to the discipline of programmingSee https://martinfowler.com/bliki/TwoHardThings.html. Whatever you want to call them, at the heart of statistical learning, machine learning and deep learning is the activity of estimation, or the approximation of functions with inputs and outputs $f:\mathcal{X}\to \mathcal{Y}$. two cultures of statisticstwo cultures of statistics

Warning

Inference is an overloaded term, and hence ambiguous without context! In generative AI the term inference refers to the forward process of evaluating a population’s distribution $p_{Y|X}(y|x)$ in order to generate samples via prediction. In classical statistics the term inference refers to the backward process of inferring a population’s distribution from observed samples

But for the classification of language modeling, how do we reconcile the goal of estimating some function $f:V^{\ast }\to V$, when GPT-2 is the probability distribution $p_{Y|X}(y|x)$ which is a function of type $p_X:\mathbb{N}\to [0,1]$? The way we reconcile such a type mismatch is to map over words from their string representations in $V$ to their index representations in $\mathbb{N}$ with a dict[str, int] instantiated by {c:i for i,c in enumerate(v)}

missing notebook: /opt/build/repo/sitp/notebooks/1.2_bigram.ipynb

However, note that the input type of a next-token classifier are sentences $(x^{(1)},x^{(2)},x^{(t)})\in V^{\ast }$. (todo: BTD)

$$V\stackrel{I}{\to }\mathbb{N}\stackrel{p}{\to }{\Delta }^{|V|}\stackrel{Q}{\to }\mathbb{N}\stackrel{I^{-1}}{\to }V$$

Let us now estimate our first language model $p:\mathbb{N}\to {\Delta }^{|V|}$

1.2.2 The Joint, Marginal, Conditional, and Chain Rule

$\hookleftarrow$ Table of Contents (§1.2 Next Token Prediction is Classification and Compression)

Let us first start by generalizing the notion of a distribution of a single random variable $p_X:\mathbb{N}\to \mathbb{R}$ which corresponds to a flat array of real-valued numbers satisfying non-negativity and normalization to that of a joint distributionjoint distribution Stanford CS109 Probability for Computer Scientists | Joint Distributions | Lecture 11, defined as a high dimensional function $p_{X_1,\dots ,X_t}$ with type $p_{X_1,\dots ,X_t}:\mathbb{R}\times \cdots \times \mathbb{R}\to [0,1]$ which corresponds to a multidimensional array. Both map the values of $t$ random variables to their probabilitiess. That is, evaluating

$$p_{X_1,\dots ,X_t}(x_1\underset{\text{and}}{\underbrace{,}}\dots ,x_t)$$

produces a single probability for random variable $X_1$ taking on the value $x_1$ and $X_2$ taking on the value $x_2$ and so on, up until the random variable $X_t$ taking on the value $x_t$. Because the joint can be a distribution of two random variables, three, up to an arbitrary t, it’s often vectorized as $p_{\mathbf{X}}(\mathbf{x})$.

In the case of (probabilistic) language modeling, the primary joint distribution of interest is the probability of a fixed t lengthed sentence occuring, which is $p_{X1,X2,\dots ,X_t}(x_1,x_2,\dots ,x_t)$ or $p(\mathbf{x})$ with $x\in {\mathbb{R}}^t$, where the $X_i$’th random variable represents the $i$’th word in the $t$-lengthed sentence. For instance, $p(x_1=\text{Hello},x_2=\text{nice},x_3=\text{to},x_4=\text{meet},x_5=\text{you})$ is the probability of the complete sentence "Hello nice to meet you" (todo, instead of distribution over words, it’s a distribtion over sentences).

Warning

Note with the random variables being modeled the way they are (that is, with order) that the joint $p(x_1,\dots ,x_t)$ is commutative. For instance, with $n:=2$, $p(x_1,x_2)=p(x_2,x_1)$. Semantically speaking, $p(x_1=\text{Hello},x_2=\text{world})=p(x_2=\text{world},x_1=\text{Hello})$ is saying that the probability of the first word being “hello” and the second word being “world” is equal to the probability of the second word being “world” and the first word being “hello”. This holds more generally for any permutation of $x_1,x_2,\dots ,x_t$.

With us defining the notion of a joint probability as a multidimensional array, it’s probably no surprise that the joint distribution $p_X(\mathbf{x})$ corresponds to a borscht.Tensor with .shape (V,...,V) where .ndim is t. That is, .shape is a t-lengthed tuple. Take a look at Figure 1.2.2 below:

Figure 1.2.2

In particular, let us consider the joint distribution of two random variables $p_{X,Y}$, whose 2 dimensional surface can be plotted. This distribution has type $p_{X,Y}:\mathbb{N}\times \mathbb{N}\to \mathbb{R}$, which corresponds to a borscht.Tensor with a .shape of (d1,d2) — a matrixmatrix supported in ${\mathbb{R}}^{d_1\times d_2}$, or in other words, a $d_1$ by $d_2$ rectangular array that maps indices in $\{0,1,\dots ,d_1\}\times \{0,1,\dots ,d_2\}$ to their values in $\mathbb{R}$. Because the array for a matrix can be visualized as a rectangle when taking the data-oriented perspective, the first dimension $d_1$ is usually referred to as the number of rowsrow, and the second dimension $d_2$ is referred to the number of columnscolumn. This naming of course breaks down for tensors with a higher number of dimensions. The code example initializes a matrix with borscht.zeros((n,m)), and plots the surface with ax.plot_surface().

import borscht
x = borscht.zeros((d1,d2))
print(x)

The joint of three random variables $p_{X,Y,Z}:\mathbb{N}\times \mathbb{N}\times \mathbb{N}\to \mathbb{R}$ is plotted as what is known as a voxel plot (it’s difficult to see 4 dimensions alltogether), which corresponds to a borscht.Tensor with a .shape of (d1,d2,d3) — a tensortensor supported in ${\mathbb{R}}^{d_1\times d_2}$, a cubic array that maps indices in $\{0,1,\dots ,d_1\}\times \{0,1,\dots ,d_2\}\times \{0,1,\dots ,d_3\}$ to their values in $\mathbb{R}$. (todo, remark on naming with tensor.) Similarly, the tensor is initialized with borscht.zeros((n,m,p)), is filled with values, and is plotted with ax.voxels(x).

import borscht
x = borscht.zeros((d1,d2,d3))
print(x)

Important

Once again, probability can be viewed as array programming!

(todo, warning on

With the notion of a joint distribution now defined, the actual joint of interest to language modeling is the probability of a sentence taking on specific words $p(W_1=v_i,W_2=v_i,\dots ,W_t=v_i)$ where $v_i\in V$, let $t:=2$ so that we can consider the joint of two random variables $p(W_1=v_i,W_2=v_j)$, which corresponds to a matrix represented by a borscht.Tensor with .shape of (V,V), in this case (5,5). In Figure 1.2.3, the surface on the left and the matrix on the right represent the same distribution $p_{W_1,W_2}(v_i,v_j)$. These values are one possible assignment to the function body $p_{W_1\times W_2}:\mathbb{N}\times \mathbb{N}\to \mathbb{R}$, or equally, one possible assignment to the borscht.Tensor. For instance, we see that $p_{W_1,W_2}(\text{to},\text{meet})=0.16$. (todo, it’s really $p(W_1=3,W_2=4)$. (todo, table is technically infinite. t is not fixed.). Note that the corresponding $v_i\in V$ is provided, even though the distribution is defined on $\mathbb{N}$.

Figure I.2.3

Note

In SITP, we will incrementally increase the expressivity of our next-token predictors or classifiers in order to model language more accurately. Throughout the book, we will maintain a visual language of color in addition to an algebraic one to indicate a consistent purpose or role regardless of what space an object is represented with. For instance, input sentences will be coded with blue whether they are strings $x\in V$ indices $i\in \mathbb{N}$ one hot vectors $\mathbf{x}\in {\mathbb{R}}^{|V|}$, or soon to be embeddings $\mathbf{x}\in {\mathbb{R}}^{|E|}$. Without loss of generality, the same holds for output tokens except coded with orange. More colors will also be introduced later on. This can be toggled at anytime by clicking the located in the top left corner next to the theme picker.

• modify left figure to highlight the 0.15 and 0.16s
• why is the marginal a uniform (0.2)?

For a single random variable $W_i$, we can see that the truth of it’s state is indeed distributed across the vocabulary $V$. For instance, scanning the blue arrow from top to bottom, we see that state of random variable $W_1$ is distributed across five words. In particular, for one of those values $i=3$ which corresponds to $\text{to}$, it’s state in of itself is distributed across $W_2$ with the orange arrow from left to right. For instance, the chance of "to meet"should be more likely than "to to". The code for this joint is the transliteration of the matrix from Figure I.2.3 into borscht.Tensor with .shapeof (5,5):

missing notebook: /opt/build/repo/sitp/notebooks/1.2_bigram.ipynb

Warning

Let us digress on two technical aspects of borscht: shape suffixes and indexing.

Digression 1: Shape Suffixes

The first place to draw your attention to is in the particular choice of naming the joint distribution as p_VV. The _VV suffix is documentation-based coding brought back from the deadSee Hungarian notation https://en.wikipedia.org/wiki/Hungarian_notation. Also see Google’s C++ Style Guide — a somewhat representative and perhaps canonical standard of most source conventions — frown on Hungarian notation. Finally see Tensor Considered Harmful (Rush 2019) and Named Tensor Notation (Chiang et al. 2021). as a reading aid to humans and fellow LLMs to denote that the borscht.Tensor has a .shape of (V, V), where in this case V evaluates to 5. We even assert such fact with assert p_VV.shape == (V, V). Such judgement can actually be enforced by typechecking dynamically at runtime with the jaxtyping Python package, particularly Float[np.ndarray, "V V"] (expand on jaxtyping and python typechecking):

missing notebook: /opt/build/repo/sitp/notebooks/1.2_bigram.ipynb

The function printp_VV takes a parameter p_VV with type Float[np.ndarray, "V V"], so the first invocation with the argument p_VV itself passes. Mutating p_VV to alias a borscht.Tensor with .shape (3,) and applying the function a second time results in a type error.

Digression 2: Indexing

Second, notice how the value of $p_{W_1,W_2}(\text{to},\text{meet})$ is retrieved with p_VV[toi, meeti]. (todo, curry by indexing into first or second):

missing notebook: /opt/build/repo/sitp/notebooks/1.2_bigram.ipynb

Note

Digression is now over.

Now, with borscht technicalities of shape suffixes and indexing aside, let us move onto the next probabilistic concept of the marginal. From this joint distribution of two random variables $p_{W_1,W_2}$, we can recover both distributions of single random variables $p_{W_1}$ and $p_{W_2}$.

Consider again that second “level” of distribution. That is, the specific state of $W_1=\text{to}$ being distributed in of itself across the first row of the matrix according to $W_2$. We can “take out” the randomness of $W_2$ by summing all the entries in the first row, and no longer have this second level or axis of state being distributed, and are simply left with $p_{W_1}(\text{to})=0.2$. In general, for any $v_i$, we have that $p_{W_1}(v_i)=\sum _{v_j\in V}{p}_{W_1,W_2}(v_i,v_j)$. Note that $v_i$ is fixed! In the context of recovering a distribution of a single random variable from a joint distribution, such single random variable distribution is referred to as the marginal distributionmarginal distribution. Transliterating this to borscht, we use borscht.Tensor.sum()borscht.Tensor.sum reduce operationreduction (expand):

missing notebook: /opt/build/repo/sitp/notebooks/1.2_bigram.ipynb

(todo: generalize for $t$ randvars) torch.sum

After generalizing the notion of a distribution into a joint distribution, and recovering marginals from said joint, we are now finally ready to backtrack to refining the notion of the conditional. $p_{W_1,W_2}(\text{to},\text{meet})=\frac{0.16}{0.2}=0.8$.

conditional distributionconditional distribution Stanford CS109 | Conditional Probability and Bayes | Lecture 4

$$\underset{\text{conditional}}{\underbrace{p_{Y|X}(y|x)}}\stackrel{△}{=}\frac{\stackrel{\text{joint}}{\overbrace{p_{Y,X}(y,x)}}}{\underset{\text{marginal}}{\underbrace{p_X(x)}}}$$

so, the conditional is defined as the joint over the marginal distribution that is being conditioned on. chain rulechain rule product ruleproduct rule (decomposing a joint of sentences with a chain of conditionals over words requires proof (ntermezzo).) with an inference loop, namely predicting a distribution over the next token, sampling and appending a token to the history and repeating.

We can further clarify our notion of language modelslanguage models to next-token predictors performing classification which then auto-regressively feeding such classifying prediction into the next evaluation via chain or product rule.

More generally, broadcastingbroadcasting broadcast divide

There is a rule you should learn at long last,
combination of tensors the task.
Dims right-aligned,
extra left 1s assigned,
match paired dimensions: Broadcast!

9 x 1 x 3 1 x 8 x 1 9 x 8 x 3

Tensor Puzzles (Sasha Rush, Marcos Treviso)

missing notebook: /opt/build/repo/sitp/notebooks/1.2_bigram.ipynb

1.2.3 Estimating Markov Biased Histograms by Counting

$\hookleftarrow$ Table of Contents (§1.2 Next Token Prediction is Classification and Compression)

We know that $p_{X|Y}(y|x)$ denotes the conditional distribution of the next word given the history of sentence thus far, we know this distribution corresponds to a borscht.Tensor with .shape (V,), and we know the function type is $p_{Y|X}:V^{\ast }\to {\Delta }^{|V|}$ (todo: ${\Delta }^{|V|}$ can be viewed as $|V|\to \mathbb{R}$) which maps sentences to distributions over next words, but how do we go about estimating the distribution $p(y|x)$?

When initially describing the distribution of next possible tokens $p(y|x^{(1)},\dots ,x^{(t)})$ with the stochasticity of probability theory, we also mentioned other phenomena like coin flips, die rolls, and card draws that are modeled stochastically as well and thus with distributions. In the case where such physical objects are deemed to be fair — that is, no one has weighted the coin, die, nor rigged the deck — we can assume that all outcomes are equally likely and subsequently calculate the probability of any event (todo, go write that exposition) is total number of outcomes in the event over the number of outcomes in the sample space.

However, it’s clear that language is not a stochastic process with equally likely outcomes. For instance words “the”, “is, “a” are more likely than “quixotic” for instance. It is for this reason that the first idea to recover or estimate $p(y|x^{(1)},\dots ,x^{(t)})$ is by simply countingcounting the number of times the event of interest occurs from some dataset? sample? and dividing it by the total number of events in the event space, thereby calculating the relative frequency estimaterelative frequency estimate:

$$\frac{\text{\# "Hello world nice to meet you"}}{\text{\# "Hello world nice to meet"}}$$

Taking the internet as our dataset, this ratio is asking out of all the sentences that start with “Hello world nice to meet”, how many times does “Hello world nice to meet you” occur? (todo, define unbiased estimator. in the limit..) While the relative frequency estimate is an unbiased estimatorunbiased estimator, the problem is that language as a phenomena continually evolves, grows, and adapts itself. The fact that natural language can be easily modified and innovated on by new generations of our culture is a feature rather than a bug. While new slang is invented every day, the classic examplar for the problem of using a relative frequency estimate is the following sentence posed by Noam ChomskySee https://en.wikipedia.org/wiki/Colorless_green_ideas_sleep_furiously:

$$\text{Colorless green ideas sleep furiously.}$$

Although such sentence now exists on the internet because of Noam Chomsky, the sentence shows one that is syntactically valid but perhaps not semantically meaningful, which, very much resembles the essence of slang. Now, since estimating $p(y|x^{(1)},\dots ,x^{(t)})$ with a relative frequency estimate is intractable, let us explore alternative ways. This motivates the need to dig deeper into the precise defintion of the conditional, what is what we explore here in 1.2.2 The Joint, Marginal, Conditional, and Chain Rule

$$\underset{\text{conditional}}{\underbrace{p_{Y|X}(y|x)}}\stackrel{△}{=}\frac{\stackrel{\text{joint}}{\overbrace{p_{Y,X}(y,x)}}}{\underset{\text{marginal}}{\underbrace{p_X(x)}}}$$

Rather than strive for perfection, we can achieve the good by introducing some biasbias, known as the Markov assumptionmarkov assumption.

• (todo, explain data d)
• counting is the maximum likelihood estimatemaximum likelihood estimate (MLE).
• simplest idea: unbiased estimator count the relative frequency (which is the MLE) $\frac{h,w}{h}$
• need chain/product rule here. but do not introduce bayes rule, for that will introduce discriminant vs generative modeling.
• intractable, so need to introduce bias. with markov assumption (truncated conditional in chain rule)

missing notebook: /opt/build/repo/sitp/notebooks/1.2_bigram.ipynb

missing notebook: /opt/build/repo/sitp/notebooks/1.2_bigram.ipynb

missing notebook: /opt/build/repo/sitp/notebooks/1.2_bigram.ipynb

1.2.4 The Scoring of Software 2.0

$\hookleftarrow$ Table of Contents (§1.2 Next Token Prediction is Classification and Compression)

Scoringscore such estimation involves probabilities of probabilities.

maximum likelihood estimation Stanford CS109 Probability for Computer Scientists | M.L.E. | Lecture 21 (Piech)

missing notebook: /opt/build/repo/sitp/notebooks/1.2_bigram.ipynb

Reinventing Entropy | Compression is Intelligence Part 1 (3Blue1Brown 2026)

1.2.5 The Evaluation of Software 2.0

$\hookleftarrow$ Table of Contents (§1.2 Next Token Prediction is Classification and Compression)

1.2.6 Summary and Further Reading

(Eisenstein 2019)’s Chapter 6. Language Models; (Jurafsky and Martin 2026)’s Chapter 3. N-gram Language Models;

(further reading for torch)

1.3 Parameterizing Classification with Discriminant Functions

$\leftarrow$ §1.2 $\uparrow$ Table of Contents §1.4 $\to$

In which we parameterize classification with a three step process of parameterization, scoring, and optimization; in particular, by directly optimizing the mean squared error of the perceptron; and understand the perceptron’s classifying by expanding our notion of vectors from vectors as lists to vectors as points, enabling not only array programming with distributions, but also with geometry

1.3.1 The Parameter Estimation of Software 2.0

$\hookleftarrow$ Table of Contents (§1.3 Parameterizing Classification with Discriminant Functions)

Whereas in §1.3 Parameterizing Classification with Discriminant Functions, we now change our approach to language modeling next-token prediction (and statistical machine learning more generally) from that of counting to parameterizingparameterization, meaning to define a function with parametersparameters $\theta$ in addition to it’s input $x$. The type of the function $f$ on the right is changed to take in a parameter $\theta$ and an input sentence $x$ in order to produce an output token $y$. However, passing in $\theta$ as an argument on every function evaluation is cumbersome since, after we find a good set of parameters $\theta$, they in fact do not changeUnlike the human brain, which is constantly being reconfigured and rewired!. Rather, what we would like is to instantiate the function with parameters once, and then evaluate such function whenever we needed to predict the next token. The solution to this is to change the type of $f$ to

$$f:\Theta \to (V^{\ast }\to V)$$

which we can partially applypartially apply by passing in $\theta$ as input, and receive the next-token predicting function $V^{\ast }\to V$ itself. In other words, we have curriedcurried the type into a higher order functionhigher order function which can return function themselves. Transliterating a toy parameterized model which multiplies an input $x$ with a parameter $\theta$ to borscht then looks like

import borscht
import borscht.nn as nn

class MyParameterizedModel(nn.Module):
  def __init__(self, theta):
    self.theta = theta

  def __call__(self, x):
    return x*self.theta

theta = 42
model = MyParameterizedModel(theta)
y = model(1)
print(y)

where $f^{\ast }\leftarrow f(\theta )$ is the first evaluation with $f^{\ast }$ denoting the next-token predictor $f:V^{\ast }\to V$, and $y\leftarrow f^{\ast }(x)$ is the second, These correspond to the Python’s class __init__ and __call__ methods respectively. When functions are parameterized, they are denoted as $f_{\theta }$ or $f(x;\theta )$ to reflect such.

Now, because the model is characterized or determined by said parameters $\theta$, such parameters must be estimatedparameter estimation, usually from the dataset $d:=\{({\mathbf{x}}^{(i)},y^{(i)}){\}}_{i=1}^n$. The way this is done is to score such parameters with an objective known as the loss functionloss function, and subsequently find the parameters that optimizeoptimization such loss function, a process denoted by $\mathrm{argmin}$. That is, the problem of finding a good model reduces to finding a good set of parameters according to the loss, which is comprised of a three step process:

todo. fri oct 2

1. risk
2. step zero?

The parameterization of functions and the subsequent parameter estimation are the primary workhorses in the function approximation of statistical learning. However with deep learning, we concern ourselves primarily with the class of models $f_{\theta }$ that are biologically inspired by the brain’s neuroanatomy. Although far from realistic, a neural network $f_{\theta }$ models the signals between the axon of one neuron to the dendrite of another with a set of weightsweights $w$ and biasesbias $b$ where $\theta :=\{w,b\}$. Such models $f_{w,b}(x)$ are referred to as deep neural networks, in which we start by training single layer networks comprised of perceptronsperceptrons Stanford CS229: Machine Learning | Lecture 5 - Perceptron and Logistic Regression. The computational circuit of a perceptron displayed on the left side of Figure 1.3.1, where a perceptron produces a real-valued scalar $y$ by taking the sum $\sum$ of $\mathbf{x}$ with parameters $w$ and $b$.

Figure 1.3.1

The perceptron $f_{w,b}$ is a weighted sum of $x_i$’s and $w_i$’s, which can be rewritten as the dot productdot product of the two vectors $x$ and $w$, and trivially transliterates to borscht.dot()

import borscht
import borscht.nn as nn

class Perceptron(nn.Module):
  def __init__(self):
    self.w = borscht.randn(d)

  def __call__(self, x: borscht.Tensor):
    return borscht.dot(self.w, x)

model = Perceptron()
y = model(x)
print(y)

where the bias term $b$ is folded by prepending $b$ and $1$ into $w_0$ and $x_0$ respectively, subsequently updating the dimension from $d$ to $d+1$. With respect to the three step process, we have completed the first one of parameterizing the architecture $f_{\theta }$. But before moving on to the second and third steps of scoring such parameters with $\ell (\theta )$ and optimizing them with $\mathrm{argmin}\mathcal{L}(\theta )$, you may have noticed that the domain of $f$ is ${\mathbb{R}}^d$ rather than sentences $x\in V^{\ast }$. While we are now familiar with high dimensional arrays, a matrix and it’s row vectors were being used to model the stochastic matrix $P\in {\mathbb{R}}^{|V|\times |V|}$ where discrete word indices were encoded into and out of such matrix with $I$. Moreover, how is it’s codomain of $\mathbb{R}$ classifying some category $y\in V$?

Warning

So how is the perceptron taking in vectors ${\mathbb{R}}^d$ directly,
and how does producing scalars $\mathbb{R}$ perform a classification?
To answer this we need to take the proverbial red pill,
let us go deeper into the rabbit hole of the borscht.Tensor.

1.3.2 Analytic Geometry as Array Programming

$\hookleftarrow$ Table of Contents (§1.3 Parameterizing Classification with Discriminant Functions)

For now, we reconcile how the real-valued output $y=f_{w,b}(x)$ is parameterizing classification with a direct decisiondecision without inferringinferring a probability distribution. We square the latter with the perceptron in §1.4 Parameterizing Classification with Probabilistic Models thereafter. The red pill needed for such reconciliation will deepen our understanding of vectors and matrices, both represented by the multidimensional array of the borscht.Tensor.

Important

The red pill is one which expands our previously conceived notion of vectors from the so-called computer scientist view of tensors as lists, to that of the geometrical view with vectorsvector as points in space. It is quite likely that you may have encountered vectors initially as geometrical points or lines, rather than as computational lists. Because of the inherent difficulty to visualize vectors $\mathbf{x}$ that have support in high dimensional spaces ${\mathbb{R}}^d$ where $d$ is greater than 3, we will start with an example of vectors in ${\mathbb{R}}^2$.

Let us consider a simplified language model setting where the goal is to perform binary classification with a context or history of length 1 on a vocabulary of length 2 where $V:=[new,york]$. And moreover, recall that we are parameterizing such next-token classification $f:V^{\ast }\to V$ with weights and biases via perceptron $f_{\theta }(\mathbf{x}):={\mathbf{w}}^{\top }\mathbf{x}$, corresponding to def __call__(self, x): borscht.dot(self.w, x).

Starting with reconciling the type mismatch of our classifier $h:V^{\ast }\to V$ with the perceptron’s domain of ${\mathbb{R}}^d$, we continue to map words from their string representation in $V$ to their positional indices in $\mathbb{N}$, however, there is an additional transformation of the input thereafter from $\mathbb{N}$ to ${\mathbb{R}}^d$. The simplest idea for such a transformation is to map or embedembedding the position index $i\in \mathbb{N}$ to a $|V|$-dimensional vector $x\in {\mathbb{R}}^{|V|}$ with all entries set to 0 except for the $i$’th dimension which is set to 1. Such a vector is referred to as a one-hot vectorone-hot vector and is known as a localist representationlocalist representation of a word, where each word’s meaning is represented by itself in isolation without a regard to the whole lanugage of English as a whole — for instance, adding two one-hot vectors is meaningless, since any vector with more than one non-zero entry is ill-defined. The final series of representational transformations of the input then looks like

where the corresponding function in borscht to map indices to their one hot vectors is borscht.nn.functional.one_hot(tensor, num_classes=-1) which takes a LongTensor with index values of shape (*) and returns a tensor of shape (*, num_classes).

import borscht
import borscht.nn.functional as F

X = borscht.Tensor([0, 1, 2, 3])
X_BV = F.one_hot(X)

When we map our words $v\in V$ to their one-hot vector representations $\mathbf{x}\in {\mathbb{R}}^2$ with the earlier constructed hashmap c2i, the point $(1,0{)}^{\top }$ for the word $\text{new}$ is located on the $x_1$-axis and the point $(1,0{)}^{\top }$ for the word $\text{york}$ on the $x_2$-axis. This is can seen in Figure 1.3.2

Figure I.3.2 On the left hand side the one hot vectors in ${\mathbb{R}}^2$ are displayed for the strings $\text{new},\text{york}\in V$. On the right hand side a decision boundary linearly separates ${\mathbb{R}}^2$ into decision regions ${\mathcal{R}}_1$ and ${\mathcal{R}}_2$.

Now that we’ve reconciled the type mismatch of our classifier $h:V^{\ast }\to V$ with the perceptron’s domain ${\mathbb{R}}^d$, let us clear our confusion for the perceptron’s codomain of $\mathbb{R}$. In terms of a desirata for the decision step — that is assigning the input $\mathbf{x}\in {\mathbb{R}}^2$ some outcome $y\in V$ — it would be nice if we could divide the plane into two decision regionsdecision regions ${\mathcal{R}}_1$ and ${\mathcal{R}}_2$ associated with the output categories of $y=\text{york}$ and $y=\text{new }$ respectively with some line acting as a decision boundarydecision boundary $l$ such that our classifier $\hat{f}$ assigned the next word $y=\text{york}$ given $(1,0{)}^{\top }$ for $\text{new}$, and perhaps $y=\text{new}$ given $(1,0{)}^{\top }$ for $\text{york}$Bear with the toy example, for although “york new” is non-sensical, it allows us to focus our attention on the mechanics of the perceptron itself..

Rather than define the decision boundary with an explicit equation for a line, another possible idea is to define it as the set of all points $\mathbf{x}$ that are perpendicularperpendicular to the weight vector $\mathbf{w}$.

length(norm) = sqrt[angle(inner/dot product)]

How do we define the notion of perpendicularity without using our eye balls?

• rotationrotation
• orthogonal projectionorthogonal projection

foobar

• distance with metricmetric
• length with normnorm
• length with Euclidean normEuclidean norm ${\ell }_2$
• length with Manhattan normManhattan norm ${\ell }_1$
• dot productdot product
• inner productinner product
• similarity with angleangle
• dot productdot product $\sum _{i=1}^{d+1}{w}_i{x}_i$ = ${\mathbf{w}}^{\top }\mathbf{x}$
• pendicularity with orthogonalityorthogonal

and the final classifier has the form $\hat{f}(\mathbf{x}):=g(\sum _{i=1}^d{w}_i{x}_i)$ where $g:\mathbb{R}\to V$ is some function

$$g(z):=\{\begin{array}{ll}\text{york}& \text{if }z>0\\ \text{new}& \text{if }z\le 0\end{array}$$

$\left[\begin{array}{c}1\\ 0\end{array}\right]\cdot \left[\begin{array}{c}1\\ 1\end{array}\right]=1\left[\begin{array}{c}0\\ 1\end{array}\right]\cdot \left[\begin{array}{c}1\\ 1\end{array}\right]=1$

$\left[\begin{array}{c}1\\ 0\end{array}\right]\cdot \left[\begin{array}{c}1\\ -1\end{array}\right]=1\left[\begin{array}{c}0\\ 1\end{array}\right]\cdot \left[\begin{array}{c}1\\ -1\end{array}\right]=-1$

Linear combinations, span, and basis vectors | Chapter 2, Essence of linear algebra (3Blue1Brown 2016)

Linear transformations and matrices | Chapter 3, Essence of linear algebra (3Blue1Brown 2016)

1.3.3 The Optimization of Scores in Software 2.0

$\hookleftarrow$ Table of Contents (§1.3 Parameterizing Classification with Discriminant Functions)

$$\mathrm{argmin}\mathcal{L}(\theta )$$

baz

• least squares lost per $x$ $\mathcal{L}(\mathbf{w}):=\sum _{i=1}^n(y^{(i)}-{\mathbf{w}}^{\top }{\mathbf{x}}^{(i)}{)}^2$
  • when d=1, this is finding the line of best fit? (with bias)
  • but then we use the line of best fit as a discriminator with $g$
• is a linear system on $X\mathbf{w}$=$\mathbf{y}$
• min $\Vert \mathbf{y}-X\mathbf{w}\Vert$

1.3.4 Evaluation with Bias Variance Tradeoff

$\hookleftarrow$ Table of Contents (§1.2 Next Token Prediction is Classification and Compression)

1.3.5 Summary and Further Reading

(Duda et al. 2001)’s Chapter 5. Linear Discriminant Functions; (Hastie et al. 2009)’s Chapter 3. Linear Methods for Regression; (Bishop 2006)’s Chapter 3. Linear Methods for Regression; (Murphy 2022)’s Chapter 11. Linear Regression

(Goodfellow et al. 2016)’s Chapter 2. Linear Algebra and Chapter 5. Machine Learning Basics; (Murphy 2022)’s Chapter 4: Statistics and Chapter 7. Linear Algebra

(Bishop 2006)’s Chapter 2. Probability Distributions; (Goodfellow et al. 2016)’s Chapter 3. Probability and Information Theory; (Murphy 2022)’s Chapter 2. Probability: Univariate Models, Chapter 3. Probability: Multivariate Models

(further reading for torch)

1.4 Parameterizing Classification with Probabilistic Models

$\hookleftarrow$ Table of Contents

In which we…

1.4.1 Probabilistic Architectures: Discriminative vs Generative

In machine learningmachine learning, there are actually three ways to estimate a classifier. Recall that the definition of the conditional distributionconditional distribution $p(y|x)$ is the joint $p(y,x)$ over the marginal $p(x)$, and decomposing such joint with the product or chain rule starting with $p(y)$ gives us a relationship of what is known as bayes rulebayes rule

$$\begin{array}{rl}\underset{\text{conditional}}{\underbrace{p_{Y|X}(y|x)}}& \stackrel{△}{=}\frac{\stackrel{\text{joint}}{\overbrace{p_{Y,X}(y,x)}}}{\underset{\text{marginal}}{\underbrace{p_X(x)}}}\\ & =\frac{\stackrel{\text{prior}}{\overbrace{p_Y(y)}}\stackrel{\text{likelihood}}{\overbrace{p_{X|Y}(x|y)}}}{\underset{\text{marginal}}{\underbrace{p_X(x)}}}\end{array}$$

Because the bigram model first produces an inferredinference conditional distribution $p(y|x^{(1)},\dots ,x^{(t)})$ which expresses it’s uncertainty of the next token followed by decidingdecision a specific token from such a distribution, the bigram classifier is an instance of a probabilistic architectureprobabilistic architecture. namely with a discriminative probabilistic modeldiscriminative probabilistic model

And so an alternative approach to modeling $p(y|x)$ directly by (expand), the conditional distribution can be modeled as $p(y|x)=\frac{p(y)p(x|y)}{p(x)}$ and is known as a generative probabilistic modelgenerative probabilistic model because it provides the opportunity to generate samples of the input space $\mathcal{X}$ with $p(x|y)$ where in the case of language modeling the input space is the set of all possible sentences $\mathcal{X}:=V^{\ast }$).

Finally, the inference stage can be skipped alltogether by simply discriminating a specific $y\in \mathcal{Y}$ in one step (in the case of language modeling, the output space is the set of all possible next tokens $\mathcal{Y}:=V$). Models that only consist of that one step are known as discriminant functionsdiscriminant functions and is where we actually begin our discussion of parametric models namely with linear regressionlinear regression corresponding to a single neuron, even though classification and regression were differentiated earlier as producing qualitative and quantitative outputs. We will revisit parameterizing the architecture of a discriminative probabilistic model with logistic regression in 1.4 Parameterizing Classification with Probabilistic Architectures.

which involves a two stage process of specification and implementation.

First, the distribution of $p(y|x)$ must first be specified as a family of distributions $\mathcal{P}$ parameterized by a set of parameters $p$, followed by modeling $p$ itself as some parameterized function class $f_{\theta }(x)$

$$\begin{array}{rl}{p}_{Y|X}(y|x)& =\stackrel{\text{step zero}}{\overbrace{\mathcal{P}(p)}}\\ & =\underset{\text{step one}}{\underbrace{\mathcal{P}(f_{\theta }(x))}}\end{array}$$

Our choice for $\mathcal{P}$ is the Gaussiangaussian distribution with a meanmean $\mu$ and variancevariance $\sigma$ where $\mathcal{P}(p):=\mathrm{Nor}(\mu ,\sigma )$. That is, we are assuming that the output words $y\in V$ are normally distributed conditioned on the input sentences $x\in V^{\ast }$.

whose parameters in turn are modeled as $f_{\theta }(x)$ which is chosen to be the perceptronperceptron, a single layer neural network of weightsweights $w$ and biasesbiases $b$ where $p:=f_{\theta =\{w,b\}}(x)$. That is, not only are we assuming that the outputs $y\in \mathcal{Y}$ are normally distributed conditioned on the outputs $x\in \mathcal{X}$, but also that the parameters to such a distribution can be modeled by a neural network $f_{\theta =\{w,b\}}(x)$.

where $f$ in turn is a function characterized or determined by it’s own set of parameters $\theta$ in addition to the input $x$, and is denoted as $f_{\theta }(x)$ or $f(x;\theta )$ to reflect such. That is, $f:\Theta \times \mathcal{X}\to \mathcal{Y}$ is some function that is curried with parameters $\theta$ and returns a map from inputs $x$ to outputs $y$, which is reflected in borscht.nn.Module where a model is first instantiated with Python’s builtin __init__, and subsequently evaluated with __call__.

With deep learning, $f_{\theta }(x)$ often takes the form of a neural network parameterized by weightsweights and biasesbiases such that $\theta :=\{w,b\}$. In §1.3 and §1.4 we will understand and train single layer networks for $f_{\theta }(x)$ to prepare us for the various network architectures explored in §3. Self-Supervised Sequence Learning with Deep Networks.

Figure 1.3.1

So, after the model family or class has been specified via parameterization, a specific model is implemented or instantiated by continuing to score the quality of the output with the likelihood (TODO) but rather as a function of parameters and moreover, optimizingoptimization said scoring function $\ell (\theta )$ and denoted by $\mathrm{argmin}$. This process is referred to as parameter estimationparameter estimation, or more specifically, maximum likelihood estimationmaximum likelihood estimation (MLE)

$$\begin{array}{rl}& \mathcal{L}(\theta ):=\stackrel{\text{step two}}{\overbrace{\ell (\theta )}}\circ {\hat{f}}_{\theta }(x)\\ \underset{\text{step three}}{\underbrace{\underset{\theta }{\mathrm{argmin}}}}& \mathcal{L}(\theta )\end{array}$$

Figure 1.3.2

1.4.1 Probabilistic Architecture with Logistic Regression (Old)

$\hookleftarrow$ Table of Contents (§1.3 Parameterizing Classification with Discriminant Functions)

In contrast, in this chapter and for deep nets more generally, we will estimate $\hat{f}:V\to V$ by parameterizingparameters probabilities $p(y^{(n+1)}|x^{(n+1)})$ rather than counting them, namely with the weightsweights $w_i$ and biasesbias $b$ of a logistic regressionlogistic regression modelPerhaps logistic classification serves as better terminology, but we will shortly see why it’s called regression which can also be viewed as a single neuronneuron Stanford CS109 Probability for Computer Scientists | Logistic Regression | Lecture 24 Stanford CS229: Machine Learning | Lecture 5 - Perceptron and Logistic Regression from a deep net. That is, some function $\hat{f}:=\text{todo}(\sigma (\sum _{i=1}^d{w}_i{x}_i))$, which we will take the remainder of §1.3.1 Architecture with Logistic Regression to better understand.

A circuit for the logistic regression model is shown on the right hand side of Figure 1.3.1. Rather than produce the orange output $p(y|\mathbf{x})$ with a stochastic matrix containing normalized relative frequency counts, the logistic regression model produces it by passing the blue input $\mathbf{x}$ into the green summing and squashing functions $\sum$ and $\sigma$ with the sum being configured — in this case weighted — with purple parameters $w$ and $b$. Intuitively speaking, by parameterizing $\hat{f}$ what we mean is that we will assume or model $\hat{f}$ as some function that has configurable parameters. You can view the logistic regression model on the right as adding some indirection between the blue input and orange output with green machinery configured by purple parameters. The problem of training a good language model then reduces to that of finding a good set of parameters, in this case $\theta :=\{\mathbf{w},b\}$. That is, the process of training is parameter estimation.

Scanning from left to right of the neuronal circuit depicted in Figure 1.3.1, we first see that the input word "to" to the conditional probability that the logistic regression is modeling is some column vector $\mathbf{x}$ with support in ${\mathbb{R}}^{|V|}$. So, starting from the string "to", we continue to map the string to it’s index $i$ in the vocabulary $V$ with some hash map $I:V\to \mathbb{N}$, but then we need an extra step that maps the index $i\in \mathbb{N}$ to some vector $\mathbf{x}\in {\mathbb{R}}^{|V|}$. If you pay closer attention to the circuit, you will see that only ${\mathbf{x}}_2:=1$ whereas all the other ${\mathbf{x}}_i:=0$. This is called a one-hot vectorone-hot vector, which semantically represents the same positional information within the vocabulary $V$ except vectorizes that information continuously with support in ${\mathbb{R}}^{|V|}$ rather than discretely in $\mathbb{N}$. Because of this, the conditional probability is denoted as $p(y|\mathbf{x})$. Like 1.2_bigram.ipynb, the first cell in 1.3_logistic_regression.ipynb performs the data processing, except with the additional step of tensorizing integer indices over to one-hot vectors with F.one_hot(torch.tensor(indices)):

missing notebook: /opt/build/repo/sitp/notebooks/1.3_logistic_regression.ipynb

missing notebook: /opt/build/repo/sitp/notebooks/1.2_bigram.ipynb

So far we in the logistic regression notebook we have the same dataset, vocab and c2i. But because logistic regression takes in continuous inputs rather than discrete, one more step must be perform in its data processing, namely mapping (tensorizing) integer indices to their one-hot vectors (local, however). We construct list[int] first and convert to tensor’s second, for tensor.append()… (todo)

missing notebook: /opt/build/repo/sitp/notebooks/1.3_logistic_regression.ipynb

Where first five ${\mathbf{x}}^{(i)}$’s are plotted as a heatmap with plt.imshow. Continuing our left to right scan of the logistic regression model depicted in Figure 1.3.1, the next aspect to the neuronal circuit is the core machinery which the one hot vectors are passed through, namely taking a weighted sum $\sum _{i=1}^{|V|}{w}_i{x}_i$ whose real-valued output is then mapped to a number between 0 and 1 with a squashing function $\sigma (\cdot )$. That is, the logistic regression assumption is $P(Y=1|X=\mathbf{x}):=\sigma ({\mathbf{w}}^{\top }\mathbf{x})$.

Figure I.2.3

$$\begin{array}{rlrl}{f}_{\text{sum}}:{\mathbb{R}}^d& \to & & \mathbb{R}\\ \mathbf{x}& \mapsto & & ({\mathbf{w}}_1{\mathbf{x}}_1+{\mathbf{w}}_2{\mathbf{x}}_2+\cdots +{\mathbf{w}}_n{\mathbf{x}}_n)+b\\ & & & =(\sum _{i=1}^d{\mathbf{w}}_i{\mathbf{x}}_i)+b\\ & & & ={\mathbf{w}}^{\top }\mathbf{x}+b\\ f_{\text{squash}}:\mathbb{R}& \to & & [0,1]\\ z& \mapsto & & \sigma (z)\\ & & & =\frac{1}{1+\exp (-z)}\end{array}$$

That is, we model $\hat{f}$ as follows:

$$\begin{array}{rlrl}{\hat{f}}_{\mathbf{w}}:{\mathbb{R}}^d& \to & & [0,1]\\ \mathbf{x}& \mapsto & & f_{\text{squash}}(f_{\text{sum}}(\mathbf{x}))\end{array}$$

One small adjustment we can make to $f$ is to fold the final bias term $b$ into the vector $\mathbf{w}$ as ${\mathbf{w}}_1$, increasing the total dimensionality so that $\mathbf{w}\in {\mathbb{R}}^{d+1}$, and adjusting $\mathbf{x}$ accordingly so that ${\mathbf{x}}_1:=1$ and $\mathbf{x}\in {\mathbb{R}}^{d+1}$. Then, the equation of our line with bias is simply parameterized as $f(\mathbf{x};\mathbf{w})$ $$\begin{array}{rl}{f}_{\text{sum}}:{\mathbb{R}}^{d+1}& \to \mathbb{R}\\ \mathbf{x}& \mapsto {\mathbf{w}}^{\top }\mathbf{x}\end{array}$$

$f(\mathbf{x};\mathbf{w},b)$

missing notebook: /opt/build/repo/sitp/notebooks/1.3_logistic_regression.ipynb

1.4.2 Optimizing Cross Entropy with Iterative Reweighted Least Squares

$\hookleftarrow$ Table of Contents (§1.3 Parameterizing Classification with Discriminant Functions)

But what is cross-entropy? | Compression is Intelligence Part 2 (3Blue1Brown 2026)

missing notebook: /opt/build/repo/sitp/notebooks/1.3_logistic_regression.ipynb

1.4.3 Evaluation with Bias Variance Tradeoff

1.4.4 Summary and Further Reading

(Eisenstein 2019)’s Chapter 2. Linear text classification; (Jurafsky and Martin 2026)’s Chapter 4. Logistic Regression

(Hastie et al. 2009)’s Chapter 4. Linear Methods for Classification; (Bishop 2006)’s Chapter 4. Linear Methods for Classification; (Goodfellow et al. 2016)’s Chapter 5. Machine Learning Basics; (Murphy 2022)’s Chapter 10. Logistic Regression and Chapter 12. Generalized Linear Models

(MacKay 2003)’s Chapter 1. Introduction to Information Theory and Chapter 2. Probability, Entropy, and Inference; (Murphy 2022)’s Chapter 6: Information Theory

(further reading for torch)

1.4.4 bogus below

1.4.4 Architecture with Linear Regression

Recall from §1.2.1 The Estimation of Software 2.0 that the regime of regression is when the function to be estimated is of the form $f:\mathcal{X}\to \mathcal{Y}$ where $\mathcal{Y}:=\mathbb{R}$.

Similarly to learning distributions $X$ from data, learning functions $f:\mathcal{X}\to \mathcal{Y}$ from data roughly follows the three step process of selecting a model for the task $T$, a performance measure $P$, and optimizing such measure on experience $E$. Let’s take a look at the linearly modeling regression with squared error loss, a problem which straddles tools from all three languages of probability theory, linear algebra, and calculus.

For the task of house price prediction, the input space $\mathcal{X}={\mathbb{R}}^d$ is the size in square feet, and the output space $\mathcal{Y}=\mathbb{R}$ is the price, meaning that the function which needs to be recovered has the type of $f:{\mathbb{R}}^d\to \mathbb{R}$. An assumption about the structure of the data needs to be made, referred to as the inductive bias. The simplest assumption to make is to assume that there exists some linear relationship between the data $\mathbf{x}$ and parameters $\mathbf{w}$ so that $f$ ends up being modeled as an inner product plus bias, parameterized as $f(\mathbf{x};\mathbf{w},b)$: $$\begin{array}{rl}f& :{\mathbb{R}}^d\to \mathbb{R}\\ f(\mathbf{x};\mathbf{w},b)& :=({\mathbf{w}}_1{\mathbf{x}}_1+{\mathbf{w}}_2{\mathbf{x}}_2+\cdots +{\mathbf{w}}_n{\mathbf{x}}_n)+b\\ & =(\sum _{i=1}^d{\mathbf{w}}_i{\mathbf{x}}_i)+b\\ & ={\mathbf{w}}^{\top }\mathbf{x}+b\end{array}$$

One small adjustment we can make to $f$ is to fold the final bias term $b$ into the vector $\mathbf{w}$ as ${\mathbf{w}}_1$, increasing the total dimensionality so that $\mathbf{w}\in {\mathbb{R}}^{d+1}$, and adjusting $\mathbf{x}$ accordingly so that ${\mathbf{x}}_1:=1$ and $\mathbf{x}\in {\mathbb{R}}^{d+1}$. Then, the equation of our line with bias is simply parameterized as $f(\mathbf{x};\mathbf{w})$ $$\begin{array}{rl}f& :{\mathbb{R}}^{d+1}\to \mathbb{R}\\ f(\mathbf{x};\mathbf{w})& :={\mathbf{w}}^{\top }\mathbf{x}\end{array}$$

Semantically speaking, each property of the input house ${\mathbf{x}}_i\in {\mathbb{R}}^d$ is weighted by some weight ${\mathbf{w}}_i$, adjusted accordingly during the learning process so that it accurately reflects the property’s influence on the output price $y\in \mathbb{R}$ when evaluating a prediction $\hat{y}$. We can evaluate $f$ on all ${\mathbf{x}}^{(i)}\in \mathcal{D}$ with a randomly intialized $\mathbf{w}$ to ensure we’ve wired everything up correctly.

However, while specifying the equation of the linear regression model on paper screen and subsequently evaluating the computation manually might have been sufficient for pre-computational mathematicians such as Gauss and Legendre, the discipline of statistical learning is entirely predicated on computational methods. Let’s continue to use PyTorch’s core n-dimensional array datastructure with torch.Tensor, this time modeling a function $f:\mathcal{X}\to \mathcal{Y}$ rather than some distribution $X$:

(todo: maybe change example to closure following torch.nn)

{{#include ../../teeny/examples/2.2-regression.py:1:4}}

in which we can update the names of our torch.Tensor variables with the ranks of their vector spacesReferred to as Shape Suffixes, described by Noam Shazeer, an engineer from Google to clarify what dimensions these operations are evaluated on:

{{#include ../../teeny/examples/2.2-regression.py:5:8}}

and finally, actually enforce them via runtime typechecking with the jaxtyping package (try swapping the return statements below to trigger a type error with the return type):

{{#include ../../teeny/examples/2.2-regression.py:9:19}}

Let’s now sanity-check that f_typechecked is wired up correctly by evaluating it on all ${\mathbf{x}}^{(i)}\in \mathcal{D}$ with random parameter w_D. (todo: use actual housing data)

{{#include ../../teeny/examples/2.2-regression.py:29:41}}

random weight vector: tensor([-0.3343,  0.0768,  2.7828, -0.1331])
expected: $500000.0, actual: $-492.46
expected: $800000.0, actual: $-690.89
expected: $250000.0, actual: $-258.08

Clearly, these outputs are unintelligible, and will only become intelligible with a “good” choice of parameter $\mathbf{w}$. Before we find such a parameter, let’s make one more modification to our function $f$ to increase it’s performance. Rather then evaluating f_typechecked on every input xi_D in the dataset zip(X_ND,Y_N) sequentially with a loop, we can evaluate all outputs with a single matrix-vector multiplication by updating $f$’s function definition to the following

$$\begin{array}{rl}{f}_{\text{batched}}& :{\mathbb{R}}^{n\times d}\to {\mathbb{R}}^n\\ f_{\text{batched}}(\mathbf{X};\mathbf{w})& :=\mathbf{Xw}\\ & =\left[\begin{array}{c}{\mathbf{x}}_1^{\top }\mathbf{w}\\ {\mathbf{x}}_2^{\top }\mathbf{w}\\ ⋮\\ {\mathbf{x}}_n^{\top }\mathbf{w}\end{array}\right]\end{array}$$

and the corresponding PyTorch updated to

{{#include ../../teeny/examples/2.2-regression.py:21:27}}

Now that we’ve selected our model for the task $T$, we can proceed with selecting our performance measure $P$ which the learner will improve on with experience $E$.

1.4.5 Loss Function with Least Squares

After the inductive bias on the family of functions has been made, the learning algorithm must find the function $\hat{f}$ with a good fit. Since artificial learning algorithms don’t have visual cortex like biological humans[], the notion of “good fit” needs to defined in a systematic fashion. This is done by selecting the parameter $\theta \in {\mathbb{R}}^m$ which maximizes the likelihood of the data $p(d;\theta )$. Returning to the linear regression inductive bias we’ve selected to model the house price data, we assume there exists noise ${ϵ}^{(i)}$ in both our model (epistemic uncertainty) and data (aleatoric uncertainty), so that $y^{(i)}={\theta }^{\top }{\mathbf{x}}^{(i)}+{ϵ}^i$ where ${ϵ}^{(i)}\sim \mathcal{N}(\mu ,{\sigma }^2)$