# Lecture 8: Sparse matrices

## Syllabus

Week 1: Python intro
Week 2: Matrices, vectors, norms, ranks
Week 3: Linear systems, eigenvectors, eigenvalues
Week 4: Singular value decomposition + test + homework seminar
Week 5: Sparse & structured matrices

## Recap of the previous lecture

• SVD, best low-rank approximations
• Applications
(Eigenfaces, collaborative filtering, integral equations, latent semantic indexing)

## Today lecture

Today we will talk about sparse matrices.

• Formats: list of lists and compressed sparse row format, relation to sparse graph
• Fast matrix-by-vector product
• Fast direct solvers for Gaussian elimination

## Why sparse matrices are important

• Describe local physical nature
• Provide storage reduction and matrix-by-vector product reduction
• Describe graphs: $a_{ij} \ne 0$ if there is an edge between vertex $i$ and vertex $j$.

## Sparse matrices are ubiqitous in PDE

Consider the simplest partial differential equation (PDE), called Poisson equation:
$$\Delta T = \frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} = f.$$

## Discretization

$$\frac{\partial^2 T}{\partial x^2} \approx \frac{T(x+h) + T(x-h) - 2T(x)}{h^2} + \mathcal{O}(h^2),$$

same for $\frac{\partial^2 T}{\partial y^2},$ and we get a linear system.
First, let us do one-dimensional case:

$$\frac{u_{i+1} + u_{i-1} - 2u_i}{h^2} = f_i,$$

or in the matrix form

$$A u = f,$$

and ($n = 5$ illustration) $$A = -\frac{1}{h^2}\begin{pmatrix} 2& -1 & 0 & 0 & 0\\ -1 & 2 & -1 & 0 &0 \\ 0 & -1 & 2& -1 & 0 \\ 0 & 0 & -1 & 2 &-1\\ 0 & 0 & 0 & -1 & 2 \end{pmatrix}$$

The matrix is triadiagonal and sparse
(and also Toeplitz: all elements on the diagonal are the same)

## Block structure in 2D

In two dimensions, we get equation of the form
$$\frac{4u_{ij} -u_{(i-1)j} - u_{(i+1)j} - u_{i(j-1)}-u_{i(j+1)}}{h^2} = f_{ij},$$ or in the Kronecker product form

$$\Delta_2 = \Delta_1 \otimes I + I \otimes \Delta_1,$$

where $\Delta_1$ is a 1D Laplacian, and $\otimes$ is a Kronecker product of matrices.
For matrices $A$ and $B$ its Kronecker product is defined as a block matrix of the form $$[a_{ij} B]$$

$$A = -\frac{1}{h^2}\begin{pmatrix} \Delta_1 + 2I & -I & 0 & 0 & 0\\ -I & \Delta_1 + 2I & -I & 0 &0 \\ 0 & -I & \Delta_1 + 2I & -I & 0 \\ 0 & 0 & -I & \Delta_1 + 2I &-I\\ 0 & 0 & 0 & -I & \Delta_1 + 2I \end{pmatrix}$$

We can create this matrix using scipy.sparse package



In [5]:

import numpy as np
import scipy as sp
import scipy.sparse
from scipy.sparse import csc_matrix
import matplotlib.pyplot as plt
%matplotlib inline
n = 20
ex = np.ones(n);
lp1 = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr');
e = sp.sparse.eye(n)
A = sp.sparse.kron(lp1, e) + sp.sparse.kron(e, lp1)
A = csc_matrix(A)
plt.spy(A, aspect='equal', marker='.', markersize=1)




Out[5]:

<matplotlib.lines.Line2D at 0x10ecb85d0>




In [3]:

#np.vstack((ex,  -2*ex, ex)).shape
print ex.shape




(20,)



We can now solve the system



In [33]:

