From matrix to insight - Fundamentals of Numerical Computation

Any two-dimensional array of numbers may be interpreted as a matrix. Whether or not this is the only point of view that matters to a particular application, it does lead to certain types of analysis. The related mathematical and computational tools are universally applicable and find diverse uses.

7.1.1Tables as matrices¶

Tables are used to represent variation of a quantity with respect to two variables. These variables may be encoded as the rows and columns of a matrix.

Example 7.1.1

A corpus is a collection of text documents. A term-document matrix has one column for each document and one row for each unique term appearing in the corpus. The $(i,j)$ entry of the matrix is the number of times term $i$ appears in document $j$ . That is, column $j$ of the matrix is a term-frequency vector quantifying all occurrences of the indexed terms. A new document could be represented by its term-frequency vector, which is then comparable to the columns of the matrix. Or, a new term could be represented by counting its appearances in all of the documents and be compared to the rows of the matrix.

It turns out that by finding the singular value decomposition of the term-document matrix, the strongest patterns within the corpus can be isolated, frequently corresponding to what we interpret as textual meaning. This is known as latent semantic analysis.

7.1.2Graphs as matrices¶

Definition 7.1.1 (Graphs and adjacency matrices)

A graph or network consists of a set $V$ of nodes and a set $E$ of edges, each of which is an ordered pair of nodes. If there is an edge $(v_i,v_j)$ , then we say that node $i$ is adjacent to node $j$ . The graph is undirected if for every edge $(v_i,v_j)$ , the pair $(v_j,v_i)$ is also an edge; otherwise the graph is directed.

The adjacency matrix of a graph with $n$ nodes $V$ and edge set $E$ is the $n\times n$ matrix whose elements are

A_{ij} = \begin{cases} 1 & \text{if $(v_i,v_j)\in E$ (i.e., node $i$ is adjacent to node $j$)},\\ 0 & \text{otherwise}. \end{cases}

(7.1.1)

Graphs are a useful way to represent the link structure of social networks, airline routes, power grids, sports teams, and web pages, to name a few examples. The natural interpretation is that the edge $(v_i,v_j)$ denotes a link from node $i$ to node $j$ , in which case we say that node $i$ is adjacent to node $j$ . One usually visualizes small graphs by drawing points for nodes and arrows or lines for the edges.

Here are some elementary results about adjacency matrices.

Example 7.1.4 (Adjacency matrix)

Julia

MATLAB

Python

Example 7.1.4

Here we create an adjacency matrix for a graph on four nodes.

A = [0 1 0 0; 1 0 0 0; 1 1 0 1; 0 1 1 0]

4×4 Matrix{Int64}:
 0  1  0  0
 1  0  0  0
 1  1  0  1
 0  1  1  0

The graphplot function makes a visual representation of this graph.

using Plots, GraphRecipes
graphplot(A, names=1:4, markersize=0.2, arrow=6)

Since this adjacency matrix is not symmetric, the edges are all directed, as indicated by the arrows. Here are the counts of all walks of length 3 in the graph:

A^3

4×4 Matrix{Int64}:
 0  1  0  0
 1  0  0  0
 3  2  0  1
 1  3  1  0

If the adjacency matrix is symmetric, the result is an undirected graph: all edges connect in both directions.

A = [0 1 1 0; 1 0 0 1; 1 0 0 0; 0 1 0 0]
graphplot(A, names=1:4, markersize=0.2)

Example 7.1.4

Here we create an adjacency matrix for a graph on four nodes.

A = [0 1 0 0; 1 0 0 0; 1 1 0 1; 0 1 1 0]

Since this adjacency matrix is not symmetric, the edges are all directed. We use digraph to create a directed graph.

G = digraph(A);
plot(G)

Here are the counts of all walks of length 3 in the graph:

A^3

If the adjacency matrix is symmetric, the result is an undirected graph: all edges connect in both directions.

A = [0 1 1 0; 1 0 0 1; 1 0 0 0; 0 1 0 0];
plot(graph(A))

A “buckyball” is an allotrope of carbon atoms with the same connection structure as a soccer ball.

plot(graph(bucky))

Example 7.1.4

Here we create an adjacency matrix for a graph on four nodes.

A = array([[0, 1, 0, 0], [1, 0, 0, 0], [1, 1, 0, 1], [0, 1, 1, 0]])
print(A)

[[0 1 0 0]
 [1 0 0 0]
 [1 1 0 1]
 [0 1 1 0]]

The networkx package has many functions for working with graphs. Here, we instruct it to create a directed graph from the adjacency matrix, then make a drawing of it.

import networkx as nx
G = nx.from_numpy_array(A, create_using=nx.DiGraph)
nx.draw(G, with_labels=True, node_color="yellow")

Here are the counts of all walks of length 3 in the graph:

print(linalg.matrix_power(A, 3))

[[0 1 0 0]
 [1 0 0 0]
 [3 2 0 1]
 [1 3 1 0]]

If the adjacency matrix is symmetric, the result is an undirected graph: all edges connect in both directions.

A = array([[0, 1, 1, 0], [1, 0, 0, 1], [1, 0, 0, 0], [0, 1, 0, 0]])
G = nx.from_numpy_array(A, create_using=nx.Graph)
nx.draw(G, with_labels=True, node_color="yellow")

