# Linear algebra in SciPy

> No one can be told what the matrix is. You have to see it for yourself.
>
> — Morpheus, *The Matrix*

Just like Chapter 4, which dealt with the Fast Fourier Transform, this chapter
will feature an elegant *method*. We
want to highlight the packages available in SciPy to do linear algebra, which forms
the basis of much scientific computing.

## Linear algebra basics

A chapter in a programming book is not really the right place to learn about
linear algebra itself, so we assume familiarity with linear algebra concepts.
At a minimum, you should know that linear algebra involves vectors (ordered
collections of numbers) and their transformations by multiplying them with
matrices (collections of vectors). If all of this sounded like gibberish to
you, you should probably pick up an introductory linear algebra textbook before
reading this. We highly recommend Gil Strang's *Linear Algebra and its
Applications*. Introductory is all you need though — we hope to convey the power
of linear algebra while keeping the operations relatively simple!

As an aside, we will break Python notation convention in order to match linear
algebra conventions: in Python, variable names should usually begin with a
lower case letter. However, in linear algebra, matrices are denoted by
a capital letter, while vectors and scalar values are lowercase. Since we're
going to be dealing with quite a few matrices and vectors, following the
linear algebra convention helps to keep them straight. Therefore, variables
that represent matrices will start with a capital letter, while vectors and
numbers will start with lowercase:

```python
import numpy as np

m, n = (5, 6)  # scalars
M = np.ones((m, n))  # a matrix
v = np.random.random((n,))  # a vector
w = M @ v  # another vector
```

In mathematical notation, the vectors would typically be written in boldface,
as in $\mathbf{v}$ and $\mathbf{w}$, while the scalars would not, as
in $m$ and $n$.  In Python code, we can't make that distinction, so we will rely instead
on context to keep scalars and vectors straight.

## Laplacian matrix of a graph

We discussed graphs in chapter 3, where we represented image regions as
nodes, connected by edges between them. But we used a rather simple method of
analysis: we *thresholded* the graph, removing all edges above some value.
Thresholding works in simple cases, but all you need
is one value to fall on the wrong side of the threshold for the approach
to fail.

As an example, suppose you are at war, and your enemy is camped just across the
river from your forces. You want to cut them off, so you decide to blow up all
the bridges between you. Intelligence suggests that you need $t$ kilograms of
TNT to blow each bridge crossing the river, but the bridges in your own
territory can withstand $t+1$ kg. You might, having just read chapter 3 of
*Elegant SciPy*, order your commandos to detonate $t$ kg of TNT on every bridge
in the region. But, if intelligence was wrong about just *one* bridge crossing
the river, and it remains standing, the enemy's army can come marching through!
Disaster!

So, in this chapter, we will explore some alternative approaches to graph
analysis, based on linear algebra. It turns out that we can think of a graph, $G$,
as an *adjacency matrix*, in which we number the nodes of the graph from $0$
to $n-1$, and place a 1 in row $i$, column $j$ of the matrix whenever there is
an edge from node $i$ to node $j$. In other words, if we call the adjacency
matrix $A$, then $A_{i, j} = 1$ if and only if the edge $(i, j)$ is in $G$. We
can then use linear algebra techniques to study this matrix, often with
striking results.

The *degree* of a node is defined as the number of edges touching it.  For
example, if a node is connected to five other nodes in a graph, its degree
is 5. (Later, we will differentiate between out-degree and in-degree, when edges
have a "from" and "to".) In matrix terms, the degree corresponds to the *sum*
of the values in a row or column.

The *Laplacian* matrix of a graph (just "the Laplacian" for short) is defined
as the *degree matrix*, $D$, which
contains the degree of each node along the diagonal and zero everywhere else,
minus the adjacency matrix $A$:

$
L = D - A
$

We definitely can't fit all of the linear algebra theory needed to understand
the properties of this matrix, but suffice it to say: it has some *great*
properties. We will exploit a couple in the following paragraphs.

First, we will look at the *eigenvectors* of $L$.
An eigenvector $v$ of a matrix $M$ is a vector that
satisfies the property $Mv = \lambda v$ for some number $\lambda$,
known as the eigenvalue.  In other words, $v$ is a special vector in
relation to $M$ because $Mv$ simply changes the size of the vector, without
changing its direction. As we will soon see, eigenvectors have many useful
properties — sometimes seeming even magical!

As an example, a 3x3 rotation matrix $R$, when multiplied by any
3-dimensional vector $p$, rotates it $30^\circ$ degrees around the z-axis.  $R$
will rotate all vectors except for those that lie *on* the z-axis.  For those,
we'll see no effect, or $Rp = p$, i.e. $Rp = \lambda p$ with
eigenvalue $\lambda = 1$.

<!-- exercise begin -->

**Exercise:** Consider the rotation matrix

$
R = \begin{bmatrix}
  \cos \theta &  -\sin \theta & 0 \\
  \sin \theta & \cos \theta & 0\\
  0 & 0 & 1\\
\end{bmatrix}
$

