Worked examples — NumPy linear algebra — np.linalg.solve, eig, svd, norm, det

This page is the exercise gym for the parent note. We do not re-teach the theory — we stress-test it. The goal: after this page there is no case — no sign, no zero, no collapse, no weird shape — that surprises you.

If any symbol below feels unfamiliar, the parent note and these prerequisites build it: Gaussian elimination and LU decomposition, Eigenvalues and diagonalization, Vector norms and metric spaces, Condition number and numerical stability, Least squares regression, NumPy arrays and broadcasting.

The scenario matrix

Every one of the five tools has normal cases and edge cases. The table below lists every class of scenario this topic can throw at you. Each worked example that follows is tagged with the cell(s) it covers. By the last example, every cell is green.

Figures accompany the geometric cells (C1, C4, C7, C8, C10).

C1 — Well-behaved square system

Look at the figure: each equation is a line in the $x$–$y$ plane. Two non-parallel lines cross at exactly one point — that crossing is the solution.

Steps.

1. Write it as $Ax=b$ with $A=\left(\begin{array}{cc}3& 2\\ 1& 2\end{array}\right)$, $b=\left(\begin{array}{c}12\\ 8\end{array}\right)$. Why this step? NumPy needs the coefficient grid $A$ and the right-hand column $b$ separated — the unknowns $x,y$ live in the vector we solve for.
2. Check the machine is invertible first: $\det A=(3)(2)-(2)(1)=6-2=4\ne 0$. Why this step? A nonzero determinant means the two column-arrows span real area (the figure's lines are not parallel), so a unique crossing exists.
3. Eliminate $x$: $R_2\to R_2-\frac{1}{3}{R}_1$ gives $2y-\frac{2}{3}y=8-4$, i.e. $\frac{4}{3}y=4\Rightarrow y=3$. Why this step? Subtracting a multiple of row 1 zeroes the first column below the pivot — exactly what Gaussian elimination and LU decomposition does inside solve.
4. Back-substitute: $x+2(3)=8\Rightarrow x=2$. Why this step? Once the system is triangular, the bottom row already gave $y$; we walk upward, plugging the known $y$ into the top row to release the last unknown $x$.

$\boxed{x=2,y=3}$

A = np.array([[3.,2.],[1.,2.]]); b = np.array([12.,8.])
np.linalg.solve(A, b)          # array([2., 3.])

C2 — Singular system ($\det =0$)

Steps.

1. Form $A=\left(\begin{array}{cc}2& 4\\ 1& 2\end{array}\right)$.
2. Determinant: $\det A=(2)(2)-(4)(1)=4-4=0$. Why this step? Row 2 is exactly half of row 1's coefficients ($2,4=2\times (1,2)$). The two column-arrows are parallel — the parallelogram has zero area, so the machine squishes the plane onto a line.
3. Because $\det =0$, np.linalg.solve raises LinAlgError: Singular matrix. Why? There is no unique crossing. Here the lines are parallel but not identical ($2x+4y=6$ says $x+2y=3$, but the second says $x+2y=5$) → no solution at all.

A = np.array([[2.,4.],[1.,2.]])
np.linalg.det(A)               # 0.0  (up to rounding)
# np.linalg.solve(A, [6.,5.])  -> raises LinAlgError: Singular matrix

C3 — Nearly singular (ill-conditioned)

Steps.

1. $\det A=(1)(1.0001)-(1)(1)=0.0001$ — nonzero, so a unique answer exists on paper. Why this step? We always test invertibility first (as in C1): a nonzero determinant promises a single crossing, so it is legitimate to ask solve for the answer.
2. Solve: np.linalg.solve returns $x=1,y=1$. Why this step? With invertibility confirmed, one solve call runs the Gaussian elimination of Gaussian elimination and LU decomposition and hands back the unique $(x,y)$ — but the tiny pivot warns us to keep our guard up.
3. Measure danger with the condition number $\kappa (A)=\Vert A\Vert \Vert A^{-1}\Vert$ (the symbol introduced at the top of the page), which asks "how much can a small error in $b$ blow up in $x$?" Why this tool and not det? A tiny det can just mean small scale, not instability. $\kappa$ is scale-independent — it directly measures error amplification (see Condition number and numerical stability).
4. Here $\kappa (A)\approx 4\times 10^4$: a change of $10^{-4}$ in $b$ can shift $x$ by order $1$. The answer is mathematically fine but numerically fragile.
5. Watch the amplification for real. Nudge $b$ by a hair to $b^{\prime }=(2,2.0002)$ — a change of only $10^{-4}$ in one entry. Re-solving gives $x=0,y=2$: the answer jumped by a whole unit. A $10^{-4}$ wobble in the input produced an $O(1)$ swing in the output — that is exactly what $\kappa \approx 4\times 10^4$ predicted. Why this step? Quoting a condition number is abstract; re-solving a perturbed system makes the fragility concrete and visible.

A = np.array([[1.,1.],[1.,1.0001]])
b  = np.array([2., 2.0001])
b2 = np.array([2., 2.0002])          # perturb b by 1e-4
np.linalg.solve(A, b)    # array([1., 1.])
np.linalg.solve(A, b2)   # array([0., 2.])  <- huge jump from a tiny nudge!
np.linalg.cond(A)        # ~ 4.0e4  -> large: distrust this solve

C4 — Negative determinant (orientation flip)

Steps.

1. $\det =(0)(0)-(1)(1)=-1$.
2. Magnitude $|-1|=1$: areas are preserved (a unit square still has area 1). Why this step? $|\det |$ is the volume-scaling factor from the parent note — here it scales by 1, no growth.
3. Sign $-1$: the transformation is a reflection. In the figure, the labelled corners $A\to B\to C$ run counter-clockwise before, and clockwise after — the orientation flipped, like a mirror. Why does sign encode this? $\det$ is a signed area; swapping the two column-arrows (which this matrix does — it swaps the axes) reverses the "handedness".

np.linalg.det(np.array([[0.,1.],[1.,0.]]))   # -1.0  (area kept, orientation flipped)

C5 — All the norms of one object (vector, Frobenius, and operator)

Steps.

1. ${\ell }_1={\sum }_i|v_i|=|-3|+|4|+|0|+|12|=3+4+0+12=19$. Why: ${\ell }_1$ is the total taxicab distance — add up every step regardless of sign.
2. ${\ell }_2=\sqrt{{\sum }_i{v}_i^2}=\sqrt{9+16+0+144}=\sqrt{169}=13$. Why: straight-line (Euclidean) length; the zero entry contributes nothing.
3. ${\ell }_{\infty }={\max }_i|v_i|=\max (3,4,0,12)=12$. Why: the single biggest component — the worst offender.
4. Ordering: ${\ell }_{\infty }\le {\ell }_2\le {\ell }_1$ always, so $12\le 13\le 19$ ✓.
5. Frobenius of $M$: flatten the whole matrix into one long vector and take its ${\ell }_2$ length: $\sqrt{1^2+(-2{)}^2+(-2{)}^2+2^2}=\sqrt{1+4+4+4}=\sqrt{13}\approx 3.606$. Why this tool? Frobenius treats the matrix as raw stored "energy" of all entries — cheap, but it ignores what the matrix does to vectors.
6. Operator (induced-2) norm of $M$: this answers a different question — "what is the biggest stretch factor $M$ applies to any unit vector?" That number is the largest singular value ${\sigma }_1$ (see C8 for singular values). For our symmetric $M$, ${\sigma }_1$ equals the largest $|\text{eigenvalue}|$. Its eigenvalues solve $(1-\lambda )(2-\lambda )-4=0\Rightarrow {\lambda }^2-3\lambda -2=0\Rightarrow \lambda =\frac{3\pm \sqrt{17}}{2}$, so ${\sigma }_1=\frac{3+\sqrt{17}}{2}\approx 3.562$. Why a separate norm? When you care about worst-case amplification of a system (stability, error bounds), Frobenius is the wrong ruler — the operator norm is the honest one. Notice the two matrix norms differ ($3.606$ vs $3.562$): they measure genuinely different things.

v = np.array([-3.,4.,0.,12.])
np.linalg.norm(v, 1)       # 19.0
np.linalg.norm(v)          # 13.0  (default ord=2)
np.linalg.norm(v, np.inf)  # 12.0
 
M = np.array([[1.,-2.],[-2.,2.]])
np.linalg.norm(M, 'fro')   # 3.6055...  'fro' is the API shorthand for the Frobenius norm
np.linalg.norm(M, 2)       # 3.5615...  ord=2 on a MATRIX = operator/induced-2 norm = largest singular value

C6 — Real distinct eigenvalues (symmetric matrix)

Steps.

1. Characteristic equation $\det (A-\lambda I)=0$, where $I$ is the identity matrix defined at the top of this page: $\det \left(\begin{array}{cc}2-\lambda & 1\\ 1& 2-\lambda \end{array}\right)=(2-\lambda {)}^2-1=0.$ Why this step? From the parent note: a fixed stretch direction exists only when $A-\lambda I$ collapses (is singular), i.e. its determinant is zero. We subtract $\lambda I$ (a scaled do-nothing machine) so the equation $Av=\lambda v$ becomes $(A-\lambda I)v=0$.
2. Solve $(2-\lambda {)}^2=1\Rightarrow 2-\lambda =\pm 1\Rightarrow \lambda =1$ or $\lambda =3$. Why this step? Taking the square root of both sides is how we undo the squaring in step 1 — and a square equals $1$ in exactly two ways ($+1$ and $-1$), which is precisely why a $2\times 2$ matrix has (up to) two eigenvalues.
3. For $\lambda =3$: $(A-3I)v=0\Rightarrow \left(\begin{array}{cc}-1& 1\\ 1& -1\end{array}\right)v=0\Rightarrow v=(1,1)$. For $\lambda =1$: $\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right)v=0\Rightarrow v=(1,-1)$. Why: we seek the direction the machine only stretches; solving $(A-\lambda I)v=0$ finds it.
4. Check perpendicular: $(1,1)\cdot (1,-1)=1-1=0$ ✓ — symmetric matrices always give orthogonal eigenvectors (Eigenvalues and diagonalization).

