glmgen vs. genlasso

[muschellij2]

@mandymejia @damondpham
It looks like genlasso is being more maintained and is on CRAN: https://github.com/glmgen/genlasso (https://cran.r-project.org/web/packages/genlasso/index.html). It also has a trendfilter function, but I'm not sure if it's the same arguments and functionality. I might ask @glmgen or @statsmaths about it and which package should be used, if it's to be on CRAN or otherwise.

[damondpham]

Thanks Dr. Muschelli! I will ask Dr. McDonald (our contact for the PCATF method that uses trendfiltering).

The results for a given lambda look the same, but one difference is that the newer function only lets you set the minimum lambda whereas the older let you specify a specific lambda value. That might not be a problem, but I'll check with Dr. McDonald.

trendfilter_comparison.pdf

[taylor-arnold]

There is a major difference between the algorithms used in the two packages. genlasso computes a path solution, a theoretically exact method that computes solutions for all values of lambda greater than a particular cutoff; glmgen uses an approximation method for computing the solution at a particular value of lambda. If you are doing one-dimensional trend filtering I would definitely recommend using genlasso; it is more stable and works for much larger datasets. It also has the ability to do cross-validation baked into the package, and this is quite fast.

The glmgen package would be recommended if you are doing two-dimensional or graph-base trend filtering. It is hard to explain why without getting into the technical details. In short, the path solution has to compute a large number of change points before getting to a non-trivial solution in the case of two-dimensional filtering; this takes a long time and runs into numerical issues.

[damondpham]

Thanks for the information @statsmaths !

I tried replacing glmgen with genlasso, but the execution time of our method on our test datasets (~150 x 5000) increased from ~10 seconds to ~15 minutes. We are trying to create a PCA-like projection where the PC scores are trend-filtered (1-D and zero order), so computing each PC requires iterating between glmgen::trendfilter/genlasso::trendfilter and orthonormalization until convergence. Therefore, there are many calls to trendfilter.

For glmgen, the call to trendfilter looks like this:

tf.kwargs <- list(x = as.vector(X %*% v), nlambda = 1, lambda = lambda, k = 0)
u <- do.call(glmgen::trendfilter, tf.kwargs)$beta

And for genlasso:

tf.kwargs <- list(y = as.vector(X %*% v), ord = 0)
tf.i <- do.call(genlasso::trendfilter, tf.kwargs)
u <- coef(tf.i, lambda)$beta

Do you have any advice for reducing the execution time of genlasso? I tried approx=TRUE but that did not work. Or, should we stick to glmgen since it's much faster? (The results still look the identical for our use case.)

[taylor-arnold]

That's interesting that they return the same result. I can think of two improvements to genlasso that would probably help. First of all, you are currently returning all the points in the solution space currently (or the first 2000, whichever is less). If you set the minlam parameter to the lambda that you want the solution for, that may significantly speed up the execution of the trend filter function.

Also, it looks like you are only using the fused lasso (ord=0). I believe it would be much better to use genlasso::fusedlasso rather than the trend filter function with a manual ord parameter. This may significantly speed up the execution as well.

[damondpham]

Thank you!

I tried both suggestions:

tf.kwargs <- list(y = as.vector(X %*% v), minlam=lambda)
tf.i <- do.call(genlasso::fusedlasso1d, tf.kwargs)
u <- coef(tf.i, lambda)$beta

It decreased the execution time to ~11 minutes (and this remained the same when I added approx=TRUE), but this is still much slower than the ~10 seconds using glmgen::genlasso.

[damondpham]

I also confirmed that it's the initial call to fusedlasso1d which is taking a while, and not the coef call.