Implementation Issues

Exercise 8.1

(a)

\[\begin{split}& \begin{bmatrix} 2 & 5 & & 6 & \\ 1 & 1 & 3 & 9 & 6\\ & & 2 & 6 & 4\\ & & & 4 & 1\\ & & -1 & -3 & -1 \end{bmatrix} \begin{matrix} \\ R_2 - 0.5 R_1\\ \\ \\ \\ \end{matrix}\\ \Rightarrow & \begin{bmatrix} 2 & 5 & & 6 & \\ 1 & -1.5 & 3 & 6 & 6\\ & & 2 & 6 & 4\\ & & & 4 & 1\\ & & -1 & -3 & -1 \end{bmatrix} \begin{matrix} \\ \\ \\ \\ R_5 - (-0.5) R_3\\ \end{matrix}\\ \Rightarrow & \begin{bmatrix} 2 & 5 & & 6 & \\ 1 & -1.5 & 3 & 6 & 6\\ & & 2 & 6 & 4\\ & & & 4 & 1\\ & & -1 & 0 & 1 \end{bmatrix}\end{split}\]

\[\begin{split}B &= \begin{bmatrix} 2 & & & & \\ 1 & -1.5 & & & \\ & & 2 & & \\ & & & 4 & \\ & & -1 & 0 & 1 \end{bmatrix} \begin{bmatrix} 2 & & & & \\ & -1.5 & & & \\ & & 2 & & \\ & & & 4 & \\ & & & & 1 \end{bmatrix}^{-1} \begin{bmatrix} 2 & 5 & & 6 & \\ & -1.5 & 3 & 6 & 6\\ & & 2 & 6 & 4\\ & & & 4 & 1\\ & & & & 1 \end{bmatrix}\\ &= \begin{bmatrix} 1 & & & & \\ 0.5 & 1 & & & \\ & & 1 & & \\ & & & 1 & \\ & & -0.5 & & 1 \end{bmatrix} \begin{bmatrix} 2 & 5 & & 6 & \\ & -1.5 & 3 & 6 & 6\\ & & 2 & 6 & 4\\ & & & 4 & 1\\ & & & & 1 \end{bmatrix}\\ &= LU\end{split}\]

(b)

Forward substitution gives

\[\begin{split}B \Delta x_\mathcal{B} &= a_j\\ LU \Delta x_\mathcal{B} &= a_j\\ Ly &= a_j\\ y &= L^{-1} a_j\\ &= L^{-1} \begin{bmatrix} 0\\ 2\\ 1\\ 3\\ 0 \end{bmatrix} &= \begin{bmatrix} 0\\ 2\\ 1\\ 3\\ 0.5 \end{bmatrix}.\end{split}\]

Backward substitution yields

\[\begin{split}U \Delta x_\mathcal{B} &= y\\ \Delta x_\mathcal{B} &= U^{-1} y\\ &= \begin{bmatrix} 2.083\\ -1.583\\ -2.375\\ 0.625\\ 0.5 \end{bmatrix}.\end{split}\]

(c)

Suppose that \(\tilde{B}\) is \(B\) (8.9) with the second column replaced by

\[a_j = \begin{bmatrix} 0 & 2 & 1 & 3 & 0 \end{bmatrix}^\top.\]

Suppose that \(\tilde{a}_j = \begin{bmatrix} 1 & 0 & -1 & 0 & 0 \end{bmatrix}^\top\) and the goal is to solve

\[\tilde{B} \Delta \tilde{x}_\mathcal{B} = BE \Delta \tilde{x}_\mathcal{B} = \tilde{a}_j.\]

Applying the forward and backward substitutions gives

\[\begin{split}B u &= \tilde{a}_j\\ LU u &= \tilde{a}_j\\ u &= \begin{bmatrix} 2.417\\ -0.9167\\ 0.125\\ 0.125\\ -0.5 \end{bmatrix}.\end{split}\]

The desired quantity is

\[\begin{split}E \Delta \tilde{x}_\mathcal{B} &= u\\ \Delta \tilde{x}_\mathcal{B} &= E^{-1} u\\ &= u - \frac{u_2}{\Delta x_2} \left( \Delta x_\mathcal{B} - e_2 \right)\\ &= \begin{bmatrix} 2.417\\ -0.9167\\ 0.125\\ 0.125\\ -0.5 \end{bmatrix} - \frac{-0.9167}{-1.583} \left( \begin{bmatrix} 2.083\\ -1.583\\ -2.375\\ 0.625\\ 0.5 \end{bmatrix} - \begin{bmatrix} 0\\ 1\\ 0\\ 0\\ 0 \end{bmatrix} \right)\\ &= \begin{bmatrix} 1.210\\ 0.5791\\ 1.500\\ -0.2369\\ -0.7895 \end{bmatrix}.\end{split}\]

