Blog

Matrices and Linear Equations

(\(\star\) Chapter 10 in the textbook)

A matrix is a rectangular array of numbers. It has \(m\) rows and \(n\) columns, and its entries \(a_{ij}\) can be real or complex: \[ \mA=\col{a_{11}&\cdots&a_{1n}\\\vdots&&\vdots\\a_{m1}&\cdots&a_{mn}},\quad n\text{ columns},\ m\text{ rows},\quad a_{ij}\in\R\text{ or }a_{ij}\in\C. \]Example \(\mA=\col{2&1\\0&4\\-5&7}\) is a real \(3\times2\) matrix, and \(\mB=\col{1+j&2\\3-j&4+2j}\) is a complex \(2\times2\) matrix.

Definition  Two matrices \(\mA=\{a_{ij}\}_{m\times n}\) and \(\mB=\{b_{ij}\}_{k\times l}\) are equal if they have the same size, \(m=k\) and \(n=l\), and all entries agree: \(a_{ij}=b_{ij}\) for \(1\leq i\leq m\) and \(1\leq j\leq n\).

Definition  Matrices of the same size are added entry by entry. For \(\mA=\{a_{ij}\}_{m\times n}\) and \(\mB=\{b_{ij}\}_{m\times n}\), \[ \mC=\mA+\mB=\{a_{ij}+b_{ij}\}_{m\times n}. \]Example Let \(\mA=\col{1&2&3\\4&5&6}\), \(\mB=\col{1&2&-5\\-1&3&7}\) and \(\mC=\col{1&2\\3&4}\). Then \[ \mA+\mB=\col{2&4&-2\\3&8&13}, \] but \(\mA+\mC\) and \(\mB+\mC\) are not defined, because \(\mC\) has a different size.

Definition  Multiplying a matrix by a scalar \(c\) multiplies every entry by \(c\): for \(\mA=\{a_{ij}\}_{m\times n}\), \(c\mA=\{ca_{ij}\}_{m\times n}\), where \(c\in\R\) (or \(c\in\C\) for complex matrices).

Example For \(\mA=\col{1&2&3\\4&5&6}\) we get \(2\mA=\col{2&4&6\\8&10&12}\) and \((-1)\mA=-\mA=\col{-1&-2&-3\\-4&-5&-6}\).

Theorem 10.2.1 Properties of matrix addition and scalar multiplication

Addition and scalar multiplication of matrices obey the same rules as for numbers. If \(\mA\), \(\mB\) and \(\mC\) are \(m\times n\) matrices, \(\mO\) is the \(m\times n\) zero matrix, and \(\alpha,\beta\) are scalars, then:

1. \(\mA+\mB=\mB+\mA\) (commutativity).

2. \((\mA+\mB)+\mC=\mA+(\mB+\mC)\) (associativity).

3. \(\mA+\mO=\mA\) (identity element of matrix addition).

The remaining properties are:

4. \(\mA+(-\mA)=\mO\) (inverse element of matrix addition).

5. \(\alpha(\beta\mA)=(\alpha\beta)\mA\) (associativity of scalar multiplication).

6. \((\alpha+\beta)\mA=\alpha\mA+\beta\mA\) (distributivity).

7. \(\alpha(\mA+\mB)=\alpha\mA+\alpha\mB\) (distributivity).

8. \(1\cdot\mA=\mA\) (identity element of scalar multiplication).

9. \(0\mA=\mO\).

10. \(\alpha\mO=\mO\).

Definition Matrix multiplication  Matrices can be multiplied when the number of columns of the first equals the number of rows of the second. If \(\mA=\{a_{ij}\}_{m\times n}\) and \(\mB=\{b_{ij}\}_{n\times p}\), their product is the \(m\times p\) matrix \(\mC=\mA\mB=\{c_{ij}\}_{m\times p}\) whose entry \(c_{ij}\) is the dot product of row \(i\) of \(\mA\) with column \(j\) of \(\mB\):

Example Let \(\mA=\col{1&2&3\\4&5&6}\), \(\mB=\col{1&2&4\\-1&0&5\\1&3&1}\) and \(\mC=\col{1&3\\2&4}\). Then \[ \mA\mB=\col{2&11&17\\5&26&47},\qquad\mC\mA=\col{13&17&21\\18&24&30} \] but \(\mA\mC\), \(\mB\mC\), \(\mB\mA\) and \(\mC\mB\) are not defined, because the sizes do not match.

Note that  even when both \(\mA\mB\) and \(\mB\mA\) exist, \(\mA\mB\) is not necessarily equal to \(\mB\mA\). The order matters. If \(\mA\mB=\mB\mA\), we say that \(\mA\) and \(\mB\) commute under multiplication.

The definition of matrix multiplication comes from linear systems. A system of \(m\) equations in \(n\) unknowns, \[ \left.\begin{aligned} a_{11}x_1+a_{12}x_2+\cdots+a_{1n}x_n&=c_1\\a_{21}x_1+a_{22}x_2+\cdots+a_{2n}x_n&=c_2\\&\vdots\\a_{m1}x_1+a_{m2}x_2+\cdots+a_{mn}x_n&=c_m \end{aligned}\right\} \] is exactly the statement that a matrix times a column of unknowns equals a column of right-hand sides: \[ \underbrace{\col{a_{11}&a_{12}&\cdots&a_{1n}\\a_{21}&a_{22}&\cdots&a_{2n}\\\vdots&\vdots&&\vdots\\a_{m1}&a_{m2}&\cdots&a_{mn}}}_{\mA\ m\times n} \underbrace{\col{x_1\\x_2\\\vdots\\x_n}}_{\xx\ n\times1}=\underbrace{\col{c_1\\c_2\\\vdots\\c_m}}_{\vc{c}\ m\times1}. \]

Definition  For a square matrix \(\mA=\{a_{ij}\}_{n\times n}\), the power \(\mA^p\) is the product of \(p\) copies of \(\mA\): \[ \mA^p=\underbrace{\mA\cdot\mA\cdots\mA}_{p\text{ factors}}. \] Definition  A square matrix \(\mA=\{a_{ij}\}_{n\times n}\) is diagonal if \(a_{ij}=0\) for \(i\ne j\), that is, if all its non-zero entries lie on the main diagonal.

Example \(\mA=\col{1&0&0\\0&5&0}\) is not diagonal, because it is not square. \(\mB=\col{1&0\\0&8}\) is a diagonal matrix.

Definition  The identity matrix is the diagonal matrix with \(a_{ii}=1\) on the diagonal.

Example \(\mI_2=\col{1&0\\0&1}\) is the \(2\times2\) identity matrix, and \(\mI_3=\col{1&0&0\\0&1&0\\0&0&1}\) is the \(3\times3\) identity matrix.

The identity matrix leaves any matrix unchanged: \(\mA_{m\times n}\mI_n=\mA_{m\times n}\) and \(\mI_m\mA_{m\times n}=\mA_{m\times n}\). This holds because in each entry of the product, only one term of the sum is non-zero, and it is multiplied by \(1\).

