4.3: LU factorization for square matrices

One system becomes two¶

The letters remind us of lower and upper. Throughout this chapter, we require ones on the diagonal of $L$ .

If $A=LU$ , then $A\vec x=\vec b$ becomes $L(U\vec x)=\vec b$ . Give the intermediate vector a name: $\vec y=U\vec x$ . Now we can work backward from $\vec b$ :

\vec x\xrightarrow{\ U\ }\vec y\xrightarrow{\ L\ }\vec b.

Although $U$ acts first on $\vec x$ , we solve the $L$ -system first. We know $\vec b$ and are recovering the vectors that led to it.

Finding the factors for the circuit¶

Return to the coefficient matrix from Chapter 4.1:

A=\begin{bmatrix}3&-1\\-1&2\end{bmatrix}.

Start with the required shapes and multiply:

\begin{bmatrix}1&0\\\ell&1\end{bmatrix} \begin{bmatrix}p&q\\0&r\end{bmatrix} =\begin{bmatrix}p&q\\\ell p&\ell q+r\end{bmatrix}.

Match this product to $A$ , one row and column at a time.

First, find the first row of $U$ . Matching the top row gives $p=3$ and $q=-1$ .

Next, find the first column of $L$ . Matching the bottom-left entry gives $3\ell=-1$ , so $\ell=-\frac13$ .

Finally, find the second row of $U$ . Matching the bottom-right entry gives

\left(-\frac13\right)(-1)+r=2\quad\Longrightarrow\quad r=\frac53.

Therefore,

\boxed{A=\underbrace{\begin{bmatrix}1&0\\-\frac13&1\end{bmatrix}}_L \underbrace{\begin{bmatrix}3&-1\\0&\frac53\end{bmatrix}}_U.}

Integer entries in $A$ can produce fractions in $L$ and $U$ . That’s fine! Multiplying the factors gives

LU=\begin{bmatrix}3&-1\\-1&\frac13+\frac53\end{bmatrix} =\begin{bmatrix}3&-1\\-1&2\end{bmatrix}=A.

Solving the circuit¶

Our original system was $A\vec V=\vec b$ , where $\vec b=\begin{bmatrix}15\\0\end{bmatrix}$ .

First solve $L\vec y=\vec b$ :

\begin{bmatrix}1&0\\-\frac13&1\end{bmatrix}\begin{bmatrix}y_1\\y_2\end{bmatrix} =\begin{bmatrix}15\\0\end{bmatrix}.

The first equation gives $y_1=15$ . The second gives $-\frac13(15)+y_2=0$ , so $y_2=5$ .

Next solve $U\vec V=\vec y$ :

\begin{bmatrix}3&-1\\0&\frac53\end{bmatrix}\begin{bmatrix}V_1\\V_2\end{bmatrix} =\begin{bmatrix}15\\5\end{bmatrix}.

Start at the bottom: $\frac53V_2=5$ , so $V_2=3$ . Then $3V_1-3=15$ , so $V_1=6$ .

\boxed{\vec V=\begin{bmatrix}6\\3\end{bmatrix}\text{ volts}.}

Does this agree with the physical model? The supply sends $(10-6)/2=2$ amperes into the first junction. Of these, $6/6=1$ ampere goes directly to ground and $(6-3)/3=1$ ampere goes to the second junction. From there, $3/3=1$ ampere flows to ground. Both junctions balance.

Reusing the same factors¶

For one small system, factoring may feel like extra work. Its value becomes clearer when we solve several systems with the same coefficient matrix.

If we change the supply voltage but keep the resistors, $A$ , $L$ , and $U$ stay the same. Only $\vec b$ changes. For a supply voltage $s$ , the circuit equations become

A\vec V=\begin{bmatrix}3s/2\\0\end{bmatrix}.

Doubling the supply from 10 to 20 volts doubles $\vec b$ . The two triangular solves then give $\vec V=\begin{bmatrix}12\\6\end{bmatrix}$ , without finding the factors again.