A = np.array([[2.,1.],[1.,2.]])
w, V = np.linalg.eig(A)     # w = [3., 1.] (order may vary)
np.allclose(A @ V[:,0], w[0]*V[:,0])   # True

C7 — Complex eigenvalues (pure rotation)

Steps.

1. Build the characteristic determinant $\det (R-\lambda I)=\det \left(\begin{array}{cc}-\lambda & -1\\ 1& -\lambda \end{array}\right)={\lambda }^2+1=0$. Why this step? Same recipe as C6: eigenvalues are the $\lambda$ that make $R-\lambda I$ singular, and singular ⇔ zero determinant. We subtract $\lambda I$ from $R$ and take the $2\times 2$ determinant $(-\lambda )(-\lambda )-(-1)(1)={\lambda }^2+1$.
2. ${\lambda }^2=-1\Rightarrow \lambda =\pm i$ (imaginary!). Why complex? Look at the figure: every real arrow $v$ gets turned 90°, so $Rv$ is never a real multiple of $v$. No real eigenvector can exist — the algebra confesses this by producing imaginary $\lambda$.
3. This is exactly the parent note's warning: eig on a rotation returns complex numbers. NumPy gives w = [0.+1.j, 0.-1.j].

R = np.array([[0.,-1.],[1.,0.]])
np.linalg.eigvals(R)   # array([0.+1.j, 0.-1.j])

