\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
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
At the moment vecpack knows how to pack scalars, vectors, matrices, arrays and cimg objects (images from the imager package).
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))
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 }
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
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 }
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
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.
You can use vecpack to express constraints as well (hiding them from the cost function). TODO: example with SPD matrices.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.