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

WHY scipy.linalg over numpy.linalg

WHY "more stable"? It's not magic — both call LAPACK. The practical differences:

scipy.linalg is always compiled against an optimized, complete LAPACK/BLAS; numpy may fall back to slower/less-complete routines depending on build.
scipy lets you pick the stable algorithm explicitly (e.g. solve via LU with partial pivoting, or via Cholesky if SPD) instead of a one-size-fits-all path.
Never compute inv(A) @ b — that's the real instability. Use solve/factorizations.

The three core factorizations

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

1. LU — for solving $Ax=b$

Formula Deriving LU from elimination (first principles)

Take $A$ . To zero out $a_{21}$ using row 1 we do $R_2 \leftarrow R_2 - m_{21}R_1$ where $m_{21}=a_{21}/a_{11}$ . This row-operation is left-multiplication by an elementary matrix $E_{21}$ . Doing all eliminations: $E_{n}\cdots E_2 E_1 \, A = U \quad\Rightarrow\quad A = (E_1^{-1}\cdots E_n^{-1})\,U = L U.$ The inverse of "subtract $m$ " is "add $m$ ", so $L$ just stores the multipliers $m_{ij}$ — no extra computation! With pivoting we swap rows first, recorded in $P$ , giving $==PA = LU==$ .

HOW we then solve $Ax=b$ : $PA x = Pb \Rightarrow LUx = Pb$ . Let $y=Ux$ .

Forward-solve $Ly = Pb$ (cheap, triangular).
Back-solve $Ux = y$ (cheap, triangular).

Worked example Solve a 2×2 with LU

$A=\begin{pmatrix}2&1\\6&5\end{pmatrix}$ , $b=\begin{pmatrix}3\\17\end{pmatrix}$ .

Pivot: $|6|>|2|$ → swap rows. $P A=\begin{pmatrix}6&5\\2&1\end{pmatrix}$ . (Why? Pivot on the largest entry to avoid dividing by something tiny.)
$m_{21}=2/6=1/3$ . $U=\begin{pmatrix}6&5\\0& 1-\tfrac13\cdot5\end{pmatrix}=\begin{pmatrix}6&5\\0&-\tfrac23\end{pmatrix}$ . (Why? eliminate the $2$ below the pivot: row2 ← row2 − ⅓·row1.)
$L=\begin{pmatrix}1&0\\1/3&1\end{pmatrix}$ , and since we swapped rows, $Pb=\begin{pmatrix}17\\3\end{pmatrix}$ .
Forward $Ly=Pb$ : $y_1=17$ ; then $\tfrac13 y_1 + y_2 = 3 \Rightarrow y_2 = 3-\tfrac13(17)=-\tfrac{8}{3}$ . (Why? substitute top-down through the triangular L.)
Back $Ux=y$ : $-\tfrac23 x_2=-\tfrac{8}{3}\Rightarrow x_2=4$ ; then $6x_1+5(4)=17\Rightarrow 6x_1=-3\Rightarrow x_1=-\tfrac12$ . (Why? substitute bottom-up through the triangular U.)
Solution $x=(-0.5,\,4)$ . Check against original $A$ : row1 $2(-0.5)+4=3$ ✓, row2 $6(-0.5)+5(4)=17$ ✓.

(Steel-manned takeaway: do arithmetic carefully and always machine-check — see the VERIFY block. In real code just call solve.)

import numpy as np
from scipy.linalg import lu, lu_factor, lu_solve, solve
A = np.array([[2.,1.],[6.,5.]]); b = np.array([3.,17.])
P, L, U = lu(A)                 # SciPy convention: P @ A = L @ U  (so A = P.T @ L @ U)
x = solve(A, b)                 # the right way to solve  -> x = [-0.5, 4.]
lu_piv = lu_factor(A)           # reuse factorization for many b's
x2 = lu_solve(lu_piv, b)

2. QR — for least squares & orthogonalization

from scipy.linalg import qr, lstsq
Q, R = qr(A, mode='economic')   # economic: Q is m×n, R is n×n
x, res, rank, sv = lstsq(A, b)  # robust least squares (uses SVD/QR internally)

3. Schur — for eigenvalues & matrix functions

from scipy.linalg import schur, eig
T, Z = schur(A)                 # A = Z @ T @ Z.conj().T
w, vl, vr = eig(A, left=True)   # scipy gives left+right eigenvectors
# generalized eigenproblem A x = λ B x  (numpy can't do this!):
w_gen = eig(A, b=np.eye(len(A)))

Recall Feynman: explain to a 12-year-old