C8 — Non-square matrix (eig would crash, SVD saves us)

Steps.

1. eig(A) is undefined — eigenvalues need a square machine that maps a space to itself. $A$ sends 2-D input to 3-D output. Why SVD instead? SVD (parent §5) works for any shape: $A=U\Sigma V^{\top }$, read right-to-left as rotate ($V^{\top }$) → stretch ($\Sigma$) → rotate ($U$).
2. Singular values are the stretches. Since $A$ scales the first axis by 3, the second by 2, and adds a zero row (no third input direction), ${\sigma }_1=3,{\sigma }_2=2$. Why: ${\sigma }_i^2$ are eigenvalues of $A^{\top }A=\left(\begin{array}{cc}9& 0\\ 0& 4\end{array}\right)$, giving $9,4$, so $\sigma =3,2$.
3. How $U$ and $V$ set the ellipse's axes. In the figure, the input unit circle maps to an ellipse. The columns of $V$ are the input directions that get stretched cleanly; the columns of $U$ are the output directions those land on — the ellipse's principal axes. Concretely: $V$'s first column $(1,0)$ is sent to ${\sigma }_1=3$ times $U$'s first column, giving the long (red) semi-axis of length 3; $V$'s second column $(0,1)$ maps to ${\sigma }_2=2$ times $U$'s second column, the short semi-axis of length 2. So $V$ says which input arrows, $\Sigma$ says how far they stretch, $U$ says where they point in output space. Why this matters: this is the exact geometric engine behind Principal Component Analysis (PCA) — the biggest $\sigma$ picks the direction of greatest spread.

A = np.array([[3.,0.],[0.,2.],[0.,0.]])
U, s, Vt = np.linalg.svd(A)
s                     # array([3., 2.])
Sig = np.zeros_like(A); np.fill_diagonal(Sig, s)
np.allclose(U @ Sig @ Vt, A)   # True

C9 — Rank-deficient matrix / compression (word problem)

