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.

2.9.1Banded matrices¶

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.

In order to exploit the savings offered by sparsity, we would need to make modifications to Function 2.4.1 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.

2.9.2Symmetric matrices¶

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.4) a bit without furnishing proof.

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

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.

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.

2.9.3Symmetric 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

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.

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.

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.

2.9.4Exercises¶

1. ✍ For each matrix, use (2.9.1) to determine whether it is diagonally dominant.

   $$\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}.$$

   (2.9.9)

2. ⌨ For each matrix, use inspection or cholesky in Julia to determine whether it is SPD.

   $$\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}.$$

   (2.9.10)

3. ✍ Show that the diagonal entries of a symmetric positive definite matrix are positive numbers. (Hint: Apply certain special cases of (2.9.5).)

4. ⌨ Using Function 2.4.1 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 Efficiency of matrix computations.)

   Test your function on the matrix with elements

   $$A_{ij} = \begin{cases} \frac{1}{i+j}, & -1 \le i-j \le 2,\\ 0 & \text{otherwise.} \end{cases}$$

   (2.9.11)

5. ⌨ 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.)

6. ✍ 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}$.)