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

This page is a practice range. The parent topic note told you what LU, QR, and Schur are. Here we fire every kind of situation at those tools — nice square systems, singular matrices, tall skinny least-squares, defective eigenvalues, symmetric-positive-definite shortcuts, and a couple of exam traps — and work each one by hand, then check it.

Before we start, four pieces of notation we lean on constantly:

The scenario matrix

Every problem this topic throws at you falls into one of these cells. The worked examples below are labelled with the cell they cover, so by the end you have hit every box.

#	Cell (case class)	Tool	What makes it tricky
C1	Square system, needs a row swap	LU	pivot on largest entry, track $P$
C2	Zero pivot forces a swap	LU	naive $a_{11} = 0$ → divide-by-zero
C3	Singular / degenerate matrix	LU	$U$ has a zero on diagonal → no unique solution
C4	Symmetric positive definite	Cholesky	skip pivoting, half the work
C5	Tall skinny (overdetermined), least squares	QR	avoid normal equations
C6	Sign/quadrant of the fit — negative slope	QR/lstsq	make sure signs come out right
C7	Defective matrix (eig would break)	Schur	eigenvectors don't span
C8	Complex-eigenvalue rotation (real Schur block)	Schur	$2 \times 2$ block on diagonal
C9	Word problem (real world)	LU	translate story → $A x = b$
C10	Exam twist: scipy's $P$ convention	LU	$P A = LU$ vs $A = P LU$ trap

Prereqs we reuse: Gaussian Elimination, Gram-Schmidt Orthogonalization, Eigenvalues and Eigenvectors, Cholesky Decomposition, Least Squares, Condition Number, LAPACK and BLAS, and the plain numpy.linalg for contrast.

C1 · LU with a genuine row swap

Worked example Solve a 3×3 that needs pivoting

$A = 144246224, b = 3610 .$ Solve $A x = b$ by hand using $P A = LU$ .

Forecast (guess before reading on): the first column is $(1, 4, 4)$ . The largest is $4$ , not the $1$ sitting on top — so we must swap. The answer should be small integers or simple fractions.

Worked steps.

Step 1 — pick the pivot in column 1. Entries are $∣1∣, ∣4∣, ∣4∣$ . Take row 2 (the first $4$ ). This swap is the permutation $P_{1}$ (identity with rows 1↔2 swapped). Why this step? Dividing by the biggest entry keeps every multiplier $\leq 1$ , which stops rounding error from exploding. $P_{1} A = 414426224, P_{1} = 010100001 .$

Step 2 — eliminate below the pivot $4$ . Multipliers $m_{21} = 1/4$ , $m_{31} = 4/4 = 1$ (each = entry-to-kill over pivot). Row2 $\leftarrow$ Row2 $- \frac{1}{4}$ Row1 $= (0, 1, 1.5)$ . Row3 $\leftarrow$ Row3 $- 1 \cdot$ Row1 $= (0, 2, 2)$ . Why this step? This is one elimination sweep — exactly Gaussian Elimination, but we store the multipliers instead of throwing them away. $400412 2 1.5 2 .$

Step 3 — pivot in column 2. Candidates below the diagonal: $∣1∣$ (row2) vs $∣2∣$ (row3). Swap them → new pivot $2$ . This second swap is $P_{2}$ (identity with rows 2↔3 swapped). Why this step? The largest-entry rule applies at every column, not just the first, so we re-check before eliminating again. $P_{2} 400412 2 1.5 2 = 400421 22 1.5, P = P_{2} P_{1} .$

Step 4 — eliminate below the second pivot. $m_{32} = 1/2$ . Row3 $\leftarrow$ Row3 $- \frac{1}{2}$ Row2 $= (0, 0, 0.5)$ . Why this step? Same knock-out move one column further right; once this is done every entry below the diagonal is zero, so what remains is upper-triangular $U$ . $U = 400420 22 0.5 .$

Step 5 — apply the total permutation $P = P_{2} P_{1}$ to $b$ . Original $b = (3, 6, 10)$ . $P_{1}$ swaps 1↔2 → $(6, 3, 10)$ ; then $P_{2}$ swaps 2↔3 → $P b = (6, 10, 3)$ . Why this step? $P A = LU$ means we must feed the solver the same re-ordered right-hand side $P b$ , or $L$ and $U$ would be paired with the wrong equations.

