scipy.linalg — more stable than numpy.linalg, lu, qr, schur: Exercises

This page is a graded workout for the parent topic. Work each problem before opening its solution. Levels climb from recognition to mastery. Every symbol is built from the ground up as we go — nothing is assumed.

Before we start, some words we will use constantly. Picture a matrix as a machine that moves arrows (vectors).

Level 1 — Recognition

Goal: can you spot which factorization/tool a situation calls for, and read scipy's output correctly?

L1.1 — Name the tool

For each task, name the single best scipy.linalg routine.

(a) Solve one system $A x = b$ for a general square $A$ . (b) Find $min_{x} ∥ A x - b ∥$ where $A$ is tall ( $m > n$ , more equations than unknowns). (c) Get eigenvalues of a possibly defective matrix (recall: one with too few independent eigenvectors to diagonalize), using only perfectly-conditioned steps. (d) Solve the generalized eigenproblem $A x = λ B x$ .

Recall Solution

(a) solve(A, b) — internally one LU factorization + two triangular solves. (b) lstsq(A, b) (QR/SVD under the hood). "Tall + minimize the leftover length" = least squares. (c) schur(A) — always exists, uses only unitary steps ( $κ = 1$ ), eigenvalues appear on the diagonal of $T$ . Diagonalization would fail here because a defective matrix has no invertible eigenvector matrix $V$ . (d) eig(A, b=B) — scipy's eig accepts a second matrix $B$ ; numpy.linalg cannot do this.

L1.2 — Read the `lu` output

scipy.linalg.lu(A) returns P, L, U. Which identity is correct?

$(i) A = P L U (ii) P A = L U (iii) A = P^{T} L U$

Recall Solution

Both (ii) and (iii) are correct — they are the same statement. SciPy's convention is $P A = LU$ . Multiplying on the left by $P^{T}$ and using $P^{T} P = I$ (a permutation matrix is orthogonal) gives $A = P^{T} LU$ . Statement (i) $A = P LU$ is the common wrong reading.

Level 2 — Application

Goal: run a factorization by hand, then confirm it is exactly what the code returns.

L2.1 — LU by hand, then solve

Let

$A = (4633), b = (1012) .$

Do partial-pivoted LU by hand, then solve $A x = b$ via forward + back substitution.

The figure below walks the two stages of this exact problem. Left panel: a bar chart comparing the first-column magnitudes $∣4∣$ (row 1) and $∣6∣$ (row 2); the taller amber bar (row 2, value $6$ ) wins and becomes the pivot, so we swap rows. Right panel: the substitution ladder — the cyan block does the forward pass $L y = P b$ from the top down ( $y_{1} = 12$ , then $y_{2} = 2$ ), and the amber block does the back pass $U x = y$ from the bottom up ( $x_{2} = 2$ , then $x_{1} = 1$ ), landing on the solution $x = (1, 2)$ .

Figure — scipy.linalg — more stable than numpy.linalg, lu, qr, schur

Recall Solution

Pivot. Compare the first-column entries $∣4∣$ and $∣6∣$ . We pivot on the largest to avoid dividing by something small (small pivots amplify rounding). Since $∣6∣ > ∣4∣$ , swap rows: $P A = (6433), P = (0110) .$ Eliminate. The multiplier is == $m_{21} = \frac{4}{6} = \frac{2}{3}$ == (the "how many pivot-rows fit into the entry below" number). Do $R_{2} \leftarrow R_{2} - \frac{2}{3} R_{1}$ :

L=\begin{pmatrix}1&0\\ \tfrac23 & 1\end{pmatrix}.$$ $L$ just **stores the multiplier** — undoing "subtract $\tfrac23$" is "add $\tfrac23$", so no new computation. **Permute $b$.** We swapped rows, so $Pb=\begin{pmatrix}12\\10\end{pmatrix}$. **Forward-solve $Ly=Pb$** (top-down through the lower triangle): $$y_1=12,\qquad \tfrac23 y_1 + y_2 = 10 \Rightarrow y_2 = 10 - \tfrac23(12)=2.$$ **Back-solve $Ux=y$** (bottom-up through the upper triangle): $$1\cdot x_2 = 2 \Rightarrow x_2 = 2,\qquad 6x_1 + 3x_2 = 12 \Rightarrow 6x_1 = 6 \Rightarrow x_1 = 1.$$ **Answer $x=(1,\,2)$.** Check in the original $A$: row 1 $= 4(1)+3(2)=10$ ✓, row 2 $=6(1)+3(2)=12$ ✓.

L2.2 — QR of a $2 \times 2$ by hand

Let

$A = (3411) .$

Use Gram–Schmidt on the columns to build $Q$ (orthonormal columns) and $R$ (upper triangular) with $A = QR$ .

Recall Solution

