Exploiting matrix structure

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

(2.9.1)#\[\begin{split} |A_{ii}| > \sum_{\substack{j=1\\ j \neq i}}^{n} |A_{ij}| \hskip 0.25in \text{for each } i=1,\ldots,n.\end{split}\]

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#

Definition 2.9.1 : Bandwidth

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.

Demo 2.9.2

Here is a small tridiagonal matrix. Note that there is one fewer element on the super- and subdiagonals than on the main diagonal.

Use fill to create an array of a given size, with each element equal to a provided value.

A = diagm( -1=>[4,3,2,1,0], 0=>[2,2,0,2,1,2], 1=>fill(-1,5) )

6×6 Matrix{Int64}:
-1   0   0   0   0
 2  -1   0   0   0
 3   0  -1   0   0
 0   2   2  -1   0
 0   0   1   1  -1
 0   0   0   0   2

We can extract the elements on any diagonal using the diag function. The main or central diagonal is numbered zero, above and to the right of that is positive, and below and to the left is negative.

The diag function extracts the elements from a specified diagonal of a matrix.

@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}}:
0   ⋅     ⋅        ⋅         ⋅    ⋅ 
0  1.0    ⋅        ⋅         ⋅    ⋅ 
0  0.75  1.0       ⋅         ⋅    ⋅ 
0  0.0   2.66667  1.0        ⋅    ⋅ 
0  0.0   0.0      0.214286  1.0   ⋅ 
0  0.0   0.0      0.0       0.0  1.0

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.

Observation 2.9.3

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.

Demo 2.9.4

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);

  2.459274 seconds (7 allocations: 763.016 MiB, 0.25% 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.003595 seconds (86 allocations: 9.920 MiB)

You can also see above that far less memory was used in the sparse case.

Symmetric matrices#

Definition 2.9.5 : Symmetric matrix

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

\[ \mathbf{A}=\mathbf{L}\mathbf{D}\mathbf{L}^T, \]

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.

Theorem 2.9.6 : Linear combination of outer products

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

(2.9.2)#\[\mathbf{A}\mathbf{D}\mathbf{B} = \sum_{k=1}^n d_k \mathbf{a}_k \mathbf{b}_k^T.\]

Let’s derive the LDL\(^T\) factorization for a small example.

Demo 2.9.7

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  -3.0   0.0   -9.0
0   0.0  -2.0   -2.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  -2.0  -2.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}:
0  -0.0  -0.0  0.0
0   1.0  -0.0  0.0
0  -0.0   1.0  0.0
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

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.

Observation 2.9.8

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

(2.9.3)#\[ \mathbf{x}^T\mathbf{A}\mathbf{x} = \sum_{i=1}^n \sum_{j=1}^n A_{ij}x_ix_j.\]

Definition 2.9.9 : Symmetric positive definite matrix

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\),

(2.9.4)#\[ \mathbf{x}^T \mathbf{A} \mathbf{x} > 0.\]

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

\[ 0 < \mathbf{x}^T\mathbf{A}\mathbf{x} = \mathbf{x}^T \mathbf{L} \mathbf{D} \mathbf{L}^T \mathbf{x} = \mathbf{z}^T \mathbf{D} \mathbf{z},\]

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

(2.9.5)#\[\begin{split} \mathbf{D} = \begin{bmatrix} D_{11} & & & \\ & D_{22} & & \\ & & \ddots & \\ & & & D_{nn} \end{bmatrix} = \begin{bmatrix} \sqrt{D_{11}} & & & \\ & \sqrt{D_{22}} & & \\ & & \ddots & \\ & & & \sqrt{D_{nn}} \end{bmatrix}^{\,2} = \bigl( \mathbf{D}^{1/2} \bigr)^2.\end{split}\]

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.

Theorem 2.9.10 : Cholesky factorization

Any SPD matrix \(\mathbf{A}\) may be factored as

\[ \mathbf{A} = \mathbf{R}^T \mathbf{R}, \]

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.

Observation 2.9.11

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.

Demo 2.9.12

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}:
0  12.0  14.0  11.0
0  10.0  14.0  14.0
0  14.0  12.0   8.0
0  14.0   8.0  18.0

Similarly, a random symmetric matrix is unlikely to be positive definite. The Cholesky algorithm always detects a non-PD matrix by quitting with an error.

The cholesky function computes a Cholesky factorization if possible, or throws an error for a non-positive-definite matrix. However, it does not check for symmetry.

cholesky(B)    # throws an error

PosDefException: matrix is not positive definite; Factorization failed.

Stacktrace:
 [1] cholesky(A::Matrix{Float64}, ::NoPivot; check::Bool)
   @ LinearAlgebra ~/Documents/julia_fast.so:-1
 [2] cholesky
   @ ~/.julia/juliaup/julia-1.11.1+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/cholesky.jl:401 [inlined]
 [3] cholesky(A::Matrix{Float64})
   @ LinearAlgebra ~/.julia/juliaup/julia-1.11.1+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/cholesky.jl:401
 [4] top-level scope
   @ In[14]:1

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}}:
 13.4536  11.6697   9.96014  12.3387
   ⋅       3.71721  4.24187  -0.804083
   ⋅        ⋅       2.60808   2.49831
   ⋅        ⋅        ⋅        3.29679

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}}:
 13.4536  11.6697   9.96014  12.3387
   ⋅       3.71721  4.24187  -0.804083
   ⋅        ⋅       2.60808   2.49831
   ⋅        ⋅        ⋅        3.29679

Here we validate the factorization:

opnorm(R'*R - B) / opnorm(B)

0.0

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 cholesky in 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 lufact given 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 Tridiagonal matrix type invokes a specialized algorithm for solving a linear system.

(a) Set n=1000 and t=0. In a loop that runs 50 times, generate a linear system via
```
A = triu( tril(rand(n,n),1), -1)
b = ones(n)
```
Using @elapsed, increment t by the time it takes to perform A\b. Print out the final value of t.

(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.