Step 6 — forward-solve $L y = P b$ with $L = 11 1/4 01 1/2 001$ (the stored multipliers, reordered by $P$ ): $y_{1} = 6$ ; $1 \cdot 6 + y_{2} = 10 \Rightarrow y_{2} = 4$ ; $\frac{1}{4} (6) + \frac{1}{2} (4) + y_{3} = 3 \Rightarrow y_{3} = 3 - 1.5 - 2 = - 0.5.$ Why this step? $L$ is lower-triangular, so each unknown $y_{i}$ depends only on ones already found — we sweep top-down, one substitution each.

Step 7 — back-solve $U x = y$ : $0.5 x_{3} = - 0.5 \Rightarrow x_{3} = - 1$ ; $2 x_{2} + 2 (- 1) = 4 \Rightarrow x_{2} = 3$ ; $4 x_{1} + 4 (3) + 2 (- 1) = 6 \Rightarrow 4 x_{1} = - 4 \Rightarrow x_{1} = - 1.$ Why this step? $U$ is upper-triangular, so now each unknown depends only on ones below it — we sweep bottom-up, mirroring step 6.

Answer: $x = (- 1, 3, - 1)$ .

Verify: plug into original $A$ . Row1: $- 1 + 6 - 2 = 3$ ✓. Row2: $- 4 + 12 - 2 = 6$ ✓. Row3: $- 4 + 18 - 4 = 10$ ✓.

from scipy.linalg import lu, solve
import numpy as np
A = np.array([[1.,2,2],[4,4,2],[4,6,4]]); b=np.array([3.,6,10])
print(solve(A,b))   # -> [-1.  3. -1.]

C2 · A zero pivot: why the naive formula divides by zero

Worked example The pivot on top is exactly

$0$ $A = (0321), b = (45) .$

Forecast (guess first): $a_{11} = 0$ . The textbook " $m_{21} = a_{21} / a_{11}$ " would be $3/0$ . A swap should save us and the answer should be clean.

Worked steps.

Step 1 — try naive elimination. $m_{21} = 3/0 = \infty$ . Undefined. This is precisely case C2. Why this step? To see the failure the parent warned about, not just be told about it.

Step 2 — swap rows (partial pivoting picks the $3$ ): the permutation $P = (0110)$ gives $P A = (3012)$ , already upper-triangular, and $P b = (5, 4)$ . Why this step? Pivoting isn't only for accuracy — here it is the only way to proceed at all, since the naive multiplier didn't exist.

Step 3 — back-solve $U x = P b$ . From row2: $2 x_{2} = 4 \Rightarrow x_{2} = 2$ . From row1: $3 x_{1} + 2 = 5 \Rightarrow x_{1} = 1.$ Why this step? $P A$ is already upper-triangular, so no elimination is left — we jump straight to bottom-up substitution.

Answer: $x = (1, 2)$ .

Verify: Row1 of original $A$ : $0 (1) + 2 (2) = 4$ ✓. Row2: $3 (1) + 1 (2) = 5$ ✓.

Takeaway: a zero pivot is not a singular matrix — it just needs a swap. Contrast with C3 next.

C3 · A truly singular matrix — no unique solution

Worked example Elimination hits a zero it

cannot swap away $A = (1224), b = (36) .$

Forecast (guess first): row 2 is exactly $2 \times$ row 1 — the rows are parallel. Expect infinitely many solutions (this particular $b$ is consistent), and $U$ will show a zero pivot with nothing to swap in.

Worked steps.

Step 1 — eliminate. $m_{21} = 2/1 = 2$ . Row2 $\leftarrow$ Row2 $- 2$ Row1 $= (0, 0)$ . Why this step? A zero entire row in $U$ means $det A = 1 \cdot 0 = 0$ → singular; this is what we wanted to expose. $U = (1020) .$

Step 2 — can pivoting rescue it? Column 2 has no nonzero entry below the diagonal to swap in. No. This is the difference from C2. Why this step? To hammer the distinction: C2's zero was cosmetic (a swap fixed it), C3's zero is structural (no swap exists).

Step 3 — consistency check on $b$ . Apply the same row-op to $b$ : $6 - 2 (3) = 0$ . The reduced system is $x_{1} + 2 x_{2} = 3$ , $0 = 0$ . Infinitely many solutions: $x = (3 - 2 t, t)$ for any $t$ . Why this step? A singular matrix can still be consistent; we must test the right-hand side to learn whether we get infinitely many solutions ( $0 = 0$ ) or none ( $0 = nonzero$ ).