Theorem 10.2.2 Exceptional Properties of matrix multiplication:

Matrix multiplication differs from multiplication of numbers in four important ways, so some familiar rules fail:

1. \(\mA\mB\ne\mB\mA\) in general. For example, with \(\mA=\col{0&1\\0&0}\) and \(\mB=\col{0&0\\1&0}\) we get \(\mA\mB=\col{1&0\\0&0}\) but \(\mB\mA=\col{0&0\\0&1}\).

2. Even if \(\mA\mB=\mA\mC\), the matrices \(\mB\) and \(\mC\) need not be equal. For example, with \(\mA=\col{1&0\\0&0}\), \(\mB=\col{1&0\\0&1}\) and \(\mC=\col{1&0\\0&5}\), both products equal \(\col{1&0\\0&0}\).

3. \(\mA\mB=\mO\) does not imply that \(\mA=\mO\) or \(\mB=\mO\). For example, \(\col{0&1\\0&0}\col{1&0\\0&0}=\mO\).

4. \(\mA^2=\mI\) does not imply \(\mA=\mI\) or \(\mA=-\mI\). For example, \(\mA=\col{1&0\\0&-1}\) has \(\mA^2=\mI\).

Theorem 10.2.3 Ordinary properties of matrix multiplication

The rules that do hold are the ones we need for algebra:

1. \((\alpha\mA)\mB=\alpha(\mA\mB)=\mA(\alpha\mB)\) (scalars can be moved around).

2. \((\mA\mB)\mC=\mA(\mB\mC)\) (associativity).

3. \((\mA+\mB)\mC=\mA\mC+\mB\mC\) (distributivity).

4. \(\mC(\mA+\mB)=\mC\mA+\mC\mB\) (distributivity).

5. \(\mA(\alpha\mB+\beta\mC)=\alpha\mA\mB+\beta\mA\mC\) (linearity).

Partitioning in Matrix Operations

A large matrix can be cut by horizontal and vertical lines into smaller matrices called blocks, and the blocks can then be treated as if they were entries: \[ \mA=\col{\mA_{11}&\mA_{12}&\cdots&\mA_{1n}\\\vdots&\vdots&&\vdots\\\mA_{m1}&\mA_{m2}&\cdots&\mA_{mn}},\quad \mB=\col{\mB_{11}&\cdots&\mB_{1q}\\\vdots&&\vdots\\\mB_{p1}&\cdots&\mB_{pq}}. \] If \(m=p\), \(n=q\) and each block \(\mA_{ij}\) has the same size as the corresponding block \(\mB_{ij}\), the sum is computed block by block: \[ \mA+\mB=\col{\mA_{11}+\mB_{11}&\cdots&\mA_{1n}+\mB_{1n}\\\vdots&&\vdots\\\mA_{m1}+\mB_{m1}&\cdots&\mA_{mn}+\mB_{mn}}. \] If \(n=p\), the product \(\mC=\mA\mB\) has blocks \[ \mC_{ij}=\sum_{k=1}^n\mA_{ik}\mB_{kj}, \] provided that each product \(\mA_{ik}\mB_{kj}\) is defined.

Example Cut two \(3\times3\) matrices into blocks as shown by the lines: \[ \mA=\left[\begin{array}{cc|c}2&4&1\\-1&3&0\\\hline5&4&6\end{array}\right]=\col{\mA_{11}&\mA_{12}\\\mA_{21}&\mA_{22}},\quad \mB=\left[\begin{array}{cc|c}0&1&3\\2&-4&-1\\\hline5&8&2\end{array}\right]=\col{\mB_{11}&\mB_{12}\\\mB_{21}&\mB_{22}}. \] The sum and the product are then computed block by block: \[ \mA+\mB=\col{\mA_{11}+\mB_{11}&\mA_{12}+\mB_{12}\\\mA_{21}+\mB_{21}&\mA_{22}+\mB_{22}} =\left[\begin{array}{cc|c}2&5&4\\1&-1&-1\\\hline10&12&8\end{array}\right]. \] \begin{gather*} \mA\mB=\col{\mA_{11}\mB_{11}+\mA_{12}\mB_{21}&\mA_{11}\mB_{12}+\mA_{12}\mB_{22}\\\mA_{21}\mB_{11}+\mA_{22}\mB_{21}&\mA_{21}\mB_{12}+\mA_{22}\mB_{22}}\\ =\left[\begin{array}{cc|c}13&-6&4\\6&-13&-6\\\hline38&37&23\end{array}\right]. \end{gather*} Partitioning is especially helpful for sparse matrices, where most of the entries are zero and many blocks can be skipped.

Definition  The transpose of \(\mA=\{a_{ij}\}_{m\times n}\) swaps its rows and columns: \[ \mA^T=\{a_{ij}^T\}_{n\times m}\quad\text{where}\quad a_{ij}^T=a_{ji}. \]Example For \(\mA=\col{1&2&3\\4&5&6}\) we get \(\mA^T=\col{1&4\\2&5\\3&6}\).

Theorem 10.3.1

Transposition behaves well with the other operations, except that it reverses the order of a product:

1. \((\mA^T)^T=\mA\).

2. \((\mA+\mB)^T=\mA^T+\mB^T\).

3. \((\alpha\mA)^T=\alpha\mA^T\).

4. \((\mA\mB)^T=\mB^T\mA^T\), and so \((\mA\mB\mC\mD)^T=\mD^T\mC^T\mB^T\mA^T\).

Proof for (iv)  Let \(\mC=\mA\mB\), so \(c_{ij}=\sum_{k=1}^n a_{ik}b_{kj}\). The \((i,j)\) entry of \(\mC^T\) is the \((j,i)\) entry of \(\mC\): \[ c_{ij}^T=c_{ji}=\sum_{k=1}^n a_{jk}b_{ki}=\sum_{k=1}^n a_{kj}^Tb_{ik}^T=\sum_{k=1}^n b_{ik}^Ta_{kj}^T \quad\Rightarrow\mC^T=\mB^T\mA^T.\quad\blacksquare \] Note that for column vectors \(\xx=\col{x_1\\\vdots\\x_n}\) and \(\yy=\col{y_1\\\vdots\\y_n}\) in \(\R^n\), the matrix product of a row with a column is a dot product:

Definition

1. If \(\mA^T=\mA\), the matrix \(\mA\) is called symmetric.

2. If \(\mA^T=-\mA\), the matrix \(\mA\) is called skew-symmetric.

Every square matrix is the sum of a symmetric part and a skew-symmetric part: \[ \mA=\tfrac12(\mA+\mA^T)+\tfrac12(\mA-\mA^T). \] The first term is symmetric, because \((\mA+\mA^T)^T=\mA^T+\mA\), and the second is skew-symmetric, because \((\mA-\mA^T)^T=\mA^T-\mA=-(\mA-\mA^T)\).

Determinants

