yixuan / admm Goto Github PK

Solving Statistical Optimization Problems Using the ADMM Algorithm

R 14.41% C++ 85.54% C 0.04%

admm's Introduction

Introduction

ADMM is an R package that utilizes the Alternating Direction Method of Multipliers (ADMM) algorithm to solve a broad range of statistical optimization problems. Presently the models that ADMM has implemented include Lasso, Elastic Net, Least Absolute Deviation and Basis Pursuit.

Installation

ADMM package is still experimental, so it has not been submitted to CRAN yet. To install ADMM, you need a C++ compiler such as g++ or clang++, and for Windows users, the Rtools software is needed (unless you can configure the toolchain by yourself).

The installation follows the typical way of R packages on Github:

library(devtools)
install_github("yixuan/ADMM")

In order to achieve high performance computation, the package uses single precision floating point type (aka float) in some part of the code, which uses single precision BLAS functions such as ssyrk(). However, the BLAS library shipped with R does not contain such functions, so for code portability, a macro called NO_FLOAT_BLAS is defined in src/Makevars to use fallback Eigen implementation.

If your R is linking to a high performance BLAS such as OpenBLAS, it is suggested to remove that macro in order to use the BLAS implementation. Moreover, the Eigen library that ADMM pakcage depends on benefits from modern CPU's vectorization support, so it is also recommended to define -msse4 or -mavx in the PKG_CXXFLAGS variable (in file src/Makevars) whenever applicable.

Models

Lasso

library(glmnet)
library(ADMM)
set.seed(123)
n <- 100
p <- 20
m <- 5
b <- matrix(c(runif(m), rep(0, p - m)))
x <- matrix(rnorm(n * p, mean = 1.2, sd = 2), n, p)
y <- 5 + x %*% b + rnorm(n)

fit <- glmnet(x, y)
out_glmnet <- coef(fit, s = exp(-2), exact = TRUE)
out_admm <- admm_lasso(x, y)$penalty(exp(-2))$fit()
out_paradmm <- admm_lasso(x, y)$penalty(exp(-2))$parallel()$fit()

data.frame(glmnet = as.numeric(out_glmnet),
           admm = as.numeric(out_admm$beta),
           paradmm = as.numeric(out_paradmm$beta))

##          glmnet         admm      paradmm
## 1   5.357410774  5.357455254  5.357429504
## 2   0.178916019  0.178915471  0.178917870
## 3   0.683606818  0.683609307  0.683610320
## 4   0.310518550  0.310507625  0.310525119
## 5   0.861034415  0.861029863  0.861012816
## 6   0.879797912  0.879794598  0.879801810
## 7   0.007854581  0.007850002  0.007853498
## 8   0.000000000  0.000000000  0.000000000
## 9   0.000000000  0.000000000  0.000000000
## 10  0.023462980  0.023467677  0.023452930
## 11  0.010952896  0.010957017  0.010950469
## 12  0.000000000  0.000000000  0.000000000
## 13 -0.003800159 -0.003811116 -0.003801103
## 14  0.000000000  0.000000000  0.000000000
## 15  0.094591923  0.094586611  0.094600648
## 16  0.000000000  0.000000000  0.000000000
## 17  0.000000000  0.000000000  0.000000000
## 18  0.000000000  0.000000000  0.000000000
## 19  0.000000000  0.000000000  0.000000000
## 20 -0.002916255 -0.002929136 -0.002919935
## 21  0.000000000  0.000000000  0.000000000

Elastic Net

fit <- glmnet(x, y, alpha = 0.5)
out_glmnet <- coef(fit, s = exp(-2), exact = TRUE)
out_admm <- admm_enet(x, y)$penalty(exp(-2), alpha = 0.5)$fit()
data.frame(glmnet = as.numeric(out_glmnet),
           admm = as.numeric(out_admm$beta))