import numpy as np
import scipy as sp
import scipy.sparse
import scipy.sparse.linalg
from scipy.sparse import csc_matrix, csr_matrix
import matplotlib.pyplot as plt
%matplotlib inline
n = 256
ex = np.ones(n);
lp1 = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr');
e = sp.sparse.eye(n)
A = sp.sparse.kron(lp1, e) + sp.sparse.kron(e, lp1)
A = csr_matrix(A) * (n + 1) ** 2
rhs = np.ones(n * n)
sol = sp.sparse.linalg.spsolve(A, rhs)



And plot the solution



In [34]:

from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
t = np.linspace(0, 1, n)
x, y = np.meshgrid(t, t)
sol = np.reshape(sol, (n, n))
ax.plot_wireframe(x, y, -sol)




Out[34]:

<mpl_toolkits.mplot3d.art3d.Line3DCollection at 0x111207610>




In [41]:

print A.shape
v256 = (np.max(-sol))
vex = (4 * v256 - v128) / 3
print vex - v256




(65536, 65536)
2.74259855287e-06



## How to store a sparse matrix

At first glance it is easy: you can use the coordinate format as

   (i, j, value)

i.e. store two integer arrays and one real array. This is not optimal.
And how to multiply a matrix by vector?

## Compressed sparse row

A matrix is stored as 3 different arrays:

ia, ja, sa

where:

• ia is an integer array of length $n+1$
• ja is an integer array of length nnz
• sa is an real-value array of length nnz
• nnz is the total number of non-zeros for the matrix

## Idea behind CSR.

• For each row $i$ we store the column number of the non-zeros (and their) values
• We stack this all together into ja and sa arrays
• We save the location of the first non-zero element in each row

## CSR helps for matrix-by-vector product as well

   for i in xrange(n):
for k in xrange(ia(i):ia(i+1)-1):
y(i) += sa(k) * x(ja(k))

Let us do a short timing test



In [44]:

import numpy as np
import scipy as sp
import scipy.sparse
import scipy.sparse.linalg
from scipy.sparse import csc_matrix, csr_matrix, coo_matrix
import matplotlib.pyplot as plt
%matplotlib inline
n = 256
ex = np.ones(n);
lp1 = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr');
e = sp.sparse.eye(n)
A = sp.sparse.kron(lp1, e) + sp.sparse.kron(e, lp1)
A = csr_matrix(A)
rhs = np.ones(n * n)
B = coo_matrix(A)
%timeit A.dot(rhs)
%timeit B.dot(rhs)




1000 loops, best of 3: 542 µs per loop
1000 loops, best of 3: 758 µs per loop



As you see, CSR is faster, and for more unstructured patterns the gain will be larger.

## Sparse matrices and efficiency

Sparse matrix give complexity reduction.
But they are not very good for parallel implementation.
And they do not give maximal efficiency due to random data access.
Typically, peak efficiency of 10-15% is considered good.
We can measure peak efficiency of an ordinary dot product:



In [45]:

import numpy as np
import time
n = 1000
a = np.random.randn(n, n)
v = np.random.randn(n)
t = time.time()
np.dot(a, v)
t = time.time() - t
print 'Efficiency:', ((2 * n ** 2)/t) / 10 ** 9, 'Gflops'




Efficiency: 2.3723438914 Gflops



And for the sparse one:



In [46]:

n = 1000
ex = np.ones(n);
a = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr');
rhs = np.random.randn(n)
t = time.time()
a.dot(rhs)
t = time.time() - t
print 'Efficiency:', (3 * n) / t / 10 ** 9, 'Gflops'




Efficiency: 0.0164913656619 Gflops



## How to solve linear systems?

How to solve linear systems with sparse matrices?
The direct methods use sparse Gaussian elimination, i.e.
they eliminate variables while trying to keep
the matrix as sparse as possible. And often, the inverse of a sparse matrix is not sparse (but $L$ and $U$ factors can be).



In [55]:

n = 5
ex = np.ones(n);
a = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr');
a = a.todense()
b = np.array(np.linalg.inv(a))
#print b
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
t = np.linspace(0, 1, n)
x, y = np.meshgrid(t, t)
ax.plot_wireframe(x, y, -b)




Out[55]:

<mpl_toolkits.mplot3d.art3d.Line3DCollection at 0x115e30610>



The L and U factors are typically better



In [56]:

p, l, u = scipy.linalg.lu(a)
print l




[[ 1.          0.          0.          0.          0.        ]
[-0.5         1.          0.          0.          0.        ]
[-0.         -0.66666667  1.          0.          0.        ]
[-0.         -0.         -0.75        1.          0.        ]
[-0.         -0.         -0.         -0.8         1.        ]]



However, for a 2D case everything is much worser:



In [58]:

n = 20
ex = np.ones(n);
lp1 = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr');
e = sp.sparse.eye(n)
A = sp.sparse.kron(lp1, e) + sp.sparse.kron(e, lp1)
A = csc_matrix(A)
T = scipy.sparse.linalg.splu(A)




In [60]:

plt.spy(T.U, marker='.', markersize=0.4)




Out[60]:

<matplotlib.lines.Line2D at 0x1140d3590>



For two-dimensional case the number of nonzeros in the L factor grows as $\mathcal{O}(N^{3/2})$.

## Sparse matrices and graph ordering

