Chapter 10 Exercises — LU and PLU Decomposition

Work these with pencil first, then verify the computational ones against numpy / scipy. Tiers: ⭐ conceptual · ⭐⭐ computation (by hand) · ⭐⭐⭐ proof (A) / coding (C) · ⭐⭐⭐⭐ application. Problems marked [code] want a short Python check; [proof] want a rigorous argument. Remember the index shift: math counts from 1 ($\ell_{21}$), numpy/scipy from 0 (L[1,0]).

⭐ Conceptual

10.1 In one sentence each, say what $L$ records and what $U$ records in the factorization $A=LU$. Why does $L$ have $1$'s on its diagonal?

10.2 True or false, with a reason: Every square matrix has a factorization $A=LU$ with $L$ unit lower-triangular. (If false, what is the correct universal statement?)

10.3 You must solve $A\mathbf{x}=\mathbf{b}$ for one fixed $A$ and $300$ different right-hand sides. Explain in two or three sentences why computing $A=LU$ once beats running full Gaussian elimination $300$ times. Roughly, what is the cost of each approach?

10.4 Why is solving a triangular system cheap ($\sim n^2$) compared to factoring ($\sim n^3$)? Answer in terms of how much work each unknown requires in back-substitution.

10.5 State precisely when pivoting is required (not merely advisable). Give the one-line reason, and give a $2\times2$ matrix that forces a swap.

10.6 A colleague writes x = np.linalg.inv(A) @ b to solve $A\mathbf{x}=\mathbf{b}$. Give the two distinct reasons (one about cost, one about accuracy) this is worse than np.linalg.solve(A, b).

10.7 In the PLU factorization $PA=LU$, what does the matrix $P$ do to $A$, and what is $P^{-1}$ in terms of $P$? Why is that last fact convenient?

⭐⭐ Computation (by hand)

10.8 Factor $A=\begin{bmatrix}1&2\\3&10\end{bmatrix}$ as $A=LU$ by hand. State the single multiplier $\ell_{21}$, write $L$ and $U$, and verify $LU=A$. [code] Confirm with scipy.linalg.lu.

10.9 Factor $A=\begin{bmatrix}2&4&-2\\4&9&-3\\-2&-3&7\end{bmatrix}$ as $A=LU$ (no swaps needed). Show every multiplier. [code] Check np.allclose(L@U, A).

10.10 Using your factors from 10.9, solve $A\mathbf{x}=\mathbf{b}$ with $\mathbf{b}=(4,8,10)$ by forward-substituting $L\mathbf{y}=\mathbf{b}$ then back-substituting $U\mathbf{x}=\mathbf{y}$. Show $\mathbf{y}$ and $\mathbf{x}$. [code] Verify against np.linalg.solve.

10.11 Solve the lower-triangular system $L\mathbf{y}=\mathbf{b}$ by forward substitution, where $L=\begin{bmatrix}1&0&0\\-2&1&0\\3&-1&1\end{bmatrix}$ and $\mathbf{b}=(2,1,4)$. Show each step.

10.12 Solve the upper-triangular system $U\mathbf{x}=\mathbf{y}$ by back substitution, where $U=\begin{bmatrix}3&-1&2\\0&4&1\\0&0&-2\end{bmatrix}$ and $\mathbf{y}=(5,9,-4)$. Show each step.

10.13 The matrix $A=\begin{bmatrix}0&1&3\\2&4&1\\1&1&2\end{bmatrix}$ has a zero in the $(1,1)$ position. Using partial pivoting, determine the permutation $P$ that swaps the first column's largest-magnitude entry to the top, write $PA$, then factor $PA=LU$ by hand. Give $P$, $L$, $U$. [code] Check np.allclose(P@A, L@U).

10.14 For $A=\begin{bmatrix}1&1&1\\2&3&5\\4&6&8\end{bmatrix}$, attempt the no-pivot factorization. At which step does a zero pivot appear? What does that tell you about $A$ (compute its determinant by any means)? Does $A=LU$ still "exist," and is it usable for solving?

10.15 Given the LU factors $L=\begin{bmatrix}1&0\\2&1\end{bmatrix}$, $U=\begin{bmatrix}5&1\\0&3\end{bmatrix}$, solve $A\mathbf{x}=\mathbf{b}$ for both $\mathbf{b}_1=(5,13)$ and $\mathbf{b}_2=(10,16)$ reusing the same factors. (This is factor-once-reuse-many in miniature.)

10.16 Build the $3\times3$ permutation matrix $P$ that sends row order $(1,2,3)\to(3,1,2)$ (row 3 to the top, then 1, then 2). Verify $P^{\mathsf{T}}P=I$ by hand, confirming $P^{-1}=P^{\mathsf{T}}$.