##          glmnet         admm
## 1   5.150556538  5.150497437
## 2   0.204543779  0.204526767
## 3   0.705652674  0.705665767
## 4   0.330650192  0.330640256
## 5   0.872594728  0.872611761
## 6   0.884433876  0.884422064
## 7   0.048044107  0.048055928
## 8   0.025072878  0.025097074
## 9   0.000000000  0.000000000
## 10  0.057804317  0.057830613
## 11  0.041853068  0.041876025
## 12 -0.004476248 -0.004499977
## 13 -0.035255637 -0.035279647
## 14  0.000000000  0.000000000
## 15  0.110919341  0.110915266
## 16  0.000000000  0.000000000
## 17  0.000000000  0.000000000
## 18  0.000000000  0.000000000
## 19  0.000000000  0.000000000
## 20 -0.021003756 -0.020984368
## 21  0.000000000  0.000000000

Least Absolute Deviation

Least Absolute Deviation (LAD) minimizes ||y - Xb||_1 instead of ||y - Xb||_2^2 (OLS), and is equivalent to median regression.

library(quantreg)
out_rq <- rq.fit(x, y)
out_admm <- admm_lad(x, y, intercept = FALSE)$fit()

data.frame(rq_br = out_rq$coefficients,
           admm = out_admm$beta[-1])

##           rq_br          admm
## 1   0.463871497  0.4630289961
## 2   0.829243353  0.8324149339
## 3   0.151432833  0.1493799430
## 4   1.074107564  1.0707590072
## 5   0.958979798  0.9569585188
## 6   0.502539859  0.5028832829
## 7   0.337640338  0.3360263689
## 8   0.209127703  0.2120946512
## 9   0.361765382  0.3630356485
## 10  0.323168985  0.3217875563
## 11 -0.002009264  0.0007319653
## 12 -0.036099511 -0.0370447075
## 13  0.328007777  0.3290499302
## 14  0.296038071  0.2992857234
## 15  0.310187867  0.3117528782
## 16  0.071713681  0.0711670377
## 17  0.166827429  0.1622600454
## 18  0.260366502  0.2580854533
## 19  0.324487629  0.3251952295
## 20  0.209758565  0.2131039214

Basis Pursuit

set.seed(123)
n <- 50
p <- 100
nsig <- 15
beta_true <- c(runif(nsig), rep(0, p - nsig))
beta_true <- sample(beta_true)

x <- matrix(rnorm(n * p), n, p)
y <- drop(x %*% beta_true)
out_admm <- admm_bp(x, y)$fit()

range(beta_true - out_admm$beta)

## [1] -0.0006052779  0.0004780069

Performance

Lasso and Elastic Net

library(microbenchmark)
library(ADMM)
library(glmnet)
# compute the full solution path, n > p
set.seed(123)
n <- 10000
p <- 1000
m <- 100
b <- matrix(c(runif(m), rep(0, p - m)))
x <- matrix(rnorm(n * p, sd = 2), n, p)
y <- x %*% b + rnorm(n)

lambdas1 = glmnet(x, y)$lambda
lambdas2 = glmnet(x, y, alpha = 0.6)$lambda

microbenchmark(
    "glmnet[lasso]" = {res1 <- glmnet(x, y)},
    "admm[lasso]"   = {res2 <- admm_lasso(x, y)$penalty(lambdas1)$fit()},
    "padmm[lasso]"  = {res3 <- admm_lasso(x, y)$penalty(lambdas1)$parallel()$fit()},
    "glmnet[enet]"  = {res4 <- glmnet(x, y, alpha = 0.6)},
    "admm[enet]"    = {res5 <- admm_enet(x, y)$penalty(lambdas2, alpha = 0.6)$fit()},
    times = 5
)

## Unit: milliseconds
##           expr      min        lq      mean    median        uq       max neval
##  glmnet[lasso] 939.0194  943.7227 1005.1232 1043.2826 1048.1939 1051.3973     5
##    admm[lasso] 320.1375  320.7603  321.6413  321.0283  321.0635  325.2172     5
##   padmm[lasso] 486.7478  510.7104  529.9721  512.5421  567.0847  572.7757     5
##   glmnet[enet] 950.9513 1048.1872 1031.4437 1049.9402 1051.2483 1056.8917     5
##     admm[enet] 286.9177  288.8735  289.0231  288.9758  289.1651  291.1834     5

