Row pivoting¶
Here is the system that “broke” LU factorization for us.
using FundamentalsNumericalComputation
A = [2 0 4 3 ; -4 5 -7 -10 ; 1 15 2 -4.5 ; -2 0 2 -13];
b = [ 4, 9, 29, 40 ];
When we use the lu function (from LinearAlgebra) with three outputs, we get the elements of the PLU factorization.
L,U,p = lu(A);
L
4×4 Array{Float64,2}:
1.0 0.0 0.0 0.0
-0.25 1.0 0.0 0.0
0.5 -0.153846 1.0 0.0
-0.5 0.153846 0.0833333 1.0
U
4×4 Array{Float64,2}:
-4.0 5.0 -7.0 -10.0
0.0 16.25 0.25 -7.0
0.0 0.0 5.53846 -9.07692
0.0 0.0 0.0 -0.166667
p
4-element Array{Int64,1}:
2
3
4
1
As you see above, the p return is a vector permutation of 1:n, rather than the permutation matrix P. We can recover the latter as follows:
P = I(4)[p,:]
4×4 SparseMatrixCSC{Bool,Int64} with 4 stored entries:
[4, 1] = 1
[1, 2] = 1
[2, 3] = 1
[3, 4] = 1
However, this is rarely necessary in practice (and the vector requires a lot less storage). We can the linear system, for example, using only p.
x = FNC.backsub( U, FNC.forwardsub(L,b[p,:]) )
4-element Array{Float64,1}:
-3.000000000000014
1.0000000000000009
4.000000000000004
-1.9999999999999973
If you call lu with just one output, it is a “factorization object”. You can access the individual parts of it using a dot syntax.
fact = lu(A)
fact.L
4×4 Array{Float64,2}:
1.0 0.0 0.0 0.0
-0.25 1.0 0.0 0.0
0.5 -0.153846 1.0 0.0
-0.5 0.153846 0.0833333 1.0
The factorization object can be used efficiently to solve linear systems by the backslash.
x = fact\b
4-element Array{Float64,1}:
-3.0000000000000138
1.000000000000001
4.000000000000004
-1.9999999999999973
The idea here is that if you have to solve many different linear systems for the same matrix, you can perform the computationally expensive factorization just once, and repeat only the much faster triangular solves for the different right-hand sides.