⭐⭐⭐ Proof (A) / Coding (C)

10.17 [proof] Prove that if $A=L_1U_1=L_2U_2$ with $L_1,L_2$ unit lower-triangular and $U_1,U_2$ upper-triangular and invertible, then $L_1=L_2$ and $U_1=U_2$. (Hint: rearrange to $L_2^{-1}L_1=U_2U_1^{-1}$ and argue both sides equal $I$.)

10.18 [proof] Prove that the product of two unit lower-triangular matrices is unit lower-triangular. (Use the definition of matrix multiplication and the locations of the zeros.)

10.19 [proof] Show that the inverse of the elementary matrix $E_{ij}$ — the identity with $-\ell$ in position $(i,j)$, $i>j$ — is the identity with $+\ell$ in position $(i,j)$. Conclude that $E_{ij}^{-1}$ is unit lower-triangular.

10.20 [code] Implement lu_decompose(A) from scratch in pure Python (no numpy in the body), returning $(L,U)$ with $A=LU$, recording each multiplier in $L$. Raise an error on a (near-)zero pivot. Verify on $5$ random $4\times4$ matrices with np.allclose(np.array(L)@np.array(U), A).

10.21 [code] Implement solve_lu(L, U, b) (forward then back substitution, pure Python). Test it against np.linalg.solve on $10$ random systems. Then time the difference between (a) calling np.linalg.solve(A, b_k) in a loop over $200$ right-hand sides and (b) factoring once with scipy.linalg.lu_factor and calling lu_solve in the loop. Report the speedup.

10.22 [code] Implement plu(A) with partial pivoting, returning $(P,L,U)$ with $PA=LU$. Be sure to swap the already-computed rows of $L$ when you swap rows of the working matrix. Verify np.allclose(P@A, L@U) on random matrices, including ones with a zero first pivot.

10.23 [proof] Suppose $A$ is invertible and $PA=LU$. Express $A^{-1}$ in terms of $P$, $L$, $U$ (and their inverses). Explain why this is not the efficient way to solve $A\mathbf{x}=\mathbf{b}$, even though it is correct.

⭐⭐⭐⭐ Application

10.24 [code] Repeated solves in a circuit. A linear circuit's node equations are $G\mathbf{v}=\mathbf{i}$ for a fixed conductance matrix $G$ (symmetric positive-definite) and a sequence of input-current vectors $\mathbf{i}_1,\dots,\mathbf{i}_{500}$ (one per timestep). Take $G$ to be a random $50\times50$ SPD matrix (e.g. M=rng.standard_normal((50,50)); G = M@M.T + 50*np.eye(50)) and $500$ random right-hand sides. (a) Solve all $500$ systems by factoring $G$ once with lu_factor and looping lu_solve. (b) Solve them by re-calling np.linalg.solve in a loop. Confirm the answers match (np.allclose) and report the runtime ratio. (c) In one sentence, connect the result to the chapter's "factor once, reuse many" argument.

10.25 [code] Operation-count reasoning. For $n=200,400,800$, time how long scipy.linalg.lu_factor(A) takes on a random $n\times n$ matrix, and how long a single lu_solve takes. (a) Verify that the factorization time grows roughly like $n^3$ (doubling $n$ should multiply the time by about $8$) while the solve time grows roughly like $n^2$ (a factor of about $4$). (b) Explain why this confirms that the per-right-hand-side cost is negligible against the one-time factorization for large batches.

10.26 [code/application] Why not invert. Construct a moderately ill-conditioned matrix (e.g. a $20\times20$ Hilbert-like matrix, from scipy.linalg import hilbert; A = hilbert(12)) and a known solution $\mathbf{x}_{\text{true}}=\mathbf{1}$, so $\mathbf{b}=A\mathbf{x}_{\text{true}}$. Solve for $\mathbf{x}$ two ways: np.linalg.solve(A,b) and np.linalg.inv(A)@b. Compare each result's error $\lVert \mathbf{x}-\mathbf{x}_{\text{true}}\rVert$ and residual $\lVert A\mathbf{x}-\mathbf{b}\rVert$. Which is more accurate, and does this support the chapter's "never invert to solve" rule? Write two sentences interpreting the numbers.

10.27 [application] Designing a solver interface. You are writing a small finite-element library. Sketch (in prose or pseudocode) the public API you would expose so that a user assembles a stiffness matrix $K$ once and then solves $K\mathbf{u}=\mathbf{f}$ for many loads $\mathbf{f}$ efficiently. Which operation is "factor" and which is "solve," and where in a typical user's loop does each belong? Reference the Common Pitfall from §10.6.