The determinant is a single number computed from a square matrix. Among other things, it tells us whether the matrix is invertible, as we will see later in this section. For a square matrix \(\mA=\{a_{ij}\}_{n\times n}\) it is defined by cofactor expansion along a row:

The minor \(M_{ik}\) is the determinant of the \((n-1)\times(n-1)\) matrix that is left when row \(i\) and column \(k\) of \(\mA\) are removed. The determinant of a \(1\times1\) or \(2\times2\) matrix is the starting point of the recursion.

Example Small determinants.

1. For a \(1\times1\) matrix, \(\det[a_{11}]=a_{11}\).

2. For a \(2\times2\) matrix, expanding along the first row gives the familiar formula \[ \det\col{a_{11}&a_{12}\\a_{21}&a_{22}}=a_{11}(-1)^{1+1}|a_{22}|+a_{12}(-1)^{1+2}|a_{21}|=a_{11}a_{22}-a_{12}a_{21}. \]

3. For a \(3\times3\) matrix, expanding along the first row reduces the problem to three \(2\times2\) determinants: \begin{align*} \det\col{4&7&-2\\0&3&2\\1&5&6} &=4(-1)^{1+1}\begin{vmatrix}3&2\\5&6\end{vmatrix}+7(-1)^{1+2}\begin{vmatrix}0&2\\1&6\end{vmatrix}\\ &\quad+(-2)(-1)^{1+3}\begin{vmatrix}0&3\\1&5\end{vmatrix}\\ &=4\cdot8-7\cdot(-2)-2\cdot(-3)=32+14+6=52. \end{align*}

Comments

1. The expansion works along any row, so it is best to choose the row with the most zero entries, since zero entries contribute nothing.

2. Cofactor expansion of large matrices is extremely time consuming: the number of operations for an \(n\times n\) matrix grows like \(N(n)\sim e\cdot n!\), which is enormous even for moderate \(n\).

Properties of Determinants

Row operations change a determinant in a simple, predictable way. This gives a much faster way to compute determinants than cofactor expansion:

• D1. Adding a multiple of one row to another does not change the determinant.

• D2. Interchanging two rows changes the sign of the determinant.

• D3. The determinant of a triangular matrix is the product of its diagonal entries.

Example  To compute the determinant of the matrix below, reduce it to triangular form using D1 and D2, and then multiply the diagonal. We first interchange the first and third rows (which flips the sign), then use \(-2R_1+R_2\to R_2\) and \(-\frac23R_2+R_3\to R_3\). \begin{align*} \det\col{0&2&-1\\4&3&5\\2&0&-4} &=-\begin{vmatrix}2&0&-4\\4&3&5\\0&2&-1\end{vmatrix} =-\begin{vmatrix}2&0&-4\\0&3&13\\0&2&-1\end{vmatrix}\\ &=-\begin{vmatrix}2&0&-4\\0&3&13\\0&0&-29/3\end{vmatrix} =-2\cdot3\cdot(-29/3)=58. \end{align*} Reducing to triangular form takes about \(N_T(n)\sim\frac{2n^3}{3}\) operations, which is far fewer than \(e\cdot n!\) when \(n\) is large.

Additional Properties of Determinants

• D4. If \(\mA\) has a row or a column of zeros, then \(\det\mA=0\).

• D5. If two rows (or two columns) are proportional, then \(\det\mA=0\).

• D6. If a row or column is a linear combination of the other rows or columns, then \(\det\mA=0\).

• D7. Multiplying every entry of one row or column by \(\alpha\) multiplies the determinant by \(\alpha\).

• D8. \(\det(\alpha\mA)=\alpha^n\det\mA\) for an \(n\times n\) matrix, because all \(n\) rows are multiplied by \(\alpha\).

• D9. If \(\mA\), \(\mB\) and \(\mC\) are identical except in one row (or column), and that row of \(\mC\) is the sum of the corresponding rows of \(\mA\) and \(\mB\), then \(\det\mC=\det\mA+\det\mB\).

• D10. \(\det(\mA^T)=\det\mA\).

• D11. In general \(\det(\mA+\mB)\ne\det\mA+\det\mB\). D9 is only a special case in which the sum of determinants does work.

• D12. \(\det(\mA\mB)=(\det\mA)(\det\mB)\).

Row and Column Spaces

A matrix can be read in two ways: as a stack of rows, or as a row of columns. Write the rows of \(\mA\) as \(\vc{r}_i^T\) (so each \(\vc{r}_i\in\R^n\)) and the columns as \(\av_i\) (so each \(\av_i\in\R^m\)): \[ \mA=\col{a_{11}&a_{12}&\cdots&a_{1n}\\a_{21}&a_{22}&\cdots&a_{2n}\\\vdots&\vdots&&\vdots\\a_{m1}&a_{m2}&\cdots&a_{mn}} =\col{\vc{r}_1^T\\\vc{r}_2^T\\\vdots\\\vc{r}_m^T},\quad\vc{r}_i=\col{a_{i1}\\a_{i2}\\\vdots\\a_{in}}\in\R^n. \] \[ \mA=\begin{bmatrix}\av_1&\cdots&\av_n\end{bmatrix}\quad\text{where }\av_i=\col{a_{1i}\\a_{2i}\\\vdots\\a_{mi}}\in\R^m. \] Definition  The row space of \(\mA\) is the span of its rows, and the column space is the span of its columns: \[ \operatorname{Row}(\mA)=\Sp\{\vc{r}_1,\ldots,\vc{r}_m\},\qquad\operatorname{Col}(\mA)=\Sp\{\av_1,\ldots,\av_n\}. \] Definition  The row rank is \(\dim(\operatorname{Row}(\mA))\) and the column rank is \(\dim(\operatorname{Col}(\mA))\).

Example Take the echelon matrix below, with rows \(\vc{r}_1,\ldots,\vc{r}_4\). \[ \mA=\col{1&-3&2&5&3\\0&0&1&-4&2\\0&0&0&1&7\\0&0&0&0&0},\qquad \vc{r}_1=\col{1\\-3\\2\\5\\3},\ \vc{r}_2=\col{0\\0\\1\\-4\\2},\ \vc{r}_3=\col{0\\0\\0\\1\\7},\ \vc{r}_4=\col{0\\0\\0\\0\\0}. \] The zero row adds nothing, so the row space is spanned by the three non-zero rows, and the row rank is \(3\): \[ \operatorname{Row}(\mA)=\Sp\{\vc{r}_1,\vc{r}_2,\vc{r}_3,\vc{r}_4\}=\Sp\{\vc{r}_1,\vc{r}_2,\vc{r}_3\}. \] The columns that contain the pivots (the leading or pivot columns) are \(\av_1,\av_3,\av_4\), and they span the column space, so the column rank is also \(3\): \[ \operatorname{Col}(\mA)=\Sp\{\av_1,\av_2,\av_3,\av_4,\av_5\}=\Sp\{\av_1,\av_3,\av_4\}. \]