Answer: singular — no unique $x$ . solve raises / warns; lstsq returns the minimum-norm solution.

Verify: determinant $= 1 \cdot 4 - 2 \cdot 2 = 0$ ✓. Pick $t = 0$ : $x = (3, 0)$ , then $A x = (3, 6) = b$ ✓; pick $t = 1$ : $x = (1, 1)$ , $A x = (3, 6) = b$ ✓ — two different valid $x$ , confirming non-uniqueness.

from scipy.linalg import lstsq
import numpy as np
A=np.array([[1.,2],[2,4]]); b=np.array([3.,6])
x,res,rank,sv = lstsq(A,b)
print(rank)          # 1  (< 2  => rank-deficient => not unique)

C4 · Symmetric positive definite → use Cholesky, skip pivoting

Worked example SPD system, the Cholesky shortcut

$A = (4223), b = (25) .$ $A$ is symmetric; both leading minors ( $4 > 0$ and $det = 8 > 0$ ) are positive → SPD.

Forecast (guess first): Cholesky Decomposition gives $A = L L^{T}$ with no pivoting needed (SPD guarantees positive pivots). Expect half the flops of LU, and a clean answer.

Worked steps.

Step 1 — factor $A = L L^{T}$ , $L$ lower-triangular. $ℓ_{11} = 4 = 2$ . $ℓ_{21} = 2/ ℓ_{11} = 1$ . $ℓ_{22} = 3 - ℓ_{21}^{2} = 2$ . Why this step? SPD ⇒ the quantity under each square root stays positive, so we never divide by tiny/zero pivots and never need $P$ . $L = (2102) .$

Step 2 — forward-solve $L y = b$ : $2 y_{1} = 2 \Rightarrow y_{1} = 1$ ; $1 (1) + 2 y_{2} = 5 \Rightarrow y_{2} = 4/ 2 = 22 .$ Why this step? Cholesky splits $A x = b$ into two triangular solves just like LU; $L$ lower-triangular ⇒ sweep top-down.

Step 3 — back-solve $L^{T} x = y$ : $2 x_{2} = 22 \Rightarrow x_{2} = 2$ ; $2 x_{1} + 1 (2) = 1 \Rightarrow x_{1} = - \frac{1}{2} .$ Why this step? The second factor $L^{T}$ is upper-triangular, so we finish with a bottom-up sweep — no separate $U$ was ever needed, which is where the flop savings come from.

Answer: $x = (- 0.5, 2)$ .

Verify: Row1: $4 (- 0.5) + 2 (2) = 2$ ✓. Row2: $2 (- 0.5) + 3 (2) = 5$ ✓.

from scipy.linalg import cho_factor, cho_solve
import numpy as np
A=np.array([[4.,2],[2,3]]); b=np.array([2.,5])
c=cho_factor(A); print(cho_solve(c,b))   # [-0.5  2. ]

C5 · Tall skinny least squares via QR

This one is geometric — look at the picture as you read. In the figure, the magenta dots are the data, the violet line is the fit we are about to derive, and the orange dashed arrows are the residuals whose total length QR minimises.

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

Worked example Fit a line

