# Solve and lu_solve¶

Here's some Python code (using the NumPy and SciPy libraries) for solving linear systems.

In :
# Imports
import numpy as np
from scipy.linalg import lu, inv, solve, lu_solve, lu_factor, norm


## Small example¶

First, let's generate a random example.

In :
np.random.seed(1)
A = np.random.randint(-5,5,(4,4))
b = np.random.randint(-5,5,4)
print(A)
print(b)

[[ 0  3  4  0]
[-5 -5 -4  2]
[ 1  4 -3 -1]
[ 0 -3 -1 -3]]
[-1  2  2  4]


Here's the simple, straightforward way to solve the system $A\vec{x}=\vec{b}$:

In :
solve(A,b)

Out:
array([-0.77777778,  0.11111111, -0.33333333, -1.33333333])

The lu command provides a way to compute the LU decomposition.

In :
P,L,U = lu(A)
print(P)
print(L)
print(U)

[[0. 1. 0. 0.]
[1. 0. 0. 0.]
[0. 0. 1. 0.]
[0. 0. 0. 1.]]
[[ 1.          0.          0.          0.        ]
[-0.          1.          0.          0.        ]
[-0.2         1.          1.          0.        ]
[-0.         -1.         -0.38461538  1.        ]]
[[-5.         -5.         -4.          2.        ]
[ 0.          3.          4.          0.        ]
[ 0.          0.         -7.8        -0.6       ]
[ 0.          0.          0.         -3.23076923]]


The matrices $P$, $L$, $U$, and $A$ should satisfy $$LU = PA.$$

In :
L.dot(U)

Out:
array([[-5., -5., -4.,  2.],
[ 0.,  3.,  4.,  0.],
[ 1.,  4., -3., -1.],
[ 0., -3., -1., -3.]])

## Checking hand computations¶

Suppose we'd like to find the permuted LU Decomposition of that same matrix $$A = \begin{bmatrix} 0 & 3 & 4 & 0 \\ -5 & -5 & -4 & 2 \\ 1 & 4 & -3 & -1 \\ 0 & -3 & -1 & -3 \\ \end{bmatrix}$$ by hand. We can illustrate the steps but still have Python do the computations as follows.

First, it's evident that we need to swap rows 1 and 2 at the beginning. That brings us to

$$P = \begin{bmatrix} 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix}.$$

So we find ourselves applying Gaussian elimination to $$PA = \begin{bmatrix} -5 & -5 & -4 & 2 \\ 0 & 3 & 4 & 0 \\ 1 & 4 & -3 & -1 \\ 0 & -3 & -1 & -3 \\ \end{bmatrix}.$$ Or, in computer speak:

In :
P = np.array([
[0,1,0,0],
[1,0,0,0],
[0,0,1,0],
[0,0,0,1]
])
PA = P.dot(A)
PA

Out:
array([[-5, -5, -4,  2],
[ 0,  3,  4,  0],
[ 1,  4, -3, -1],
[ 0, -3, -1, -3]])

The first step in Gaussian elimination should be $\frac{1}{5}R_1 + R_3$. We can accomplish this via multiplication with the following elementary matrix:

In :
E1 = np.array([
[1,0,0,0],
[0,1,0,0],
[1/5,0,1,0],
[0,0,0,1]
])
E1.dot(PA)

Out:
array([[-5. , -5. , -4. ,  2. ],
[ 0. ,  3. ,  4. ,  0. ],
[ 0. ,  3. , -3.8, -0.6],
[ 0. , -3. , -1. , -3. ]])

Our next two steps are $-R2+R3$ and $R_2+R4$, which we can combine into a single step as follows:

In :
E32 = np.array([
[1, 0, 0, 0],
[0, 1, 0, 0],
[0,-1, 1, 0],
[0, 1, 0, 1]
])
E32.dot(E1).dot(PA)

Out:
array([[-5. , -5. , -4. ,  2. ],
[ 0. ,  3. ,  4. ,  0. ],
[ 0. ,  0. , -7.8, -0.6],
[ 0. ,  0. ,  3. , -3. ]])

Finally, we do $\frac{3}{7.8}R_3+R_4$:

In :
E4 = np.array([
[1, 0,  0,    0],
[0, 1,  0,    0],
[0, 0,  1,    0],
[0, 0, 3/7.8, 1]
])
U = E4.dot(E32.dot(E1).dot(PA))
U

Out:
array([[-5.        , -5.        , -4.        ,  2.        ],
[ 0.        ,  3.        ,  4.        ,  0.        ],
[ 0.        ,  0.        , -7.8       , -0.6       ],
[ 0.        ,  0.        ,  0.        , -3.23076923]])

There's $U$!

We could compute $L$ as the inverse of the product of the $E$s:

In :
L = inv(E4.dot(E32).dot(E1))
print(L)

[[ 1.          0.          0.         -0.        ]
[ 0.          1.          0.         -0.        ]
[-0.2         1.          1.         -0.        ]
[ 0.         -1.         -0.38461538  1.        ]]


Note how easy it is, though, to see $L$ in the row operations that we used in the first place.

## Timing bigger examples¶

Let's see how to use the LU decomposition to solve large systems and examine the time savings.

In :
n = 1000
A = np.random.randint(-50,50,(n,n))
b = np.random.randint(-50,50,n)

In :
%%timeit
solve(A,b)

17.7 ms ± 1.32 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)


SciPy provides an lu_solve command that allows us to solve a system expressed in terms of a compactified form of the LU decomposition. The compactified form is generated with lu_factor.

In :
%%timeit
lu_factor(A)

10.2 ms ± 487 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In :
LU,P = lu_factor(A)

In :
%%timeit
lu_solve((LU,P),b)

268 µs ± 5.44 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)


## Using the inverse¶

It takes much longer to compute an inverse than it does to compute an LU factorization - or even just to solve $A\vec{x} = \vec{b}$ directly in the first place.

In :
%%timeit
inv(A)

27.6 ms ± 364 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)


Once we have the inverse, though, it's pretty fast to use it to solve systems.

In :
AInv = inv(A)

In :
%%timeit
AInv.dot(b)

132 µs ± 12.4 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)


However, using the inverse is usually not as accurate as Gaussian elimination.

In :
sol1 = AInv.dot(b)

In :
norm(A.dot(sol1) - b)

Out:
1.2979574249688793e-09
In :
sol2 = lu_solve((LU,P),b)

In :
norm(A.dot(sol2) - b)

Out:
7.05154876210495e-10