The Numerical Linear Algebra of Large Language Models

Abdelkader Baggag; Yousef Saad

Figure 1: Major milestones in the evolution of artificial intelligence. The timeline highlights selected developments that shaped the field, from early perceptrons and expert systems to deep neural networks, Transformers, and the current era of large language models and generative AI. It is intended as a schematic historical overview rather than an exhaustive chronology.

Figure 2: Growth of parameter sizes in large language models. The plot shows the rapid increase in parameter counts across representative models over time, from early Transformer-era systems to recent frontier models. The horizontal reference lines and the indicated “data wall” emphasize that simple parameter scaling cannot continue indefinitely and may eventually be constrained by data availability and practical training considerations.

Figure 3: Order-of-magnitude argument for recent AI progress, adapted from [ 13 ] . The figure illustrates the claim that advances in effective compute, algorithmic efficiency, and related “unhobbling” factors may combine multiplicatively, leading to rapid capability gains over relatively short time scales. It is included as a schematic argument about scaling trends, not as a precise predictive model.

Figure 4: Main ingredients of the recent AI wave. The schematic emphasizes the interaction among three mutually reinforcing factors: improved hardware, the availability of large amounts of data, and advances in methods and algorithms. Together, these ingredients explain much of the rapid progress of modern deep learning and large language models.

Figure 5: Components of artificial intelligence: the big picture. The schematic situates large language models within the broader AI landscape and highlights the interplay among core ingredients such as data, optimization, numerical linear algebra, hardware, and software frameworks. Its purpose is to emphasize that modern AI systems emerge from the interaction of several mathematical, algorithmic, and computational components rather than from a single isolated idea.

Figure 6: A simple multilayer perceptron with three hidden layers. The input vector

\bm{x}

is propagated through successive affine transformations and nonlinear activation functions to produce the final output

\phi(\bm{x})

. The figure is intended to illustrate the layered composition that underlies feed-forward neural networks and serves here mainly to establish notation.

Figure 7: A simple multilayer perceptron with two hidden layers. The figure illustrates the layer structure used in the classification example and, together with Table 1 , helps explain the dimensions of the weight matrices and the total number of trainable parameters.

	Layer 0 (=in) to 1	Layer 1 to 2	Layer 2 to 3 (=out)
$\text{size}(W_{l})$	(400 $\times$ 100)	(100 $\times$ 100)	(100 $\times$ 10)
$\text{size}(b_{l})$	(100 $\times$ 1)	(100 $\times$ 1)	(10 $\times$ 1)
Total Params	40,100	10,100	1,010

Table 1: Dimensions and parameter counts for a simple multilayer perceptron. The table lists the sizes of the weight matrices and bias vectors for each layer in the digit-classification example, together with the corresponding number of trainable parameters. It illustrates how the total parameter count is determined directly by the layer widths. The total number of parameters is

51{,}210

Figure 8: A simple computational graph for composite-function evaluation. The leaf nodes

x,y

, and

z

are primary variables; the intermediate nodes

a(x,y,z),b(x,y,z)

and

c(x,y,z)

are derived quantities; and the root node is the target function. The figure illustrates how forward evaluation computes intermediate quantities first, while reverse-mode differentiation (back-propagation) propagates derivatives backward from

f

to the leaves via the chain rule.

Figure 9: Local back-propagation step at an intermediate node in a computational graph. Assuming the partial derivatives of the objective

f

with respect to the downstream nodes

a_{j},a_{\ell}

, and

a_{m}

are already known, the derivative with respect to the current node

a_{k}

is obtained by applying the chain rule and summing the contributions along all outgoing paths. The figure illustrates the recursive local structure that underlies reverse-mode automatic differentiation.

Figure 10: A local computational-graph view of one neural-network layer. The figure shows how the pre-activation

s_{\ell}=W_{\ell}^{T}z_{\ell-1}+b_{\ell}

is formed from the previous-layer representation

z_{\ell-1}

, the weight matrix

W_{\ell}

, and the bias

b_{\ell}

and is then mapped by the activation function

\sigma

to produce

z_{\ell}

. It illustrates how each layer of a neural network can be interpreted as a small computational subgraph whose variables participate in forward evaluation and back-propagation. The circled nodes represent input or parameter variables, whereas the rectangular nodes represent computational operations or derived quantities.

Figure 11: Backward propagation in a multilayer perceptron. The figure schematically illustrates the reverse flow of derivatives through the layered computational graph associated with equations ( 2.6.9 )–( 2.6.12 ). Starting from the loss