Proof  If \(\mA\) is row equivalent to \(\mB\), then \(\mB\) is obtained from \(\mA\) by a finite sequence of elementary row operations. Every elementary row operation replaces a row by a linear combination of rows, so every row of \(\mB\) is a linear combination of the rows of \(\mA\). Hence \(\operatorname{Row}(\mB)\subset\operatorname{Row}(\mA)\). Row operations can be undone by row operations, so \(\mA\) can also be obtained from \(\mB\), and the same argument gives \(\operatorname{Row}(\mA)\subset\operatorname{Row}(\mB)\). Therefore \(\operatorname{Row}(\mA)=\operatorname{Row}(\mB)\). \(\blacksquare\)

Proof  Let \(\ev_1,\ev_2,\ldots,\ev_k\) be the non-zero rows of an echelon matrix \(\mE\), and let the first non-zero entry of \(\ev_i\) be in column \(p_i\). Since \(\mE\) is in echelon form, these positions move to the right as we go down: \(p_1<p_2<\cdots<p_k\). The picture shows the pattern, with the first non-zero entry of each row marked.

Suppose, for contradiction, that these rows are linearly dependent. Then there are coefficients \(\alpha_1,\ldots,\alpha_k\), not all zero, such that \[ \alpha_1\ev_1+\cdots+\alpha_k\ev_k=\nul. \] Look at column \(p_1\) of this equation: \[ \alpha_1e_{1p_1}+\alpha_2e_{2p_1}+\cdots+\alpha_ke_{kp_1}=0. \] All the rows below the first have a zero in column \(p_1\), that is, \(e_{2p_1}=e_{3p_1}=\cdots=e_{kp_1}=0\), while \(e_{1p_1}\ne0\). Hence \(\alpha_1=0\).

Now look at column \(p_2\). The first row’s coefficient is already zero, and the rows below the second have zeros in this column, so \(\alpha_2e_{2p_2}=0\) and \(\alpha_2=0\).

Continuing like this, we find \(\alpha_i=0\) for \(1\leq i\leq k\). This contradicts the assumption that not all the coefficients are zero, so the non-zero rows are linearly independent. \(\blacksquare\)

So the row space and the row rank of any matrix can be found by reducing it to echelon form with elementary row operations and counting the non-zero rows.

Now we turn to the column space and the column rank.

Proof  Consider the solutions of \(\mA\xx=\nul\). Elementary row operations turn this into \(\mE\xx=\nul\) without changing the solution set, so any dependence among the columns of \(\mA\) is also a dependence among the columns of \(\mE\), and vice versa. The proof continues on the next page.

Writing this out in terms of the columns \(\av_i\) of \(\mA\) and the columns \(\ev_i\) of \(\mE\), we get \[ \begin{bmatrix}\av_1&\cdots&\av_n\end{bmatrix}\col{\alpha_1\\\vdots\\\alpha_n}=\nul \Longleftrightarrow\begin{bmatrix}\ev_1&\cdots&\ev_n\end{bmatrix}\col{\alpha_1\\\vdots\\\alpha_n}=\nul. \] So the columns of \(\mA\) and of \(\mE\) satisfy exactly the same dependence relations.

The pivot columns of the echelon matrix \(\mE\) span its column space and are linearly independent. (To see independence, suppose a combination of the pivot columns is zero. In the last pivot row only the last pivot column has a non-zero entry, so its coefficient must be zero; then the same argument applies to the previous pivot row, and so on, which is the same argument as in the theorem about the rows of an echelon matrix.) Denote the pivot columns by \(\{\ev_{p_1},\ldots,\ev_{p_r}\}\). Every other column \(\ev_k\) of \(\mE\) is a combination of these, so for every non-pivot column \(k\), \[ \ev_k=\sum_{i=1}^r c_{ki}\ev_{p_i}\Rightarrow\ev_k-\sum_{i=1}^r c_{ki}\ev_{p_i}=\nul. \] Let \(\vc{c}_k\) be the vector that has the entry \(c_{ki}\) in position \(p_i\) and zeros elsewhere. Then the relation above says \[ \mE\left(\col{0\\\vdots\\1\\\vdots\\0}-\vc{c}_k\right)=\nul, \] where the \(1\) is in the \(k^{\mathrm{th}}\) position. By the equivalence of the two solution sets, \[ \mA\left(\col{0\\\vdots\\1\\\vdots\\0}-\vc{c}_k\right)=\nul. \] This says that if a column \(\ev_k\) of \(\mE\) is a combination of the pivot columns of \(\mE\), then the matching column \(\av_k\) of \(\mA\) is the same combination of the matching columns \(\{\av_{p_1},\ldots,\av_{p_r}\}\) of \(\mA\). Therefore \[ \Sp\{\av_{p_1},\ldots,\av_{p_r}\}=\Sp\{\av_1,\ldots,\av_n\}. \] These columns are also linearly independent. We show this by contradiction. Suppose some non-zero coefficients solve \[ \begin{bmatrix}\av_{p_1}&\av_{p_2}&\cdots&\av_{p_r}\end{bmatrix}\col{\beta_{p_1}\\\vdots\\\beta_{p_r}}=\nul. \] Here \(p_i\) is the index of the \(i^{\mathrm{th}}\) pivot column. This means \(\mA\zz=\nul\) for the vector \(\zz\) with \[ z_i=\begin{cases}0&\text{if }i\notin\{p_1,\ldots,p_r\}\\\beta_i&\text{if }i\in\{p_1,\ldots,p_r\}.\end{cases} \] Since \(\mA\xx=\nul\) and \(\mE\xx=\nul\) have the same solutions, also \(\mE\zz=\nul\), that is, \[ \begin{bmatrix}\ev_{p_1}&\cdots&\ev_{p_r}\end{bmatrix}\col{\beta_{p_1}\\\vdots\\\beta_{p_r}}=\nul, \] so \(\{\ev_{p_1},\ldots,\ev_{p_r}\}\) would be LD. This is a contradiction, because pivot columns are linearly independent. Thus \(\{\av_{p_1},\ldots,\av_{p_r}\}\) is LI and spans \(\Sp\{\av_1,\ldots,\av_n\}\), so it is a basis for \(\operatorname{Col}(\mA)\).

Hence \(\dim(\operatorname{Col}(\mA))=r=\dim(\operatorname{Col}(\mE))\). \(\blacksquare\)

We can now compare rows and columns. This is one of the central facts of linear algebra.

Proof  Let \(\mE\) be an echelon form of \(\mA\). We have shown two things:

1. \(\operatorname{Row}(\mA)=\operatorname{Row}(\mE)\), so the row rank of \(\mA\) is the number of pivots of \(\mE\).

2. The column rank of \(\mA\) equals the column rank of \(\mE\), which is also the number of pivots of \(\mE\).

Both ranks equal the number of pivots, so they are equal. \(\blacksquare\)

We call this common number the rank of \(\mA\) and write \(r(\mA)\).

