CODEWITHSUNDEEP
pythonnumpymathscipylinear-algebra
Junhong Xu
Junhong Xu

Reputation: 145

The necessity of LU decomposition (using numpy as an example)

I am trying to understand the necessity of LU decomposition using numpy and scipy libraries. From what I understand is that we want to solve Ax = b, we first factorize A into two triangular matrices L and U then solve LUx = b by solving Ly = b then Ux = y. By solving triangular matrices, we can reduce time compare to Gaussian Elimination.

So, I tired this idea in python using numpy and scipy.

I first construct A and b using toy examples:

A = np.array([[2, 1, 0, 5], [1, 2, 1, 2], [0, 1, 2, 4], [1, 3, 6, 4.5]])
b = np.array([9, 10, -2, 3])

Then first solve this toy example in np.solve

%timeit np.linalg.solve(A, b )

The time is

9.76 µs ± 782 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Then I use factorization to solve this system:

lu, piv = linalg.lu_factor(A)
%timeit linalg.lu_solve((lu, piv), b)

I saw the output is

18.8 µs ± 213 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

, which is quiet slow compared to np.solve.

So, my question is, why np.solve is faster than linalg.lu_factor? My guess is that numpy.solve does not use Gaussian Elimination to solve the equations? A little confused with the result here.

Edit

Now, I use a much larger matrix to do the experiment (10000 x 10000).

Here is the result: for np.linalg.solve

8.64 s ± 180 ms per loop (mean ± std. dev. of 7 runs, 1 loop each);

for scipy.linalg.lu_solve

121 ms ± 3.79 ms per loop (mean ± std. dev. of 7 runs, 10 loops each).

For lu_solve, I only count the time on solving, the decomposition part is not counted. It is now much faster!

Upvotes: 0

Views: 2193

Answers (2)

Rory Daulton
Rory Daulton

Reputation: 22534

Here is a partial answer, since I dispute one of your premises.

You write that "LU solve should be faster than Gaussian Elimination." You seem to misunderstand the purpose of LU decomposition. If you are solving just one such problem (Ax=b where matrix A and vector b are given), LU decomp is no faster than Gaussian elimination. Indeed, the decomposition's algorithm is very similar to the elimination and is no faster.

The advantage of LU decomposition comes when you are given matrix A and you want to solve the equation Ax=b for multiple different given vectors b. Gaussian elimination needs to start over from scratch, and each solution will take the same amount of time. In LU decomposition you can store the resulting matrices L and U from the first calculation, and that greatly speeds up the solutions to the succeeding equations that use different vectors b.

You can read more about this at the section in Numerical Recipes in C about LU Decomposition and Its Applications.

Upvotes: 4

Warren Weckesser
Warren Weckesser

Reputation: 114781

Look at the docstring for numpy.linalg.solve. It says in the "Notes" section "The solutions are computed using LAPACK routine _gesv". (The underscore is a place-holder for a character that corresponds to a data type. For example, dgesv uses double precision.)

The documentation for dgesv explains that it uses the LU decomposition. So you are more-or-less replicating the calculation, but you are doing more steps in Python, so your code is slower.

Upvotes: 1

Related Questions