# Sparse Matrices¶

%matplotlib inline
import numpy as np
import pandas as pd
from scipy import sparse
import scipy.sparse.linalg as spla
import matplotlib.pyplot as plt
import seaborn as sns

sns.set_context('notebook', font_scale=1.5)


## Creating a sparse matrix¶

There are many applications in which we deal with matrices that are mostly zeros. For example, a matrix representing social networks is very sparse - there are 7 billion people, but most people are only connected to a few hundred or thousand others directly. Storing such a social network as a sparse rather than dense matrix will offer orders of magnitude reductions in memory requirements and corresponding speed-ups in computation.

### Coordinate format¶

The simplest sparse matrix format is built from the coordinates and values of the non-zero entries.

#### From dense matrix¶

A = np.random.poisson(0.2, (5,15)) * np.random.randint(0, 10, (5, 15))
A

array([[ 0,  5,  0,  0,  8,  0,  4,  0,  0,  0,  7,  0,  0,  0,  0],
[ 0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
[ 0,  0,  0,  8,  0,  0,  2,  2,  4,  0,  0,  0,  7,  0,  0],
[ 0,  0,  7,  0,  9,  0,  0,  0,  0,  1, 12,  0,  0,  0,  0],
[ 9,  2,  0,  2,  0,  0,  0,  0,  4,  0,  0,  0,  8,  0,  0]])

rows, cols = np.nonzero(A)
vals = A[rows, cols]

vals

array([ 5,  8,  4,  7,  8,  2,  2,  4,  7,  7,  9,  1, 12,  9,  2,  2,  4,
8])

rows

array([0, 0, 0, 0, 2, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4])

cols

array([ 1,  4,  6, 10,  3,  6,  7,  8, 12,  2,  4,  9, 10,  0,  1,  3,  8,
12])

X1 = sparse.coo_matrix(A)
X1

<5x15 sparse matrix of type '<class 'numpy.int64'>'
with 18 stored elements in COOrdinate format>

print(X1)

(0, 1)    5
(0, 4)    8
(0, 6)    4
(0, 10)   7
(2, 3)    8
(2, 6)    2
(2, 7)    2
(2, 8)    4
(2, 12)   7
(3, 2)    7
(3, 4)    9
(3, 9)    1
(3, 10)   12
(4, 0)    9
(4, 1)    2
(4, 3)    2
(4, 8)    4
(4, 12)   8


#### From coordinates¶

Note that the (values, (rows, cols)) argument is a single tuple.

X2 = sparse.coo_matrix((vals, (rows, cols)))
X2

<5x13 sparse matrix of type '<class 'numpy.int64'>'
with 18 stored elements in COOrdinate format>

print(X2)

(0, 1)    5
(0, 4)    8
(0, 6)    4
(0, 10)   7
(2, 3)    8
(2, 6)    2
(2, 7)    2
(2, 8)    4
(2, 12)   7
(3, 2)    7
(3, 4)    9
(3, 9)    1
(3, 10)   12
(4, 0)    9
(4, 1)    2
(4, 3)    2
(4, 8)    4
(4, 12)   8


#### Convert back to dense matrix¶

X2.todense()

matrix([[ 0,  5,  0,  0,  8,  0,  4,  0,  0,  0,  7,  0,  0],
[ 0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
[ 0,  0,  0,  8,  0,  0,  2,  2,  4,  0,  0,  0,  7],
[ 0,  0,  7,  0,  9,  0,  0,  0,  0,  1, 12,  0,  0],
[ 9,  2,  0,  2,  0,  0,  0,  0,  4,  0,  0,  0,  8]])


### Compressed Sparse Row and Column formats¶

When we have 2 or more repeated entries in the rows or cols, we can remove the redundancy by indicating the location of the first occurrence of a value and its increment instead of the full coordinates. These are known as CSR or CSC formats.

np.vstack([rows, cols])

array([[ 0,  0,  0,  0,  2,  2,  2,  2,  2,  3,  3,  3,  3,  4,  4,  4,  4,
4],
[ 1,  4,  6, 10,  3,  6,  7,  8, 12,  2,  4,  9, 10,  0,  1,  3,  8,
12]])

indptr = np.r_[np.searchsorted(rows, np.unique(rows)), len(rows)]
indptr

array([ 0,  4,  9, 13, 18])

X3 = sparse.csr_matrix((vals, cols, indptr))
X3

<4x13 sparse matrix of type '<class 'numpy.int64'>'
with 18 stored elements in Compressed Sparse Row format>

X3.todense()

matrix([[ 0,  5,  0,  0,  8,  0,  4,  0,  0,  0,  7,  0,  0],
[ 0,  0,  0,  8,  0,  0,  2,  2,  4,  0,  0,  0,  7],
[ 0,  0,  7,  0,  9,  0,  0,  0,  0,  1, 12,  0,  0],
[ 9,  2,  0,  2,  0,  0,  0,  0,  4,  0,  0,  0,  8]])


#### Casting from COO format¶

Because the coordinate format is more intuitive, it is often more convenient to first create a COO matrix then cast to CSR or CSC form.

X4 = X2.tocsr()

X4

<5x13 sparse matrix of type '<class 'numpy.int64'>'
with 18 stored elements in Compressed Sparse Row format>


### COO summation convention¶

When entries are repeated in a COO matrix, they are summed. This provides a quick way to construct confusion matrices for evaluation of multi-class classification algorithms.

rows = np.r_[np.zeros(4), np.ones(4)]
cols = np.repeat([0,1], 4)
vals = np.arange(8)

rows

array([ 0.,  0.,  0.,  0.,  1.,  1.,  1.,  1.])

cols

array([0, 0, 0, 0, 1, 1, 1, 1])

vals

array([0, 1, 2, 3, 4, 5, 6, 7])

X5 = sparse.csr_matrix((vals, (rows, cols)))

print(X5)

(0, 0)    6
(1, 1)    22


## Application: Confusion matrix¶

### Creating a 2 by 2 confusion matrix¶

obs = np.random.randint(0, 2, 100)
pred = np.random.randint(0, 2, 100)
vals = np.ones(100).astype('int')

pred

array([1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0,
1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1,
1, 0, 0, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 0,
1, 0, 1, 0, 1, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1,
0, 0, 1, 1, 0, 1, 0, 1])

vals.shape, obs.shape , pred.shape

((100,), (100,), (100,))

X6 = sparse.coo_matrix((vals, (pred, obs)))

X6.todense()

matrix([[21, 22],
[34, 23]])


#### Creating an $$n$$ by $$n$$ confusion matrix¶

For classifications with a large number of classes (e.g. image segmentation), the savings are even more dramatic.

from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.neighbors import KNeighborsClassifier

iris = datasets.load_iris()

knn = KNeighborsClassifier()
X_train, X_test, y_train, y_test = train_test_split(iris.data, iris.target,
test_size=0.5, random_state=42)

pred = knn.fit(X_train, y_train).predict(X_test)

X7 = sparse.coo_matrix((np.ones(len(pred)).astype('int'), (pred, y_test)))
pd.DataFrame(X7.todense(), index=iris.target_names, columns=iris.target_names)

setosa versicolor virginica
setosa 29 0 0
versicolor 0 23 4
virginica 0 0 19

## Application: PageRank¶

SciPy provides efficient routines for solving large sparse systems as for dense matrices. We will illustrate by calculating the page rank for airports using data from the Bureau of Transportation Statisitcs.

data = pd.read_csv('data/airports.csv', usecols=[0,1])

data.shape

(445827, 2)

data.head()

ORIGIN_AIRPORT_ID DEST_AIRPORT_ID
0 10135 10397
1 10135 10397
2 10135 10397
3 10135 10397
4 10135 10397
lookup = pd.read_csv('data/names.csv', index_col=0)

lookup.shape

(6404, 1)

lookup.head()

Description
Code
10001 Afognak Lake, AK: Afognak Lake Airport
10003 Granite Mountain, AK: Bear Creek Mining Strip
10004 Lik, AK: Lik Mining Camp
10005 Little Squaw, AK: Little Squaw Airport
10006 Kizhuyak, AK: Kizhuyak Bay
import networkx as nx


### Construct the sparse adjacency matrix¶

g = nx.from_pandas_dataframe(data, source='ORIGIN_AIRPORT_ID', target='DEST_AIRPORT_ID')

airports = np.array(g.nodes())


### Construct the transition matrix¶

out_degrees = np.ravel(adj_matrix.sum(axis=1))
diag_matrix = sparse.diags(1 / out_degrees).tocsr()


### Modify the transition matrix with a damping factor¶

The PageRank algorithm assumes that every node can be reached from every other node. To guard against case where a node has out-degree 0, we allow every node a small random chance of transitioning to any other node using a damping factor $$d$$. Then we solve the linear system to find the pagerank score $$r$$.

$r = (I - dM)^{-1}\frac{1-d}{N}\mathbb{1}$

or equivalently in the $$Ax = b$$ format

$(I - dM)r = \frac{1-d}{N}\mathbb{1}$
n = len(airports)
d = 0.85
I = sparse.eye(n, format='csc')
A = I - d * M
b = (1-d) / n * np.ones(n) # so the sum of all page ranks is 1

A.todense()

matrix([[ 1.,  0.,  0., ...,  0.,  0.,  0.],
[ 0.,  1.,  0., ...,  0.,  0.,  0.],
[ 0.,  0.,  1., ...,  0.,  0.,  0.],
...,
[ 0.,  0.,  0., ...,  1.,  0.,  0.],
[ 0.,  0.,  0., ...,  0.,  1.,  0.],
[ 0.,  0.,  0., ...,  0.,  0.,  1.]])

from scipy.sparse.linalg import spsolve

r =  spsolve(A, b)
r.sum()

0.99999999999999978

idx = np.argsort(r)

top10 = idx[-10:][::-1]
bot10 = idx[:10]

df = lookup.loc[airports[top10]]
df['degree'] = out_degrees[top10]
df['pagerank']= r[top10]
df

Description degree pagerank
Code
10397 Atlanta, GA: Hartsfield-Jackson Atlanta Intern... 158 0.043286
13930 Chicago, IL: Chicago O'Hare International 139 0.033956
11292 Denver, CO: Denver International 129 0.031434
11298 Dallas/Fort Worth, TX: Dallas/Fort Worth Inter... 108 0.027596
13487 Minneapolis, MN: Minneapolis-St Paul Internati... 108 0.027511
12266 Houston, TX: George Bush Intercontinental/Houston 110 0.025967
11433 Detroit, MI: Detroit Metro Wayne County 100 0.024738
14869 Salt Lake City, UT: Salt Lake City International 78 0.019298
14771 San Francisco, CA: San Francisco International 76 0.017820
14107 Phoenix, AZ: Phoenix Sky Harbor International 79 0.017000
df = lookup.loc[airports[bot10]]
df['degree'] = out_degrees[bot10]
df['pagerank']= r[bot10]
df

Description degree pagerank
Code
12265 Niagara Falls, NY: Niagara Falls International 1 0.000693
14025 Plattsburgh, NY: Plattsburgh International 1 0.000693
11695 Flagstaff, AZ: Flagstaff Pulliam 1 0.000693
16218 Yuma, AZ: Yuma MCAS/Yuma International 1 0.000693
14905 Santa Maria, CA: Santa Maria Public/Capt. G. A... 1 0.000710
13964 North Bend/Coos Bay, OR: Southwest Oregon Regi... 1 0.000710
10157 Arcata/Eureka, CA: Arcata 1 0.000710
14487 Redding, CA: Redding Municipal 1 0.000710
12177 Hobbs, NM: Lea County Regional 1 0.000711
11049 College Station/Bryan, TX: Easterwood Field 1 0.000711

### Visualize the airport connections graph and label the top and bottom 5 airports by pagerank¶

import warnings

labels = {airports[i]: lookup.loc[airports[i]].str.split(':').str[0].values[0]
for i in np.r_[top10[:5], bot10[:5]]}

with warnings.catch_warnings():
warnings.simplefilter('ignore')
nx.draw(g, pos=nx.spring_layout(g), labels=labels,
node_color='blue', font_color='red', alpha=0.5,
node_size=np.clip(5000*r, 1, 5000*r), width=0.1)