Exploiting matrix structure
Contents
2.9. Exploiting matrix structure#
A common situation in computation is that a problem has certain properties or structure that can be used to get a faster or more accurate solution. There are many properties of a matrix that can affect LU factorization. For example, an \(n \times n\) matrix \(A\) is diagonally dominant if
It turns out that a diagonally dominant matrix is guaranteed to be invertible, and row pivoting is not required for stability.
We next consider three important types of matrices that cause the LU factorization to be specialized in some important way.
Banded matrices#
A matrix \(\mathbf{A}\) has upper bandwidth \(b_u\) if \(j-i > b_u\) implies \(A_{ij}=0\), and lower bandwidth \(b_\ell\) if \(i-j > b_\ell\) implies \(A_{ij}=0\). We say the total bandwidth is \(b_u+b_\ell+1\). When \(b_u=b_\ell=1\), we have the important case of a tridiagonal matrix.
A = diagm( -1=>[4,3,2,1,0], 0=>[2,2,0,2,1,2], 1=>fill(-1,5) )
6×6 Matrix{Int64}:
2 -1 0 0 0 0
4 2 -1 0 0 0
0 3 0 -1 0 0
0 0 2 2 -1 0
0 0 0 1 1 -1
0 0 0 0 0 2
@show diag_main = diag(A);
@show diag_minusone = diag(A,-1);
diag_main = diag(A) = [2, 2, 0, 2, 1, 2]
diag_minusone = diag(A, -1) = [4, 3, 2, 1, 0]
The lower and upper bandwidths of \(\mathbf{A}\) are repeated in the factors from the unpivoted LU factorization.
L,U = FNC.lufact(A)
L
6×6 LowerTriangular{Float64, Matrix{Float64}}:
1.0 ⋅ ⋅ ⋅ ⋅ ⋅
2.0 1.0 ⋅ ⋅ ⋅ ⋅
0.0 0.75 1.0 ⋅ ⋅ ⋅
0.0 0.0 2.66667 1.0 ⋅ ⋅
0.0 0.0 0.0 0.214286 1.0 ⋅
0.0 0.0 0.0 0.0 0.0 1.0
U
6×6 UpperTriangular{Float64, Matrix{Float64}}:
2.0 -1.0 0.0 0.0 0.0 0.0
⋅ 4.0 -1.0 0.0 0.0 0.0
⋅ ⋅ 0.75 -1.0 0.0 0.0
⋅ ⋅ ⋅ 4.66667 -1.0 0.0
⋅ ⋅ ⋅ ⋅ 1.21429 -1.0
⋅ ⋅ ⋅ ⋅ ⋅ 2.0
If row pivoting is not used, the \(\mathbf{L}\) and \(\mathbf{U}\) factors preserve the lower and upper bandwidths of \(\mathbf{A}\). This fact implies computational savings in both the factorization and the triangular substitutions because the zeros appear predictably and we can skip multiplication and addition with them.
The number of flops needed by LU factorization without pivoting is \(O(b_u b_\ell n)\) when the upper and lower bandwidths are \(b_u\) and \(b_\ell\).
In order to exploit the savings offered by sparsity, we would need to make modifications to Function 2.4.6 and the triangular substitution routines. Alternatively, we can get Julia to take advantage of the structure automatically by converting the matrix into a special type called sparse. Sparse matrices are covered in more detail in Chapters 7 and 8.
If we use an ordinary or dense matrix, then there’s no way to exploit a banded structure such as tridiagonality.
n = 10000
A = diagm(0=>1:n, 1=>n-1:-1:1, -1=>ones(n-1))
lu(rand(3,3)) # force compilation
@time lu(A);
8.652754 seconds (5 allocations: 763.016 MiB, 0.23% gc time)
If instead we construct a proper sparse matrix, though, the speedup can be dramatic.
A = spdiagm(0=>1:n, 1=>n-1:-1:1, -1=>ones(n-1))
lu(A); # compile for sparse case
@time lu(A);
0.005833 seconds (65 allocations: 9.766 MiB)
You can also see above that far less memory was used in the sparse case.
Symmetric matrices#
A square matrix \(\mathbf{A}\) satisfying \(\mathbf{A}^T = \mathbf{A}\) is called symmetric.
Symmetric matrices arise frequently in applications because many types of interactions, such as gravitation and social-network befriending, are inherently symmetric. Symmetry in linear algebra simplifies many properties and algorithms. As a rule of thumb, if your matrix has symmetry, you want to exploit and preserve it.
In \(\mathbf{A}=\mathbf{L}\mathbf{U}\) we arbitrarily required the diagonal elements of \(\mathbf{L}\), but not \(\mathbf{U}\), to be one. That breaks symmetry, so we need to modify the goal to
where \(\mathbf{L}\) is unit lower triangular and \(\mathbf{D}\) is diagonal. To find an algorithm for this factorization, we begin by generalizing (2.4.2) a bit without furnishing proof.
Let \(\mathbf{D}\) be an \(n\times n\) diagonal matrix with diagonal elements \(d_1,d_2,\ldots,d_n\), and suppose \(\mathbf{A}\) and \(\mathbf{B}\) are \(n\times n\) as well. Write the columns of \(\mathbf{A}\) as \(\mathbf{a}_1,\dots,\mathbf{a}_n\) and the rows of \(\mathbf{B}\) as \(\mathbf{b}_1^T,\dots,\mathbf{b}_n^T\). Then
Let’s derive the LDL\(^T\) factorization for a small example.
We begin with a symmetric \(\mathbf{A}\).
A₁ = [ 2 4 4 2
4 5 8 -5
4 8 6 2
2 -5 2 -26 ];
We won’t use pivoting, so the pivot element is at position (1,1). This will become the first element on the diagonal of \(\mathbf{D}\). Then we divide by that pivot to get the first column of \(\mathbf{L}\).
L = diagm(ones(4))
d = zeros(4)
d[1] = A₁[1,1]
L[:,1] = A₁[:,1]/d[1]
A₂ = A₁ - d[1]*L[:,1]*L[:,1]'
4×4 Matrix{Float64}:
0.0 0.0 0.0 0.0
0.0 -3.0 0.0 -9.0
0.0 0.0 -2.0 -2.0
0.0 -9.0 -2.0 -28.0
We are now set up the same way for the submatrix in rows and columns 2–4.
d[2] = A₂[2,2]
L[:,2] = A₂[:,2]/d[2]
A₃ = A₂ - d[2]*L[:,2]*L[:,2]'
4×4 Matrix{Float64}:
0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0
0.0 0.0 -2.0 -2.0
0.0 0.0 -2.0 -1.0
We continue working our way down the diagonal.
d[3] = A₃[3,3]
L[:,3] = A₃[:,3]/d[3]
A₄ = A₃ - d[3]*L[:,3]*L[:,3]'
d[4] = A₄[4,4]
@show d;
L
d = [2.0, -3.0, -2.0, 1.0]
4×4 Matrix{Float64}:
1.0 -0.0 -0.0 0.0
2.0 1.0 -0.0 0.0
2.0 -0.0 1.0 0.0
1.0 3.0 1.0 1.0
We have arrived at the targeted factorization.
A₁ - L*diagm(d)*L'
4×4 Matrix{Float64}:
0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0
In practice we don’t actually have to carry out any arithmetic in the upper triangle of \(\mathbf{A}\) as we work, since the operations are always the mirror image of those in the lower triangle. As a result, it can be shown that LDL\(^T\) factorization takes about half as much work as the standard LU.
LDL\(^T\) factorization on an \(n \times n\) symmetric matrix, when successful, takes \(\sim \frac{1}{3}n^3\) flops.
Just as pivoting is necessary to stabilize LU factorization, the LDL\(^T\) factorization without pivoting may be unstable or even fail to exist. We won’t go into the details, because our interest is in specializing the factorization to matrices that also possess another important property.
Symmetric positive definite matrices#
Suppose that \(\mathbf{A}\) is \(n\times n\) and \(\mathbf{x}\in\mathbb{R}^n\). Observe that \(\mathbf{x}^T\mathbf{A}\mathbf{x}\) is the product of \(1\times n\), \(n\times n\), and \(n\times 1\) matrices, so it is a scalar, sometimes referred to as a quadratic form. It can be expressed as
A real \(n\times n\) matrix \(\mathbf{A}\) is called a symmetric positive definite matrix (or SPD matrix) if it is symmetric and, for all nonzero \(\mathbf{x}\in\mathbb{R}^n\),
The definiteness property is usually difficult to check directly from the definition. There are some equivalent conditions, though. For instance, a symmetric matrix is positive definite if and only if its eigenvalues are all real positive numbers. SPD matrices have important properties and appear in applications in which the definiteness is known for theoretical reasons.
Let us consider what definiteness means to the LDL\(^T\) factorization. We compute
where \(\mathbf{z}=\mathbf{L}^T \mathbf{x}\). Note that since \(\mathbf{L}\) is unit lower triangular, it is nonsingular, so \(\mathbf{x}=\mathbf{L}^{-T}\mathbf{z}\). By taking \(\mathbf{z}=\mathbf{e}_k\) for \(k=1,\ldots,n\), we can read the equalities from right to left and conclude that \(D_{kk}>0\) for all \(k\). That permits us to write a kind of square root formula:1
Now we have \(\mathbf{A}=\mathbf{L}\mathbf{D}^{1/2}\mathbf{D}^{1/2}\mathbf{L}^T= \mathbf{R}^T \mathbf{R}\), where \(\mathbf{R} =\mathbf{D}^{1/2}\mathbf{L}^T\) is an upper triangular matrix whose diagonal entries are positive.
Any SPD matrix \(\mathbf{A}\) may be factored as
where \(\mathbf{R}\) is an upper triangular matrix with positive diagonal elements. This is called the Cholesky factorization.
While the unpivoted LDL\(^T\) factorization is not stable and not even always possible, in the SPD case one can prove that pivoting is not necessary for the existence nor the stability of the Cholesky factorization.
Cholesky factorization of an \(n \times n\) SPD matrix takes \(\sim \frac{1}{3}n^3\) flops.
The speed and stability of the Cholesky factorization make it the top choice for solving linear systems with SPD matrices. As a side benefit, the Cholesky algorithm fails (in the form of an imaginary square root or division by zero) if and only if the matrix \(\mathbf{A}\) is not positive definite. This is often the best way to test the definiteness of a symmetric matrix about which nothing else is known.
A randomly chosen matrix is extremely unlikely to be symmetric. However, there is a simple way to symmetrize one.
A = rand(1.0:9.0,4,4)
B = A + A'
4×4 Matrix{Float64}:
8.0 14.0 9.0 8.0
14.0 12.0 8.0 10.0
9.0 8.0 6.0 14.0
8.0 10.0 14.0 14.0
cholesky(B) # throws an error
PosDefException: matrix is not positive definite; Cholesky factorization failed.
Stacktrace:
[1] checkpositivedefinite
@ /Applications/julia-rosetta/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/factorization.jl:18 [inlined]
[2] cholesky!(A::Hermitian{Float64, Matrix{Float64}}, ::Val{false}; check::Bool)
@ LinearAlgebra /Applications/julia-rosetta/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/cholesky.jl:266
[3] cholesky!(A::Matrix{Float64}, ::Val{false}; check::Bool)
@ LinearAlgebra /Applications/julia-rosetta/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/cholesky.jl:298
[4] #cholesky#143
@ /Applications/julia-rosetta/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/cholesky.jl:394 [inlined]
[5] cholesky (repeats 2 times)
@ /Applications/julia-rosetta/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/cholesky.jl:394 [inlined]
[6] top-level scope
@ In[14]:1
[7] eval
@ ./boot.jl:373 [inlined]
[8] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1196
It’s not hard to manufacture an SPD matrix to try out the Cholesky factorization.
B = A'*A
cf = cholesky(B)
Cholesky{Float64, Matrix{Float64}}
U factor:
4×4 UpperTriangular{Float64, Matrix{Float64}}:
11.0905 10.0085 7.57402 8.02486
⋅ 6.76973 5.78976 6.30497
⋅ ⋅ 6.86388 0.686969
⋅ ⋅ ⋅ 2.31886
What’s returned is a factorization object. Another step is required to extract the factor as a matrix:
R = cf.U
4×4 UpperTriangular{Float64, Matrix{Float64}}:
11.0905 10.0085 7.57402 8.02486
⋅ 6.76973 5.78976 6.30497
⋅ ⋅ 6.86388 0.686969
⋅ ⋅ ⋅ 2.31886
Here we validate the factorization:
opnorm(R'*R - B) / opnorm(B)
7.37383175345327e-17
Exercises#
✍ For each matrix, use (2.9.1) to determine whether it is diagonally dominant.
\[\begin{split}\mathbf{A} = \begin{bmatrix} 3 & 1 & 0 & 1 \\ 0 & -2 & 0 & 1 \\ -1 & 0 & 4 & -1 \\ 0 & 0 & 0 & 6 \end{bmatrix}, \quad \mathbf{B} = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{bmatrix}, \quad \mathbf{C} = \begin{bmatrix} 2 & -1 & 0 & 0 \\ -1 & 2 & -1 & 0 \\ 0 & -1 & 2 & -1 \\ 0 & 0 & -1 & 2 \end{bmatrix}.\end{split}\]⌨ For each matrix, use inspection or
choleskyin Julia to determine whether it is SPD.\[\begin{split}\mathbf{A} = \begin{bmatrix} 1 & 0 & -1 \\ 0 & 4 & 5 \\ -1 & 5 & 10 \end{bmatrix}, \qquad \mathbf{B} = \begin{bmatrix} 1 & 0 & 1 \\ 0 & 4 & 5 \\ -1 & 5 & 10 \end{bmatrix}, \qquad \mathbf{C} = \begin{bmatrix} 1 & 0 & 1 \\ 0 & 4 & 5 \\ 1 & 5 & 1 \end{bmatrix}.\end{split}\]✍ Show that the diagonal entries of a symmetric positive definite matrix are positive numbers. (Hint: Apply certain special cases of (2.9.4).)
⌨ Using Function 2.4.6 as a guide, write a function
function luband(A,upper,lower)
that accepts upper and lower bandwidth values and returns LU factors (without pivoting) in a way that avoids doing arithmetic using the locations that are known to stay zero. (Hint: Refer to the more efficient form of
lufactgiven in Section 2.5.)Test your function on the matrix with elements
\[\begin{split} A_{ij} = \begin{cases} \frac{1}{i+j}, & -1 \le i-j \le 2,\\ 0 & \text{otherwise.} \end{cases} \end{split}\]⌨ The
Tridiagonalmatrix type invokes a specialized algorithm for solving a linear system.(a) Set
n=1000andt=0. In a loop that runs 50 times, generate a linear system viaA = triu( tril(rand(n,n),1), -1) b = ones(n)
Using
@elapsed, incrementtby the time it takes to performA\b. Print out the final value oft.(b) Repeat the experiment of part (a), but generate the matrix via
A = Tridiagonal(rand(n,n))
What is the ratio of running times for part (a) and (b)?
(c) Now perform the experiment of part (b) for \(n=1000,1200,1400,\ldots,3000\), keeping the total time for each value of \(n\) in a vector. Plot running time as a function of \(n\) on a log-log scale. Is the time most like \(O(n)\), \(O(n^2)\), or \(O(n^3)\)? (If the answer is unclear, try increasing the number of solves per value of \(n\) to 100 or more.)
✍ Prove that if \(\mathbf{A}\) is any real invertible square matrix, then \(\mathbf{A}^T\mathbf{A}\) is SPD. (Hint: First show that \(\mathbf{x}^T\mathbf{A}^T\mathbf{A}\mathbf{x} \ge 0\) for all \(\mathbf{x}\). Then explain why zero is ruled out if \(\mathbf{x}\neq \boldsymbol{0}\).)
- 1
Except for this diagonal, positive definite case, it’s not trivial to define the square root of a matrix, so don’t generalize the notation used here.