# difference of results
diffs = matrix(0, 3, 2)
rownames(diffs) = c("glmnet-admm [lasso]", "glmnet-padmm[lasso]", "glmnet-admm [enet]")
colnames(diffs) = c("min", "max")
diffs[1, ] = range(coef(res1) - res2$beta)
diffs[2, ] = range(coef(res1) - res3$beta)
diffs[3, ] = range(coef(res4) - res5$beta)
diffs

##                               min          max
## glmnet-admm [lasso] -0.0002873333 7.259293e-05
## glmnet-padmm[lasso] -0.0005554722 7.382258e-05
## glmnet-admm [enet]  -0.0002195360 8.176991e-05

# p > n
set.seed(123)
n <- 1000
p <- 2000
m <- 100
b <- matrix(c(runif(m), rep(0, p - m)))
x <- matrix(rnorm(n * p, sd = 2), n, p)
y <- x %*% b + rnorm(n)

lambdas1 = glmnet(x, y)$lambda
lambdas2 = glmnet(x, y, alpha = 0.6)$lambda

microbenchmark(
    "glmnet[lasso]" = {res1 <- glmnet(x, y)},
    "admm[lasso]"   = {res2 <- admm_lasso(x, y)$penalty(lambdas1)$fit()},
    "padmm[lasso]"  = {res3 <- admm_lasso(x, y)$penalty(lambdas1)$parallel()$fit()},
    "glmnet[enet]"  = {res4 <- glmnet(x, y, alpha = 0.6)},
    "admm[enet]"    = {res5 <- admm_enet(x, y)$penalty(lambdas2, alpha = 0.6)$fit()},
    times = 5
)

## Unit: milliseconds
##           expr       min        lq      mean    median        uq       max neval
##  glmnet[lasso]  197.9279  198.2767  199.1576  199.3774  200.0559  200.1499     5
##    admm[lasso]  230.6332  237.1298  245.8944  247.4240  250.6257  263.6596     5
##   padmm[lasso] 5159.2198 5170.7308 5513.3797 5345.6322 5426.8846 6464.4313     5
##   glmnet[enet]  195.7102  196.7432  197.3977  197.7577  198.3421  198.4355     5
##     admm[enet]  225.8397  239.6536  247.2756  249.9269  252.2080  268.7499     5

# difference of results
diffs[1, ] = range(coef(res1) - res2$beta)
diffs[2, ] = range(coef(res1) - res3$beta)
diffs[3, ] = range(coef(res4) - res5$beta)
diffs

##                              min         max
## glmnet-admm [lasso] -0.001518947 0.002055109
## glmnet-padmm[lasso] -0.001898237 0.002052009
## glmnet-admm [enet]  -0.001615556 0.001948477

LAD

library(ADMM)
library(quantreg)

set.seed(123)
n <- 1000
p <- 500
b <- runif(p)
x <- matrix(rnorm(n * p, sd = 2), n, p)
y <- x %*% b + rnorm(n)

microbenchmark(
    "quantreg[br]" = {res1 <- rq.fit(x, y)},
    "quantreg[fn]" = {res2 <- rq.fit(x, y, method = "fn")},
    "admm"         = {res3 <- admm_lad(x, y, intercept = FALSE)$fit()},
    times = 5
)

## Unit: milliseconds
##          expr        min         lq       mean     median         uq
##  quantreg[br] 2420.34862 2422.79215 2493.76695 2426.54150 2429.49592
##  quantreg[fn]  451.87866  452.94572  454.71406  453.56378  455.81717
##          admm   50.55354   51.17922   51.69386   51.59099   52.32357
##         max neval
##  2769.65653     5
##   459.36498     5
##    52.82196     5

