Convex.jl is a package in Julia for disciplined convex programming (DCP). Similar to CVX in MATLAB, Convex.jl makes it easy to solve general convex problems (for more details, see here). In 2011, I wrote the R package CVXfromR so that I could use CVX from R. I had assumed then that someone would very soon write a proper R version of CVX. In 2017, this finally happened – see the R package cvxr!
We begin by generating some data in R:
n <- 50 p <- 10 x <- matrix(rnorm(n * p), n, p) beta <- rnorm(p) y <- x %*% beta + 0.1 * rnorm(n)
Now suppose we would like to perform the Lasso using Convex.jl. In R, we begin by loading the convexjulia R package and identifying the location of Julia on our machine:
library(convexjulia) julia.call <- "full/path/to/julia"
For example, the path on OS X might be something like "/Applications/Julia-0.3.3.app/Contents/Resources/julia/bin/julia". We are now ready to send off an optimization problem to Convex.jl. The idea is to define a DCP problem called pr. In this example, we solve the lasso in “bound form” with an unpenalized intercept:
pr.def <- paste("pr = minimize(square(norm(y - b0 - x * b)))", "pr.constraints += norm(b, 1) <= s", sep = ";") lasso <- callconvex(opt.vars = "b = Variable(p); b0 = Variable(1)", pr.def = pr.def, const.vars = list(x = x, y = y, s = 70, p = p), opt.var.names = c("b", "b0"), julia.call = julia.call)
The output contains the optimal value, the status, the number of seconds elapsed for the call to Convex.jl (not counting the time opening Julia and transferring data between R and Julia), the code that was run in Julia, and an optimal point:
> names(lasso) [1] "optval" "status" "time_elapsed" "command" "b" [6] "b0"
Suppose we wish to solve this problem for a sequence of values of the bound parameter s. The function callconvex.varyparams makes this possible in a single call to Julia:
svalues <- seq(100, 65, length = 5) lasso <- callconvex.varyparams(opt.vars = "b = Variable(p); b0 = Variable(1)", pr.def = pr.def, const.vars = list(x = x, y = y, p = p), vary.param = list(s = svalues), opt.var.names = c("b", "b0"), julia.call = julia.call)
The returned value b is a p by 5 matrix, with each column corresponding to a value of the bound parameter s. Likewise, b0, optval, and status are vectors of length 5.
The latest version of the package is available on github.
My approach is very low-tech in the sense that it actually just forms a very long string that is passed to Julia through a call to the command line. It passes data from R to julia and then back to R by saving temporary files.