Since \(\tilde{B} \Delta \tilde{x}_\mathcal{B}\) gives \(\tilde{a}_j\), the solution is correct.

(d)

Recall from (8.5) that

\[\tilde{B} = B + (a_j - a_i) e_i^\top\]

where \(a_j\) is associated with the entering variable \(x_j\) and \(a_i\) is associated with the leaving variable \(x_i\).

Suppose that \(\tilde{B}\) is \(B\) (8.9) with the second column replaced by

\[\begin{split}a_j = \begin{matrix} 1\\ 2\\ 3\\ 4\\ 5 \end{matrix} \begin{bmatrix} 0\\ 2\\ 1\\ 3\\ 0 \end{bmatrix} \Rightarrow \tilde{B} = \begin{bmatrix} 2 & & & 6 & \\ 1 & 2 & 3 & 9 & 6\\ & 1 & 2 & 6 & 4\\ & 3 & & 4 & 1\\ & & -1 & -3 & -1 \end{bmatrix}.\end{split}\]

The goal is to solve

\[\tilde{B} \Delta \tilde{x}_\mathcal{B} = L \tilde{E}^{-1} \tilde{U} \Delta \tilde{x}_\mathcal{B} = \tilde{a}_j\]

where

\[\begin{split}\tilde{a}_j = \begin{matrix} 1\\ 2\\ 3\\ 4\\ 5 \end{matrix} \begin{bmatrix} 1\\ 0\\ -1\\ 0\\ 0 \end{bmatrix}.\end{split}\]

Since the decomposition of \(B = LU\) did not use any permutations, the LU-factorization is updated as

\[\begin{split}L^{-1} \tilde{B} = \begin{bmatrix} 2 & & & 6 & \\ & 2 & 3 & 6 & 6\\ & 1 & 2 & 6 & 4\\ & 3 & & 4 & 1\\ & 0.5 & & & 1 \end{bmatrix}\end{split}\]

without permuting the columns and rows. To remove the spike and convert the matrix back into upper-triangular form, define a permutation matrix

\[\begin{split}P = \begin{bmatrix} 1 & & & & \\ & & 1 & & \\ & & & 1 & \\ & & & & 1\\ & 1 & & & \end{bmatrix}\end{split}\]

and

\[\begin{split}\tilde{E} = \begin{bmatrix} 1 & & & & \\ & 1 & & & \\ & & 1 & & \\ & & & 1 & \\ & & & -3 & 1 \end{bmatrix} \begin{bmatrix} 1 & & & & \\ & 1 & & & \\ & & 1 & & \\ & & & 1 & \\ & & 3 & & 2 \end{bmatrix} \begin{bmatrix} 1 & & & & \\ & 1 & & & \\ & & 1 & & \\ & & & 1 & \\ & -3 & & & 2 \end{bmatrix} = \begin{bmatrix} 1 & & & & \\ & 1 & & & \\ & & 1 & & \\ & & & 1 & \\ & -6 & 3 & -3 & 4 \end{bmatrix}\end{split}\]

so that

\[\begin{split}\tilde{E} P L^{-1} \tilde{B} P^\top = \begin{bmatrix} 2 & & 6 & & \\ & 2 & 6 & 4 & 1\\ & & 4 & 1 & 3\\ & & & 1 & 0.5\\ & & & & 9.5 \end{bmatrix} = \tilde{U}.\end{split}\]

Forward substitution gives

\[\begin{split}L \left( P^\top y \right) = \tilde{a}_j \quad \Rightarrow \quad y = \begin{bmatrix} 1\\ -1\\ 0\\ -0.5\\ -0.5 \end{bmatrix}.\end{split}\]

Given \(y\), the next step is to solve

\[\begin{split}\tilde{E}^{-1} z = y \quad \Rightarrow \quad z = \tilde{E} y = \begin{bmatrix} 1\\ -1\\ 0\\ -0.5\\ 5.5 \end{bmatrix}.\end{split}\]

Finally, backward substitution yields

\[\begin{split}\tilde{U} \left( P \Delta \tilde{x}_\mathcal{B} \right) = z \quad \Rightarrow \quad \Delta \tilde{x}_\mathcal{B} = \begin{bmatrix} 1.21\\ 0.58\\ 1.50\\ -0.24\\ -0.79 \end{bmatrix}.\end{split}\]

Exercise 8.2