# difference of results
range(res1$coefficients - res3$beta[-1])

## [1] -0.006989109  0.006061505

set.seed(123)
n <- 5000
p <- 1000
b <- runif(p)
x <- matrix(rnorm(n * p, sd = 2), n, p)
y <- x %*% b + rnorm(n)

microbenchmark(
    "quantreg[fn]" = {res1 <- rq.fit(x, y, method = "fn")},
    "admm"         = {res2 <- admm_lad(x, y, intercept = FALSE)$fit()},
    times = 5
)

## Unit: seconds
##          expr      min       lq     mean   median       uq      max neval
##  quantreg[fn] 6.156911 6.231430 7.686811 6.280464 9.861437 9.903813     5
##          admm 2.184311 2.187431 2.193200 2.189173 2.202042 2.203044     5

# difference of results
range(res1$coefficients - res2$beta[-1])

## [1] -0.003577610  0.004135838

Basis Pursuit

set.seed(123)
n <- 1000
p <- 2000
nsig <- 100
beta_true <- c(runif(nsig), rep(0, p - nsig))
beta_true <- sample(beta_true)
x <- matrix(rnorm(n * p), n, p)
y <- drop(x %*% beta_true)

system.time(out_admm <- admm_bp(x, y)$fit())

##    user  system elapsed
##   0.996   0.169   0.292

range(beta_true - out_admm$beta)

## [1] -0.001267782  0.002108828

set.seed(123)
n <- 1000
p <- 10000
nsig <- 200
beta_true <- c(runif(nsig), rep(0, p - nsig))
beta_true <- sample(beta_true)
x <- matrix(rnorm(n * p), n, p)
y <- drop(x %*% beta_true)

system.time(out_admm <- admm_bp(x, y)$fit())

##    user  system elapsed
##  19.315   0.573   4.969

range(beta_true - out_admm$beta)

## [1] -0.1575968  0.3361001

admm's People

Contributors

Stargazers

Watchers

admm's Issues

Fail to install on Centos

I am trying to install ADMM on my centos server, it fails too.

Here is the bug.

* installing *source* package ?.DMM?....
** libs
g++ -I/usr/local/R/lib64/R/include -DNDEBUG  -I/usr/local/include -I"/usr/local/R/lib64/R/library/Rcpp/include" -I"/usr/local/R/lib64/R/library/RcppEigen/include"  -fopenmp -fpic  -g -O2  -c BP.cpp -o BP.o
g++ -I/usr/local/R/lib64/R/include -DNDEBUG  -I/usr/local/include -I"/usr/local/R/lib64/R/library/Rcpp/include" -I"/usr/local/R/lib64/R/library/RcppEigen/include"  -fopenmp -fpic  -g -O2  -c Enet.cpp -o Enet.o
In file included from Spectra/SymEigsSolver.h:20:0,
                 from ADMMLassoTall.h:6,
                 from ADMMEnet.h:4,
                 from Enet.cpp:3:
Spectra/LinAlg/TridiagEigen.h: In instantiation of ?.oid Spectra::TridiagEigen<Scalar>::compute(Spectra::TridiagEigen<Scalar>::ConstGenericMatrix&) [with Scalar = float; Spectra::TridiagEigen<Scalar>::ConstGenericMatrix = const Eigen::Ref<const Eigen::Matrix<float, -1, -1> >; typename Eigen::internal::conditional<const Eigen::Matrix<LhsScalar, -1, -1, 0>::IsVectorAtCompileTime, Eigen::InnerStride<1>, Eigen::OuterStride<> >::type = Eigen::OuterStride<>]?.
Spectra/LinAlg/TridiagEigen.h:46:20:   required from ?.pectra::TridiagEigen<Scalar>::TridiagEigen(Spectra::TridiagEigen<Scalar>::ConstGenericMatrix&) [with Scalar = float; Spectra::TridiagEigen<Scalar>::ConstGenericMatrix = const Eigen::Ref<const Eigen::Matrix<float, -1, -1> >; typename Eigen::internal::conditional<const Eigen::Matrix<LhsScalar, -1, -1, 0>::IsVectorAtCompileTime, Eigen::InnerStride<1>, Eigen::OuterStride<> >::type = Eigen::OuterStride<>]?
Spectra/SymEigsSolver.h:308:42:   required from ?.oid Spectra::SymEigsSolver<Scalar, SelectionRule, OpType>::retrieve_ritzpair() [with Scalar = float; int SelectionRule = 3; OpType = MatOpSymLower<float>]?
Spectra/SymEigsSolver.h:514:27:   required from ?.nt Spectra::SymEigsSolver<Scalar, SelectionRule, OpType>::compute(int, Scalar, int) [with Scalar = float; int SelectionRule = 3; OpType = MatOpSymLower<float>]?
ADMMLassoTall.h:199:33:   required from here
Spectra/LinAlg/TridiagEigen.h:93:107: error: no matching function for call to ?.ridiagonal_qr_step(float*&, float*&, int&, int&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1> >::Scalar*, int&)?
             Eigen::internal::tridiagonal_qr_step<Eigen::ColMajor>(maind, subd, start, end, evecs.data(), n);
                                                                                                           ^
Spectra/LinAlg/TridiagEigen.h:93:107: note: candidate is:
In file included from /usr/local/R/lib64/R/library/RcppEigen/include/Eigen/Eigenvalues:31:0,
                 from /usr/local/R/lib64/R/library/RcppEigen/include/Eigen/Dense:7,
                 from /usr/local/R/lib64/R/library/RcppEigen/include/RcppEigenForward.h:30,
                 from /usr/local/R/lib64/R/library/RcppEigen/include/RcppEigen.h:25,
                 from FADMMBase.h:4,
                 from ADMMLassoTall.h:4,
                 from ADMMEnet.h:4,
                 from Enet.cpp:3:
/usr/local/R/lib64/R/library/RcppEigen/include/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h:740:13: note: template<class RealScalar, class Scalar, class Index> void Eigen::internal::tridiagonal_qr_step(RealScalar*, RealScalar*, Index, Index, Scalar*, Index)
 static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n)
             ^
/usr/local/R/lib64/R/library/RcppEigen/include/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h:740:13: note:   template argument deduction/substitution failed:
make: *** [Enet.o] Error 1
ERROR: compilation failed for package ?.DMM?
* removing ?.usr/local/R/lib64/R/library/ADMM?
Error: Command failed (1)

Nonnegative constraints on fitted coefficients?

Would it be possible by any chance to also allow constraints on the fitted coefficients to be specified (e.g. nonnegativity), as one can do in glmnet using options lower.limits and upper.limits ?

Questions about the sparse lasso soluction from the implemented ADMM

I have a question about the estimated sparse LASSO coefficients. Based on Dr. Boyd MATLAB code in standford the estimated coefficients are small, but not zero.

The readme file shows that the implemented ADMM method has sparse-patten estimated coefficients, similar to that from the glmnet. I was curious how could you achieve that? I am not familiar with cpp so am not able to find the trick.

Thanks

Use sparse input matrix?

Would it also be possible by any chance to allow for a sparse input matrix (of class "sparseMatrix"), as one can do in glmnet?

Installation problem

Dear Sir,

I have problem with ADMM installation. I am using latest R-3.3.0 .
This is what I have finally using the install_github("yixuan/ADMM") command
ParLasso.o:ParLasso.cpp:(.text$_ZN17PADMMLasso_Master4initEdd[_ZN17PADMMLasso_Master4initEdd]+0x12b): undefined reference to ssyrk_' ParLasso.o:ParLasso.cpp:(.text$_ZN17PADMMLasso_Master4initEdd[__ZN17PADMMLasso_Master4initEdd]+0x38e): undefined reference to ssyrk'
collect2.exe: error: ld returned 1 exit status
no DLL was created
ERROR: compilation failed for package 'ADMM'