Steps.

1. Notice column 2 $=2\times$ column 1. So all data lies on one line — effectively 1-D. Why care? This is the heart of Principal Component Analysis (PCA) and compression: redundant directions carry zero singular value.
2. Take the SVD. $A^{\top }A=\left(\begin{array}{cc}14& 28\\ 28& 56\end{array}\right)$, whose determinant is $14\cdot 56-28^2=784-784=0$ → one eigenvalue is $0$.
3. Eigenvalues of $A^{\top }A$: trace $=70$, so they are $70$ and $0$. Singular values ${\sigma }_1=\sqrt{70}\approx 8.367$, ${\sigma }_2=0$. Why: a zero singular value = a direction with no stretch = a collapsed dimension → the matrix has rank 1.
4. Answer: exactly one independent pattern. We could store this data with a single number per day (times the fixed ratio) — that's compression.

A = np.array([[1.,2.],[2.,4.],[3.,6.]])
U, s, Vt = np.linalg.svd(A)
s                          # array([8.3666, 0.]) -> rank 1
np.linalg.matrix_rank(A)   # 1

C10 — Non-square / over-determined system (lstsq, not solve)

Steps.

1. Write one equation per point: $\{\begin{array}{l}m(0)+c=1\\ m(1)+c=2\\ m(2)+c=2\end{array}$, i.e. $A\left(\begin{array}{c}m\\ c\end{array}\right)=b$ with $A=\left(\begin{array}{cc}0& 1\\ 1& 1\\ 2& 1\end{array}\right)$, $b=\left(\begin{array}{c}1\\ 2\\ 2\end{array}\right)$. Why this step? Each data point is a demand on the line; stacking them makes $A$ tall (3 rows, 2 columns) — a non-square, over-determined system.
2. np.linalg.solve refuses this: it requires a square $A$. There is no exact $(m,c)$ hitting all three points (they are not collinear), so no true inverse exists. Why this step? solve answers "undo the machine exactly"; here no exact undo exists, so we must change the question.
3. Change the question to least squares: find $(m,c)$ minimizing the total squared miss $\Vert Ax-b{\Vert }_2^2$ — the tool for this is np.linalg.lstsq (Least squares regression). Why this tool and not solve? When you cannot satisfy every equation, the next-best is to be as close as possible in Euclidean distance; that is precisely what least squares minimizes.
4. The minimizer solves the normal equations $A^{\top }Ax=A^{\top }b$: $A^{\top }A=\left(\begin{array}{cc}5& 3\\ 3& 3\end{array}\right)$, $A^{\top }b=\left(\begin{array}{c}6\\ 5\end{array}\right)$. Solving this square $2\times 2$ system: $\det =5\cdot 3-3\cdot 3=6$, giving $m=\frac{1}{2},c=\frac{7}{6}\approx 1.167$. Why this step? $A^{\top }A$ is small and square — safe for solve — and its solution is provably the closest fit. In the figure, the red line threads between the points, minimizing the vertical gaps.

$\boxed{m=\frac{1}{2},c=\frac{7}{6}}$

A = np.array([[0.,1.],[1.,1.],[2.,1.]]); b = np.array([1.,2.,2.])
# np.linalg.solve(A, b)  -> LinAlgError: Last 2 dims must be square
sol, *_ = np.linalg.lstsq(A, b, rcond=None)   # array([0.5, 1.1666...])

C11 — Exam twist: pick the right tool, then verify

Steps.

1. One number, two answers. Compute $\det A=(4)(3)-(3)(6)=12-18=-6$. Why this step? $\det \ne 0$ ⇒ unique solution for every $b$ (answer a); $|\det |=6$ ⇒ areas scale by 6, and the minus sign ⇒ orientation flips (answer b). One determinant settles both sub-questions.
2. Now solve for $b=(1,0)$ with solve, never inv: $\{\begin{array}{l}4x+3y=1\\ 6x+3y=0\end{array}$. Subtract row 1 from row 2: $(6-4)x=-1\Rightarrow x=-\frac{1}{2}$. Then $6(-\frac{1}{2})+3y=0\Rightarrow y=1$. Why this order? Confirm invertibility (step 1) before asking for a solution — no point solving a collapsed system; and elimination gives the answer directly.

$\boxed{\det A=-6,x=-\frac{1}{2},y=1}$

A = np.array([[4.,3.],[6.,3.]]); b = np.array([1.,0.])
np.linalg.det(A)          # -6.0
np.linalg.solve(A, b)     # array([-0.5, 1.])

Quick self-test