The Matrix R package is a fantastic tool, allowing the user to construct structured matrices of different types (triangular, symmetric, sparse, etc). A common set of operations are offered for all types (\(+\), \(-\), `solve`

, …), meaning differently typed matrices can all be treated the same. Importantly, the method of storage for each matrix type takes advantage of its structure, as do the operations on these matrices.

Below I construct 4 \(d\times d\) matrices of different types: general dense, positive semidefinite, sparse and symmetric sparse.

```
library(Matrix)
d <- 1000
A <- Matrix(rnorm(d*d), nrow=d, ncol=d) # Dense matrix
Asym <- as(crossprod(A), "dppMatrix") # Symmetric, positive semidefinite matrix
Asp <- drop0(A, tol=2) # Sparse matrix
Aspsym <- drop0(Asym, tol=50) # Sparse symmetric matrix
```

`crossprod`

is from the `Matrix`

package; this works similarly to the `base::crossprod`

function but returns a symmetric matrix class. I have used `as( , "dppMatrix")`

to cast this to a *packed* positive-semidefinite matrix, so only one triangle is stored. `drop0`

is used to remove any elements with absolute value below `tol`

. This code produces `Asp`

and `Aspsym`

which are 4.6% and 11.5% sparse respectively.

The classes of the matrices reflect their structure.

` vapply(list(A, Asym, Asp, Aspsym), class, character(1))`

`## [1] "dgeMatrix" "dppMatrix" "dgCMatrix" "dsCMatrix"`

The leading `d`

in each case shows it is a matrix of doubles. The following two letters give us the structure: `ge`

tells us it is a *general* matrix; `pp`

tells us that the matrix is positive-semidefinite in packed storage (only one triangle stored); `gC`

is for a general, column-compressed sparse matrix, and `sC`

indicates a symmetric column-compressed sparse matrix. Only `A`

stores all \(d^2\) elements, so has the largest memory footprint.

A | Asym | Asp | Aspsym |
---|---|---|---|

8 MB | 4.01 MB | 553 kB | 704 kB |

My favourite feature is the very general `Matrix::solve`

method. If we have matrices \(\mathbf{A}\) and \(\mathbf{B}\), `X <- solve(A, B)`

solves \(\mathbf{AX} = \mathbf{B}\). This works very similarly to `base::solve`

, but solves for \(\mathbf{X}\) using a decomposition appropriate for the type of \(\mathbf{A}\). Furthermore, whenever a decomposition is performed, this is saved in the `@factors`

slot of the original matrix. This decomposition is then reused whenever a decomposition would be performed, such as solving a system, or calculating a determinant or condition number.

```
B <- Matrix(rnorm(d*10), ncol=10)
X1 <- solve(A, B)
X2 <- solve(Asym, B)
X3 <- solve(Asp, B)
X4 <- solve(Aspsym, B)
```

It is interesting to inspect the saved decompositions, which I do below. For more detail, consult the documentation on the Matrix CRAN page.

` str(A@factors)`

```
## List of 1
## $ LU:Formal class 'denseLU' [package "Matrix"] with 4 slots
## .. ..@ x : num [1:1000000] 3.8103 0.4652 0.1406 -0.1984 0.0461 ...
## .. ..@ perm : int [1:1000] 495 442 460 694 765 540 108 67 216 600 ...
## .. ..@ Dimnames:List of 2
## .. .. ..$ : NULL
## .. .. ..$ : NULL
## .. ..@ Dim : int [1:2] 1000 1000
```

We can see that for the general dense matrix, an LU decomposition was performed, which is fairly standard. Looking at the decomposition for the symmetric dense matrix, we see it has type `pCholesky`

: a packed Cholesky decomposition (does not store the triangle of zeroes). This is good behaviour; we do not need to worry about the most numerically stable method of solving our system, since `Matrix`

figures it out for us.

` str(Asym@factors)`

```
## List of 1
## $ pCholesky:Formal class 'pCholesky' [package "Matrix"] with 5 slots
## .. ..@ x : num [1:500500] 32.713 0.216 32.874 1.603 0.716 ...
## .. ..@ Dim : int [1:2] 1000 1000
## .. ..@ Dimnames:List of 2
## .. .. ..$ : NULL
## .. .. ..$ : NULL
## .. ..@ uplo : chr "U"
## .. ..@ diag : chr "N"
```

More interesting decompositions appear in the sparse matrices. Starting with the symmetric matrix, we wee a decomposition of type `sparseLU`

.

` str(Asp@factors)`

```
## List of 1
## $ LU:Formal class 'sparseLU' [package "Matrix"] with 5 slots
## .. ..@ L :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
## .. .. .. ..@ i : int [1:462576] 0 7 22 27 46 49 51 120 136 158 ...
## .. .. .. ..@ p : int [1:1001] 0 60 115 167 226 273 308 371 416 463 ...
## .. .. .. ..@ Dim : int [1:2] 1000 1000
## .. .. .. ..@ Dimnames:List of 2
## .. .. .. .. ..$ : NULL
## .. .. .. .. ..$ : NULL
## .. .. .. ..@ x : num [1:462576] 1 0.617 0.53 -0.581 -0.589 ...
## .. .. .. ..@ uplo : chr "L"
## .. .. .. ..@ diag : chr "N"
## .. ..@ U :Formal class 'dtCMatrix' [package "Matrix"] with 7 slots
## .. .. .. ..@ i : int [1:469581] 0 1 2 3 4 5 6 7 8 8 ...
## .. .. .. ..@ p : int [1:1001] 0 1 2 3 4 5 6 7 8 9 ...
## .. .. .. ..@ Dim : int [1:2] 1000 1000
## .. .. .. ..@ Dimnames:List of 2
## .. .. .. .. ..$ : NULL
## .. .. .. .. ..$ : NULL
## .. .. .. ..@ x : num [1:469581] 3.81 3.64 -3.54 -3.21 -3.67 ...
## .. .. .. ..@ uplo : chr "U"
## .. .. .. ..@ diag : chr "N"
## .. ..@ p : int [1:1000] 494 294 459 679 261 903 107 484 215 202 ...
## .. ..@ q : int [1:1000] 0 1 2 3 4 5 6 7 8 9 ...
## .. ..@ Dim: int [1:2] 1000 1000
```

This is similar to the regular LU decomposition \(\mathbf{A} = \mathbf{LU}\), but has two new pivot matrices \(\mathbf{P}\), \(\mathbf{Q}\). The sparse LU decomposes \(\mathbf{A} = \mathbf{P}^\top\mathbf{LUQ}\), where the pivots are chosen such that \(\mathbf{L}\) and \(\mathbf{U}\) are sparse. This can be seen in the slots `@L`

, `@U`

. Such decomposition-sparsity is neccesary when dealing with sparse matrices; if we required a dense decomposition of the same size as the original matrix, then we would lose all benefit of using a sparse matrix in the first place.

Finally, the class of the decomposition for `Aspsym`

is `dCHMsimpl`

. According to the Matrix reference manual, this inherits from a more general class of ‘’CHOLMOD-based Cholesky factorizations of symmetric, sparse, compressed, column-oriented matrices’’.

In overview, I really like the matrix package because it becomes possible to write a single method that can deal with many matrix types without detailed knowledge of large numbers of matrix decompositions. Even without matrix inversions, I have experienced significant speed increases just from storing symmetric matrices as their appropriate `Matrix`

class, since any addition/multiplication only needs to consider a single triangle. Methods `crossprod`

and `symmpart`

are also lifesavers when constructing symmetric matrices, as it means avoiding going through large general dense matrices during construction.