Call the columns $a_{1} = (4 3)$ , $a_{2} = (1 1)$ . Gram–Schmidt turns these two arrows into two perpendicular unit arrows. Step 1 — normalise $a_{1}$ . Its length (norm) is $∥ a_{1} ∥ = 3^{2} + 4^{2} = 5$ . So $q_{1} = \frac{1}{5} (4 3) = (0.8 0.6) .$ The first $R$ entry is $r_{11} = ∥ a_{1} ∥ = 5$ (how much of $q_{1}$ builds $a_{1}$ ). Step 2 — remove the $q_{1}$ -part of $a_{2}$ . The overlap is the dot product $r_{12} = q_{1} \cdot a_{2} = 0.6 (1) + 0.8 (1) = 1.4$ . The leftover perpendicular piece: $a_{2} - r_{12} q_{1} = (1 1) - 1.4 (0.8 0.6) = (1 - 1.12 1 - 0.84) = (- 0.12 0.16) .$ Its length is $0.1 6^{2} + 0.1 2^{2} = 0.04 = 0.2 = r_{22}$ . Normalise: $q_{2} = \frac{1}{0.2} (- 0.12 0.16) = (- 0.6 0.8) .$ Assemble. $Q = (0.6 0.8 0.8 - 0.6), R = (50 1.4 0.2) .$ Check $Q^{T} Q = I$ (columns perpendicular unit) and $QR = A$ . (SciPy may flip a sign on a whole $q_{i}$ column and the matching $R$ row — both are valid QRs since $q$ and $- q$ are both unit and perpendicular.)

Level 3 — Analysis

Goal: reason about why one method is more accurate than another.

L3.1 — Why not the normal equations?

To solve least squares $min_{x} ∥ A x - b ∥$ , a classic route is the normal equations $A^{T} A x = A^{T} b$ . QR instead reduces the problem to $R x = Q^{T} b$ . Using the condition number $κ (A)$ (a number $\geq 1$ measuring how much $A$ can amplify input error — see Condition Number), explain in one sentence why QR is safer, and give the numeric blow-up factor when $κ (A) = 1 0^{6}$ .

Recall Solution

Forming $A^{T} A$ squares the condition number: $κ (A^{T} A) = κ (A)^{2}$ . QR never forms $A^{T} A$ ; it works with $A$ directly, so it keeps $κ (A)$ . With $κ (A) = 1 0^{6}$ , the normal equations behave like a system of condition $κ (A)^{2} = 1 0^{12}$ — you lose about $12$ decimal digits instead of $6$ . On standard double precision ( $\approx 16$ digits) the normal-equations answer can be garbage, while QR still keeps $\approx 10$ good digits.

L3.2 — Pivot choice changes the error

Solve $(1 0^{- 16} 1 11) (x _{2} x _{1}) = (2 1)$ two ways: (a) no pivoting (use the tiny $1 0^{- 16}$ as pivot), (b) partial pivoting (swap so $1$ is the pivot). Which one is safe, and why?

Recall Solution

True answer: $x_{1} = 1/ (1 - 1 0^{- 16}) \approx 1.0000 \dots$ , $x_{2} = 1 - 1 0^{- 16} x_{1} \approx 1$ . So $x \approx (1, 1)$ . (a) No pivot. Multiplier $m = 1/1 0^{- 16} = 1 0^{16}$ . Row 2 becomes $1 - 1 0^{16} \cdot 1 = - 1 0^{16}$ (the original $1$ is swallowed by rounding — it is $16$ orders smaller). Back-substitution then divides by this and you recover $x_{1}$ with total loss of the small digits — a catastrophic error. (b) Partial pivoting. Swap first: pivot is $1$ , multiplier $m = 1 0^{- 16}$ , tiny and harmless. Elimination gives a well-behaved $U$ and the computed answer matches $(1, 1)$ to full precision. Moral: partial pivoting picks the largest-magnitude available pivot precisely so the multiplier $m \leq 1$ and nothing gets swamped. This is exactly why scipy.linalg.solve/lu pivot by default — the same lesson as Gaussian Elimination.

Level 4 — Synthesis

Goal: combine factorizations, or connect them across problems.

L4.1 — Reuse one factorization for many right-hand sides

You must solve $A x = b_{k}$ for the same $A = (4633)$ but three different vectors $b_{1} = (12 10)$ , $b_{2} = (9 7)$ , $b_{3} = (3 1)$ . Why is lu_factor + lu_solve better than calling solve three times, and what are the three answers?

Recall Solution

