Skip to main content

Subsection 5.6.1 Additional homework

In this chapter, we discussed how the LU factorization (with pivoting) can be used to solve \(A x = y \text{.}\) Why don't we instead discuss how to compute the inverse of the matrix \(A\) and compute \(x = A^{-1} y \text{?}\) Through a sequence of exercises, we illustrate why one should (almost) never compute the inverse of a matrix.

Homework 5.6.1.1.

Let \(A \in \mathbb C^{m \times m} \) be nonsingular and \(B \) its inverse. We know that \(A B = I \) and hence

\begin{equation*} A \left( \begin{array}{c | c | c} b_0 \amp \cdots \amp b_{m-1} \end{array} \right) = \left( \begin{array}{c | c | c} e_0 \amp \cdots \amp e_{m-1} \end{array} \right), \end{equation*}

where \(e_j \) can be thought of as the standard basis vector indexed with \(j \) or the column of \(I \) indexed with \(j \text{.}\)

  1. Justify the following algorithm for computing \(B \text{:}\)

    \begin{equation*} \begin{array}{l} {\bf for~} j = 0, \ldots , m-1 \\ ~~~ \mbox{Compute the LU factorization with pivoting }: P( p ) A = L U \\ ~~~ \mbox{Solve } L z = P( p ) e_j \\ ~~~ \mbox{Solve } U b_j = z \\ {\bf endfor} \end{array} \end{equation*}

  2. What is the cost, in flops, of the above algorithm?

  3. How can we reduce the cost in the most obvious way and what is the cost of this better algorithm?

  4. If we want to solve \(A x = y \) we can now instead compute \(x = B y \text{.}\) What is the cost of this multiplication and how does this cost compare with the cost of computing it via the LU factorization, once the LU factorization has already been computed:

    \begin{equation*} \begin{array}{l} ~~~ \mbox{Solve } L z = P( p ) y \\ ~~~ \mbox{Solve } U x = z \end{array} \end{equation*}

What do we conclude about the wisdom of computing the inverse?

Homework 5.6.1.2.

Let \(L \) be a unit lower triangular matrix. Partition

\begin{equation*} L = \left( \begin{array}{c | c} 1 \amp 0 \\ \hline l_{21} \amp L_{22} \end{array} \right). \end{equation*}

  1. Show that

    \begin{equation*} L^{-1} = \left( \begin{array}{c | c} 1 \amp 0 \\ \hline - L_{22}^{-1} l_{21} \amp L_{22}^{-1} \end{array} \right). \end{equation*}

  2. Use the insight from the last part to complete the following algorithm for computing the inverse of a unit lower triangular matrix:

    \begin{equation*} \newcommand{\routinename}{[ L ] = inv( L )} \newcommand{\guard}{ n( L_{TL} ) \lt n( L ) } \newcommand{\partitionings}{ L \rightarrow \FlaTwoByTwo{L_{TL}}{L_{TR}} {L_{BL}}{L_{BR}} } \newcommand{\partitionsizes}{ L_{TL} {\rm ~is~} 0 \times 0 } \newcommand{\repartitionings}{ \FlaTwoByTwo{L_{TL}}{L_{TR}} {L_{BL}}{L_{BR}} \rightarrow \FlaThreeByThreeBR{L_{00}}{l_{01}}{L_{02}} {l_{10}^T}{\lambda_{11}}{l_{12}^T} {L_{20}}{l_{21}}{L_{22}} } \newcommand{\moveboundaries}{ \FlaTwoByTwo{L_{TL}}{L_{TR}} {L_{BL}}{L_{BR}} \leftarrow \FlaThreeByThreeTL{L_{00}}{l_{01}}{L_{02}} {l_{10}^T}{\lambda_{11}}{l_{12}^T} {L_{20}}{l_{21}}{L_{22}} } \newcommand{\update}{ \begin{array}{l} ~~~~~~~~~~~~~~~~~~\\ l_{21} := ~~~~~~~~~~~~~~~~~~~~~~~~~\\ ~~~~~~~~~~~~~~~~~~ \end{array} } \FlaAlgorithm \end{equation*}

  3. The correct algorithm in the last part will avoid inverting matrices and will require, approximately, \(\frac{1}{3} m^3 \) flops. Analyze the cost of your algorithm.

Homework 5.6.1.3.

LINPACK, the first software package for computing various operations related to solving (dense) linear systems, includes routines for inverting a matrix. When a survey was conducted to see what routines were in practice most frequently used, to the dismay of the developers, it was discovered that routine for inverting matrices was among them. To solve \(A x = y \) users were inverting \(A\) and then computing \(x = A^{-1} y \text{.}\) For this reason, the successor to LINPACK, LAPACK, does not even include a routine for inverting a matrix. Instead, if a user wants to compute the inverse, the user must go through the steps.

\begin{equation*} \begin{array}{l} ~~~ \mbox{Compute the LU factorization with pivoting }: P( p ) A = L U \\ ~~~ \mbox{Invert } L, \mbox{ overwriting } L \mbox{ with the result} \\ ~~~ \mbox{Solve } U X = L \mbox{ for } X\\ ~~~ \mbox{Compute } A^{-1} := X P( p ) \mbox{ (permuting the columns of } X \mbox{)} \end{array} \end{equation*}

  1. Justify that the described steps compute \(A^{-1} \text{.}\)

  2. Propose an algorithm for computing \(X \) that solves \(U X = L \text{.}\) Be sure to take advantage of the triangular structure of \(U \) and \(L \text{.}\)

  3. Analyze the cost of the algorithm in the last part of this question. If you did it right, it should require, approximately, \(m^3 \) operations.

  4. What is the total cost of inverting the matrix?

login