For large systems, finding the factors is usually the expensive part. Once we have them, each new right-hand side requires only two triangular solves.

Activity 1

The heat-flow system from Chapter 4.1 is

\begin{bmatrix}2&-1\\-1&2\end{bmatrix}\begin{bmatrix}T_1\\T_2\end{bmatrix} =\begin{bmatrix}0\\60\end{bmatrix}.

Find its LU factorization by matching entries. Find the first row of $U$ , then the unknown entry of $L$ , then the final entry of $U$ .
Solve $L\vec y=\vec b$ , then $U\vec T=\vec y$ .
Check the temperatures in the original heat-balance equations. Do they make physical sense?
Use the same factors when the right endpoint is held at $90^\circ\mathrm C$ instead.

Solution

1. Write

\begin{bmatrix}2&-1\\-1&2\end{bmatrix}=\begin{bmatrix}1&0\\\ell&1\end{bmatrix}\begin{bmatrix}p&q\\0&r\end{bmatrix}.

The first row gives $p=2$ , $q=-1$ . The bottom-left entry gives $2\ell=-1$ , so $\ell=-\frac12$ . The bottom-right entry gives $\frac12+r=2$ , so $r=\frac32$ . Therefore

L=\begin{bmatrix}1&0\\-\frac12&1\end{bmatrix},\qquad U=\begin{bmatrix}2&-1\\0&\frac32\end{bmatrix}.

2. Solving $L\vec y=\begin{bmatrix}0\\60\end{bmatrix}$ gives $y_1=0$ and $-\frac12y_1+y_2=60$ , so $\vec y=\begin{bmatrix}0\\60\end{bmatrix}$ . Then $U\vec T=\vec y$ gives

\frac32T_2=60\quad\Longrightarrow\quad T_2=40,\qquad2T_1-40=0\quad\Longrightarrow\quad T_1=20.

3. In the original heat-balance equations, $40-20=20-0$ and $60-40=40-20$ . Each segment carries 20 watts toward the left. The temperatures $0,20,40,60$ increase evenly from the cold endpoint to the hot endpoint, as expected.

4. With $\vec b=\begin{bmatrix}0\\90\end{bmatrix}$ , the same $L$ gives $\vec y=\begin{bmatrix}0\\90\end{bmatrix}$ . Then $\frac32T_2=90$ gives $T_2=60$ , and $2T_1-T_2=0$ gives $T_1=30$ . The temperatures are now $0,30,60,90$ degrees Celsius.

A three-by-three example¶

The same approach works for larger matrices. Let’s factor

A=\begin{bmatrix}4&-2&2\\2&2&0\\-4&8&1\end{bmatrix}.

Write the required shapes:

L=\begin{bmatrix}1&0&0\\\ell_{21}&1&0\\\ell_{31}&\ell_{32}&1\end{bmatrix},\qquad U=\begin{bmatrix}u_{11}&u_{12}&u_{13}\\0&u_{22}&u_{23}\\0&0&u_{33}\end{bmatrix}.

Step 1: First row of $U$ . The first row of $LU$ equals the first row of $U$ , so

u_{11}=4,\qquad u_{12}=-2,\qquad u_{13}=2.

Step 2: First column of $L$ . Matching entries $(2,1)$ and $(3,1)$ gives

4\ell_{21}=2,\qquad4\ell_{31}=-4, \quad\Longrightarrow\quad\ell_{21}=\frac12,\quad\ell_{31}=-1.

Step 3: Second row of $U$ . Matching entries $(2,2)$ and $(2,3)$ gives

\frac12(-2)+u_{22}=2,\qquad\frac12(2)+u_{23}=0, \quad\Longrightarrow\quad u_{22}=3,\quad u_{23}=-1.

Step 4: Second column of $L$ . Matching entry $(3,2)$ gives

(-1)(-2)+3\ell_{32}=8\quad\Longrightarrow\quad\ell_{32}=2.

Step 5: Final entry of $U$ . Matching entry $(3,3)$ gives