When $R$ is multiplied with a 3-dimensional column-vector $p =
\left[ x\, y\, z \right]^T$, the resulting vector $R p$ is rotated
by $\theta$ degrees around the z-axis.

1. For $\theta = 45^\circ$, verify (by testing on a few arbitrary
   vectors) that $R$ rotates these vectors around the z axis.
   Remember that matrix multiplication in Python is denoted with `@`.

2. What does the matrix $S = RR$ do? Verify this in Python.

3. Verify that multiplying by $R$ leaves the vector
   $\left[ 0\, 0\, 1\right]^T$ unchanged.  In other words, $R p = 1
   p$, which means $p$ is an eigenvector of $R$ with eigenvalue 1.

4. Use `np.linalg.eig` to find the eigenvalues and eigenvectors of $R$, and
   verify that $\left[0, 0, 1\right]^T$ is indeed among them, and that it
   corresponds to the eigenvalue 1.

<!-- solution begin -->

**Solution:**

Part 1:

```python
import numpy as np

theta = np.deg2rad(45)
R = np.array([[np.cos(theta), -np.sin(theta), 0],
              [np.sin(theta),  np.cos(theta), 0],
              [            0,              0, 1]])

print("R times the x-axis:", R @ [1, 0, 0])
print("R times the y-axis:", R @ [0, 1, 0])
print("R times a 45 degree vector:", R @ [1, 1, 0])
```

Part 2:

Since multiplying a vector by $R$ rotates it 45 degrees, multiplying the result
by $R$ again should result in the original vector being rotated 90 degrees.
Matrix multiplication is associative, which means that $R(Rv) = (RR)v$, so
$S = RR$ should rotate vectors by 90 degrees around the z-axis.

```python
S = R @ R
S @ [1, 0, 0]
```

Part 3:

```python
print("R @ z-axis:", R @ [0, 0, 1])
```

R rotates both the x and y axes, but not the z-axis.

Part 4:

Looking at the documentation of `eig`, we see that it returns two values:
a 1D array of eigenvalues, and a 2D array where each column contains the
eigenvector corresponding to each eigenvalue.

```python
np.linalg.eig(R)
```

In addition to some complex-valued eigenvalues and vectors, we see the value 1
associated with the vector $\left[0, 0, 1\right]^T$.

<!-- solution end -->

<!-- exercise end -->

Back to the Laplacian. A common problem in network analysis is visualization.
How do you draw nodes and edges in such a way that you don't get a complete
mess such as this one?

<img src="https://upload.wikimedia.org/wikipedia/commons/9/90/Visualization_of_wiki_structure_using_prefuse_visualization_package.png"/>
<!-- caption text="Visualization of wikipedia structure. Created by Chris Davis and released under CC-BY-SA-3.0 (https://commons.wikimedia.org/wiki/GNU_Free_Documentation_License)." -->


