Here I will introduce **matpow**, a package to flexibly and conveniently compute matrix powers. But even if you are not interested in matrices, I think many of you will find that this post contains much general material on R that you’ll find useful. Indeed, most of this post will be about general R issues, not so much about matrices *per se*. So, bear with me.

**Why matrix powers? **

Sadly, most university beginning linear algebra courses say very little about why the material is important. Worse, if some motivation is presented, it tends to be physics-oriented. Indeed, the current trend in U.S. universities is to fold the introductory linear algebra curriculum into the differential equations course.

Most readers of this blog know better, as they are (to various extents) aware of the major roles that matrices play in statistics, e.g. in linear models and principal components analysis. But did you know that matrices play a major role in the analysis of social networks? We’re not in Physicsland anymore, Toto.

Matrix powers are useful here. For instance, suppose we wish to determine whether a network (or *graph*) is *connected, *meaning that every node leads to every other node in one or more steps. (The famous Six Degrees of Separation notion comes from this context.) It turns out that can be determined by taking powers of the adjacency matrix A of the network, whose (i,j) element is 1 if there is a one-step link from i to j, 0 otherwise. (For technical reasons, we will set the diagonal of A to all 1s.) If some power of A has all its elements nonzero, then the graph is connected; it’s disconnected if this situation is never reached. (One need calculate only up to the power n-1, where n is the number of rows and columns of the matrix.) Moreover, eigenanalysis of A yields other important information about the network, and one way to compute eigenvectors of a matrix involves finding its powers too.

**But why have a package to compute the powers?**

Isn’t this already straightforward in R? For instance, the third power of A is computed via the code

` m %*% m %*% m`

and higher powers can be coded with a simple loop. But we needed to be much more flexible and versatile than this in developing **matpow**.

**1. Generality: ** Our **matpow** package accommodates various classes of matrices.

There are many of these in the R world. In addition to the basic **“matrix”** class, there is also the **“Matrix”** class (from the **Matrix** package), the **“big.matrix”** class (from **bigmemory**) and so on. Syntax can vary; with **bigmemory**, for instance, one must use brackets in assignment, e.g.

` m1[,]`

And most important, not all of these use **%*%** as their multiplication operator. For example, in **gputools **one uses the function **gpuMatMult()** for this purpose.

**2. Accommodation of parallel computation:** Today’s matrices can be huge, so parallel computation would be nice. Our **matpow** package does not itself include facilities for parallel computation, but it is designed to integrate well with external parallel methods, as mentioned above for **gputools**, an R package that interfaces to GPUs.

**3. Capability to include callback functions,** to perform application-specific operations after each new power is calculated.

For instance, consider the social network example again. As we compute higher and higher powers of A, we want to check for an “all non zeroes” condition after each iteration; if the condition is found to hold, there is no point in doing any more iterations. We can have our callback function check for this.

And there is more. The (i,j) element in the r^{th }power of A can be shown to be the number of ways to go from node i to node j in r steps. If that element is now positive but was 0 in power r-1 of A, we now know that the shortest path from i to j consists of r steps. So, we can have our callback function watch for changes from zero to nonzero, and record them, and thus have the overall internode distance matrix in the end.

**A few quick examples:**

Here is a quick tour, including of the **squaring** option. There we square the input matrix, then square the result and so on, thus requiring only a logarithmic number of steps to reach our desired power. This is much faster if we don’t need the intermediate powers.

`> m <- rbind(1:2,3:4)`

> ev <- matpow(m,16) # find m to 16th power

> ev$prod1 # here it is

[,1] [,2]

[1,] 115007491351 1.67615e+11

[2,] 251422553235 3.66430e+11

> ev$i # last iteration was 15th

[1] 15

> ev <- matpow(m,16,squaring=TRUE)

> ev$prod1 # same result

[,1] [,2]

[1,] 115007491351 1.67615e+11

[2,] 251422553235 3.66430e+11

> ev$i # but with only 4 iterations

[1] 4

`# test network connectivity`