(-1)(2)+2(-1)+u_{33}=1\quad\Longrightarrow\quad u_{33}=5.

Thus,

\boxed{L=\begin{bmatrix}1&0&0\\\frac12&1&0\\-1&2&1\end{bmatrix},\qquad U=\begin{bmatrix}4&-2&2\\0&3&-1\\0&0&5\end{bmatrix}.}

Always check the product:

LU=\begin{bmatrix}4&-2&2\\2&-1+3&1-1\\-4&2+6&-2-2+5\end{bmatrix} =\begin{bmatrix}4&-2&2\\2&2&0\\-4&8&1\end{bmatrix}=A.

The diagonal entries of $U$ are called pivots in this factorization procedure. To find entries of $L$ below a pivot, we divide by that pivot. In the next chapter, we’ll see what can happen when a pivot is zero.

Activity 2

Find an LU factorization of $\begin{bmatrix}3&6&-3\\-4&-5&10\\5&8&-4\end{bmatrix}$ without reordering its rows. Check your product.
Use the worked example’s factors for $A=\begin{bmatrix}4&-2&2\\2&2&0\\-4&8&1\end{bmatrix}$ to solve $A\vec x=\begin{bmatrix}6\\6\\15\end{bmatrix}$ .
Explain why the intermediate vector $\vec y$ is generally different from the solution $\vec x$ .

Solution

1. For $A=\begin{bmatrix}3&6&-3\\-4&-5&10\\5&8&-4\end{bmatrix}$ , the first row of $U$ is $\begin{bmatrix}3&6&-3\end{bmatrix}$ . Matching the first column gives $\ell_{21}=-\frac43$ and $\ell_{31}=\frac53$ .

Next, $(-\frac43)6+u_{22}=-5$ gives $u_{22}=3$ , and $(-\frac43)(-3)+u_{23}=10$ gives $u_{23}=6$ . Then $(\frac53)6+3\ell_{32}=8$ gives $\ell_{32}=-\frac23$ . Finally, $(\frac53)(-3)+(-\frac23)6+u_{33}=-4$ gives $u_{33}=5$ .

Thus

L=\begin{bmatrix}1&0&0\\-\frac43&1&0\\\frac53&-\frac23&1\end{bmatrix},\qquad U=\begin{bmatrix}3&6&-3\\0&3&6\\0&0&5\end{bmatrix}.

The product checks:

LU=\begin{bmatrix}3&6&-3\\-4&-8+3&4+6\\5&10-2&-5-4+5\end{bmatrix}=A.

2. For the different matrix $A=\begin{bmatrix}4&-2&2\\2&2&0\\-4&8&1\end{bmatrix}$ , the factors are

L=\begin{bmatrix}1&0&0\\\frac12&1&0\\-1&2&1\end{bmatrix},\qquad U=\begin{bmatrix}4&-2&2\\0&3&-1\\0&0&5\end{bmatrix}.

Solve $L\vec y=\begin{bmatrix}6\\6\\15\end{bmatrix}$ : $y_1=6$ , $\frac12(6)+y_2=6$ gives $y_2=3$ , and $-6+2(3)+y_3=15$ gives $y_3=15$ .

Then $U\vec x=\begin{bmatrix}6\\3\\15\end{bmatrix}$ gives $x_3=3$ , $3x_2-3=3$ gives $x_2=2$ , and $4x_1-4+6=6$ gives $x_1=1$ .

\boxed{\vec x=\begin{bmatrix}1\\2\\3\end{bmatrix}.}

Indeed, $A\vec x=\begin{bmatrix}4-4+6\\2+4\\-4+16+3\end{bmatrix}=\begin{bmatrix}6\\6\\15\end{bmatrix}$ .

3. The intermediate vector satisfies $\vec y=U\vec x$ , so it is the result of applying $U$ to the solution. It need not equal $\vec x$ . Here $\vec y=\begin{bmatrix}6\\3\\15\end{bmatrix}$ while $\vec x=\begin{bmatrix}1\\2\\3\end{bmatrix}$ .