Solve and lu_solve

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

In [1]:
# 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 [2]:
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 [3]:
solve(A,b)
Out[3]:
array([-0.77777778,  0.11111111, -0.33333333, -1.33333333])

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

In [4]:
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 [5]:
L.dot(U)
Out[5]:
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 [6]:
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[6]:
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 [7]:
E1 = np.array([
    [1,0,0,0],
    [0,1,0,0],
    [1/5,0,1,0],
    [0,0,0,1]
])
E1.dot(PA)
Out[7]:
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 [8]:
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[8]:
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 [9]:
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[9]:
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 [10]:
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 [11]:
n = 1000
A = np.random.randint(-50,50,(n,n))
b = np.random.randint(-50,50,n)
In [12]:
%%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 [13]:
%%timeit
lu_factor(A)

10.2 ms ± 487 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
In [14]:
LU,P = lu_factor(A)
In [15]:
%%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 [16]:
%%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 [17]:
AInv = inv(A)
In [18]:
%%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 [19]:
sol1 = AInv.dot(b)
In [20]:
norm(A.dot(sol1) - b)
Out[20]:
1.2979574249688793e-09
In [21]:
sol2 = lu_solve((LU,P),b)
In [22]:
norm(A.dot(sol2) - b)
Out[22]:
7.05154876210495e-10