Example How many linearly independent vectors are there in \[ \left\{\col{2\\2\\0\\2},\col{1\\4\\3\\1},\col{-3\\-2\\1\\-3},\col{4\\5\\3\\-2}\right\}\ ? \] Put the vectors in the columns of a matrix \(\mA\) and reduce it to echelon form: \begin{gather*} \col{2&1&-3&4\\2&4&-2&5\\0&3&1&3\\2&1&-3&-2} \xrightarrow{R_1/2\to R_1}\col{1&1/2&-3/2&2\\2&4&-2&5\\0&3&1&3\\2&1&-3&-2}\\[2pt] \xrightarrow[ R_4-2R_1\to R_4]{R_2-2R_1\to R_2} \col{1&1/2&-3/2&2\\0&3&1&1\\0&3&1&3\\0&0&0&-6}\\ \xrightarrow{R_3-R_2\to R_3}\col{1&1/2&-3/2&2\\0&3&1&1\\0&0&0&2\\0&0&0&-6} \xrightarrow{R_4+3R_3\to R_4}\col{1&1/2&-3/2&2\\0&3&1&1\\0&0&0&2\\0&0&0&0}. \end{gather*} There are three pivots, in columns 1, 2 and 4, so three of the vectors are linearly independent (the first, second and fourth). The third vector is a combination of the first two: it equals \(-\tfrac53\) of the first plus \(\tfrac13\) of the second.

Example Now consider a system of four equations in six unknowns: \[ \col{1&-1&1&3&0&2\\1&0&3&3&-1&6\\2&-1&2&1&-1&7\\1&0&5&8&-1&7}\col{x_1\\x_2\\x_3\\x_4\\x_5\\x_6}=\col{4\\3\\9\\1}. \] Reduce the augmented matrix to echelon form: \[ \col{1&-1&1&3&0&2&4\\1&0&3&3&-1&6&3\\2&-1&2&1&-1&7&9\\1&0&5&8&-1&7&1} \xrightarrow[\text{Row Operations}]{\text{Elementary}} \col{1&-1&1&3&0&2&4\\0&1&2&0&-1&4&-1\\0&0&2&5&0&1&-2\\0&0&0&0&0&0&0}. \] Here \(x_1,x_2,x_3\) are the pivot variables and \(x_4,x_5,x_6\) are the free variables.

Observations

1. The rank of \(\mA\) and the rank of the augmented matrix \([\mA:\vc{c}]\) are the same, so the system is consistent.

2. There are three free variables, so the solutions form a family with three parameters.

We give the free variables the names \(x_6=\alpha\), \(x_5=\beta\) and \(x_4=\gamma\), and solve for the pivot variables by back substitution: \begin{align*} x_3&=-1-\tfrac12\alpha-\tfrac52\gamma\\ x_2&=1-3\alpha+\beta+5\gamma\\ x_1&=6-\tfrac92\alpha+\beta+\tfrac92\gamma. \end{align*} In vector form this is a constant vector plus a combination of three direction vectors: \[ \xx=\underbrace{\col{6\\1\\-1\\0\\0\\0}}_{\xx_0}+\alpha\underbrace{\col{-9/2\\-3\\-1/2\\0\\0\\1}}_{\xx_1} +\beta\underbrace{\col{1\\1\\0\\0\\1\\0}}_{\xx_2}+\gamma\underbrace{\col{9/2\\5\\-5/2\\1\\0\\0}}_{\xx_3}. \] Observe that:

1. \(\xx_0\) is a particular solution of \(\mA\xx=\vc{c}\) (take \(\alpha=\beta=\gamma=0\)).

2. \(\xx_1\), \(\xx_2\) and \(\xx_3\) solve the homogeneous equation \(\mA\xx=\nul\). Indeed, for every \(\alpha,\beta,\gamma\), \begin{align*} \mA\xx&=\mA(\xx_0+\alpha\xx_1+\beta\xx_2+\gamma\xx_3)\\ &=\underbrace{\mA\xx_0}_{\vc{c}}+\alpha\mA\xx_1+\beta\mA\xx_2+\gamma\mA\xx_3=\vc{c}, \end{align*} so \(\alpha\mA\xx_1+\beta\mA\xx_2+\gamma\mA\xx_3=\nul\) for all \(\alpha,\beta,\gamma\), and therefore \(\mA\xx_i=\nul\) for \(1\leq i\leq3\).

3. The matrix of these three vectors has rank \(3\) (its last three rows contain an identity matrix), so \(\xx_1,\xx_2,\xx_3\) are linearly independent: \[ \begin{bmatrix}\xx_1&\xx_2&\xx_3\end{bmatrix}=\col{-9/2&1&9/2\\-3&1&5\\-1/2&0&-5/2\\0&0&1\\0&1&0\\1&0&0}. \]

This example illustrates a general fact. For \(\mA\xx=\vc{c}\) with \(\mA\in\R^{m\times n}\) (or \(\C^{m\times n}\)), writing \(r(\mA)\) for the rank:

(iii) If the system is consistent and \(n>r(\mA)\), every solution can be written as

where \(\xx_p\) is one particular solution of \(\mA\xx=\vc{c}\) and \(\xx_n\) solves the homogeneous system \(\mA\xx_n=\nul\).

Definition  The set of all solutions of \(\mA\xx=\nul\) is called the null space of \(\mA\).

The number of pivots and the size of the null space always add up to the number of unknowns:

