5.4.11 · D5Scientific Computing (Python)
Question bank — scipy.linalg — more stable than numpy.linalg, lu, qr, schur
True or false — justify
scipy.linalg is faster than numpy.linalg because it uses a totally different, secret algorithm.
False. Both call LAPACK; the real difference is scipy is always linked to a complete optimized LAPACK/BLAS and lets you choose the stable path — see LAPACK and BLAS.
Computing inv(A) @ b gives the same numerical answer as solve(A, b).
False. Forming costs ~3× more flops and each entry of the inverse carries its own rounding error, so
solve (one LU + back-substitution) is both faster and more accurate.In scipy.linalg.lu the returned P satisfies A = P @ L @ U.
False. SciPy's convention is , i.e.
A = P.T @ L @ U. Getting this backwards silently corrupts your reconstruction.The L in requires storing extra numbers we had to compute separately.
False. just stores the elimination multipliers we already computed while making — zero extra arithmetic, because "undo subtract " is simply "add ".
Orthogonal matrices preserve the length of every vector they act on.
Every square matrix has a Schur decomposition .
True. Unlike diagonalization, Schur always exists because it only needs a unitary , which you can always build — even for defective matrices.
The normal-equations approach is just as stable as QR for least squares.
False. Forming squares the Condition Number (), amplifying error; QR keeps — see Gram-Schmidt Orthogonalization.
A matrix that is symmetric positive-definite should still be solved with plain LU.
False (wasteful). For SPD matrices use Cholesky Decomposition (
cho_factor) — it's ~2× cheaper and needs no pivoting because the pivots are guaranteed positive.Spot the error
P, L, U = lu(A); assert np.allclose(P @ L @ U, A) — why does the assert fail?
The reconstruction should be
P @ (L @ U) only if A = P @ L @ U, but scipy gives , so the correct check is np.allclose(L @ U, P @ A) or np.allclose(P.T @ L @ U, A).x = inv(A) @ b "to be safe against a bad solve".
This is the opposite of safe — it's the classic instability trap; explicit inversion amplifies rounding. Use
solve(A, b) or lu_solve(lu_factor(A), b).Someone computes eig(A) for a defective matrix to get its diagonalization .
For a defective matrix is singular (not enough independent eigenvectors), so doesn't exist and the result is garbage — use schur instead, which always works.
Q, R = qr(A) then x = solve(R, b) for least squares.
The right-hand side must be projected first: least squares gives , not . Forgetting solves the wrong problem.
Reusing lu(A)'s P,L,U to solve for 100 different b by re-factorizing each time.
Wasteful: factorize once with
lu_factor(A), then call lu_solve(lu_piv, b_i) per right-hand side — the expensive elimination is shared across all .Calling numpy.linalg.eig(A, B) for the generalized problem .
numpy's
eig takes only one argument; the generalized eigenproblem needs scipy.linalg.eig(A, b=B) — this is a feature numpy lacks entirely.Using qr(A) (full mode) then multiplying by the full Q in a tall-skinny least-squares problem.
For the full is a huge matrix; use
mode='economic' so is and is — same answer, far less memory.Why questions
Why do we swap rows (partial pivoting) before eliminating in LU?
To avoid dividing by a tiny (or zero) pivot, which would blow up the multipliers and destroy accuracy — pivoting picks the largest available entry, see Gaussian Elimination.
Why does QR "keep" the condition number while normal equations square it?
QR only multiplies by orthogonal (perfectly conditioned, ) matrices, so the sensitivity of is untouched; forming literally multiplies by itself, squaring how much error can be amplified.
Why does Schur use a unitary rather than a general invertible matrix?
Unitary matrices have condition number exactly , so applying never amplifies rounding error — a general in diagonalization can be arbitrarily ill-conditioned.
Why is solving two triangular systems (forward + back) considered "cheap"?
Triangular solves need no elimination — each unknown is found by direct substitution using already-solved unknowns, costing versus for the factorization itself.
Why does scipy.linalg expose more routines than numpy.linalg?
numpy exposes only the convenience subset (
inv, solve, eig, svd); scipy wraps the full LAPACK interface — lu, qr, schur, cho_factor, banded solvers, generalized eig.Why does lstsq return rank alongside the solution?
The rank tells you whether actually had full column rank; a deficient rank means the least-squares solution isn't unique, which you'd otherwise never notice — relevant to Least Squares.
Edge cases
What does LU do when a pivot is exactly zero and no row swap can fix it?
The matrix is singular — LU cannot proceed (division by zero), signalling has no unique solution;
solve will raise or warn rather than return nonsense.For a matrix, what are , , and ?
, , — the degenerate case still fits trivially, a good sanity check for the convention.
If is already orthogonal, what is its QR factorization?
and (up to sign conventions on the diagonal), since 's columns are already orthonormal — Gram–Schmidt has nothing left to do.
What does the real Schur form look like for a matrix with complex-conjugate eigenvalues?
is block upper-triangular with a block on the diagonal for each complex pair, so everything stays real — you never need complex arithmetic for a real .
For a symmetric matrix, are the eigenvalues from Schur guaranteed real?
Yes — a real symmetric matrix has a Schur form that is actually diagonal with real entries, so all its Eigenvalues and Eigenvectors are real, no blocks appear.
What happens to as approaches singular, and why should you care before calling solve?
, meaning tiny input errors get amplified without bound; a huge Condition Number warns you the computed may be meaningless even if
solve returns a number.Recall One-line summary of the traps
Almost every trap here reduces to one of three rules ::: (1) never form the inverse — factor and substitute, (2) prefer orthogonal/unitary operations because , (3) match the factorization to the matrix's structure (LU / QR / Schur / Cholesky).