$y = a + b x$ through 3 points Data: $(x, y) = (0, 1), (1, 2), (2, 2)$ . Three equations, two unknowns $(a, b)$ → overdetermined. We collect the right-hand values into the vector $r_{rhs} = (1, 2, 2)$ (this plays the role of $b$ ; we rename it so it doesn't clash with the slope $b$ ): $A = 111012, r_{rhs} = 122 .$

Forecast (guess first): the three points don't lie on one line, so there's no exact solution — we want the line closest in vertical distance (the orange arrows in the figure). Expect slope near $+ 0.5$ , intercept near $1.17$ .

Worked steps.

Step 1 — why QR, not the normal equations? $A^{T} A$ squares the Condition Number ( $κ (A^{T} A) = κ (A)^{2}$ ). QR keeps $κ (A)$ . See the parent's Least Squares argument. Why this step? Choosing the tool up front is the whole point of the scenario matrix — QR trades a tiny bit of work for a big stability win.

Step 2 — orthonormalize the columns (Gram-Schmidt Orthogonalization). Column 1 is $(1, 1, 1)$ , length $3$ → $q_{1} = \frac{1}{3} (1, 1, 1)$ . Column 2 is $(0, 1, 2)$ ; subtract its $q_{1}$ -part: $q_{1}^{T} (0, 1, 2) = \frac{3}{3} = 3$ , projection $= 3 q_{1} = (1, 1, 1)$ . Residual $(0, 1, 2) - (1, 1, 1) = (- 1, 0, 1)$ , length $2$ → $q_{2} = \frac{1}{2} (- 1, 0, 1)$ . Why this step? Orthonormal $Q$ preserves length, which is exactly what turns the messy fit into a clean triangular solve in step 4.

Step 3 — form $R = Q^{T} A$ : $R = (3032)$ . Why this step? $R$ records how each original column of $A$ is rebuilt from the orthonormal basis $Q$ ; being upper-triangular is what makes back-substitution possible.

Step 4 — compute $Q^{T} r_{rhs}$ : $q_{1}^{T} r_{rhs} = \frac{1 + 2 + 2}{3} = \frac{5}{3}$ ; $q_{2}^{T} r_{rhs} = \frac{- 1 + 0 + 2}{2} = \frac{1}{2}$ . Why this step? Because $Q$ preserves length, $min ∥ A x - r_{rhs} ∥$ becomes $min ∥ R x - Q^{T} r_{rhs} ∥$ — so we need $Q^{T} r_{rhs}$ as the new right-hand side.

Step 5 — back-solve $R x = Q^{T} r_{rhs}$ : $2 b = \frac{1}{2} \Rightarrow b = \frac{1}{2}$ ; $3 a + 3 (\frac{1}{2}) = \frac{5}{3} \Rightarrow a = \frac{5}{3} - \frac{1}{2} = \frac{7}{6} .$ Why this step? $R$ is upper-triangular, so a single bottom-up sweep hands us the slope then the intercept — no normal equations, no squared condition number.

Answer: slope $b = 0.5$ , intercept $a = \frac{7}{6} \approx 1.1667$ (matching the violet line in the figure).

Verify: the fitted residual vector should be orthogonal to both columns of $A$ (that's the defining property of least squares). Fitted $\overset{y}{^} = (7/6, 5/3, 13/6)$ , residual $r = (- 1/6, 1/3, - 1/6)$ ; $1^{T} r = 0$ ✓ and $x^{T} r = 0 (- 1/6) + 1 (1/3) + 2 (- 1/6) = 0$ ✓.

C6 · Sign check — a negative slope

Worked example Same machinery, downward line

Data: $(0, 3), (1, 1), (2, - 1)$ . These lie exactly on $y = 3 - 2 x$ . The right-hand vector is $r_{rhs} = (3, 1, - 1)$ .

Forecast (guess first): since the points are perfectly collinear, least squares should return the exact line with zero residual and a clearly negative slope $b = - 2$ .

Worked steps.

Step 1 — reuse $Q, R$ from C5. Same $A = 111012$ , only the right-hand side changed to $r_{rhs} = (3, 1, - 1)$ . Why this step? The basis $Q$ and the triangular $R$ depend only on $A$ , so changing the data $r_{rhs}$ costs us nothing but two dot products.

Step 2 — compute $Q^{T} r_{rhs}$ : $q_{1}^{T} r_{rhs} = \frac{3 + 1 - 1}{3} = \frac{3}{3} = 3$ ; $q_{2}^{T} r_{rhs} = \frac{- 3 + 0 - 1}{2} = \frac{- 4}{2} = - 22 .$ Why this step? Same reasoning as C5 step 4 — project the new data onto the orthonormal basis to get the transformed right-hand side. The negative second component already signals a downward slope.

Step 3 — back-solve $R x = Q^{T} r_{rhs}$ : $2 b = - 22 \Rightarrow b = - 2$ ; $3 a + 3 (- 2) = 3 \Rightarrow a = 1 + 2 = 3.$ Why this step? Bottom-up triangular solve as before; the sign of $b$ falls straight out of the negative right-hand component, confirming the fit really does slope down.

Answer: $a = 3, b = - 2$ — the exact line, negative slope confirmed.

Verify: residual $= A x - r_{rhs}$ : $(3, 1, - 1) - (3, 1, - 1) = (0, 0, 0)$ , norm $0$ ✓.

C7 · Defective matrix — why `eig` breaks and Schur doesn't

Worked example A matrix whose eigenvectors don't span

$A = (2012) .$ Both eigenvalues equal $2$ , but there is only one independent eigenvector — a defective matrix.

Forecast (guess first): you cannot write $A = V Λ V^{- 1}$ (no invertible $V$ ). Diagonalization via Eigenvalues and Eigenvectors should fail, but Schur $A = Z T Z^{H}$ still exists.

Worked steps.

Step 1 — eigenvalues. $det (A - λ I) = (2 - λ)^{2} = 0 \Rightarrow λ = 2$ (a double root). Why this step? The eigenvalues go on the diagonal of both $Λ$ and the Schur $T$ , so we need them first.

Step 2 — eigenvectors. $(A - 2 I) v = (0010) v = 0 \Rightarrow v = (1, 0)$ only. Just one direction → the eigenvector matrix $V$ has a repeated/zero column → $V$ is singular → $A = V Λ V^{- 1}$ is impossible (we can't invert $V$ ). Why this step? This is exactly the case the parent said kills diagonalization: $Λ$ would be fine, but $V^{- 1}$ doesn't exist.

Step 3 — Schur decomposition. $A$ is already upper-triangular, so take $Z = I$ (unitary: $Z^{H} Z = I$ ), $T = A$ . The eigenvalues $2, 2$ sit on $T$ 's diagonal. Schur exists trivially. Why this step? Schur only needs a unitary $Z$ (condition number $1$ ), never an invertible eigenvector matrix — so it survives defectiveness where diagonalization dies.

Answer: eigenvalues ${2, 2}$ ; diagonalizable? No. Schur form? Yes, $Z = I, T = A$ .

Verify: $det (A - 2 I) = 0$ ✓ and $rank (A - 2 I) = 1$ (so only $2 - 1 = 1$ eigenvector) ✓.

C8 · Real Schur with a complex-eigenvalue $2 \times 2$ block

This is a rotation — see the figure: solid arrows are input vectors $v$ , dashed arrows are their images $A v$ , and no solid/dashed pair is parallel.

Worked example A pure rotation matrix

$A = (01 - 1 0) (rotate by 9 0^{\circ}) .$

Forecast (guess first): rotations spin every real vector, so no real vector maps to a multiple of itself — eigenvalues must be complex ( $\pm i$ ). Expect real Schur to keep $A$ as one $2 \times 2$ block.

Worked steps.

Step 1 — eigenvalues. $det (A - λ I) = λ^{2} + 1 = 0 \Rightarrow λ = \pm i .$ Why this step? Complex eigenvalues are the signature of rotation — and, as the figure shows, no solid input lines up with its dashed image, so no real eigenvector can exist.

Step 2 — real Schur form. Because the eigenvalues are a complex-conjugate pair, real Schur cannot split them into two real $1 \times 1$ diagonal entries; instead it leaves a $2 \times 2$ block on the diagonal — this is called quasi-triangular (block upper-triangular). Here $A$ itself is already that single block, so $Z = I$ , $T = A$ . Why this step? Staying in real arithmetic keeps everything perfectly conditioned ( $κ (Z) = 1$ ) — no complex round-off. The block is still upper-triangular in the block sense (nothing below it), just not entry-by-entry triangular over the reals.

Answer: eigenvalues $\pm i$ ; real (quasi-triangular) Schur = one $2 \times 2$ block, essentially $A$ .

Verify: $det (A - λ I) = λ^{2} + 1$ ; roots $\pm i$ satisfy $i^{2} + 1 = 0$ ✓. Also $A^{T} A = I$ (orthogonal/length-preserving) ✓.

C9 · Word problem — mixing two solutions

Worked example Real-world: a chemistry blend

A lab needs $7$ litres of acid and $4$ litres of base. Solution X is $50%$ acid / $50%$ base. Solution Y is $80%$ acid / $20%$ base. How many litres $x$ of X and $y$ of Y?

Forecast (guess first): two ingredients, two requirements → a clean $2 \times 2$ system. Expect a few litres of each.

Worked steps.

Step 1 — translate the story to $A v = c$ . Acid: $0.5 x + 0.8 y = 7$ . Base: $0.5 x + 0.2 y = 4$ . Why this step? The parent's whole point: real questions become $A x = b$ , then LU (Gaussian elimination) solves them mechanically. $A = (0.5 0.5 0.8 0.2), c = (74), v = (x y) .$

Step 2 — eliminate. Pivot is $0.5$ (top-left, nonzero — no swap needed), multiplier $m_{21} = 0.5/0.5 = 1$ : Row2 $\leftarrow$ Row2 $-$ Row1 $\Rightarrow (0, - 0.6 ∣ - 3)$ . Why this step? One elimination sweep clears the entry below the pivot, leaving a triangular system ready for back-substitution.

Step 3 — back-solve. $- 0.6 y = - 3 \Rightarrow y = 5$ ; $0.5 x + 0.8 (5) = 7 \Rightarrow 0.5 x = 3 \Rightarrow x = 6.$ Why this step? The system is now upper-triangular, so a bottom-up sweep reads off $y$ then $x$ .

Answer: $x = 6$ L of X, $y = 5$ L of Y.

Verify (units + plug-in): total acid $= 0.5 (6) + 0.8 (5) = 3 + 4 = 7$ L ✓; total base $= 0.5 (6) + 0.2 (5) = 3 + 1 = 4$ L ✓. Total volume $6 + 5 = 11$ L, sensible.

C10 · Exam twist — scipy's permutation convention

Worked example The

$P A = LU$ vs $A = P LU$ trap You call P, L, U = scipy.linalg.lu(A) on $A = (2615)$ (the parent's example). Which reconstruction is correct: A = P @ L @ U or A = P.T @ L @ U?

Forecast (guess first): the parent flagged this. scipy returns $P$ such that $A = P @ L @ U$ (its docstring differs from the " $P A = LU$ " textbook), so the other form should fail — but only when $P \neq = P^{T}$ .

Worked steps.

Step 1 — reproduce the factorization. Largest of $∣2∣, ∣6∣$ is $6$ , so rows swap. $m_{21} = 2/6 = 1/3$ , $U = (60 5 - 2/3)$ , $L = (1 1/3 01)$ . Why this step? You must know the factorization to know what $P$ has to undo.

Step 2 — read scipy's convention carefully. scipy's returned $P$ satisfies $A = P L U$ (it already pre-multiplies $LU$ back into the original row order). So the textbook $P A = LU$ corresponds to scipy's $P^{T}$ . Careful: here the permutation is a single row-swap, so $P = P^{T}$ (a swap is its own inverse) — this $2 \times 2$ example cannot distinguish the two forms. Why this step? Mixing conventions silently gives a wrong but same-shaped matrix; a symmetric swap hides the bug, so we need a case where $P \neq = P^{T}$ to expose it.

Step 3 — build a distinguishing example. Take $A = 002100010$ , whose pivoting produces a cyclic permutation $P$ (a three-way rotation of rows, $P \neq = P^{T}$ ). Now $A = P @ L @ U$ holds but $A = P^{T} @ L @ U$ does not. Why this step? Only a non-symmetric permutation forces the two conventions to disagree, turning the silent trap into a visible failure.

Answer: with scipy, use $A = P @ L @ U$ (equivalently $P^{T} @ A = L @ U$ ).

Verify: for the $2 \times 2$ swap $P = P^{T}$ , so both forms coincide (see VERIFY). For the $3 \times 3$ cyclic matrix, $∥ A - P LU ∥ = 0$ while $∥ A - P^{T} LU ∥ \neq = 0$ — that is the case that exposes the trap.

from scipy.linalg import lu
import numpy as np
# 2x2: single swap, P == P.T, both forms agree (does NOT expose the trap)
A2=np.array([[2.,1],[6,5]]); P,L,U=lu(A2)
print(np.allclose(A2, P@L@U), np.allclose(A2, P.T@L@U))   # True True
 
# 3x3 cyclic: P != P.T, only the scipy convention works
A3=np.array([[0.,1,0],[0,0,1],[2,0,0]]); P,L,U=lu(A3)
print(np.allclose(A3, P@L@U))     # True  <- scipy convention
print(np.allclose(A3, P.T@L@U))   # False <- textbook PA=LU form differs here

Recall Scenario checklist — can you name the tool for each?

Row swap needed ::: LU with partial pivoting (track $P$ ). Pivot is 0 but matrix is fine ::: swap rows, still LU. Whole row of $U$ is 0 ::: singular → no unique solution (use lstsq for min-norm). Symmetric positive definite ::: Cholesky $A = L L^{T}$ , no pivoting. More equations than unknowns ::: QR least squares (avoid normal equations). Eigenvectors don't span (defective) ::: Schur, not diagonalization. Complex eigenvalues in a real matrix ::: real (quasi-triangular) Schur $2 \times 2$ block. scipy lu reconstruction ::: $A = P @ L @ U$ , not $P . T @ L @ U$ .