Let \(P_{ij}\) denote a permutation matrix that swaps rows \(i\) and \(j\) when post-multiplied by a matrix; \(P_{ij}\) swaps columns \(i\) and \(j\) when pre-multiplied by a matrix. By inspection, \(P_{ij}^2 = P_{ij} P_{ij} = I\).

\[\begin{split}& P_{14} \begin{bmatrix} 2 & 5 & & 6 & \\ 1 & 1 & 3 & 9 & 6\\ & & 2 & 6 & 4\\ & & & 4 & 1\\ & & -1 & -3 & -1 \end{bmatrix} P_{15} &= \begin{bmatrix} 1 & & & 4 & \\ 6 & 1 & 3 & 9 & 1\\ 4 & & 2 & 6 & \\ & 5 & & 6 & 2\\ -1 & & -1 & -3 & \end{bmatrix} \begin{matrix} \\ R_2 - 6 R_1\\ R_3 - 4 R_1\\ \\ R_5 - (-R_1) \end{matrix}\\ \Rightarrow & P_{23} \begin{bmatrix} 1 & & & 4 & \\ 6 & 1 & 3 & -15 & 1\\ 4 & & 2 & -10 & \\ & 5 & & 6 & 2\\ -1 & & -1 & 1 & \end{bmatrix} P_{23} &= \begin{bmatrix} 1 & & & 4 & \\ 4 & 2 & & -10 & \\ 6 & 3 & 1 & -15 & 1\\ & & 5 & 6 & 2\\ -1 & -1 & & 1 & \end{bmatrix} \begin{matrix} \\ \\ R_3 - \frac{3}{2} R_2\\ \\ R_5 - \left( -\frac{1}{2} R_2 \right) \end{matrix}\\ \Rightarrow & P_{35} \begin{bmatrix} 1 & & & 4 & \\ 4 & 2 & & -10 & \\ 6 & 3 & 1 & & 1\\ & & 5 & 6 & 2\\ -1 & -1 & & -4 & \end{bmatrix} P_{34} &= \begin{bmatrix} 1 & & 4 & & \\ 4 & 2 & -10 & & \\ -1 & -1 & -4 & & \\ & & 6 & 5 & 2\\ 6 & 3 & & 1 & 1 \end{bmatrix} \begin{matrix} \\ \\ \\ R_4 - (-\frac{6}{4} R_3)\\ \\ \end{matrix}\\ \Rightarrow & &= \begin{bmatrix} 1 & & 4 & & \\ 4 & 2 & -10 & & \\ -1 & -1 & -4 & & \\ & & 6 & 5 & 2\\ 6 & 3 & & 1 & 1 \end{bmatrix} \begin{matrix} \\ \\ \\ \\ R_5 - \frac{1}{5} R_4 \end{matrix}\\ \Rightarrow & \begin{bmatrix} 1 & & 4 & & \\ 4 & 2 & -10 & & \\ -1 & -1 & -4 & & \\ & & 6 & 5 & 2\\ 6 & 3 & & 1 & 0.6 \end{bmatrix}\end{split}\]

\[\begin{split}\begin{aligned} L = \begin{bmatrix} 1 & & & & \\ 4 & 2 & & & \\ -1 & -1 & -4 & & \\ & & 6 & 5 & \\ 6 & 3 & & 1 & 0.6 \end{bmatrix} \begin{bmatrix} 1 & & & & \\ & 2 & & & \\ & & -4 & & \\ & & & 5 & \\ & & & & 0.6 \end{bmatrix}^{-1} \end{aligned} \quad \text{and} \quad \begin{aligned} U = \begin{bmatrix} 1 & & 4 & & \\ & 2 & -10 & & \\ & & -4 & & \\ & & & 5 & 2\\ & & & & 0.6 \end{bmatrix} \end{aligned}\end{split}\]

Note that the row indices of \(L\) and column indices of \(U\) have been permuted by \(P_R = P_{35} P_{23} P_{14}\) and \(P_C = P_{15} P_{23} P_{34}\) respectively. These need to be accounted for when applying the forward and backward substitutions to solve (8.9). The original linear system can be reconstructed as \(B = P_R^{-1} LU P_C^{-1} = P_R^\top LU P_C^\top\).

The forward substitution now gives

\[\begin{split}B \Delta x_\mathcal{B} &= a_j\\ L y &= P_R a_j\\ y &= L^{-1} P_R a_j.\end{split}\]

The backward substitution now yields