removing 'C:/Users/Alex Beylin/Documents/R/win-library/3.3/ADMM'
Error: Command failed (1)

I am very interesting in your package.
Please advise.

Regards

Alex Beylin

API not very friendly to programming

The current design of ADMM package funtions adopts the fluent style of programming, that is, the operations are done in a pipeline.

This can be friendly to interactive usage of model fitting functions, easier to read and write, rather than confused by so many parameters in one call.

However, the downside is obvious too. It is not very friendly to programming purposes. If I'm creating a general framework of modeling which supports a number of existing modeling packages (e.g., glmnet, xgboost, randomForest), it would be not easy to forward arguments to ADMM functions.

It would be appreciated if a simple argument-based interface is provided.

Thanks again, great work!

Fail to install on Mac

I have attempted to seek any possible solutions in StackOverflow，however, I am failed to find a good way to fix these bugs. As the error said:

clang: error: no such file or directory: '/usr/local/lib/libfontconfig.a'
clang: error: no such file or directory: '/usr/local/lib/libreadline.a'

some C++ problem occur and many said the problem is Xcode Path, but I have installed Xcode and command line.

xcode-select -p
/Applications/Xcode.app/Contents/Developer

I need your help.Here is some other useful information which may help you.

My sessionInfo()

> sessionInfo()
R version 3.2.2 (2015-08-14)
Platform: x86_64-apple-darwin14.5.0 (64-bit)
Running under: OS X 10.11.3 (El Capitan)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
 [1] httr_1.0.0     R6_2.1.1       magrittr_1.5   tools_3.2.2    curl_0.9.3    
 [6] memoise_0.2.1  stringi_0.5-5  stringr_1.0.0  digest_0.6.8   devtools_1.9.1

R.home()
[1] "/Library/Frameworks/R.framework/Resources"

uname

Darwin harryzhu.local 15.3.0 Darwin Kernel Version 15.3.0: Thu Dec 10 18:40:58 PST 2015; root:xnu-3248.30.4~1/RELEASE_X86_64 x86_64

Demo bug from http://statr.me/admm/

When I run this command, I find a bug.

fit3 = mod3$fit()
Error in .Call("admm_dantzig", .self$x, .self$y, .self$lambda, .self$nlambda,  : 
  "admm_dantzig" not available for .Call() for package "ADMM"

Fit GLMs with other distributions?

Nice package! Would it also be possible to use the same approach to fit other families/distributions, as one can do in glmnet, such as multivariate gaussian, multinomial, binomial or Poisson?

Adding interfaces for platforms other than R

Hi,
I would like to develop additional interfaces for your library.
I'm not experienced in R packages so can you give me some directions?

Where should I look first in order to develop interfaces.

Recommend Projects

React

A declarative, efficient, and flexible JavaScript library for building user interfaces.
Vue.js

🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
Typescript

TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
TensorFlow

An Open Source Machine Learning Framework for Everyone
Django

The Web framework for perfectionists with deadlines.
Laravel

A PHP framework for web artisans
D3

Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

javascript

JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
web

Some thing interesting about web. New door for the world.
server

A server is a program made to process requests and deliver data to clients.
Machine learning

Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Visualization

Some thing interesting about visualization, use data art
Game

Some thing interesting about game, make everyone happy.

Recommend Org

Facebook

We are working to build community through open source technology. NB: members must have two-factor auth.
Microsoft

Open source projects and samples from Microsoft.
Google

Google ❤️ Open Source for everyone.
Alibaba

Alibaba Open Source for everyone
D3

Data-Driven Documents codes.
Tencent

China tencent open source team.

yixuan / admm Goto Github PK

admm's Introduction

Introduction

Installation

Models

Lasso

Elastic Net

Least Absolute Deviation

Basis Pursuit

Performance

Lasso and Elastic Net

LAD

Basis Pursuit

admm's People

Contributors

Stargazers

Watchers

Forkers

admm's Issues

Recommend Projects

Recommend Topics

Recommend Org