Sparsity and structure - Fundamentals of Numerical Computation

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_{ij})$ for all the nonzero $(i,j)$ locations.

8.1.1Computing with sparse matrices¶

Most graphs with real applications have many fewer edges than the maximum possible $n^2$ for $n$ nodes. Accordingly, their adjacency matrices have mostly zero elements and should be represented sparsely. Julia functions to deal with sparse matrices are found in the SparseArrays package in the standard library.

Example 8.1.1 (Sparsity)

Julia

MATLAB

Python

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

@load "roswell.jld2" A;      # file is on the book's website

We may define the density of $\mathbf{A}$ as the number of nonzeros divided by the total number of entries.

Use nnz to count the number of nonzeros in a sparse matrix.

m, n = size(A)
@show density = nnz(A) / (m * n);

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
end

F * x;
@elapsed for i in 1:300
    F * x
end

Example 8.1.1

load roswelladj
a = whos('A')

We may define the density of $\mathbf{A}$ as the number of nonzeros divided by the total number of entries.

Use nnz to count the number of nonzeros in a sparse matrix.

sz = size(A);  n = sz(1);
density = nnz(A) / prod(sz)

F = full(A);
f = whos('F');
storage_ratio = f.bytes / a.bytes

Matrix-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 = toc

tic, for i = 1:200, F*x; end
dense_time = toc

However, 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 = toc

Example 8.1.1

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 $\mathbf{A}$ as the number of nonzeros divided by the total number of entries.

m, n = A.shape
print(f"density is {A.nnz / (m * n):.3%}")

We can compare the storage space needed for the sparse $\mathbf{A}$ with the space needed for its dense / full counterpart.

F = A.todense()
print(f"{A.data.nbytes/1e6:.3f} MB for sparse form, {F.nbytes/1e6:.3f} 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)

Julia

MATLAB

Python

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.

@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, the entries of $\mathbf{A}^k$ give the number of walks of length $k$ between pairs of nodes, as with “k degrees of separation” within a social network. As $k$ grows, the density of $\mathbf{A}^k$ also grows.

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
plt

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.

load smallworld.mat
G = graph(A);
plot(G, nodecolor='r')

Because each node connects to relatively few others, the adjacency matrix is quite sparse.

spy(A)

clf
tiledlayout(2, 2)
for k = [2, 3, 4, 6]
    nexttile
    spy(A^k)
    title(sprintf("A^{%d}", k))
end

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.

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

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.2Banded matrices¶

A particularly important type of sparse matrix is a banded matrix. Recall from Exploiting matrix structure that $\mathbf{A}$ has upper bandwidth $p$ if $j-i > p$ implies $A_{ij}=0$ , and lower bandwidth $q$ if $i-j > q$ implies $A_{ij}=0$ . We say the total bandwidth is $p+q+1$ . Banded matrices appear naturally in many applications where each element interacts directly with only a few neighbors.

Without pivoting, an LU factorization preserves bandwidth, but pivoting can change or destroy bandedness.

Example 8.1.3 (Banded matrices)

Julia

MATLAB

Python

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

Without pivoting, the LU factors have the same lower and upper bandwidth as the original matrix.

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

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.

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

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.3Systems 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)$ needed in the general case.

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)

Julia

MATLAB

Python

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 finds a small number eigenvalues meeting some criterion. First, we ask for the 5 of largest (complex) magnitude using which=:LM.

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

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)

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

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

@show sparse_err;
@show dense_err;

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, 5

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

x = 1 ./ (1:n)';  
b = A * x;
tic, sparse_err = norm(x - A\b), sparse_time = toc

F = full(A);
tic, dense_err = norm(x - F\b), dense_time = toc

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

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)

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}")

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}")

8.1.4Exercises¶

⌨ Use spdiagm to build the $50\times 50$ matrices
$\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}.$
(8.1.1)
For each matrix, use spy and an inspection of the $5\times 5$ submatrices in the corners to verify the correctness of your matrices.
⌨ This problem requires the matrix used in Demo 8.1.2; you can load it via
```
url = "https://tobydriscoll.net/fnc-julia/_static/resources/smallworld.jld2"
datafile = download(url)
@load datafile A
```
(a) Find the density of $\mathbf{A}$ (number of nonzeros divided by total number of elements), $\mathbf{A}^2$ , $\mathbf{A}^4$ , and $\mathbf{A}^8$ . (You should find that it increases with the power of $A$ .)
(b) The LU factors tend to at least partially retain sparsity. Find the density of the $\mathbf{L}$ and $\mathbf{U}$ factors of $\mathbf{A}$ using the built-in lu. (See Demo 2.9.1 for usage of lu in the sparse case.)
(c) Repeat part (b) for the QR factorization using qr.
⌨ One use of adjacency matrices is to analyze the links between members of a collection. Obtain the adjacency matrix $\mathbf{A}$ from Demo 8.1.1 via the following:
```
url = "https://tobydriscoll.net/fnc-julia/_static/resources/roswell.jld2"
datafile = download(url)
@load datafile A
```
The matrix catalogs the links between web sites related to the town of Roswell, NM, with $A_{ij}=1$ if and only if site $i$ links to site $j$ .
(a) Verify numerically that the matrix does not include any links from a site to itself.
(b) Verify numerically that $\mathbf{A}$ is not symmetric. (Thus, its graph is a directed one.)
(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 2790² 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 $\mathbf{L}=\mathbf{D}-\mathbf{A}$ , where $\mathbf{A}$ is the adjacency matrix and $\mathbf{D}$ is the degree matrix, a diagonal matrix with diagonal entries $d_{jj}=\sum_{i=1}^n a_{ij}$ .
Follow the directions in Exercise 3 to obtain an adjacency matrix $\mathbf{A}$ . Then find the five eigenvalues of $\mathbf{L}$ having largest magnitude.
⌨ See Exercise 7.1.5 for instructions on loading a matrix $\mathbf{A}$ that contains information about the appearances of 392,400 actors in 127,823 movies, as given by the Internet Movie Database. Specifically, $A_{ij}=1$ if actor $j$ appeared in movie $i$ , and all other elements are zero.
(a) What is the maximum number of actors appearing in any one movie?
(b) How many actors appeared in exactly three movies?
(c) Define $\mathbf{C}=\mathbf{A}^T\mathbf{A}$ . How many nonzero entries does $\mathbf{C}$ have? What is the interpretation of $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$ is a real parameter. Let $n=50$ .
(a) Let $k=1$ . What is the size of $\mathbf{A}$ ? What is its density?
(b) Still with $k=1$ , use eigs to find the four largest and four smallest (in magnitude) eigenvalues of $\mathbf{A}$ . (See Demo 8.1.4 for examples.)
(c) The eigenvalues are all real. Find a value of $k$ so that $\mathbf{A}$ has exactly three negative eigenvalues.

Preface

Krylov methods in linear algebra

Preface

Power iteration