Sparsity and structure Tobin A. Driscoll Richard J. Braun
Very large matrices cannot be stored all within primary memory of a computer unless they are sparse . A sparse matrix has structural zeros , meaning entries that are known to be exactly zero. For instance, the adjacency matrix of a graph has zeros where there are no links in the graph. To store and operate with a sparse matrix efficiently, it is not represented as an array of all of its values. There is a variety of sparse formats available; for the most part, you can imagine that the matrix is stored as triples ( i , j , A i j ) (i,j,A_{ij}) ( i , j , A ij ) ( i , j ) (i,j) ( i , j )
8.1.1 Computing with sparse matrices ¶ Most graphs with real applications have many fewer edges than the maximum possible n 2 n^2 n 2 n n n
Example 8.1.1 Here we load the adjacency matrix of a graph with 2790 nodes. Each node is a web page referring to Roswell, NM, and the edges represent links between web pages. (Credit goes to Panayiotis Tsaparas and the University of Toronto for making this data public.)
using SparseArrays, JLD2
@load "roswell.jld2" A; # file is on the book's website
We may define the density of A \mathbf{A} A
Tip
Use nnz to count the number of nonzeros in a sparse matrix.
m, n = size(A)
@show density = nnz(A) / (m * n);density = nnz(A) / (m * n) = 0.0010902994565845762 The computer memory consumed by any variable can be discovered using summarysize. We can use it to compare the space needed for the sparse representation to its dense counterpart, that is, the space needed to store all the elements, whether zero or not.
F = Matrix(A)
Base.summarysize(F) / Base.summarysize(A)As you can see, the storage savings are dramatic. Matrix-vector products are also much faster using the sparse form because operations with structural zeros are skipped.
x = randn(n)
A * x; # make sure * is loaded and compiled
@elapsed for i in 1:300
A * x
endF * x;
@elapsed for i in 1:300
F * x
endExample 8.1.1 Here we load the adjacency matrix of a graph with 2790 nodes. Each node is a web page referring to Roswell, NM, and the edges represent links between web pages. (Credit goes to Panayiotis Tsaparas and the University of Toronto for making this data public.)
load roswelladj
a = whos('A')We may define the density of A \mathbf{A} A
Tip
Use nnz to count the number of nonzeros in a sparse matrix.
sz = size(A); n = sz(1);
density = nnz(A) / prod(sz)The computer memory consumed by any variable can be discovered using whos. We can use it to compare the space needed for the sparse representation to its dense counterpart, that is, the space needed to store all the elements, whether zero or not.
F = full(A);
f = whos('F');
storage_ratio = f.bytes / a.bytesMatrix-vector products are also much faster using the sparse form because operations with structural zeros are skipped.
x = randn(n,1);
tic, for i = 1:200, A*x; end
sparse_time = toctic, for i = 1:200, F*x; end
dense_time = tocHowever, the sparse storage format in MATLAB is column-oriented. Operations on rows may take a lot longer than similar ones on columns.
v = A(:, 1000);
tic, for i = 1:n, A(:, i) = v; end
column_time = toc
r = v';
tic, for i = 1:n, A(i, :) = r; end
row_time = tocExample 8.1.1 Here we load the adjacency matrix of a graph with 2790 nodes. Each node is a web page referring to Roswell, NM, and the edges represent links between web pages. (Credit goes to Panayiotis Tsaparas and the University of Toronto for making this data public.)
import scipy.sparse as sp
from scipy.io import loadmat
vars = loadmat("roswelladj.mat") # get from the book's website
A = sp.csr_matrix(vars["A"])
We may define the density of A \mathbf{A} A
m, n = A.shape
print(f"density is {A.nnz / (m * n):.3%}")We can compare the storage space needed for the sparse A \mathbf{A} A
F = A.todense()
print(f"{A.data.nbytes/1e6:.3f} MB for sparse form, {F.nbytes/1e6:.3f} MB for dense form")0.068 MB for sparse form, 62.273 MB for dense form
Matrix-vector products are also much faster using the sparse form because operations with structural zeros are skipped.
from timeit import default_timer as timer
x = random.randn(n)
start = timer()
for i in range(1000):
A @ x
print(f"sparse time: {timer() - start:.4g} sec")start = timer()
for i in range(1000):
F @ x
print(f"dense time: {timer() - start:.4g} sec")Arithmetic operations such as +, -, *, and ^ respect and exploit sparsity if the matrix operands are sparse. However, matrix operations may substantially decrease the amount of sparsity, a phenomenon known as fill-in .
Example 8.1.2 (Fill-in of a sparse matrix)
Example 8.1.2 Here is the adjacency matrix of a graph representing a small-world network, featuring connections to neighbors and a small number of distant contacts.
using GraphRecipes
@load "smallworld.jld2" A
graphplot(A, linealpha=0.5)Because each node connects to relatively few others, the adjacency matrix is quite sparse.
spy(A, title="Nonzero locations", m=2, color=:blues)By Theorem 7.1.1 A k \mathbf{A}^k A k k k k k degrees of separation” within a social network. As k k k A k \mathbf{A}^k A k
plt = plot(layout=(1, 3), legend=:none, size=(600, 240))
for k in 2:4
spy!(A^k;
subplot=k - 1, color=:blues,
title=latexstring("\\mathbf{A}^$k"))
end
pltExample 8.1.2 Here is the adjacency matrix of a graph representing a small-world network, featuring connections to neighbors and a small number of distant contacts.
load smallworld.mat
G = graph(A);
plot(G, nodecolor='r')Because each node connects to relatively few others, the adjacency matrix is quite sparse.
By Theorem 7.1.1 A k \mathbf{A}^k A k k k k k degrees of separation” within a social network. As k k k A k \mathbf{A}^k A k
clf
tiledlayout(2, 2)
for k = [2, 3, 4, 6]
nexttile
spy(A^k)
title(sprintf("A^{%d}", k))
endExample 8.1.2 Here is the adjacency matrix of a graph representing a small-world network, featuring connections to neighbors and a small number of distant contacts.
import networkx as nx
wsg = nx.watts_strogatz_graph(200, 4, 0.02)
Because each node connects to relatively few others, the adjacency matrix is quite sparse.
A = nx.adjacency_matrix(wsg)
spy(A)
title("Adjacency matrix $A$");By Theorem 7.1.1 A k \mathbf{A}^k A k k k k k degrees of separation” within a social network. As k k k A k \mathbf{A}^k A k
Tip
While A**6 is valid syntax here, it means elementwise power, not matrix power.
from scipy.sparse.linalg import matrix_power
spy(matrix_power(A, 6))
title(("$A^6$"));8.1.2 Banded matrices ¶ A particularly important type of sparse matrix is a banded matrix. Recall from Exploiting matrix structure that A \mathbf{A} A upper bandwidth p p p j − i > p j-i > p j − i > p A i j = 0 A_{ij}=0 A ij = 0 lower bandwidth q q q i − j > q i-j > q i − j > q A i j = 0 A_{ij}=0 A ij = 0 bandwidth is p + q + 1 p+q+1 p + q + 1
Without pivoting, an LU factorization preserves bandwidth, but pivoting can change or destroy bandedness.
Example 8.1.3 (Banded matrices)
Example 8.1.3 The spdiagm function creates a sparse matrix given its diagonal elements. 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.
n = 50;
A = spdiagm(-3 => fill(n, n - 3),
0 => ones(n),
1 => -(1:n-1),
5 => fill(0.1, n - 5))
Matrix(A[1:7, 1:7])7×7 Matrix{Float64}:
1.0 -1.0 0.0 0.0 0.0 0.1 0.0
0.0 1.0 -2.0 0.0 0.0 0.0 0.1
0.0 0.0 1.0 -3.0 0.0 0.0 0.0
50.0 0.0 0.0 1.0 -4.0 0.0 0.0
0.0 50.0 0.0 0.0 1.0 -5.0 0.0
0.0 0.0 50.0 0.0 0.0 1.0 -6.0
0.0 0.0 0.0 50.0 0.0 0.0 1.0
Without pivoting, the LU factors have the same lower and upper bandwidth as the original matrix.
Tip
The sparse function converts any matrix to sparse form. But it’s usually better to construct a sparse matrix directly, as the standard form might not fit in memory.
L, U = FNC.lufact(A)
plot(layout=2)
spy!(sparse(L), m=2, subplot=1, title=L"\mathbf{L}", color=:blues)
spy!(sparse(U), m=2, subplot=2, title=L"\mathbf{U}", color=:blues)However, if we introduce row pivoting, bandedness may be expanded or destroyed.
fact = lu(A)
plot(layout=2)
spy!(sparse(fact.L), m=2, subplot=1, title=L"\mathbf{L}", color=:blues)
spy!(sparse(fact.U), m=2, subplot=2, title=L"\mathbf{U}", color=:blues)Example 8.1.3 The spdiags function creates a sparse matrix given its diagonal elements. 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.
n = 50;
n = 50;
% Put constant values on 3 diagonals
A = spdiags([n, 1, 0.1], [-3, 0, 5], n, n);
% Put other values on 1st superdiagonal
A = spdiags(-(0:n-1)', 1, A);
full(A(1:7, 1:7))Without pivoting, the LU factors have the same lower and upper bandwidth as the original matrix.
Tip
The sparse function converts any matrix to sparse form. But it’s usually better to construct a sparse matrix directly, as the standard form might not fit in memory.
[L, U] = lufact(A);
clf
subplot(1, 2, 1), spy(L), title('L')
subplot(1, 2, 2), spy(U), title(('U'));However, if we introduce row pivoting, bandedness may be expanded or destroyed.
[L, U, p] = plufact(A);
subplot(1, 2, 1), spy(L(p, :)), title('L')
subplot(1, 2, 2), spy(U), title(('U'));Example 8.1.3 The scipi.sparse.diags function creates a sparse matrix given its diagonal elements and the diagonal indexes to put them on. 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.
n = 50
data = [n * ones(n-3), ones(n), linspace(-1, 1-n, n-1)]
offsets = [-3, 0, 1] # 3rd below, main, 1st above
A = sp.diags(data, offsets, format="lil")
print(A[:7, :7].todense())[[ 1. -1. 0. 0. 0. 0. 0.]
[ 0. 1. -2. 0. 0. 0. 0.]
[ 0. 0. 1. -3. 0. 0. 0.]
[50. 0. 0. 1. -4. 0. 0.]
[ 0. 50. 0. 0. 1. -5. 0.]
[ 0. 0. 50. 0. 0. 1. -6.]
[ 0. 0. 0. 50. 0. 0. 1.]]
Without pivoting, the LU factors have the same lower and upper bandwidth as the original matrix.
L, U = FNC.lufact(A.todense())
subplot(1, 2, 1), spy(L)
subplot(1, 2, 2), spy(U);However, if we introduce row pivoting, bandedness may be expanded or destroyed.
L, U, p = FNC.plufact(A.todense())
subplot(1, 2, 1), spy(L[p, :])
subplot(1, 2, 2), spy(U)8.1.3 Systems and eigenvalues ¶ If given a sparse matrix, the backslash operator will automatically try a form of sparse-aware Cholesky or pivoted LU factorization. Depending on the sparsity pattern of the matrix, the time taken to solve the linear system may be well below the O ( n 3 ) O(n^3) O ( n 3 )
For very large matrices, it’s unlikely that you will want to find all of its eigenvalues and eigenvectors. In Krylov subspaces we describe some of the math behind an algorithm that can find a selected number of eigenvalues of largest magnitude, lying to the extreme left or right, or nearest a given complex number.
Example 8.1.4 (Eigenvalues of sparse matrices)
Example 8.1.4 The following generates a random sparse matrix with prescribed eigenvalues.
n = 4000
density = 4e-4
λ = @. 1 + 1 / (1:n) # exact eigenvalues
A = FNC.sprandsym(n, density, λ);
The eigs function from Arpack finds a small number eigenvalues meeting some criterion. First, we ask for the 5 of largest (complex) magnitude using which=:LM.
using Arpack
λmax, V = eigs(A, nev=5, which=:LM) # Largest Magnitude
fmt = ft_printf("%20.15f")
pretty_table([λmax λ[1:5]], header=["found", "exact"], formatters=fmt)┌──────────────────────┬──────────────────────┐
│ found │ exact │
├──────────────────────┼──────────────────────┤
│ 1.999999999999999 │ 2.000000000000000 │
│ 1.500000000000001 │ 1.500000000000000 │
│ 1.333333333333337 │ 1.333333333333333 │
│ 1.250000000000000 │ 1.250000000000000 │
│ 1.199999999999997 │ 1.200000000000000 │
└──────────────────────┴──────────────────────┘
Now we find the 5 closest to the value 1 in the complex plane, via sigma=1.
λ1, V = eigs(A, nev=5, sigma=1) # closest to sigma
data = [λ1 λ[end:-1:end-4]]
pretty_table(data, header=["found", "exact"], formatters=fmt)┌──────────────────────┬──────────────────────┐
│ found │ exact │
├──────────────────────┼──────────────────────┤
│ 1.000250000000000 │ 1.000250000000000 │
│ 1.000250062515629 │ 1.000250062515629 │
│ 1.000250125062531 │ 1.000250125062531 │
│ 1.000250187640731 │ 1.000250187640731 │
│ 1.000250250250250 │ 1.000250250250250 │
└──────────────────────┴──────────────────────┘
The time needed to solve a sparse linear system is not easy to predict unless you have some more information about the matrix. But it will typically be orders of magnitude faster than the dense version of the same problem.
x = @. 1 / (1:n);
b = A * x;
norm(x - A \ b); # force compilation
t = @elapsed sparse_err = norm(x - A \ b)
println("Time for sparse solve: $t")Time for sparse solve: 0.000542542
D = Matrix(A) # convert to regular matrix
norm(x - D \ b);
t = @elapsed dense_err = norm(x - D \ b)
println("Time for dense solve: $t")Time for dense solve: 0.118054875
@show sparse_err;
@show dense_err;sparse_err = 2.2855717272351584e-17
dense_err = 2.2392539063344847e-16
Example 8.1.4 The following generates a random sparse matrix with prescribed eigenvalues.
n = 4000;
density = 4e-4;
lambda = 1 ./ (1:n);
A = sprandsym(n, density, lambda);
clf, spy(A)
title('Sparse symmetric matrix')The eigs function finds a small number eigenvalues meeting some criterion. First, we ask for the 5 of largest (complex) magnitude.
[V, D] = eigs(A, 5); % largest magnitude
1 ./ diag(D) % should be 1, 2, 3, 4, 5Now we find the 4 closest to the value 0.03 in the complex plane.
[V, D] = eigs(A, 4, 0.03); % closest to 0.03
diag(D)The time needed to solve a sparse linear system is not easy to predict unless you have some more information about the matrix. But it will typically be orders of magnitude faster than the dense version of the same problem.
x = 1 ./ (1:n)';
b = A * x;
tic, sparse_err = norm(x - A\b), sparse_time = tocF = full(A);
tic, dense_err = norm(x - F\b), dense_time = tocExample 8.1.4 The following generates a random sparse matrix with prescribed eigenvalues.
n = 4000
density = 4e-4
ev = 1 / arange(1, n + 1)
A = FNC.sprandsym(n, density, eigvals=ev)
print(f"density is {A.nnz / prod(A.shape):.3%}")The eigs function finds a small number eigenvalues meeting some criterion. First, we ask for the 5 of largest (complex) magnitude using which="LM".
from scipy.sparse.linalg import eigs
ev, V = eigs(A, k=5, which="LM") # largest magnitude
print(1 / ev)[1.+0.j 2.+0.j 3.+0.j 4.+0.j 5.+0.j]
Now we find the 4 closest to the value 1 in the complex plane, via sigma=1.
from scipy.sparse.linalg import eigs
ev, V = eigs(A, k=4, sigma=0.03) # closest to sigma
print(ev)[0.03030303+0.j 0.02941176+0.j 0.03125 +0.j 0.02857143+0.j]
The time needed to solve a sparse linear system is not easy to predict unless you have some more information about the matrix. But it will typically be orders of magnitude faster than the dense version of the same problem.
from scipy.sparse.linalg import spsolve
x = 1 / arange(1, n + 1)
b = A @ x
start = timer()
xx = spsolve(A, b)
print(f"sparse time: {timer() - start:.3g} sec")
print(f"residual: {norm(b - A @ xx, 2):.1e}")sparse time: 0.0037 sec
residual: 1.1e-16
from numpy.linalg import solve
F = A.todense()
start = timer()
xx = solve(F, b)
print(f"dense time: {timer() - start:.3g} sec")
print(f"residual: {norm(b - A @ xx, 2):.1e}")dense time: 0.863 sec
residual: 1.1e-16
8.1.4 Exercises ¶ ⌨ Use spdiagm to build the 50 × 50 50\times 50 50 × 50
A = [ − 2 1 1 − 2 1 ⋱ ⋱ ⋱ 1 − 2 1 1 − 2 ] , B = [ − 2 1 1 1 − 2 1 ⋱ ⋱ ⋱ 1 − 2 1 1 1 − 2 ] . \mathbf{A} =
\begin{bmatrix}
-2 & 1 & & & \\
1 & -2 & 1 & & \\
& \ddots & \ddots & \ddots & \\
& & 1 & -2 & 1 \\
& & & 1 & -2
\end{bmatrix}, \qquad
\mathbf{B} =
\begin{bmatrix}
-2 & 1 & & & 1 \\
1 & -2 & 1 & & \\
& \ddots & \ddots & \ddots & \\
& & 1 & -2 & 1 \\
1 & & & 1 & -2
\end{bmatrix}. A = ⎣ ⎡ − 2 1 1 − 2 ⋱ 1 ⋱ 1 ⋱ − 2 1 1 − 2 ⎦ ⎤ , B = ⎣ ⎡ − 2 1 1 1 − 2 ⋱ 1 ⋱ 1 ⋱ − 2 1 1 1 − 2 ⎦ ⎤ . For each matrix, use spy and an inspection of the 5 × 5 5\times 5 5 × 5
⌨ This problem requires the matrix used in Example 8.1.2
Download smallworld.mat by clicking the link and saving (you may need to fix the file name).
using MAT
A = matread("smallworld.mat")["A"]Download smallworld.mat by clicking the link and saving (you may need to fix the file name).
Download smallworld.mtx by clicking the link and saving (you may need to fix the file name).
import scipy.io as spio
A = spio.mmread("smallworld.mtx")(a) Find the density of A \mathbf{A} A A 2 \mathbf{A}^2 A 2 A 4 \mathbf{A}^4 A 4 A 8 \mathbf{A}^8 A 8 A \mathbf{A} A
(b) The LU factors tend to at least partially retain sparsity. Find the density of the L \mathbf{L} L U \mathbf{U} U A \mathbf{A} A lufact (Algorithm 2.4.1
(c) Repeat part (b) for the QR factorization using qrfact (Algorithm 3.4.1
⌨ One use of adjacency matrices is to analyze the links between members of a collection. Obtain the adjacency matrix A \mathbf{A} A Example 8.1.1
Download roswell.mat by clicking the link and saving (you may need to fix the file name).
using MAT
A = matread("roswell.mat")["A"]Download roswell.mat by clicking the link and saving (you may need to fix the file name).
Download roswell.mtx by clicking the link and saving (you may need to fix the file name).
import scipy.io as spio
A = spio.mmread("roswell.mtx")The matrix catalogs the links between web sites related to the town of Roswell, NM, with A i j = 1 A_{ij}=1 A ij = 1 i i i j j j
(a) Verify numerically that the matrix does not include any links from a site to itself.
(b) Verify numerically that A \mathbf{A} A
(c) How many sites in the group are not pointed to by any other sites in the group?
(d) Which site points to the most other sites?
(e) Which site is pointed to the most by the other sites? This is a crude way to establish the most important site.
(f) There are 27902 possible ways to connect ordered pairs of sites. What fraction of these pairs is connected by a walk of links that is no greater than three in length?
⌨ The graph Laplacian matrix is L = D − A \mathbf{L}=\mathbf{D}-\mathbf{A} L = D − A A \mathbf{A} A D \mathbf{D} D degree matrix , a diagonal matrix with diagonal entries d j j = ∑ i = 1 n a i j d_{jj}=\sum_{i=1}^n a_{ij} d jj = ∑ i = 1 n a ij
Follow the directions in Exercise 8.1.3 A \mathbf{A} A L \mathbf{L} L
⌨ See Exercise 7.1.5 A \mathbf{A} A Internet Movie Database . Specifically, A i j = 1 A_{ij}=1 A ij = 1 j j j i i i
(a) What is the maximum number of actors appearing in any one movie?
(b) How many actors appeared in exactly three movies?
(c) Define C = A T A \mathbf{C}=\mathbf{A}^T\mathbf{A} C = A T A C \mathbf{C} C C i j C_{ij} C ij
⌨ A matrix that arises from the Helmholtz equation for wave propagation can be specified using
A = FNC.poisson(n) - k^2*I;where k k k n = 50 n=50 n = 50
(a) Let k = 1 k=1 k = 1 A \mathbf{A} A
(b) Still with k = 1 k=1 k = 1 eigs to find the four largest and four smallest (in magnitude) eigenvalues of A \mathbf{A} A Example 8.1.4
(c) The eigenvalues are all real. Find a value of k k k A \mathbf{A} A