\eta

, sensitivities are propagated backward through the intermediate states

s_{3},s_{2},s_{1}

in order to compute gradients with respect to the parameters

(W_{3},b_{3}),(W_{2},b_{2})

, and

(W_{1},b_{1})

.

Figure 12: Recovering an underlying function from noisy data. The figure illustrates the basic problem used to motivate attention: given a cloud of noisy training points

(x_{i},y_{i})

, estimate the value of an unknown function at a new query point. In the numerical-analysis viewpoint this suggests interpolation or least-squares approximation, whereas in the machine-learning viewpoint it motivates weighted averaging based on relevance, i.e., attention.

Figure 13: Illustration of Nadaraya–Watson kernel regression. Given a query point

\bm{q}

, the prediction is formed as a weighted average of the observed values

{\bm{v}}_{i}

, with weights that depend on the similarity between

\bm{q}

and the keys

{\bm{k}}_{i}

. This provides a useful bridge to attention, where output representations are likewise obtained by similarity-based weighted averaging.

Figure 14: Nadaraya–Watson kernel regression on periodic data. The figure compares an attention-based kernel-regression estimate with a classical low-degree approximation for an oscillatory target function. In this setting, the similarity-based weighted average captures the periodic behavior more accurately than a simple cubic approximation, illustrating how attention can recover useful structure directly from data when the functional form is not known in advance.

Figure 15: High-level view of autoregressive language modeling. Given a context

u_{\,<\,t}=(u_{1},\dots,u_{t-1})

, the model assigns a probability distribution to the next token

u_{t}

, and text generation proceeds by repeating this prediction step sequentially.

Figure 16: Schematic pipeline of a decoder-only large language model. Tokens are mapped to embeddings with positional structure, propagated through repeated Transformer blocks, and finally projected to logits for next-token prediction.

Figure 17: Sinusoidal positional encoding coordinates across dimensions. Each subplot shows the two coordinates

{\bm{P}}(i,2j)=\sin\left(i/10000^{2j/d_{\text{model}}}\right)

and

{\bm{P}}(i,2j+1)=\cos\left(i/10000^{2j/d_{\text{model}}}\right)

as functions of token position

i

for a representative pair of dimensions. Low-index dimensions vary rapidly with position, while higher-index dimensions vary more slowly, so the encoding provides a multiscale representation of sequence order.

Figure 18:

\operatorname{Pre-LayerNorm}

Transformer block at layer

\ell

. The block applies layer normalization before the attention and MLP sublayers, and writes both outputs additively into the same residual stream. This makes the block an identity-plus-correction map composed of token mixing through attention and channel mixing through the MLP.

Figure 19: Encoder-decoder Transformer architecture. The encoder computes source-side contextual representations, while the decoder combines masked self-attention and encoder-decoder attention to generate target-side outputs.

Figure 20: Decoder-only Transformer architecture. Repeated masked self-attention and feed-forward updates act along a shared residual stream, producing contextual representations used for autoregressive prediction.

Figure 21:

\operatorname{Post-LayerNorm}

and

\operatorname{Pre-LayerNorm}

Transformer blocks. In the

\operatorname{Post-LayerNorm}

form, normalization is applied after each residual addition, whereas in the

\operatorname{Pre-LayerNorm}

form, it is applied before each sublayer. The comparison highlights that

\operatorname{Pre-LayerNorm}

preserves an exact identity path along the residual stream, which is one reason it is favored in modern LLMs.

Figure 22: Generic residual unit in Post-LayerNorm and Pre-LayerNorm form. In both cases the sublayer contributes a learned correction

\mathcal{F}

to the incoming state, but the location of normalization differs: Post-LayerNorm applies normalization after the residual sum, whereas Pre-LayerNorm applies it before the sublayer.

		GPT-3 175B	Llama-3 405B
Vocabulary Size	$n_{\texttt{vocab}}\equiv\|\mathcal{V}\|$	$50,257$	$128,256$
Layers	$n_{\texttt{blocks}}$	96	126
Model Dimension	$d_{\texttt{model}}$	$12,288$	$16,384$
MLP Dimension	$d_{\texttt{mlp}}$	$49,152$	$53,248$
Attention Heads	$n_{\texttt{heads}}$	96	128
Key/Value Heads			8

