Speed of LU decomposition vs traditional Ax=b

Question

Learning about LU decompositions, and the textbook I'm using claimed that once the initial matrices L and U have been computed (assuming no pivoting matrix P is needed), it is much faster/more efficient (aka less flops) to solve Ly = b and then Ux = y than Ax = b from scratch. However, when I ran it on my system with a random A ⊆ R^(10 x 10) matrix and b ⊆ R^(10 x 1), it ended up being much slower (12.49 seconds vs 6.17 seconds).

Here's my code:

import numpy as np
from scipy.linalg import lu
from numpy.random import rand
from numpy.linalg import solve
from timeit import timeit as duration

A, b = rand(500, 500), rand(500, 1)
P,L,U = lu(A)

t_vanilla = duration("solve(A,b)", setup="from __main__ import solve, A, b")
print(t_vanilla)
t_decomps = duration("solve(U, solve(L,b))", setup="from __main__ import solve, L, U, b")
print(t_decomps)

Any ideas on why this might be the case? This is my first time using the timeit library, so I may well have messed something up haha. I'm running this on Mac OSX, python 3.8.

Also, I noticed that scipy.solve was WAY slower (over double the time) than numpy.solve for some reason ... any intuition as to why?

Thanks!!!

user_10692380 · Answer

When running solve(U, solve(L,b)), you just give to LAPACK gesv routine behind np.linalg.solve two distinct linear systems that it blindly solves without accounting for the triangular structure of L and U. This is done by performing respectively for each a LU factorization followed by forward, Ly = b , and backward, Ux = y, substitution. That is why your timings seem to indicate roughly twice the time as a single solve call.

Now according to your textbook, and more generally, the numerical efficiency that comes from already having the LU decomposition of matrix A, is that you can simply reuse it if you have multiple right-hand side b in order to avoid a costly and redundant factorization step. In order to do so, you only need to solve the two triangular systems Ly = b by forward substitution then Ux = y using backward substitution. Let me time these as follows:

import numpy as np
from scipy.linalg import lu_factor, lu_solve, lu

A, b = np.random.rand(500, 500), np.random.rand(500, 1)
P,L,U = lu(A)
lu_A = lu_factor(A)

# LU + forward backward
%timeit np.linalg.solve(A, b)
>>> 1.4 ms ± 10.6 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

# 2 x (LU + forward backward)
%timeit np.linalg.solve(U, np.linalg.solve(L, b))
>>> 2.39 ms ± 85.4 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

# forward backward
%timeit lu_solve(lu_A, b)
>>> 69.8 µs ± 6.48 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

and it is indeed far more efficient if you have multiple b to solve for.

Note that I used scipy.linalg.lu_factor, which stores a compact form of the LU factorization, followed by scipy.linalg.lu_solve, which solves the system with the decomposition given by lu_factor.

Note also that reusing the LU factorization is useful when these multiple b are not available at the same time. On the contrary if they are readily available, you could simply make a single solve call with horizontally stacked b with equivalent performance as the solver itself will reuse the LU factorization computed for the first one for every remaining b.

Concerning the performances issues with scipy it's rather hard to say without more informations about numpy/scipy (np.show_config()).

Hope it is clearer now.

Speed of LU decomposition vs traditional Ax=b

1 Answers