Imagine a tangled machine that twists arrows around. It's hard to understand the tangle directly. So we break the machine into simple pieces: one piece only stretches (triangular), one piece only spins without stretching (orthogonal/unitary). Once it's in simple pieces, answering questions like "which arrow comes out as itself?" (eigenvalues) or "undo the machine to get back the original arrow" (solve $Ax=b$ ) becomes easy. scipy.linalg is the box of tools that does these clean breaks carefully so tiny rounding mistakes don't blow up.

Why use scipy.linalg over numpy.linalg?

Full LAPACK interface (lu, qr, schur, cho_factor, generalized eig), always linked to optimized LAPACK/BLAS, and lets you pick the stable algorithm explicitly.

What is the LU factorization with pivoting?

PA=LU

— P permutation (row swaps for stability), L unit-lower-triangular (elimination multipliers), U upper triangular.

How do you solve Ax=b once you have PA=LU?

Forward-solve

Ly=Pb

, then back-solve

Ux=y

Why never compute inv(A)@b?

~3× more flops and amplifies rounding error; use solve()/factorizations instead.

What convention does scipy.linalg.lu return?

P, L, U with P @ A = L @ U (equivalently A = P.T @ L @ U), NOT A = P @ L @ U.

Define QR factorization.

A=QR

with Q orthogonal (

Q^TQ=I

) and R upper triangular.

Why does QR solve least squares stably?

Orthogonal Q preserves length, so min reduces to

Rx=Q^Tb

; avoids normal equations which square the condition number.

What sits on the diagonal of T in the Schur form

A=ZTZ^H

The eigenvalues of A.

Why prefer Schur over diagonalization numerically?

Schur always exists and uses unitary Z (condition number 1), even for defective/non-diagonalizable matrices.

Which scipy feature does numpy.linalg lack for eigenproblems?

The generalized eigenproblem

Ax=\lambda Bx

via eig(A, b=B).

Relationship between κ(A) and the normal equations?

Normal equations have condition number

\kappa(A)^2

, making them less accurate than QR.

Connections

numpy.linalg — the convenience subset; same LAPACK backend
Gaussian Elimination — the algorithm LU records
Gram-Schmidt Orthogonalization — the idea behind QR
Eigenvalues and Eigenvectors — what Schur reveals
Condition Number — why algorithm choice matters
Least Squares — solved stably by QR
Cholesky Decomposition — the SPD special case of LU
LAPACK and BLAS — the Fortran muscle underneath

Concept Map

Hinglish (regional understanding)

Intuition Hinglish mein samjho

Dekho, matrix ek "machine" hai jo vectors ko transform karti hai. Seedha-seedha is machine ke saath kaam karna mushkil hota hai, isliye hum use simple pieces me todte hain — yahi "factorization" kehlata hai. scipy.linalg me yeh saare tools available hain, aur numpy.linalg se zyada complete hai: isme lu, qr, schur, generalized eigenvalues sab milte hain.

Teen main factorizations yaad rakho. LU ( $PA=LU$ ) = Gaussian elimination ko record karna, jise hum $Ax=b$ solve karne ke liye use karte hain — pehle $Ly=Pb$ (forward), phir $Ux=y$ (back). QR ( $A=QR$ ) me Q orthogonal hota hai (length preserve karta hai), isliye least-squares problem $\|Ax-b\|$ minimize karna easy ho jaata hai bina normal equations ( $A^TA$ ) use kiye — kyunki normal equations condition number ko square kar dete hain, accuracy gir jaati hai. Schur ( $A=ZTZ^H$ ) me eigenvalues T ke diagonal pe baith jaate hain, aur yeh hamesha exist karta hai even jab matrix diagonalize nahi hoti — kyunki Z unitary hai (perfectly conditioned).

Sabse important practical baat: kabhi inv(A)@b mat likho. Yeh slow bhi hai aur rounding error ko badhata hai. Hamesha solve(A,b) use karo. Note: scipy.linalg.lu(A) ka convention hai $P @ A = L @ U$ (yaani $A = P^T @ L @ U$ ), galti se $A=P@L@U$ mat samajhna. Mnemonic: "LU Solves, QR Squares, Schur Spectra." Bas yeh ek line yaad rahi toh poora chapter ka 80% cover ho gaya.

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

WHY scipy.linalg over numpy.linalg

The three core factorizations

1. LU — for solving $Ax=b$

2. QR — for least squares & orthogonalization

3. Schur — for eigenvalues & matrix functions

Connections

Concept Map

Hinglish (regional understanding)

Go deeper — visual, from zero

Test yourself — Scientific Computing (Python)

Connections