Table 2: Comparison of key architectural hyperparameters for GPT-3 175B and Llama-3 405B. The table shows that, although both models follow the decoder-only Transformer template, they differ substantially in vocabulary size, depth, model width, feed-forward dimension, head configuration, and activation design.

	GPT-3 175B		# parameters
$\operatorname{MultiHead}$	$(3d_{\texttt{model}}d_{\texttt{head}}+d_{\texttt{head}}d_{\texttt{model}})\times n_{\texttt{heads}}\penalty\ \times n_{\texttt{blocks}}$	=	$57,982,058,496$
$\operatorname{MLP}$	$(2d_{\texttt{mlp}}d_{\texttt{model}}+(d_{\texttt{mlp}}+d_{\texttt{model}}))\penalty\ \times n_{\texttt{blocks}}$	=	$115,970,015,232$
$\bm{W}_{E},\bm{W}_{U}$	$(2d_{\texttt{model}}n_{\texttt{vocab}})$	=	$1,235,116,032$
$\bm{\gamma},\bm{\beta}$	$2\times(2d_{\texttt{model}})\times n_{\texttt{blocks}}$	=	$4,718,592$
	Total	=	$175,191,908,352$

Table 3: Parameter breakdown of GPT-3 175B by major architectural component. The table decomposes the total parameter count into attention, MLP, embedding/unembedding, and normalization contributions, illustrating that the dense MLP blocks dominate the parameter budget.

	Llama-3 405B		# parameters
$\operatorname{MultiHead},{\bm{W}}_{K},{\bm{W}}_{V}$	$2\times d_{\texttt{head}}d_{\texttt{model}}\times n_{\texttt{blocks}}$	=	$4,227,858,432$
$\operatorname{MultiHead},{\bm{W}}_{Q},{\bm{W}}_{O}$	$2\times d_{\texttt{model}}d_{\texttt{model}}\times n_{\texttt{blocks}}$	=	$67,645,734,912$
$\operatorname{MLP}$	$(3d_{\texttt{mlp}}d_{\texttt{model}}+(2d_{\texttt{mlp}}+d_{\texttt{model}}))\penalty\ \times n_{\texttt{blocks}}$	=	$329,788,440,576$
${\bm{W}}_{E},{\bm{W}}_{U}$	$(2d_{\texttt{model}}n_{\texttt{vocab}})$	=	$4,202,692,608$
$\bm{\gamma},\bm{\beta}$	$2\times(2d_{\texttt{model}})\times n_{\texttt{blocks}}$	=	$8,257,536$
	Total	=	$405,872,984,064$

Table 4: Parameter breakdown of Llama-3 405B by major architectural component. The table decomposes the total parameter count into attention, MLP, embedding/unembedding, and normalization contributions, illustrating that the gated MLP blocks account for the dominant share of parameters in this model.

Figure 23: Training and inference in a decoder-only Transformer. During training, the model processes an entire token sequence in parallel and adjusts its weights so that the output probabilities match the target next tokens at all positions. During inference, by contrast, generation is autoregressive: the predicted next token is fed back as input, and the process is repeated sequentially.

		GPT-2	GPT-3	GPT-4	Llama-3 405B
Model Dimension	$d_{\texttt{model}}$	768	$12,288$	$12,288$	$16,384$

Table 5: Representative model dimensions illustrating the scale of modern embedding spaces. The table lists the model dimension

d_{\text{model}}

for GPT-2, GPT-3, GPT-4, and Llama-3 405B, showing that contemporary large language models operate in ambient spaces of dimension on the order of

10^{4}

. Such dimensions help explain why many randomly chosen directions are nearly orthogonal, as illustrated in Figure 24 .

Figure 24: Concentration of pairwise angles in high dimension. The figure shows the distribution of angles between 10,000 random vectors in a 12,288-dimensional embedding space. The sharp concentration around 90°, together with the narrow 88°–92° band, illustrates the near-orthogonality phenomenon that characterizes random directions in high-dimensional spaces.

Component $i$	$\sigma_{i}$ (Energy)	Cumulative
1	0.40	0.40
2	0.20	0.60
3	0.10	0.70
4	0.08	0.78
5	0.06	0.84
6-10	0.03 (each)	0.99

Table 6: Illustrative low-rank concentration of gradient energy. The table shows a normalized singular-value distribution

\sum_{i}\sigma_{i}=1

, for which the cumulative energy exceeds 90% within roughly 10 components. Despite the ambient parameter dimension

D=110\times 10^{6}

, most of the gradient energy is concentrated in a very low-dimensional subspace, illustrating the motivation behind low-rank gradient methods.