The representation of a graph by its adjacency matrix opens up the possibility of many kinds of analysis of the graph. One might ask whether the nodes admit a natural partition into clusters, for example. Or one might ask to rank the nodes in order of importance to the network as determined by some objective criteria—an application made famous by Google’s PageRank algorithm, and one which is mathematically stated as an eigenvalue problem.

7.1.3Images as matrices¶

Computers typically represent images as rectangular arrays of pixels, each of which is colored according to numerical values for red (R), green (G), and blue (B) components of white light. Most often, these are given as integers in the range from zero (no color) to 255 (full color). Thus, an image that is $m$ -by- $n$ pixels can be stored as an $m$ -by- $n$ -by-3 array of integer values. In Julia, we can work with an $m\times n$ matrix of 3-vectors representing entire colors.

Example 7.1.5 (Images as matrices)

Julia

MATLAB

Python

Example 7.1.5

The Images package has many functions for image manipulation, and TestImages has some standard images to play with.

using Images, TestImages
img = testimage("mandrill")

The variable img is a matrix.

size(img)

(512, 512)

However, its entries are colors, not numbers.

img[100, 10]

You can use eltype to find out the type of the elements of any array.

eltype(img)

RGB{N0f8}

It’s possible to extract matrices of red, green, and blue intensities, scaled from 0 to 1.

R, G, B = red.(img), green.(img), blue.(img);
@show minB, maxB = extrema(B);

(minB, maxB) = extrema(B) = (0.0N0f8, 1.0N0f8)

Or we can convert the pixels to gray, each pixel again scaled from 0 to 1.

Gray.(img)

In order to do our usual operations, we need to tell Julia that we want to interpret the elements of the image matrix as floating-point values.

A = Float64.(Gray.(img))
A[1:4, 1:5]

4×5 Matrix{Float64}:
 0.568627  0.219608  0.192157  0.34902   0.537255
 0.454902  0.396078  0.156863  0.262745  0.352941
 0.301961  0.447059  0.180392  0.180392  0.384314
 0.278431  0.533333  0.372549  0.188235  0.32549

We can use Gray to reinterpret a matrix of floating-point values as grayscale pixels.

Gray.(reverse(A, dims=1))

Example 7.1.5

MATLAB ships with a few test images to play with.

A = imread('peppers.png');
color_size = size(A)

Use imshow to display the image.

imshow(A)

The image has three layers or channels for red, green, and blue. We can deal with each layer as a matrix, or (as below) convert it to a single matrix indicating shades of gray from black (0) to white (255). Either way, we have to explicitly convert the entries to floating-point values rather than integers.

A = im2gray(A);   % collapse from 3 dimensions to 2
gray_size = size(A)
imshow(A)

Before we can do any numerical computation, we need to convert the image to a matrix of floating-point numbers.

A = double(A);

Example 7.1.5

We will use a test image from the well-known scikit-image package.

from skimage import data as testimages
img = getattr(testimages, "coffee")()
imshow(img)

The variable img is a matrix.

size(img)

720000

However, its entries are colors, not numbers.

print(f"image has shape {img.shape}")
print(f"first pixel has value {img[0, 0]}")

image has shape (400, 600, 3)
first pixel has value [21 13  8]

The three values at each pixel are for intensities of red, green, and blue. We can convert each of those layers into an ordinary matrix of values between 0 and 255, which is maximum intensity.

R = img[:, :, 0]
print("upper left corner of the red plane is:")
print(R[:5, :5])
print(f"red channel values range from {R.min()} to {R.max()}")

upper left corner of the red plane is:
[[21 21 20 21 21]
 [21 21 20 22 20]
 [21 23 20 21 21]
 [21 21 21 21 20]
 [21 21 21 21 21]]
red channel values range from 0 to 255

It may also be convenient to convert the image to grayscale, which has just one layer of values from zero (black) to one (white).

from skimage.color import rgb2gray
A = rgb2gray(img)
A[:5, :5]
print("upper left corner of grayscale:")
print(A[:5, :5])
print(f"gray values range from {A.min()} to {A.max()}")

upper left corner of grayscale:
[[0.05623333 0.05651608 0.04978902 0.05708157 0.05903882]
 [0.05595059 0.05651608 0.05792275 0.05987216 0.05511725]
 [0.05875608 0.05846549 0.05848824 0.05651608 0.06212706]
 [0.05566784 0.05651608 0.05623333 0.05679882 0.05174627]
 [0.0553851  0.05595059 0.05510235 0.05595059 0.05623333]]
gray values range from 0.0002827450980392157 to 1.0

imshow(A, cmap='gray')
axis('off');

Some changes we make to the grayscale matrix are easy to interpret visually.

imshow(flipud(A), cmap='gray')
axis('off');

Representation of an image as a matrix allows us to describe some common image operations in terms of linear algebra. For example, in Singular value decomposition we will use the singular value decomposition to compress the information, and in Matrix-free iterations we will see how to apply and remove blurring effects.

7.1.4Exercises¶

Exercise 7.1.5

⌨ For this problem you need to download and import data:

Julia

MATLAB

Python

Download actors.mat by clicking the link and saving (you may need to fix the file name).

using MAT
A = matread("actors.mat")["A"]

Based on data provided by the Self-Organized Networks Database at the University of Notre Dame, the matrix A contains information about the appearances of 392,400 actors in 127,823 movies, as given by the Internet Movie Database. It has $A_{ij}=1$ if actor $j$ appeared in movie $i$ and zero elements elsewhere.

(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}$ ?

Preface

7. Matrix analysis

Preface

Eigenvalue decomposition