Why factor once: the expensive part of solving is the LU factorization itself (about $\frac{2}{3} n^{3}$ operations). Each right-hand side then costs only two triangular solves ( $\approx n^{2}$ ). Calling solve(A, b_k) three times redoes the factorization every time. lu_factor(A) once, then lu_solve per $b_{k}$ , pays the cubic cost once. Answers (reuse the $L, U, P$ from L2.1, where we already found $A^{- 1}$ 's action):

$b_{1} = (10, 12)$ → $x = (1, 2)$ (done in L2.1).
$b_{2} = (7, 9)$ : $P b_{2} = (9, 7)$ ; $y_{1} = 9, y_{2} = 7 - \frac{2}{3} (9) = 1$ ; back: $x_{2} = 1, 6 x_{1} + 3 (1) = 9 \Rightarrow x_{1} = 1$ . So $x = (1, 1)$ .
$b_{3} = (1, 3)$ : $P b_{3} = (3, 1)$ ; $y_{1} = 3, y_{2} = 1 - \frac{2}{3} (3) = - 1$ ; back: $x_{2} = - 1, 6 x_{1} + 3 (- 1) = 3 \Rightarrow x_{1} = 1$ . So $x = (1, - 1)$ .

from scipy.linalg import lu_factor, lu_solve
lup = lu_factor(A)                 # factor ONCE
xs  = [lu_solve(lup, bk) for bk in (b1, b2, b3)]

L4.2 — Schur to read eigenvalues off the diagonal

For $A = (2013)$ , the matrix is already upper triangular. Explain why its Schur form is trivial, read off the eigenvalues, and confirm they match eig. Why is Schur a safer engine than diagonalization here?

Recall Solution

Already-triangular case. The Schur decomposition is $A = Z T Z^{H}$ with $Z$ unitary and $T$ upper triangular. If $A$ is already upper triangular, one valid choice is $Z = I$ (the identity is unitary) and $T = A$ . The eigenvalues sit on the diagonal of $T$ , so they are the diagonal of $A$ : $λ_{1} = 2, λ_{2} = 3$ . Check: the characteristic polynomial is $det (A - λ I) = (2 - λ) (3 - λ) = 0$ , giving $λ = 2, 3$ . eig(A) returns exactly [2., 3.]. See Eigenvalues and Eigenvectors. Why Schur is safer: diagonalization $A = V Λ V^{- 1}$ needs $V$ invertible and can have a huge $κ (V)$ (blows up for near-defective matrices). Schur uses only the unitary $Z$ with $κ (Z) = 1$ — no ill-conditioned inverse ever appears, and it always exists. That is why expm/sqrtm are built on it.

Level 5 — Mastery

Goal: pick and justify the whole pipeline, including degenerate cases.

L5.1 — Choose the right tool under constraints

For each scenario, name the routine AND one sentence of justification.

(a) $A$ is symmetric with all eigenvalues positive (SPD), size $5000$ , and you must solve it fast. (b) $A$ is $2000 \times 40$ , and you want $min_{x} ∥ A x - b ∥$ plus the numerical rank of $A$ . (c) You need $e^{A}$ (matrix exponential) for a non-symmetric $A$ . (d) $A$ is singular (columns dependent) and you call solve(A,b) — what happens, what should you do?

Recall Solution

(a) cho_factor / cho_solve — SPD matrices admit Cholesky $A = L L^{T}$ , which needs no pivoting and costs half of LU. Fastest and stable. (Cholesky Decomposition.) (b) lstsq(A, b) — it returns x, residual, rank, singular_values; the SVD-based driver gives a reliable numerical rank, which plain QR does not report cleanly. (Least Squares.) (c) expm(A) — built on the Schur decomposition (unitary, well-conditioned), so it stays accurate even when $A$ is not diagonalizable. (d) solve(A, b) on an exactly singular $A$ raises LinAlgError ("singular matrix"); on a merely near-singular $A$ it returns a wildly amplified, wrong answer with no warning. Fix: switch to lstsq(A, b), which returns the minimum-norm least-squares solution and reports the numerical rank, so it does not choke on rank deficiency and tells you the matrix is degenerate.

L5.2 — Degenerate LU: a zero pivot

Try LU on $A = (0321)$ without pivoting, then with partial pivoting. Show why the naive path breaks and the pivoted path succeeds.

Recall Solution

No pivoting. The would-be first pivot is $a_{11} = 0$ . The multiplier $m_{21} = a_{21} / a_{11} = 3/0$ is undefined — division by zero. Naive LU simply cannot start. This is the degenerate case that pivoting exists to handle. Partial pivoting. Compare column-1 entries $∣0∣$ and $∣3∣$ ; the largest is $3$ , so swap rows: $P A = (3012), P = (0110) .$ Now this is already upper triangular: $m_{21} = 0/3 = 0$ , so $L = (1001) = I$ and $U = (3012)$ . Factorization succeeds. This is exactly what scipy.linalg.lu does automatically — the permutation $P$ guarantees a usable pivot whenever the matrix is nonsingular. If a column were entirely zero (truly singular), no swap could help — and that is precisely when you switch to lstsq.