Figure 25: Pre-LN Transformer layer as residual-stream accumulation. The input residual stream

{\bm{x}}_{\ell-1}

is updated by adding the attention-block contribution

{\bm{a}}_{\ell}

and the MLP-block contribution

{\bm{m}}_{\ell}

, giving

{\bm{x}}_{\ell}={\bm{x}}_{\ell-1}+{\bm{a}}_{\ell}+{\bm{m}}_{\ell}

. Across layers, the residual stream can be expanded as

{\bm{x}}_{\ell}={\bm{x}}_{0}+\sum_{i=1}^{\ell}{\bm{a}}_{i}+\sum_{i=1}^{\ell}{\bm{m}}_{i}

, showing how attention and MLP updates accumulate additively through the network.

Figure 26: Empirical power-law scaling in large language models. Test loss decreases

L(N,D)

predictably as model size

N

, dataset size

D

, and compute budget

C

increase, following approximate power-law relationships. The schematic summarizes the scaling-law perspective that performance improves smoothly with scale, with diminishing returns captured by the exponents in

L(N,D)

defined in ( 3.16.1 )

Figure 27: Teacher–student knowledge distillation. A pretrained teacher model produces logits and softened output probabilities for the training data, and a smaller student model is trained to match these outputs while also fitting the original targets. The distillation loss transfers information from the teacher’s predictive distribution to the student, allowing the student to imitate the teacher’s behavior with a more compact model.

Figure 28: Forward step in LoRa. A frozen pretrained weight matrix

\bm{W}_{0}

is augmented by a trainable low-rank update

{\bm{A}}{\bm{B}}

, where

{\bm{A}}\in{\mathbb{R}}^{d\times r}

,

{\bm{B}}\in{\mathbb{R}}^{r\times d_{v}}

, and

r\ll\min(d,d_{v})

. Instead of updating the full weight matrix, LoRa keeps

\bm{W}_{0}

fixed and learns only the small-rank factors

\bm{A}

, and

\bm{B}

, thereby reducing the number of trainable parameters while preserving the original forward structure.

Figure 29: Illustration of linearized attention. If the softmax is temporarily omitted, the attention map takes the form

{\bm{A}}_{\ell}^{(\text{lin})}({\bm{X}}_{\ell-1})={\bm{Q}}_{\ell}({\bm{K}}_{\ell}^{T}{\bm{V}}_{\ell})/{\sqrt{d_{k}}}

, which can be evaluated by first forming the smaller matrix

{\bm{K}}_{\ell}^{T}{\bm{V}}_{\ell}\in{\mathbb{R}}^{d\times d_{v}}

and then multiplying on the left by

{\bm{Q}}_{\ell}

. This reordering avoids the explicit formation of the

n\times n

score matrix and highlights the computational motivation behind so-called linear Transformers.

Figure 30: Gradient descent versus natural gradient descent on a two-dimensional Gaussian model with direction-dependent geometry. The left panel shows the decay of the objective

L(w)

over the iterations, while the right panel displays the corresponding trajectories on the loss surface. Standard gradient descent is slowed by the distorted geometry of the parametrization, whereas natural gradient descent accounts for this geometry through the Fisher metric and converges much more rapidly. Standard (Cauchy) steepest descent is shown by blue plus markers, while natural gradient descent is shown by red circles.

	$\displaystyle L_{t}$	$\displaystyle=L_{t-1}+G_{t}G_{t}^{T}$		(5.4.1)
	$\displaystyle R_{t}$	$\displaystyle=R_{t-1}+G_{t}^{T}G_{t}$		(5.4.2)

Algorithm 5 Shampoo Preconditioning – matrix case

The Numerical Linear Algebra of Large Language Models

Abstract

Figures & tables

Advanced Linear Algebra with Applications - Part I (Numerical linear algebra for PDEs, machine learning, and data assimilation)

Data-Driven Dynamic Algorithm Dispatch with Large Language Models

RL4RLA: Teaching ML to Discover Randomized Linear Algebra Algorithms Through Curriculum Design and Graph-Based Search

Abstract

Explore similar work

Advanced Linear Algebra with Applications - Part I (Numerical linear algebra for PDEs, machine learning, and data assimilation)

Data-Driven Dynamic Algorithm Dispatch with Large Language Models

RL4RLA: Teaching ML to Discover Randomized Linear Algebra Algorithms Through Curriculum Design and Graph-Based Search