One way is to put nodes that share many edges close together. It turns out
that this can be done by using the second-smallest eigenvalue of the Laplacian
matrix, and its corresponding eigenvector, which is so important it has its
own name: the
[Fiedler vector](https://en.wikipedia.org/wiki/Algebraic_connectivity#The_Fiedler_vector).

Let's use a minimal network to illustrate this. We start by creating the
adjacency matrix:

```python
import numpy as np
A = np.array([[0, 1, 1, 0, 0, 0],
              [1, 0, 1, 0, 0, 0],
              [1, 1, 0, 1, 0, 0],
              [0, 0, 1, 0, 1, 1],
              [0, 0, 0, 1, 0, 1],
              [0, 0, 0, 1, 1, 0]], dtype=float)
```

We can use NetworkX to draw this network. First, we initialize matplotlib as
usual:

```python
# Make plots appear inline, set custom plotting style
%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('style/elegant.mplstyle')
```

Now, we can plot it:

```python
import networkx as nx
g = nx.from_numpy_matrix(A)
layout = nx.spring_layout(g, pos=nx.circular_layout(g))
nx.draw(g, pos=layout,
        with_labels=True, node_color='white')
```
<!-- caption text="A simple network plotted with NetworkX" -->


You can see that the nodes fall naturally into two groups, 0, 1, 2 and 3, 4, 5.
Can the Fiedler vector tell us this? First, we must compute the degree matrix
and the Laplacian. We first get the degrees by summing along either axis of $A$.
(Either axis works because $A$ is symmetric.)

```python
d = np.sum(A, axis=0)
print(d)
```

We then put those degrees into a diagonal matrix of the same shape
as A, the *degree matrix*. We can use the `np.diag` function to do this:

```python
D = np.diag(d)
print(D)
```

Finally, we get the Laplacian from the definition:

```python
L = D - A
print(L)
```

Because $L$ is symmetric, we can use the `np.linalg.eigh` function to compute
the eigenvalues and eigenvectors:

```python
val, Vec = np.linalg.eigh(L)
```

You can verify that the values returned satisfy the definition of eigenvalues
and eigenvectors. For example, one of the eigenvalues is 3:

```python
np.any(np.isclose(val, 3))
```

And we can check that multiplying the matrix $L$ by the corresponding eigenvector
does indeed multiply the vector by 3:

```python
idx_lambda3 = np.argmin(np.abs(val - 3))
v3 = Vec[:, idx_lambda3]

print(v3)
print(L @ v3)
```

As mentioned above, the Fiedler vector is the vector corresponding to the
second-smallest eigenvalue of $L$. Sorting the eigenvalues tells us which one
is the second-smallest:

```python
plt.plot(np.sort(val), linestyle='-', marker='o');
```
<!-- caption text="Eigenvalues of $L$" -->

It's the first non-zero eigenvalue, close to 0.4. The Fiedler vector is the
corresponding eigenvector:

```python
f = Vec[:, np.argsort(val)[1]]
plt.plot(f, linestyle='-', marker='o');
```
<!-- caption text="Fiedler vector of $L$" -->

It's pretty remarkable: by looking at the *sign* of the elements of the Fiedler
vector, we can separate the nodes into the two groups we identified in the
drawing!

```python
colors = ['orange' if eigv > 0 else 'gray' for eigv in f]
nx.draw(g, pos=layout, with_labels=True, node_color=colors)
```
<!-- caption text="Nodes colored by their sign in the Fiedler vector of $L$" -->

## Laplacians with brain data

Let's demonstrate this process in a real-world example by laying out the brain cells in a worm, as shown in
[Figure 2](http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1001066)
from the
[Varshney *et al* paper](http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1001066)
that we introduced in Chapter 3. (Information on
how to do this is in the
[supplementary material](http://journals.plos.org/ploscompbiol/article/asset?unique&id=info:doi/10.1371/journal.pcbi.1001066.s001)
for the paper.) To obtain their
layout of the worm brain neurons, they used a related matrix, the
*degree-normalized Laplacian*.

Because the order of the neurons is important in this analysis, we will use a
preprocessed dataset, rather than clutter this chapter with data cleaning. We
got the original data from [Lav Varshney's GitHub](https://github.com/lrvarshney/elegans),
and the processed data is in our `data/` directory.

First, let's load the data. There are four components:
- the network of chemical synapses, through which a *pre-synaptic neuron*
  sends a chemical signal to a *post-synaptic* neuron,
- the gap junction network, which contains direct electrical contacts between
  neurons),
- the neuron IDs (names), and
- the three neuron types:
  - *sensory neurons*, those that detect signals coming from the outside world,
    encoded as 0;
  - *motor neurons*, those that activate muscles, enabling the worm to move,
    encoded as 2; and
  - *interneurons*, the neurons in between, which enable complex signal processing
    to occur between sensory neurons and motor neurons, encoded as 1.

```python
import numpy as np
Chem = np.load('data/chem-network.npy')
Gap = np.load('data/gap-network.npy')
neuron_ids = np.load('data/neurons.npy')
neuron_types = np.load('data/neuron-types.npy')
```

We then simplify the network, adding the two kinds of connections together,
and removing the directionality of the network by taking the average of
in-connections and out-connections of neurons. This seems a bit like cheating
but, since we are only looking for the *layout* of the neurons on a graph, we
only care about *whether* neurons are connected, not in which direction.
We are going to call the resulting matrix the *connectivity* matrix, $C$, which
is just a different kind of adjacency matrix.

```python
A = Chem + Gap
C = (A + A.T) / 2
```

To get the Laplacian matrix $L$, we need the degree matrix $D$, which contains
the degree of node i at position [i, i], and zeros everywhere else.

```python
degrees = np.sum(C, axis=0)
D = np.diag(degrees)
```

Now, we can get the Laplacian just like before:

```python
L = D - C
```

The vertical coordinates in Fig 2 are given by arranging nodes such that, on
average, neurons are as close as possible to "just above" their downstream
neighbors. Varshney _et al_ call this measure "processing depth," and it's
obtained by solving a linear equation involving the Laplacian. We use
`scipy.linalg.pinv`, the
[pseudoinverse](https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_pseudoinverse),
to solve it:

```python
from scipy import linalg
b = np.sum(C * np.sign(A - A.T), axis=1)
z = linalg.pinv(L) @ b
```

(Note the use of the `@` symbol, which was introduced in Python 3.5 to denote
matrix multiplication. As we noted in the preface and in Chapter 5, in previous
versions of Python, you would need to use the function `np.dot`.)

In order to obtain the degree-normalized Laplacian, $Q$, we need the inverse
square root of the $D$ matrix:

```python
Dinv2 = np.diag(1 / np.sqrt(degrees))
Q = Dinv2 @ L @ Dinv2
```

Finally, we are able to extract the $x$ coordinates of the neurons to ensure that
highly-connected neurons remain close: the eigenvector of $Q$ corresponding to
its second-smallest eigenvalue, normalized by the degrees:

```python
val, Vec = linalg.eig(Q)
```

Note from the documentation of `numpy.linalg.eig`:

> "The eigenvalues are not necessarily ordered."

Although the documentation in SciPy's `eig` lacks this warning, it remains true
in this case. We must therefore sort the eigenvalues and the corresponding
eigenvector columns ourselves:

```python
smallest_first = np.argsort(val)
val = val[smallest_first]
Vec = Vec[:, smallest_first]
```

Now we can find the eigenvector we need to compute the affinity coordinates:

```python
x = Dinv2 @ Vec[:, 1]
```

(The reasons for using this vector are too long to explain here, but appear in
the paper's supplementary material, linked above. The short version is that
choosing this vector minimizes the total length of the links between neurons.)

There is one small kink that we must address before proceeding: eigenvectors
are only defined up to a multiplicative constant. This follows simply from the
definition of an eigenvector: suppose $v$ is an eigenvector of the matrix $M$,
with corresponding eigenvalue $\lambda$. Then $\alpha v$ is also an eigenvector
of $M$ for any scalar number $\alpha$,
because $Mv = \lambda v$ implies $M(\alpha v) = \lambda (\alpha v)$.
So, it is
arbitrary whether a software package returns $v$ or $-v$ when asked for the
eigenvectors of $M$. In order to make sure we reproduce the layout from the
Varshney *et al.* paper, we must make sure that the vector is pointing in the
same direction as theirs, rather than the opposite direction. We do this by
choosing an arbitrary neuron from their Figure 2, and checking the sign of `x`
at that position. We then reverse it if it doesn't match its sign in Figure 2
of the paper.

```python
vc2_index = np.argwhere(neuron_ids == 'VC02')
if x[vc2_index] < 0:
    x = -x
```

Now it's just a matter of drawing the nodes and the edges. We color them
according to the type stored in `neuron_types`, using the appealing and
functional "colorblind"
[colorbrewer palette](https://chrisalbon.com/python/data_visualization/seaborn_color_palettes/):

```python
from matplotlib.colors import ListedColormap
from matplotlib.collections import LineCollection


def plot_connectome(x_coords, y_coords, conn_matrix, *,
                    labels=(), types=None, type_names=('',),
                    xlabel='', ylabel=''):
    """Plot neurons as points connected by lines.

    Neurons can have different types (up to 6 distinct colors).

    Parameters
    ----------
    x_coords, y_coords : array of float, shape (N,)
        The x-coordinates and y-coordinates of the neurons.
    conn_matrix : array or sparse matrix of float, shape (N, N)
        The connectivity matrix, with non-zero entry (i, j) if and only
        if node i and node j are connected.
    labels : array-like of string, shape (N,), optional
        The names of the nodes.
    types : array of int, shape (N,), optional
        The type (e.g. sensory neuron, interneuron) of each node.
    type_names : array-like of string, optional
        The name of each value of `types`. For example, if a 0 in
        `types` means "sensory neuron", then `type_names[0]` should
        be "sensory neuron".
    xlabel, ylabel : str, optional
        Labels for the axes.
    """
    if types is None:
        types = np.zeros(x_coords.shape, dtype=int)
    ntypes = len(np.unique(types))
    colors = plt.rcParams['axes.prop_cycle'][:ntypes].by_key()['color']
    cmap = ListedColormap(colors)

    fig, ax = plt.subplots()

    # plot neuron locations:
    for neuron_type in range(ntypes):
        plotting = (types == neuron_type)
        pts = ax.scatter(x_coords[plotting], y_coords[plotting],
                         c=cmap(neuron_type), s=4, zorder=1)
        pts.set_label(type_names[neuron_type])

    # add text labels:
    for x, y, label in zip(x_coords, y_coords, labels):
        ax.text(x, y, '   ' + label,
                verticalalignment='center', fontsize=3, zorder=2)

    # plot edges
    pre, post = np.nonzero(conn_matrix)
    links = np.array([[x_coords[pre], x_coords[post]],
                      [y_coords[pre], y_coords[post]]]).T
    ax.add_collection(LineCollection(links, color='lightgray',
                                     lw=0.3, alpha=0.5, zorder=0))

    ax.legend(scatterpoints=3, fontsize=6)

    ax.set_xlabel(xlabel, fontsize=8)
    ax.set_ylabel(ylabel, fontsize=8)

    plt.show()
```

Now, let's use that function to plot the neurons:


```python
plot_connectome(x, z, C, labels=neuron_ids, types=neuron_types,
                type_names=['sensory neurons', 'interneurons',
                            'motor neurons'],
                xlabel='Affinity eigenvector 1', ylabel='Processing depth')
```
<!-- caption text="Spectral layout of the neurons of a nematode worm" -->

There you are: a worm brain!
As discussed in the original paper, you can see the top-down processing from
sensory neurons to motor neurons through a network of interneurons. You can
also see two distinct groups of motor neurons: these correspond to the neck
(left) and body (right) body segments of the worm.

<!-- exercise begin -->

**Exercise:** How do you modify the above code to show the affinity view in
Figure 2B from the paper?

<!-- solution begin -->
**Solution:** In the affinity view, instead of using the processing depth on the y-axis,
we use the normalized third eigenvector of Q, just like we did with x. (And we
invert it if necessary, just like we did with x!)

```python
y = Dinv2 @ Vec[:, 2]
asjl_index = np.argwhere(neuron_ids == 'ASJL')
if y[asjl_index] < 0:
    y = -y

plot_connectome(x, y, C, labels=neuron_ids, types=neuron_types,
                type_names=['sensory neurons', 'interneurons',
                            'motor neurons'],
                xlabel='Affinity eigenvector 1',
                ylabel='Affinity eigenvector 2')
```
<!-- caption text="Spectral layout of the neurons of a nematode worm, using two
spectral dimensions" -->

<!-- solution end -->

<!-- exercise end -->

<!-- exercise begin -->

### Challenge: linear algebra with sparse matrices

The above code uses numpy arrays to hold the matrix and perform
the necessary computations. Because we are using a small graph of fewer than 300
nodes, this is feasible. However, for larger graphs, it would fail.

For example, one might want to analyse the relationships between libraries
listed on the Python Package Index, or PyPI, which contains over one hundred thousand packages.
Holding the Laplacian matrix for this graph would take 
up $8 \left(100 \times 10^3\right)^2 = 8 \times 10^{10}$ bytes, or 80GB,
of RAM. If you add to that the adjacency, symmetric adjacency, pseudoinverse,
and, say, two temporary matrices used during calculations, you climb up to
480GB, beyond the reach of most desktop computers.

"Ha!", some of you might be thinking. "Ha! My desktop has 512GB of RAM! It would
make short work of this so-called 'large' graph!"

Perhaps. But you might also want to analyze the Association for Computing
Machinery (ACM) citation graph, a network of over two million scholarly works
and references. *That* Laplacian would take up 32 terabytes of RAM.

However, we know that the dependency and reference graphs are *sparse*:
packages usually depend on just a few other packages, not on the whole of PyPI.
And papers and books usually only reference a few others, too. So we can hold
the above matrices using the sparse data structures from `scipy.sparse` (see
Chapter 5), and use the linear algebra functions in `scipy.sparse.linalg` to
compute the values we need.

Try to explore the documentation in `scipy.sparse.linalg` to come up with a
sparse version of the above computation.

Hint: the pseudoinverse of a sparse matrix is, in general, not sparse, so you
can't use it here. Similarly, you can't get all the eigenvectors of a sparse
matrix, because they would together make up a dense matrix.

You'll find parts of the solution below (and of course in the solutions
chapter), but we highly recommend that you try it out on your own.

<!-- solution begin -->

### Challenge accepted

For the purposes of this challenge, we are going to use the small connectome
above, because it's easier to visualise what is going on. In later parts of the
challenge we'll use these techniques to analyze larger networks.

First, we start with the adjacency matrix, A, in a sparse matrix format, in
this case, CSR, which is the most common format for linear algebra. We'll
append `s` to the names of all the matrices to indicate that they are sparse.

```python
from scipy import sparse
import scipy.sparse.linalg

As = sparse.csr_matrix(A)
```

We can create our connectivity matrix in the same way:

```python
Cs = (As + As.T) / 2
```

In order to get the degrees matrix, we can use the "diags" sparse format, which
stores diagonal and off-diagonal matrices.

```python
degrees = np.ravel(Cs.sum(axis=0))
Ds = sparse.diags(degrees)
```

Getting the Laplacian is straightforward:

```python
Ls = Ds - Cs
```

Now we want to get the processing depth. Remember that getting the
pseudoinverse of the Laplacian matrix is out of the question, because it will
be a dense matrix (the inverse of a sparse matrix is not generally sparse
itself). However, we were actually using the pseudoinverse to compute a
vector $z$ that would satisfy $L z = b$,
where $b = C \odot \textrm{sign}\left(A - A^T\right) \mathbf{1}$.
(You can see this in the supplementary material for Varshney *et al*.) With
dense matrices, we can simply use $z = L^+b$. With sparse ones, though, we can
use one of the *solvers* (see sidebox, "Solvers") in `sparse.linalg.isolve` to get the `z` vector after
providing `L` and `b`, no inversion required!

```python
b = Cs.multiply((As - As.T).sign()).sum(axis=1)
z, error = sparse.linalg.isolve.cg(Ls, b, maxiter=10000)
```

Finally, we must find the eigenvectors of $Q$, the degree-normalized Laplacian,
corresponding to its second and third smallest eigenvalues.

You might recall from Chapter 5 that the numerical data in sparse matrices is
in the `.data` attribute. We use that to invert the degrees matrix:

```python
Dsinv2 = Ds.copy()
Dsinv2.data = 1 / np.sqrt(Ds.data)
```

Finally, we use SciPy's sparse linear algebra functions to find the desired
eigenvectors. The $Q$ matrix is symmetric, so we can use the `eigsh` function,
specialized for symmetric matrices, to compute them. We use the `which` keyword
argument to specify that we want the eigenvectors corresponding to the smallest
eigenvalues, and `k` to specify that we need the 3 smallest:

```python

Qs = Dsinv2 @ Ls @ Dsinv2
vals, Vecs = sparse.linalg.eigsh(Qs, k=3, which='SM')
sorted_indices = np.argsort(vals)
Vecs = Vecs[:, sorted_indices]
```

Finally, we normalize the eigenvectors to get the x and y coordinates
(and flip these if necessary):

```python
_dsinv, x, y = (Dsinv2 @ Vecs).T
if x[vc2_index] < 0:
    x = -x
if y[asjl_index] < 0:
    y = -y
```

(Note that the eigenvector corresponding to the smallest eigenvalue is always a
vector of all ones, which we're not interested in.)
We can now reproduce the above plots!

```python
plot_connectome(x, z, C, labels=neuron_ids, types=neuron_types,
                type_names=['sensory neurons', 'interneurons',
                            'motor neurons'],
                xlabel='Affinity eigenvector 1', ylabel='Processing depth')

plot_connectome(x, y, C, labels=neuron_ids, types=neuron_types,
                type_names=['sensory neurons', 'interneurons',
                            'motor neurons'],
                xlabel='Affinity eigenvector 1',
                ylabel='Affinity eigenvector 2')
```
<!-- caption text="Spectral layout of a nematode brain, computed using sparse
matrices" -->

<!-- solution end -->

<!-- exercise end -->

> **Solvers {.callout}**
>
> SciPy has several sparse iterative solvers available, and it is not always
> obvious which to use.  Unfortunately, that question also has no easy answer:
> different algorithms have different strengths in terms of speed of
> convergence, stability, accuracy, and memory use (amongst others).  It is also
> not possible to predict, by looking at the input data, which algorithm will
> perform best.
> 
> Here is a rough guideline for choosing an iterative solver:
> 
> - If A, the input matrix, is symmetric and positive definite, use the
>   Conjugate Gradient solver `cg`.  If A is symmetric, but
>   near-singular or indefinite, try the Minimum Residual iteration
>   method `minres`.
> 
> - For non-symmetric systems, try the Biconjugate Gradient Stabilized
>   method, `bicgstab`.  The Conjugate Gradient Squared method, `cgs`,
>   is a bit faster, but has more erratic convergence.
> 
> - If you need to solve many similar systems, use the LGMRES algorithm `lgmres`.
> 
> - If A is not square, use the least squares algorithm `lsmr`.
> 
> For further reading, see
> 
> - **How Fast are Nonsymmetric Matrix Iterations?**,
>   Noël M. Nachtigal, Satish C. Reddy, and Lloyd N. Trefethen
>   SIAM Journal on Matrix Analysis and Applications 1992 13:3, 778-795.
> 
> - **Survey of recent Krylov methods**, Jack Dongarra,
>   http://www.netlib.org/linalg/html_templates/node50.html


## PageRank: linear algebra for reputation and importance

Another application of linear algebra and eigenvectors is Google's PageRank
algorithm, which is punnily named both for webpages and for one of its
co-founders, Larry Page.

To rank webpages by importance, you might count
how many other webpages link to it. After all, if everyone is linking to a
particular page, it must be good, right? But this metric is
easily gamed: to make your own webpage rise in the rankings, just
create as many other webpages as you can and have them all link to your
original page.

The key insight that drove Google's early success was that important webpages
are not linked to by just many webpages, but by *important*
webpages. And how do we know that those other pages are important? Because
they themselves are linked to by important pages. And so on.

This recursive definition implies that page
importance can be measured by an eigenvector
of the so-called *transition matrix*, which contains the links
between webpages. Suppose you have your vector of importance $\boldsymbol{r}$,
and your matrix of links $M$. You don't know $\boldsymbol{r}$ yet, but you
do know that the importance of a page is proportional to the sum of
importances of the pages that link to it: $\boldsymbol{r} = \alpha M \boldsymbol{r}$,
or $M \boldsymbol{r} = \lambda \boldsymbol{r}$, for $\lambda = 1/\alpha$. That's
just the definition of an eigenvalue!

By ensuring some special properties are satisfied by the transition matrix, we
can further determine that the required eigenvalue is 1, and that it is the
largest eigenvalue of $M$.

The transition matrix imagines a web
surfer, often named Webster, randomly clicking a link from each webpage he
visits, and then asks, what's the probability that he ends up at any given
page? This probability is called the PageRank.

Since Google's rise, researchers have been applying PageRank to all sorts of
networks. We'll use an example by Stefano Allesina and Mercedes Pascual,
which they
[published](http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1000494)
in PLoS Computational Biology. They thought to apply the method in ecological
*food webs*, networks that link species to those that they eat.

Naively, if you wanted to see how critical a species was for an ecosystem, you
would look at how many species eat it. If it's many, and that species
disappeared, then all its "dependent" species might disappear with it. In
network parlance, you could say that its *in-degree* determines its ecological
importance.

Could PageRank be a better measure of importance for an ecosystem?

Professor Allesina kindly provided us with a few food webs to play around
with. We've saved one of these, from the St Marks National Wildlife Refuge in
Florida, in the Graph Markup Language format. The web was
[described](http://www.sciencedirect.com/science/article/pii/S0304380099000228)
in 1999 by Robert R. Christian and Joseph J. Luczovich. In the dataset, a
node $i$ has an edge to node $j$ if species $i$ eats species $j$.

We'll start by loading in the data, which NetworkX knows how to read trivially:

```python
import networkx as nx

stmarks = nx.read_gml('data/stmarks.gml')
```

Next, we get the sparse matrix corresponding to the graph. Because a matrix
only holds numerical information, we need to maintain a separate list of
package names corresponding to the matrix rows/columns:

```python
species = np.array(list(stmarks.nodes()))  # array for multi-indexing
Adj = nx.to_scipy_sparse_matrix(stmarks, dtype=np.float64)
```

From the adjacency matrix, we can derive a *transition probability* matrix,
where every edge is replaced by a *probability* of 1 over the number of
outgoing edges from that species. In the food web, it might make more sense
to call this a lunch probability matrix.

The total number of species in our matrix is going to be used a lot, so let's
call it $n$:

```python
n = len(species)
```

Next, we need the degrees, and, in particular, the *diagonal matrix* containing
the inverse of the out-degrees of each node on the diagonal:

```python
np.seterr(divide='ignore')  # ignore division-by-zero errors
from scipy import sparse

degrees = np.ravel(Adj.sum(axis=1))
Deginv = sparse.diags(1 / degrees).tocsr()
```

```python
Trans = (Deginv @ Adj).T
```

Normally, the PageRank score would simply be the first eigenvector of the
transition matrix. If we call the transition matrix $M$ and the vector of
PageRank values $r$, we have:

$$
\boldsymbol{r} = M\boldsymbol{r}
$$

But the `np.seterr` call above is a clue that it's not quite
so simple. The PageRank approach only works when the
transition matrix is a *column-stochastic* matrix, in which every
column sums to 1. Additionally, every page must be reachable
from every other page, even if the path to reach it is very long.

In our food web, this causes problems, because the bottom of the food chain,
what the authors call *detritus* (basically sea sludge), doesn't actually *eat*
anything (the Circle of Life notwithstanding), so you can't reach other species
from it.

> *Young Simba:* But, Dad, don't we eat the antelope?
>
> *Mufasa:* Yes, Simba, but let me explain. When we die, our bodies become the
> grass, and the antelope eat the grass. And so we are all connected in the
> great Circle of Life.
>
> — *The Lion King*

To deal with this, the PageRank algorithm uses a so-called "damping
factor", usually taken to be 0.85. This means that 85% of the time, the
algorithm follows a link at random, but for the other 15%, it randomly jumps to
any arbitrary page. It's as if every page had a low probability link to every
other page. Or, in our case, it's as if shrimp, on rare occasions, ate sharks.
It might seem non-sensical but bear with us! It is, in fact, the mathematical
representation of the Circle of Life. We'll set it to 0.85, but actually it
doesn't really matter for this analysis: the results are similar for a large
range of possible damping factors.

If we call the damping factor $d$, then the modified PageRank equation is:

$$
\boldsymbol{r} = dM\boldsymbol{r} + \frac{1-d}{n} \boldsymbol{1}
$$

and

$$
(\boldsymbol{I} - dM)\boldsymbol{r} = \frac{1-d}{n} \boldsymbol{1}
$$

We can solve this equation using `scipy.sparse.linalg`'s direct
solver, `spsolve`. Depending on the structure and size of a linear algebra
problem, though, it might be more efficient to use an iterative solver. See
the [`scipy.sparse.linalg` documentation](http://docs.scipy.org/doc/scipy/reference/sparse.linalg.html#solving-linear-problems)
for more information on this.

```python
from scipy.sparse.linalg import spsolve

damping = 0.85
beta = 1 - damping

I = sparse.eye(n, format='csc')  # Same sparse format as Trans

pagerank = spsolve(I - damping * Trans,
                   np.full(n, beta / n))
```

We now have the "foodrank" of the St. Marks food web!

So how does a species' foodrank compare to the number of other species eating
it?

```python
def pagerank_plot(in_degrees, pageranks, names, *,
                  annotations=[], **figkwargs):
    """Plot node pagerank against in-degree, with hand-picked node names."""

    fig, ax = plt.subplots(**figkwargs)
    ax.scatter(in_degrees, pageranks, c=[0.835, 0.369, 0], lw=0)
    for name, indeg, pr in zip(names, in_degrees, pageranks):
        if name in annotations:
            text = ax.text(indeg + 0.1, pr, name)

    ax.set_ylim(0, np.max(pageranks) * 1.1)
    ax.set_xlim(-1, np.max(in_degrees) * 1.1)
    ax.set_ylabel('PageRank')
    ax.set_xlabel('In-degree')
```

We now draw the plot. Having explored the dataset before writing this, we have
pre-labeled some interesting nodes in the plot:

```python
interesting = ['detritus', 'phytoplankton', 'benthic algae', 'micro-epiphytes',
               'microfauna', 'zooplankton', 'predatory shrimps', 'meiofauna',
               'gulls']
in_degrees = np.ravel(Adj.sum(axis=0))
pagerank_plot(in_degrees, pagerank, species, annotations=interesting)
```

Sea sludge ("detritus") is the most important element both by number of
species feeding on it (15) and by PageRank (>0.003). But the second most
important element is *not* benthic algae, which feeds 13 other species, but
rather phytoplankton, which feeds just 7! That's because other *important*
species feed on it. On the bottom left, we've got sea gulls, who, we can now
confirm, do bugger-all for the ecosystem. Those vicious *predatory shrimps*
(we're not making this up) support the same number of species as phytoplankton,
but they are less essential species, so they end up with a lower foodrank.

Although we won't do it here, Allesina and Pascual go on to model the
ecological impact of species extinction, and indeed find that PageRank
predicts ecological importance better than in-degree.

Before we wrap up though, we'll note that PageRank can be computed several
different ways. One way, complementary to what we did above, is called the
*power method*, and it's quite, well, powerful! It stems from the
[Perron-Frobenius theorem](https://en.wikipedia.org/wiki/Perron%E2%80%93Frobenius_theorem),
which states, among other things, that a stochastic matrix has 1 as an
eigenvalue, and that this is its *largest* eigenvalue. (The corresponding
eigenvector is the PageRank vector.) What this means is that, whenever we
multiply *any* vector by $M$, its component pointing towards this major
eigenvector stays the same, while *all other components shrink* by a
multiplicative factor. The consequence is that if we multiply some random
starting vector by $M$ repeatedly, we should eventually get the PageRank
vector!

SciPy makes this very efficient with its sparse matrix module:

```python
def power(Trans, damping=0.85, max_iter=10**5):
    n = Trans.shape[0]
    r0 = np.full(n, 1/n)
    r = r0
    for _iter_num in range(max_iter):
        rnext = damping * Trans @ r + (1 - damping) / n
        if np.allclose(rnext, r):
            break
        r = rnext
    return r
```

<!-- exercise begin -->

**Exercise:** In the above iteration, note that `Trans` is *not*
column-stochastic, so the $r$ vector gets shrunk at each iteration. In order to
make the matrix stochastic, we have to replace every zero-column by a column of
all $1/n$. This is too expensive, but computing the iteration is cheaper. How
can you modify the code above to ensure that $r$ remains a probability vector
throughout?

<!-- solution begin -->

**Solution:** In order to have a stochastic matrix, all columns of the
transition matrix must sum to 1. This is not satisfied when a species isn't
eaten by any others: that column will consist of all zeroes. Replacing all
those columns by $1/n \boldsymbol{1}$, however, would be expensive.

The key is to realise that *every row* will contribute the *same amount* to the
multiplication of the transition matrix by the current probability vector. That
is to say, adding these columns will add a single value to the result of the
iteration multiplication. What value? $1/n$ times the elements of $r$ that
correspond to a dangling node. This can be expressed as a dot-product of a
vector containing $1/n$ for positions corresponding to dangling nodes, and zero
elswhere, with the vector $r$ for the current iteration.

```python
def power2(Trans, damping=0.85, max_iter=10**5):
    n = Trans.shape[0]
    dangling = (1/n) * np.ravel(Trans.sum(axis=0) == 0)
    r0 = np.full(n, 1/n)
    r = r0
    for _ in range(max_iter):
        rnext = (damping * (Trans @ r + dangling @ r) +
                 (1 - damping) / n)
        if np.allclose(rnext, r):
            return rnext
        else:
            r = rnext
    return r
```

Try this out manually for a few iterations. Notice that if you start with a
stochastic vector (a vector whose elements all sum to 1), the next vector will
still be a stochastic vector. Thus, the output PageRank from this function will
be a true probability vector, and the values will represent the probability
that we end up at a particular species when following links in the food chain.

<!-- solution end -->

<!-- exercise end -->


<!-- exercise begin -->

**Exercise:** Verify that these three methods all give the same ranking for the
nodes. `numpy.corrcoef` might be a useful function for this.

<!-- solution begin -->

**Solution:** `np.corrcoef` gives the Pearson correlation coefficient between
all pairs of a list of vectors. This coefficient will be equal to 1 if and only
if two vectors are scalar multiples of each other. Therefore, a correlation
coefficient of 1 is sufficient to show that the above methods produce the same
ranking.

```python
pagerank_power = power(Trans)
pagerank_power2 = power2(Trans)
np.corrcoef([pagerank, pagerank_power, pagerank_power2])
```

<!-- solution end -->

<!-- exercise end -->

## Concluding remarks

The field of linear algebra is far too broad to adequately cover in a chapter,
but this chapter gave you a glimpse into its power, and of
the way Python, NumPy, and SciPy make its elegant algorithms accessible.
