In [ ]:
```python shell
import numpy as np
def strassen(A: np.ndarray, B: np.ndarray, cutoff: int = 128) -> np.ndarray:
assert A.shape[1] == B.shape[0]
n = max(A.shape + B.shape)
# pad to next power of two (square) for brevity
k = 1 << (n - 1).bit_length()
Ap = np.zeros((k, k), dtype=A.dtype)
Ap[: A.shape[0], : A.shape[1]] = A
Bp = np.zeros((k, k), dtype=B.dtype)
Bp[: B.shape[0], : B.shape[1]] = B
def _stras(X, Y):
m = X.shape[0]
if m <= cutoff:
return X @ Y
h = m // 2
X11, X12, X21, X22 = X[:h, :h], X[:h, h:], X[h:, :h], X[h:, h:]
Y11, Y12, Y21, Y22 = Y[:h, :h], Y[:h, h:], Y[h:, :h], Y[h:, h:]
M1 = _stras(X11 + X22, Y11 + Y22)
M2 = _stras(X21 + X22, Y11)
M3 = _stras(X11, Y12 - Y22)
M4 = _stras(X22, Y21 - Y11)
M5 = _stras(X11 + X12, Y22)
M6 = _stras(X21 - X11, Y11 + Y12)
M7 = _stras(X12 - X22, Y21 + Y22)
Z11 = M1 + M4 - M5 + M7
Z12 = M3 + M5
Z21 = M2 + M4
Z22 = M1 - M2 + M3 + M6
Z = np.empty_like(X)
Z[:h, :h], Z[:h, h:], Z[h:, :h], Z[h:, h:] = Z11, Z12, Z21, Z22
return Z
Cp = _stras(Ap, Bp)
return Cp[: A.shape[0], : B.shape[1]]
A = np.random.randn(200, 150)
B = np.random.randn(150, 200)
# Check precision
assert np.allclose(strassen(A, B), A @ B)
strassen(A, B)
```
> \[!note\] when to use strassen
>
> Use only for very large dense square-ish matrices, good arithmetic intensity, and when you can tolerate a modest increase in rounding error. Otherwise, high-quality blocked GEMM is typically faster and more stable.
> \[!warning\] Common pitfalls
>
> - Forming $A\otimes B$ explicitly explodes memory (sizes multiply). Use the vec/Kronecker identities.
> - Confusing $A\circ B$ with $AB$: elementwise vs. composition. The Frobenius pairing relates them via $\langle A\circ B, C\rangle=\langle A, B\circ C\rangle$ but they are different operations.
> - Ignoring data layout: mismatched row/column-major assumptions or unaligned packing can erase SIMD/cache benefits.
### Chain order and blocking
- Matrix-chain ordering affects cost; dynamic programming finds the cheapest parenthesization for a fixed chain of shapes.
- Blocking (tiling) yields cache-friendly performance; BLAS uses blocked algorithms under the hood.
### Worked examples
- Vec/Kronecker: let $A=\begin{bmatrix}1&2\\0&1\end{bmatrix}$, $B=\begin{bmatrix}2&0\\0&3\end{bmatrix}$, $X=\begin{bmatrix}1&4\\2&5\end{bmatrix}$. Then $AXB=\begin{bmatrix}10&42\\4&15\end{bmatrix}$ so $\operatorname{vec}(AXB)=[10,\,4,\,42,\,15]^\top$. Also $(B^\top\!\otimes A)\operatorname{vec}(X)=[10,\,4,\,42,\,15]^\top$.
- Chain order cost: $A(10\!\times\!100)$, $B(100\!\times\!5)$, $C(5\!\times\!50)$. $(AB)C$: $10\cdot100\cdot5+10\cdot5\cdot50=7{,}500$ mults; $A(BC)$: $100\cdot5\cdot50+10\cdot100\cdot50=75{,}000$.
- Parenthesization matters!!
## linear transformations
> \[!abstract\] Linear maps unify algebra and geometry
>
> A map $T:V\to W$ is linear if $T(u+v)=T(u)+T(v)$ and $T(\alpha v)=\alpha T(v)$. With bases $B_V,B_W$, $T$ has a matrix $[T]_{B_W\leftarrow B_V}$ so that $[T x]_{B_W}=[T]_{B_W\leftarrow B_V}[x]_{B_V}$.
### Structure and examples
- Kernel and image: $\ker T=\{x:T(x)=0\}$, $\mathrm{im}\,T=\{T(x)\}$ with $\dim\ker T+\dim\mathrm{im}\,T=\dim V$.
- Projections and reflections: $P=QQ^\top$ for orthonormal $Q$ (projection onto $\mathrm{span}(Q)$); reflections $R=I-2uu^\top$ for unit $u$.
- Rotations: in $\mathbb{R}^2$, $R_\theta=\begin{bmatrix}\cos\theta&-\sin\theta\\ \sin\theta&\cos\theta\end{bmatrix}$; in $\mathbb{R}^n$, orthogonal matrices with $R^\top R=I$.
- Similarity: changing basis in $V$ gives $[T]'=S^{-1}[T]S$; eigenvalues are basis-invariant.
> \[!note\] Linearization and Jacobian
>
> A differentiable $f:\mathbb{R}^n\to\mathbb{R}^m$ is locally $f(x_0+h)\approx f(x_0)+J_f(x_0)h$, where $J_f$ is the Jacobian; see [[thoughts/Vector calculus#Jacobian matrix|Jacobian]].
### Small projection example
Project onto $\mathrm{span}\{a\}$ with $a\ne0$: $P=\dfrac{a a^\top}{a^\top a}$, so $P^2=P$ and $P^\top=P$.
### Isomorphisms, adjoints, invariant subspaces
- Isomorphism: $T$ is invertible (an isomorphism) iff $\ker T=\{0\}$ and $\dim\mathrm{im}\,T=\dim V$ (square), equivalently $[T]$ has full rank.
- Adjoint (Euclidean case): $T^\dagger$ satisfies $\langle Tx,y\rangle=\langle x,T^\dagger y\rangle$; with the standard inner product, $[T^\dagger]=[T]^\top$. Symmetric operators ($T=T^\dagger$) admit orthonormal eigenbases; see [[lectures/411/notes#eigenvectors|eigenvectors]].
- Invariant subspace: $S\subseteq V$ with $T(S)\subseteq S$. Decomposing $V$ into invariant subspaces block-diagonalizes $[T]$ in a suitable basis (Jordan/Schur for general, orthogonal diagonalization for symmetric).
### Affine maps and change of basis
- Affine map: $f(x)=Ax+b$ is linear plus translation; linearity properties apply to the $Ax$ part. Solution sets to $Ax=b$ are affine translates of $\ker A$.
- Change of basis in domain/codomain: with bases $B_V,B_W$ and new bases $\tilde B_V,\tilde B_W$ related by $S_V,S_W$, the matrix changes by $[T]_{\tilde B_W\leftarrow \tilde B_V}=S_W^{-1}[T]_{B_W\leftarrow B_V} S_V$.
### Operator norm and exponential
- Operator 2-norm: $\|T\|=\sup_{\|x\|=1}\|Tx\|$ equals the top singular value of $[T]$.
- Matrix exponential: $e^{At}=\sum_{k\ge0}\frac{(At)^k}{k!}$; if $A=V\Lambda V^{-1}$ then $e^{At}=V e^{\Lambda t} V^{-1}$. Solves $\dot x=Ax$ (LTI systems), with stability governed by eigenvalues.
## eigenvectors
> \[!abstract\] Eigenvalues and eigenvectors
>
> For square $A\in\mathbb{R}^{n\times n}$, a nonzero $v$ and scalar $\lambda$ satisfy $Av=\lambda v$. The set of all eigenvectors with eigenvalue $\lambda$ (plus $0$) is the eigenspace $\mathcal{E}_\lambda=\ker(A-\lambda I)$.
> Left-eigenvectors satisfy $w^\top A = \lambda\, w^\top$ (equivalently, $A^\top w=\lambda w$); for symmetric/Hermitian $A$, left and right eigenvectors coincide up to transpose/conjugation.
### Diagonalization and multiplicities
- Algebraic multiplicity: multiplicity of $\lambda$ as a root of $\det(A-\lambda I)=0$.
- Geometric multiplicity: $\dim\ker(A-\lambda I)\le$ algebraic multiplicity.
- Diagonalizable iff the eigenvectors span $\mathbb R^n$ (equivalently, sum of geometric multiplicities is $n$), then $A=V\Lambda V^{-1}$.
### Spectral theorem (real symmetric/Hermitian)
- If $A=A^\top$ (or $A=A^{\dagger}$), then there exists an orthonormal basis of eigenvectors: $A=Q\,\Lambda\,Q^\top$ with $Q$ orthogonal and $\Lambda$ real diagonal.
- Eigenvectors for distinct eigenvalues are orthogonal. Quadratic form simplifies: $x^\top A x=\sum_i \lambda_i\,\alpha_i^2$ when $x=\sum_i \alpha_i q_i$.
- Positive (semi)definiteness: $A\succ 0$ ($\succeq 0$) iff all eigenvalues are $>$ ($\ge$) 0.
> \[!tip\] Power iteration
>
> For symmetric $A$ with $|\lambda_1|> |\lambda_2|\ge\cdots$, $x_{k+1}=\frac{Ax_k}{\|Ax_k\|}$ converges to the top eigenvector for almost any $x_0$ with a nonzero component on it; the Rayleigh quotient gives $\lambda_1$.
### Intuition first (geometry)
- What stays straight: an eigenvector is a direction $v$ that $A$ does not turn, scaling it only by $\lambda$ (and possibly flipping it if $\lambda<0$).
- Symmetric picture: for $A=A^\top$, the unit circle becomes an ellipse. The ellipse’s axes are the eigenvectors; their radii are $|\lambda_i|$.
- Quadratic form view: on the unit circle, $x^\top A x$ measures “how much $A$ prefers direction $x$”; the most preferred directions are the eigenvectors with largest/smallest eigenvalues.
- Dynamics: repeated application $A^k$ scales components along each eigenvector by $\lambda^k$. The largest $|\lambda|$ dominates long-run behavior; see [[lectures/411/notes#linear transformations|linear transformations]].