\[\begin{split}U \Delta x'_\mathcal{B} &= y\\ \Delta x'_\mathcal{B} &= U^{-1} y\\ P_{34} P_{23} P_{15} \Delta x_\mathcal{B} &= U^{-1} y\\ \Delta x_\mathcal{B} &= P_C U^{-1} y.\end{split}\]

Exercise 8.3

Given a permutation \(\pi \colon \{1, \ldots, m\} \mapsto \{1, \ldots, m\}\), a permutation matrix can be defined as

\[\begin{split}P_\pi = \begin{bmatrix} \mathbf{e}_{\pi(1)}^\top\\ \vdots\\ \mathbf{e}_{\pi(m)}^\top \end{bmatrix}\end{split}\]

where \(\mathbf{e}_{\pi(j)} = \mathbf{e}_{\pi_j}\) denotes a standard basis vector.

Let \(\pi' \colon \{1, \ldots, m\} \mapsto \{1, \ldots, m\}\) denote a mapping for the columns of \(P_\pi\) such that

\[\begin{split}\mathbf{e}_{\pi'(j)} = \mathbf{e}_{\pi'_j} = \begin{bmatrix} {\mathbf{e}_{\pi_1}}_j^\top\\ \vdots\\ {\mathbf{e}_{\pi_m}}_j^\top \end{bmatrix}\end{split}\]

and \(\pi'_j\) is the row index of column \(j\) that contains the one.

(a)

Let \(B\) be an \(m \times m\) matrix. Pre-multiplying \(B\) by \(P_\pi\)

\[\begin{split}P_\pi B = \begin{bmatrix} \mathbf{e}_{\pi_1}^\top B_{\colon 1} & \cdots & \mathbf{e}_{\pi_1}^\top B_{\colon m}\\ \vdots & \ddots & \vdots\\ \mathbf{e}_{\pi_m}^\top B_{\colon 1} & \cdots & \mathbf{e}_{\pi_m}^\top B_{\colon m} \end{bmatrix} = \begin{bmatrix} B_{\pi_1 1} & \cdots & B_{\pi_1 m}\\ \vdots & \ddots & \vdots\\ B_{\pi_m 1} & \cdots & B_{\pi_m m} \end{bmatrix}\end{split}\]

and post-multiplying \(B\) by \(P_\pi\)

\[\begin{split}B P_\pi = \begin{bmatrix} B_{1 \colon} \mathbf{e}_{\pi'_1} & \cdots & B_{1 \colon} \mathbf{e}_{\pi'_m}^\top\\ \vdots & \ddots & \vdots\\ B_{m \colon} \mathbf{e}_{\pi'_1} & \cdots & B_{m \colon} \mathbf{e}_{\pi'_m}^\top \end{bmatrix} = \begin{bmatrix} B_{1 \pi'_1} & \cdots & B_{1 \pi'_m}\\ \vdots & \ddots & \vdots\\ B_{m \pi'_1} & \cdots & B_{m \pi'_m} \end{bmatrix}\end{split}\]

will yield the same permutation when \(\pi\) is an identity permutation or the elements of the swapped rows and columns have corresponding values. Otherwise, post-multiplying by \(P_\pi^\top\) is the only way to permute the column indices in the same order as row indices.

(b)

As derived in (a), every permutation of the rows of a matrix \(B\) can be represented as a mapping \(\pi\).

(c)

Recall that \(\mathbf{e}_i^\top \mathbf{e}_j = 1\) if and only if \(i = j\); otherwise, the scalar product is zero.

\[\begin{split}P_\pi P_\pi^\top = \begin{bmatrix} \mathbf{e}_{\pi_1}^\top\\ \vdots\\ \mathbf{e}_{\pi_m}^\top \end{bmatrix} \begin{bmatrix} \mathbf{e}_{\pi_1} & \cdots & \mathbf{e}_{\pi_m} \end{bmatrix} = \begin{bmatrix} \mathbf{e}_{\pi_1}^\top \mathbf{e}_{\pi_1} & \cdots & \mathbf{e}_{\pi_1}^\top \mathbf{e}_{\pi_m}\\ \vdots & \ddots & \vdots\\ \mathbf{e}_{\pi_m}^\top \mathbf{e}_{\pi_1} & \cdots & \mathbf{e}_{\pi_m}^\top \mathbf{e}_{\pi_m}\\ \end{bmatrix} = \mathbf{I} \quad \Rightarrow \quad P_\pi^{-1} = P_\pi^\top.\end{split}\]

Exercise 8.4

Given \(B = LU\), solving \(B^\top x = U^\top L^\top x = b\) can be done in two steps.

1. Use forward substitution on \(U^\top y = b\).

2. Follow up with backward substitution on \(L^\top x = y\).