Proof  Consider \(\mA\xx=\nul\) and reduce \(\mA\) to an echelon matrix \(\mE\) to find the pivots and the free variables. The rank \(r(\mA)\) is the number of pivot variables, and \(\dim\operatorname{Null}(\mA)\) is the number of free variables, since each free variable contributes one independent direction. Every one of the \(n\) unknowns is either a pivot variable or a free variable, so \[\underbrace{\#\text{pivots}}_{r(\mA)}+\underbrace{\#\text{free variables}}_{\dim(\operatorname{Null}(\mA))}=n.\quad\blacksquare\]

Inverse Matrix; Cramer’s Rule, Factorization

To solve the scalar equation \(ax=c\) we divide by \(a\), which is the same as multiplying by \(a^{-1}\), as long as \(a\ne0\): \[\underbrace{a^{-1}a}_{1}x=a^{-1}c\Rightarrow x=a^{-1}c.\] We would like to do the same for \(\mA\xx=\vc{c}\): multiply both sides by a suitable matrix \(\mA^{-1}\) such that \[\underbrace{\mA^{-1}\mA}_{\text{we want this to be }\mI}\xx=\mA^{-1}\vc{c}.\] But we have to be careful. For example, the matrix below has a left inverse, a matrix that multiplies it to give the identity: \[ \col{-1&-2\\1&1\\0&1}\col{x_1\\x_2}=\col{1\\1\\1},\quad\text{and note that }\col{1&2&0\\1&1&2}\col{-1&-2\\1&1\\0&1}=\col{1&0\\0&1}. \]

Here is what goes wrong. We might multiply both sides of the system by the matrix \(\col{1&2&0\\1&1&2}\): \[ \col{1&0\\0&1}\col{x_1\\x_2}=\col{1&2&0\\1&1&2}\col{1\\1\\1}=\col{3\\4} \] and conclude that \(x_1=3\), \(x_2=4\) is a solution. But it is not. The new system is not equivalent to the original one: the original system has no solution at all, and multiplying by a matrix can destroy that information. A matrix that is not square, or not of full rank, cannot serve as an inverse in this way, as the next theorem shows.

Proof  Suppose such an \(\mA^{-1}\) existed for a matrix \(\mA\) that is not a full-rank square matrix. We show that this leads to a contradiction. There are two cases to examine.

1. Suppose \(\mA\xx=\vc{c}\) is inconsistent for some \(\vc{c}\), so its solution set is empty. But \((\mA^{-1}\mA)\xx=\mA^{-1}\vc{c}\) reads \(\xx=\mA^{-1}\vc{c}\) and always has exactly one solution. The two systems would then have different solution sets. So \(\mA\xx=\vc{c}\) must be consistent for all \(\vc{c}\), which means that \(\mA\) has full row rank, \(r(\mA_{m\times n})=m\).

2. Suppose \(\mA\xx=\vc{c}\) has more than one solution for some \(\vc{c}\). Again \((\mA^{-1}\mA)\xx=\mA^{-1}\vc{c}\) has the single solution \(\xx=\mA^{-1}\vc{c}\), so the solution sets would differ. Hence \(\mA\xx=\vc{c}\) can never have more than one solution, which means that \(\mA\) has full column rank, \(r(\mA_{m\times n})=n\).

Combining the two cases gives \(r(\mA_{m\times n})=m=n\), so \(\mA\) is square and of full rank. This contradicts the assumption. \(\blacksquare\)

The converse is also true: a full-rank square matrix has an inverse.

Proof  If \(\mA^{-1}\mA=\mI\), then transposing gives \(\mA^T(\mA^{-1})^T=\mI\). We show that this transposed equation has a unique solution. Let \((\mA^{-1})^T=\mB\); if we can find \(\mB\), then \(\mA^{-1}=\mB^T\). The equation \(\mA^T\mB=\mI\) can be written column by column, where \(\ev_i\) is the \(i^{\mathrm{th}}\) column of the identity (a \(1\) in position \(i\), zeros elsewhere): \[ \mA^T\begin{bmatrix}\bv_1&\bv_2&\cdots&\bv_n\end{bmatrix}=\begin{bmatrix}\ev_1&\ev_2&\cdots&\ev_n\end{bmatrix},\quad\ev_i=\col{0\\\vdots\\0\\1\\0\\\vdots\\0}. \] This is \(n\) separate linear systems: \[ \mA^T\bv_i=\ev_i,\qquad1\leq i\leq n. \] Since \(\mA\) is square and of full rank, so is \(\mA^T\). Therefore each system \(\mA^T\bv_i=\ev_i\) is consistent and has a unique solution. So there is exactly one \(\mB\) with \(\mA^T\mB=\mI\). Then \(\mA^{-1}=\mB^T\) exists, because \(\mB^T\mA=(\mA^T\mB)^T=\mI\), and it is unique. \(\blacksquare\)

Definition  A full-rank square matrix is called invertible (or non-singular), and its unique inverse is written \(\mA^{-1}\).

So far we only know \(\mA^{-1}\mA=\mI\). The next theorem shows that the inverse also works from the other side.

Proof  We have shown that \(\mA^{-1}\mA=\mI\), so it remains to show that \(\mA\mA^{-1}\) is also \(\mI\). Multiply \(\mA^{-1}\mA=\mI\) on the left by \(\mA\): \[ \mA^{-1}\mA=\mI\Rightarrow\mA\mA^{-1}\mA=\mA\Rightarrow(\mA\mA^{-1})\mA=\mA. \]

So \((\mA\mA^{-1})\mA\xx=\mA\xx\) for every \(\xx\). Since \(\mA\) has full rank, the vectors \(\mA\xx\) fill all of \(\R^n\), so \[ (\mA\mA^{-1})\zz=\zz\quad\text{for every }\zz\text{ in }\R^n. \] Hence \((\mA\mA^{-1}-\mI)\zz=\nul\) for every \(\zz\). Taking \(\zz=\ev_1,\ldots,\ev_n\) shows that all the entries of \(\mA\mA^{-1}-\mI\) are zero. Thus \(\mA\mA^{-1}=\mI\). \(\blacksquare\)

Definition  Negative powers are powers of the inverse: \(\mA^{-n}\equiv(\mA^{-1})^n=\underbrace{(\mA^{-1})\cdots(\mA^{-1})}_{n\text{ times}}\).

Theorem Algebra of inverse matrices:

Inverses combine with the other operations as follows:

Proof

1. \(\mA^{-1}\mA=\mI\) and \(\mA\mA^{-1}=\mI\), so \(\mA\) is the inverse of \(\mA^{-1}\).

2. This follows from (c) and the definition of the power: \[ \mA^n=\underbrace{\mA\mA\cdots\mA}_{n\text{ times}}\Rightarrow(\mA^n)^{-1} =[\mA\mA\cdots\mA]^{-1}=\underbrace{(\mA^{-1})(\mA^{-1})\cdots(\mA^{-1})}_{n\text{ times}}\equiv\mA^{-n}=(\mA^{-1})^n. \]

3. We check both products with the claimed inverse \(\mB^{-1}\mA^{-1}\): \begin{align*} \mA\mB(\mB^{-1}\mA^{-1})&=\mA(\mB\mB^{-1})\mA^{-1}=\mA\mI\mA^{-1}=\mA\mA^{-1}=\mI,\\ (\mB^{-1}\mA^{-1})\mA\mB&=\mB^{-1}(\mA^{-1}\mA)\mB=\mB^{-1}\mI\mB=\mB^{-1}\mB=\mI. \end{align*}

4. \(\mA^T(\mA^{-1})^T=(\mA^{-1}\mA)^T=\mI^T=\mI\) and \((\mA^{-1})^T\mA^T=(\mA\mA^{-1})^T=\mI^T=\mI\).

Theorem  Inverses solve linear systems directly.

Proof  \(\mA(\mA^{-1}\vc{c})=(\mA\mA^{-1})\vc{c}=\vc{c}\), so \(\mA^{-1}\vc{c}\) is a solution. It is the only one, because \(\mA\) has full column rank (it is invertible), so all \(n\) columns are pivot columns. \(\blacksquare\)

Determinants and Invertibility

The determinant gives a quick test for invertibility.

Example For \(\mA=\col{1&2\\3&4}\) we have \(\det\mA=4-6=-2\ne0\), so \(\mA\) is invertible.

Proof  Let \(\mE\) be the echelon form obtained from \(\mA\) by Gauss elimination. Row interchanges and adding multiples of one row to another change the determinant at most by a sign, so \(\det\mA=\pm\det\mE\). In particular, \(\det\mA\ne0\) exactly when \(\det\mE\ne0\). Since \(\mE\) is triangular, \(\det\mE\) is the product of its diagonal entries, the pivots. Now we prove both directions.

1. If \(\mA\) is invertible, it has full rank, so \(\mE\) has full rank: there are \(n\) pivots, all non-zero. Then \(\det(\mE)\ne0\), and so \(\det(\mA)\ne0\).

2. If \(\det\mA\ne0\), then \(\det(\mE)\ne0\), so all the pivots are non-zero and \(\mE\) has full rank. Hence \(\mA\) has full rank, and \(\mA\) is invertible.

Theorem Cramer’s Rule

For a square system with \(\det\mA\ne0\) there is an explicit formula for each unknown, in terms of determinants. Consider the \(n\times n\) system \(\mA\xx=\bv\) with \(\mA=\begin{bmatrix}\av_1&\cdots&\av_n\end{bmatrix}\) written by its columns.

Proof  Since \(\det\mA\ne0\), the matrix \(\mA\) is invertible, so the system has a unique solution \(\xx=\mA^{-1}\bv\). We now find the entries of \(\xx\).

Let \(\xx\) be the unique solution. Writing the product \(\mA\xx\) as a combination of the columns of \(\mA\) shows that this combination equals \(\bv\): \[ \mA\xx=\begin{bmatrix}\av_1&\cdots&\av_n\end{bmatrix}\col{x_1\\\vdots\\x_n} =\underbrace{\av_1x_1+\cdots+\av_nx_n}_{=\ \bv}=\sum_{k=1}^n\av_kx_k. \] Now replace the \(i^{\mathrm{th}}\) column of \(\mA\) by \(\bv\), and substitute this combination for \(\bv\). Because the determinant is additive and homogeneous in a single column (D9 and D7), it splits into a sum of \(n\) determinants. In every term with \(k\ne i\) the column \(\av_k\) appears twice, so that determinant is zero (D5). Only the term \(k=i\) survives, and it is \(x_i\det\mA\): \begin{align*} &\det\begin{bmatrix}\av_1&\cdots&\av_{i-1}&\bv&\av_{i+1}&\cdots&\av_n\end{bmatrix}\\ &=\det\begin{bmatrix}\av_1&\cdots&\av_{i-1}&\sum_{k=1}^n\av_kx_k&\av_{i+1}&\cdots&\av_n\end{bmatrix}\\ &=\sum_{k=1}^n x_k\underbrace{\det\begin{bmatrix}\av_1&\cdots&\av_{i-1}&\av_k&\av_{i+1}&\cdots&\av_n\end{bmatrix}}_{0\text{ if }k\ne i;\ \det\mA\text{ if }k=i}. \end{align*} Dividing by \(\det\mA\) gives the formula: \[ \Rightarrow\det\begin{bmatrix}\av_1&\cdots&\av_{i-1}&\bv&\av_{i+1}&\cdots&\av_n\end{bmatrix}=x_i\det\mA \quad\Rightarrow\quad x_i=\frac{\det\begin{bmatrix}\av_1&\cdots&\av_{i-1}&\bv&\av_{i+1}&\cdots&\av_n\end{bmatrix}}{\det\mA}.\quad\blacksquare \] Cramer’s rule gives a closed-form expression for \(\xx\), which is useful in analytical investigations.

Example Use Cramer’s rule to solve \[ \begin{aligned}x_1+4x_2+5x_3&=2\\4x_1+2x_2+5x_3&=3\\-3x_1+3x_2-x_3&=1.\end{aligned} \] First compute \(\det\mA\) by reducing to triangular form: \[ \det\mA=\begin{vmatrix}1&4&5\\4&2&5\\-3&3&-1\end{vmatrix}=\det\col{1&4&5\\0&-14&-15\\0&15&14}=-14\cdot14+15\cdot15=29. \] Then each unknown is a determinant, with one column replaced by the right-hand side, divided by \(29\): \[ x_1=\frac{\det\col{2&4&5\\3&2&5\\1&3&-1}}{29}=\frac{33}{29},\qquad x_2=\frac{\det\col{1&2&5\\4&3&5\\-3&1&-1}}{29}=\frac{35}{29},\qquad x_3=\frac{\det\col{1&4&2\\4&2&3\\-3&3&1}}{29}=-\frac{23}{29}. \]

Inverses and Adjoint Matrix

The inverse of a square, full-rank matrix \(\mA\) can be found by solving \(\mA\mX=\mI\). Written column by column, this means solving \(n\) linear systems, one for each column of the identity: \[ \mA\mX=\mI\quad\text{or}\quad\begin{bmatrix}\av_1&\cdots&\av_n\end{bmatrix}\begin{bmatrix}\xx_1&\cdots&\xx_n\end{bmatrix}=\begin{bmatrix}\ev_1&\cdots&\ev_n\end{bmatrix}, \qquad\mA\xx_j=\ev_j,\quad1\leq j\leq n. \]

By Cramer’s rule, the \(i^{\mathrm{th}}\) entry of \(\xx_j\) is a determinant in which column \(i\) is replaced by \(\ev_j\). Expanding that determinant along the replaced column leaves a single cofactor, \(C_{ji}\): \[ x_{ij}=\frac{\left|\av_1\ \cdots\ \av_{i-1}\ \ev_j\ \av_{i+1}\ \cdots\ \av_n\right|}{\det\mA}=\frac{C_{ji}}{\det\mA}, \] where \(C_{ji}\) is the \((j,i)\) cofactor of \(\mA\). Collecting these entries gives the inverse as the transposed matrix of cofactors, called the adjoint (or adjugate) of \(\mA\), divided by the determinant:

Example Take again \(\mA=\col{1&4&5\\4&2&5\\-3&3&-1}\) with \(\det\mA=29\). We solve \(\mA\begin{bmatrix}\xx_1&\xx_2&\xx_3\end{bmatrix}=\begin{bmatrix}\ev_1&\ev_2&\ev_3\end{bmatrix}\), where \[ \xx_1=\col{x_{11}\\x_{21}\\x_{31}},\quad\xx_2=\col{x_{12}\\x_{22}\\x_{32}},\quad\xx_3=\col{x_{13}\\x_{23}\\x_{33}}. \] By Cramer’s rule, the entries of the first column \(\xx_1\) are \begin{align*} x_{11}&=\frac{\det\col{1&4&5\\0&2&5\\0&3&-1}}{\det\mA}=\frac{C_{11}}{\det\mA},\\ x_{21}&=\frac{\det\col{1&1&5\\4&0&5\\-3&0&-1}}{\det\mA}=\frac{C_{12}}{\det\mA},\\ x_{31}&=\frac{\det\col{1&4&1\\4&2&0\\-3&3&0}}{\det\mA}=\frac{C_{13}}{\det\mA}, \end{align*} and the other columns follow in the same way. Together, the inverse is the transposed cofactor matrix over \(\det\mA\): \begin{align*} \mX=\mA^{-1}&=\frac1{\det\mA}\col{C_{11}&C_{21}&C_{31}\\C_{12}&C_{22}&C_{32}\\C_{13}&C_{23}&C_{33}}\\ &=\frac1{\det\mA}\underbrace{\col{C_{11}&C_{12}&C_{13}\\C_{21}&C_{22}&C_{23}\\C_{31}&C_{32}&C_{33}}^T}_{\operatorname{adj}\mA}. \end{align*} A more practical technique for manual computation of \(\mA^{-1}\): the adjoint formula needs many determinants, which is slow by hand. A much faster way is to solve all \(n\) systems \(\mA\mA^{-1}=\mI\) at once.

Let \(\mA^{-1}=\begin{bmatrix}\xx_1&\cdots&\xx_n\end{bmatrix}\), where the columns \(\xx_i\) are the unknowns. Then \[ \mA\mA^{-1}=\mA\begin{bmatrix}\xx_1&\cdots&\xx_n\end{bmatrix}=\mI=\begin{bmatrix}\ev_1&\cdots&\ev_n\end{bmatrix}, \] where \(\ev_i\) are the standard basis vectors. So \(\mA\xx_i=\ev_i\) for each \(i\). These could be solved one by one, but they can all be solved simultaneously by putting the identity next to \(\mA\) and row reducing:

Example Find the inverse of \(\mA=\col{1&3&0\\-2&3&1\\0&1&1}\). \begin{gather*} [\mA\ \mI]=\aug{1&3&0&1&0&0\\-2&3&1&0&1&0\\0&1&1&0&0&1}\\ \xrightarrow{R_2\to R_2+2R_1}\aug{1&3&0&1&0&0\\0&9&1&2&1&0\\0&1&1&0&0&1} \xrightarrow{R_2\to R_2/9}\aug{1&3&0&1&0&0\\0&1&1/9&2/9&1/9&0\\0&1&1&0&0&1}\\ \xrightarrow{R_3\to R_3-R_2}\aug{1&3&0&1&0&0\\0&1&1/9&2/9&1/9&0\\0&0&8/9&-2/9&-1/9&1} \xrightarrow{R_3\to R_3(9/8)}\aug{1&3&0&1&0&0\\0&1&1/9&2/9&1/9&0\\0&0&1&-1/4&-1/8&9/8}\\ \xrightarrow{R_2\to R_2-\frac19R_3}\aug{1&3&0&1&0&0\\0&1&0&1/4&1/8&-1/8\\0&0&1&-1/4&-1/8&9/8}\\ \xrightarrow{R_1\to R_1-3R_2}\left[\begin{array}{ccc|ccc}1&0&0&1/4&-3/8&3/8\\0&1&0&1/4&1/8&-1/8\\0&0&1&-1/4&-1/8&9/8\end{array}\right]. \end{gather*} The three columns on the right of the bar are \(\mA^{-1}\).

LU-Factorization

LU factorization is another way to solve \(\mA\xx=\vc{c}\) for a square matrix \(\mA\). The idea is to write \(\mA\) as a product of two simple matrices, one lower triangular and one upper triangular. If Gauss elimination of \(\mA\) needs no row interchanges, then

If row interchanges are needed, one instead factorizes \(\mP\mA=\mL\mU\) with a permutation matrix \(\mP\).

The factorization is not unique. To see why, count unknowns. A \(3\times3\) matrix \(\mA\) has \(9\) entries. The two triangular factors have \(6\) parameters each: \[ \mA=\underbrace{\col{l_{11}&0&0\\l_{21}&l_{22}&0\\l_{31}&l_{32}&l_{33}}}_{6\text{ parameters}}\underbrace{\col{u_{11}&u_{12}&u_{13}\\0&u_{22}&u_{23}\\0&0&u_{33}}}_{6\text{ parameters}}. \] That is \(12\) unknowns for only \(9\) equations (for an \(n\times n\) matrix, \(n^2+n\) unknowns). In Doolittle’s method we fix the freedom by setting \(l_{11}=l_{22}=l_{33}=1\) and solving for the remaining \(9\) parameters. For an \(n\times n\) matrix, \(l_{11}=\cdots=l_{nn}=1\) and there are \(n^2\) remaining parameters.

Example Factorize \(\mA=\col{2&-3&3\\6&-8&7\\-2&6&-1}\) with Doolittle’s method. Multiplying out the factors gives \begin{gather*} \underbrace{\col{2&-3&3\\6&-8&7\\-2&6&-1}}_{\mA} =\col{1&0&0\\l_{21}&1&0\\l_{31}&l_{32}&1}\col{u_{11}&u_{12}&u_{13}\\0&u_{22}&u_{23}\\0&0&u_{33}}\\ =\col{u_{11}&u_{12}&u_{13}\\l_{21}u_{11}&l_{21}u_{12}+u_{22}&l_{21}u_{13}+u_{23}\\l_{31}u_{11}&l_{31}u_{12}+l_{32}u_{22}&l_{31}u_{13}+l_{32}u_{23}+u_{33}}. \end{gather*} Comparing entries with \(\mA\), row by row, each unknown can be found from the ones already known: \begin{align*} 2&=u_{11},\quad-3=u_{12},\quad3=u_{13},\\ 6&=l_{21}u_{11}\Rightarrow l_{21}=6/2=3;\quad-8=l_{21}u_{12}+u_{22}\\ &\Rightarrow u_{22}=-8-3\cdot(-3)=1,\\ 7&=l_{21}u_{13}+u_{23}\Rightarrow u_{23}=7-3\cdot3=-2,\\ -2&=l_{31}u_{11}\Rightarrow l_{31}=-1;\quad6=l_{31}u_{12}+l_{32}u_{22}\Rightarrow l_{32}=3,\\ -1&=l_{31}u_{13}+l_{32}u_{23}+u_{33}\Rightarrow u_{33}=-1-(-1)\cdot3-3\cdot(-2)=8. \end{align*} \[\Rightarrow\mA=\col{1&0&0\\3&1&0\\-1&3&1}\col{2&-3&3\\0&1&-2\\0&0&8}.\] Once \(\mA=\mL\mU\) is known, the system \(\mA\xx=\vc{c}\) is solved in two easy steps. Name the product \(\mU\xx=\yy\). First solve \(\mL\yy=\vc{c}\) for \(\yy\) by forward substitution, then solve \(\mU\xx=\yy\) for \(\xx\) by backward substitution: \[ \mL\underbrace{\mU\xx}_{\yy}=\vc{c}\Rightarrow\mL\yy=\vc{c}\ \text{(forward substitution)},\qquad\mU\xx=\yy\ \text{(backward substitution)}. \] This needs fewer operations than Gauss–Jordan: about \(n^3/3\) multiplications for the factorization, versus about \(n^3/2\) for Gauss–Jordan, and the factors \(\mL\) and \(\mU\) can be reused for every new right-hand side \(\vc{c}\).