In dahtah/vecpack: Helps Turn Objects and Lists into Vectors, to Facilitate Interaction with Optimisers and Samplers

\newcommand{\bm}[1]{\boldsymbol{\mathbf{#1}}}

library(vecpack)
library(purrr)

vecpack solves a problem commonly encountered when working with optimisers (like optim) and MCMC samplers: they expect to work on a function f that takes a vector as input. This is fine when optimising over vectors, but sometimes f is really a function of a matrix and a vector, or a matrix and an image and a scalar, etc. You could concatenate all these arguments into a single vector, but that always involves lots of boilerplate and book-keeping, and it's easy to make errors. vecpack saves you from having to write most of the boilerplate.

Packing a list of values into a single vector is easy enough:

l <- list(a= matrix(1:4,2,2),b = 0.4)
sapply(l,as.vector) %>% do.call(c,.)

In vecpack, that functionality is provided by the vpack function:

vpack(l)

Going the other way is more complicated: e.g., the vector 1,2,3,4 could be a 2x2 matrix followed by a scalar, or one scalar followed by a 2x2 matrix, or two vectors of length 2 and 3, etc. To invert the packing, we need to know the structure, and in vecpack the structure is inferred from an example:

vunpack <- gen.vunpack(l)
vunpack(1:5)
vpack(l) %>% vunpack

Using vecpack with optim

The first example is somewhat silly but hopefully you'll see where I'm going with this: let's say we're doing linear regression. The data are pairs (x,y), our model is $$ E(y) = a x + b $$ and the corresponding log-likelihood is "cfun", below:

x <- seq(0,1,l=20)
y <- 3*x - 2 +rnorm(length(x))
cfun <- function(a,b) sum((y-(a*x+b))^2)

We'll use optim and vecpack to optimise cfun:

guess <- list(a=0,b=1)
up <- gen.vunpack(guess) #Unpacking function
cfun.lifted <- lift_up(cfun,up) #Now cfun accepts a vector argument
opt <- optim(vpack(guess),cfun.lifted)
up(opt$par)

vecpack includes an interface to optim called vpoptim, which removes all of the boilerplate: just provide a cost function, a guess and you're set.

opt <- vpoptim(guess,cfun)
opt$par

Here's a better example: we seek to find a low-rank + diagonal decomposition of a matrix, i.e.

$$ \bm{X} \approx \bm{A}\bm{A}^t + \sigma^2 \bm{I} $$

X <- cov(USArrests)
cfun <- function(A,sigma2)
{
    lrank <- tcrossprod(A,A) + sigma2*diag(nrow(A)) 
    sum((X-lrank)^2)
}
guess <- list(A=matrix(0,4,2),sigma2=1)
opt <- vpoptim(guess,cfun)
opt$par

Which objects can vecpack pack?

At the moment vecpack knows how to pack scalars, vectors, matrices, arrays and cimg objects (images from the imager package).

Computing numerical gradients

There's an interface to the function "grad" from the numDeriv package, that works along the same lines as vpoptim:

#gradient of a function of two arguments
vpgrad(function(x,y) x+y,list(x=1,y=3))
#gradient of a matrix function, returned as a matrix
vpgrad(det,diag(2))

Extending vecpack

vecpack is extensible: you can define your own ways of packing objects into vectors. For example, suppose we want to pack tri-diagonal matrices into vectors. All we need to do is create an S3 class for tridiagonal matrices, and define just three generics: as.vector and reshapefun.

First, some utility functions to get and set values:

#Extract or set values above diagonal
utri <- function(M)
{
    M[row(M) == col(M) - 1]
}

`utri<-` <- function(M,value)
{
    M[row(M) == col(M) - 1] <- value
    M
}

#Extract or set values below
ltri <- function(M)
{
    M[row(M) == col(M) + 1]
}

`ltri<-` <- function(M,value)
{
    M[row(M) == col(M) + 1] <- value
    M
}

First generic: as.vector

as.vector should take a (tri-diagonal) matrix and pack it into a vector:

as.vector.trid <- function(M,...)
{
    c(ltri(M),diag(M),utri(M))
}

trid <- function(M)
{
    class(M) <- c('trid','matrix')
    M
}

diag(3) %>% trid %>% as.vector

Second generic: length

length should return the length of the tridiagonal matrix when packed as a vector: a tridiagonal matrix has $n + 2*(n-1) = 3n - 2$ degrees of freedom.

length.trid <- function(V)
{
    n <- nrow(V)
    3*n-2
}

Third generic: reshapefun

And next we need to define "reshapefun". Reshaping turns the vector back into the original object. "reshapefun.trid" should return a reshaping function appropriate for a given matrix. Here's reshapefun.matrix, for instance:

reshapefun.matrix <- function(x,...) function(k) matrix(k,nrow(x),ncol(x))
rs <- matrix(0,2,2) %>% reshapefun
rs(1:4)

We define reshapefun.trid along the same lines:

reshapefun.trid <- function(x,...)
{
    n <- nrow(x)
    #The following function returns a matrix of the right size
    function(k)
    {
        M <- matrix(0,n,n)
        ltri(M) <- k[1:(n-1)]
        diag(M) <- k[n+0:(n-1)]
        utri(M) <- k[(2*n):length(k)]
        M
    }
}
V <- diag(3) %>% trid
rs <- reshapefun(V)
as.vector(V) %>% rs

We're now good to go with vpack:

l <- list(V=V,a=3)
up <- gen.vunpack(l)
vpack(l) %>% up

We can now use it to find, e.g., a tri-diagonal approximate inverse:

#The cost function encourages A*X = Id
X <- rnorm(16) %>% matrix(4,4)
cfun <- function(A)
{
    P <- A%*%X
    sum(abs(eigen(P)$val- 1))
}
A <- diag(4) %>% trid
opt <- vpoptim(A,cfun,method="BFGS",control=list(trace=T,maxit=1e5))
opt$par

Making your own wrapper

If you'd like to wrap another optimiser/MCMC sampler, it should be easy enough. The source for vpoptim is just a few lines:

vpoptim <- function(guess,cfun,...)
{
    up <- gen.vunpack(guess) 
    cfun.lifted <- lift_up(cfun,up) 
    opt <- optim(vpack(guess),cfun.lifted,...)
    opt$par <- up(opt$par)
    opt
}

If you've made another interface, consider contributing it to the vecpack package.

Other uses for vecpack

You can use vecpack to express constraints as well (hiding them from the cost function). TODO: example with SPD matrices.

dahtah/vecpack documentation built on May 14, 2019, 3:27 p.m.