> m <- rbind(c(1,0,0,1),c(1,0,1,1),c(0,1,0,0),c(0,0,1,1))

> ev <- matpow(m,callback=cgraph,mindist=T)

> ev$connected

[1] TRUE

> ev$dists

[,1] [,2] [,3] [,4]

[1,] 1 3 2 1

[2,] 1 1 1 1

[3,] 2 1 1 2

[4,] 3 2 1 1

**How can we deal with different matrix classes/multiplication syntaxes?**

As noted, one problem we needed to deal with in developing **matpow** was how to accommodate diverse matrix classes and multiplication syntaxes. We solved that problem by using R’s **eval()** and **parse()** functions. Consider the following toy example:

` > x <- 28`

> s <- "x <- 16"

> eval(parse(text=s))

> x

[1] 16

Of course, here it is just a very roundabout way to set **x** to 16. But for **matpow**, it’s just what we need. We want our multiplication command in string form to be something like **prod <- m %*% m** in the **“matrix”** case, **prod[,] <- m[,] %*% m[,] **in the **“big.matrix”** (**bigmemory**) case, and so on.

The (built-in or user-supplied) argument **genmulcmd** (“generate multiplication command”) to **matpow() **handles this. For instance, here is the built-in one for **gputools**:

```
genmulcmd.gputools <- function (a, b, c)
paste(c, " <- gpuMatMult(", a, ",", b, ")")
```

We then can plug the string into **eval()** and **parse()**. In this manner, we can switch matrix types by changing not even a full line of code, just an argument to **matpow()**.

**How can we implement callbacks?**

Recall our example **matpow()** call:

```
matpow(m,callback=cgraph,mindist=T)
```

Here we are specifying that there is a callback function (the default is NULL), and it is to be **cgraph()**. This happens to be built-in for the package, but it could be user-written.

The key issue is how data will be communicated between **matpow()** and the callback–in BOTH directions, i.e. we need the callback to write information back to **matpow()**. This is not easy in R, which as a *functional language* tries to avoid *side effects*.

Consider for example this:

` > x <- sample(1:20,8)`

> x

[1] 5 12 9 3 2 6 16 18

> sort(x)

[1] 2 3 5 6 9 12 16 18

> x

[1] 5 12 9 3 2 6 16 18

The point is that **x** didn’t change, i.e. we had no side effect to the argument **x**. If we do want **x** to change, we must write

`> x <- sort(x)`

But environments are different. The **matpow()** function maintains an environment **ev**, which it passes to the callback as an argument. Though **ev** is like an R list, and its components are accessed via the familiar **$** operator, the difference is that the callback can change those components (a side effect), thus communicate information back to **matpow()**.

For instance, look at this code from **cgraph()**:

```
if (all(prd > 0)) {
ev$connected <- TRUE
ev$stop <- TRUE
}
```

We’ve discovered that the network is connected, so we set the proper two components of **ev**. Back in **matpow()**, the code will sense that **ev$stop** is TRUE, and cease the iteration process. Since **matpow()** uses **ev** as its return value, the user then can see that the network is connected.

**Other issues:**

Again, keep in mind that **matpow()** will use whatever form of matrix multiplication you give it. If you give it a parallel form of multiplication, then **matpow()**‘s computation will be parallel. Or if your R is configured with a multicore-capable BLAS, such as OpenBLAS, you again will be computing in parallel.

Though we have gotten good performance in this context with **gputools**, it should be noted that **gpuMatMult()** returns the product to the R caller. This causes unnecessary copying, which on GPUs can sap speed. Code that maintains the product on the GPU from iteration to iteration would be best.

**Where to get matpow():**

We plan submission to CRAN, but for now download it from here. Unzip the package, then run R CMD INSTALL on the directory that is created as a result of unzipping.

Just curious: why should it be “worse” if an application is physics-oriented? Matrices are useful in physics, and without physics, you would not even have R (physics being absolutely essential for computers and all).

I don’t disagree with anything you said. I hope you also agree with my point that it’s wrong for university curricula to give students the impression that the ONLY use for matrices is in physics.