The number of non-zeros in LU decomposition has a deep connection to the graph theory.
(I.e., there is an edge between $(i, j)$ if $a_{ij} \ne 0$.



In [61]:

import networkx as nx
n = 10
ex = np.ones(n);
lp1 = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr');
e = sp.sparse.eye(n)
A = sp.sparse.kron(lp1, e) + sp.sparse.kron(e, lp1)
A = csc_matrix(A)
G = nx.Graph(A)
nx.draw(G, pos=nx.spring_layout(G), node_size=10)






## Strategies for elimination

The reordering that minimizes the fill-in is important, so we can you graph theory to find one.

• Minimum degree ordering - order by the degree of the vertex
• Cuthill–McKee algorithm (and reverse Cuthill-McKee) -- order for a small bandwidth
• Nested dissection: split the graph into two with minimal number of vertices on the separator


In [65]:

import networkx as nx
from networkx.utils import reverse_cuthill_mckee_ordering, cuthill_mckee_ordering
n = 10
ex = np.ones(n);
lp1 = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr');
e = sp.sparse.eye(n)
A = sp.sparse.kron(lp1, e) + sp.sparse.kron(e, lp1)
A = csc_matrix(A)
G = nx.Graph(A)
#rcm = list(reverse_cuthill_mckee_ordering(G))
rcm = list(reverse_cuthill_mckee_ordering(G))
A1 = A[rcm, :][:, rcm]
plt.spy(A1, marker='.', markersize=3)
#p, L, U = scipy.linalg.lu(A1.todense())
#plt.spy(L, marker='.', markersize=0.8)
#nx.draw(G, pos=nx.spring_layout(G), node_size=10)




Out[65]:

<matplotlib.lines.Line2D at 0x116f0d8d0>



The most efficient one is nested dissection algorithm, which finds the ordering recursively by computing the separator. No good implementation was found by me at the moment.

## Florida sparse matrix collection

Florida sparse matrix collection which contains all sorts of matrices for different applications. It also allows for finding test matrices as well! We can have a look.



In [66]:

from IPython.display import HTML
HTML('<iframe src=http://yifanhu.net/GALLERY/GRAPHS/search.html width=700 height=450></iframe>')




Out[66]:



## Test some

Let us check some sparse matrix (and its LU).



In [67]:

fname = 'crystm02.mat'
!wget http://www.cise.ufl.edu/research/sparse/mat/Boeing/$fname   --2014-10-02 10:32:43-- http://www.cise.ufl.edu/research/sparse/mat/Boeing/crystm02.mat Resolving www.cise.ufl.edu... 128.227.205.222 Connecting to www.cise.ufl.edu|128.227.205.222|:80... connected. HTTP request sent, awaiting response... 200 OK Length: 797624 (779K) [text/plain] Saving to: ‘crystm02.mat’ 100%[======================================>] 797,624 362KB/s in 2.2s 2014-10-02 10:32:46 (362 KB/s) - ‘crystm02.mat’ saved [797624/797624]   In [68]: from scipy.io import loadmat import scipy.sparse q = loadmat(fname) print q mat = q['Problem']['A'][0, 0] T = scipy.sparse.linalg.splu(mat)   {'Problem': array([[ ([u'FEM CRYSTAL FREE VIBRATION MASS MATRIX'], <13965x13965 sparse matrix of type '<type 'numpy.float64'>' with 322905 stored elements in Compressed Sparse Column format>, [u'Boeing/crystm02'], [[354]], [u'1995'], [u'R. Grimes'], [u'T. Davis'], [u'materials problem'])]], dtype=[('title', 'O'), ('A', 'O'), ('name', 'O'), ('id', 'O'), ('date', 'O'), ('author', 'O'), ('ed', 'O'), ('kind', 'O')]), '__version__': '1.0', '__header__': 'MATLAB 5.0 MAT-file, Platform: GLNXA64, Created on: Sat Sep 6 04:12:47 2008', '__globals__': []}   In [69]: #Compute its LU %matplotlib inline import matplotlib.pyplot as plt plt.spy(T.L, markersize=0.1)   Out[69]: <matplotlib.lines.Line2D at 0x1142977d0>  ## Take home message • CSR format for storage • Sparse matrices & graphs ordering • Ordering is important for LU fill-in ## Next lecture • Structured matrices: the idea • Circulant matrices and the Fast Fourier Transform • Toeplitz matrices and shift-invariant structure • Displacement structure, solution of linear systems with Toeplitz matrices • Toeplitz and low-rank matrices ##### Questions?  In [134]: from IPython.core.display import HTML def css_styling(): styles = open("./styles/custom.css", "r").read() return HTML(styles) css_styling()   Out[134]: @font-face { font-family: "Computer Modern"; src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf'); } div.cell{ /*width:80%;*/ /*margin-left:auto !important; margin-right:auto;*/ } h1 { font-family: 'Alegreya Sans', sans-serif; } h2 { font-family: 'Fenix', serif; } h3{ font-family: 'Fenix', serif; margin-top:12px; margin-bottom: 3px; } h4{ font-family: 'Fenix', serif; } h5 { font-family: 'Alegreya Sans', sans-serif; } div.text_cell_render{ font-family: 'Alegreya Sans',Computer Modern, "Helvetica Neue", Arial, Helvetica, Geneva, sans-serif; line-height: 1.2; font-size: 160%; /*width:70%;*/ /*margin-left:auto;*/ margin-right:auto; } .CodeMirror{ font-family: "Source Code Pro"; font-size: 90%; } /* .prompt{ display: None; }*/ .text_cell_render h1 { font-weight: 200; font-size: 50pt; line-height: 110%; color:#CD2305; margin-bottom: 0.5em; margin-top: 0.5em; display: block; } .text_cell_render h5 { font-weight: 300; font-size: 16pt; color: #CD2305; font-style: italic; margin-bottom: .5em; margin-top: 0.5em; display: block; } .warning{ color: rgb( 240, 20, 20 ) } MathJax.Hub.Config({ TeX: { extensions: ["AMSmath.js"] }, tex2jax: { inlineMath: [ ['$','\$'], ["\$","\$"] ],
displayMath: [ ['$$','$$'], ["\$","\$"] ]
},
displayAlign: 'center', // Change this to 'center' to center equations.
"HTML-CSS": {
styles: {'.MathJax_Display': {"margin": 4}}
}
});