Distributed Statistical Model Fitting with ADMM Algorithm

• What is ADMM?
• Distributed Computing with ADMM
• Implementation in Scala/Spark
• Computing on Real Data

• Full name: Alternating Direction Method of Multipliers
• An algorithm for constrained convex optimization problems
• Highlighted features:
  • Provides easy-to-apply solutions to many statistical models
  • Can be easily incorporated into the distributed computing framework
  • A good friend of statistical models

• ADMM is used to solve the following optimization problem:

minimizef(θ)+g(z)subject toAθ+Bz=c

• θ and z are vectors to be optimized
• A, B and c are known coefficient matrices/vectors

• f and g are convex functions

• Looks bizarre? Let's see two familiar examples

minf(θ)+g(z)s.t.Aθ+Bz=c

• Solve minβ12n‖y−Xβ‖22+λ‖β‖1
• f(β)=12n‖y−Xβ‖22
• g(z)=λ‖z‖1
• With constraint β−z=0

minf(θ)+g(z)s.t.Aθ+Bz=c

• Solve minβ‖y−Xβ‖1
• f(β)=0
• g(z)=‖z‖1
• With constraint Xβ+z=y

θk+1:=argminθ(f(θ)+ρ2‖Aθ+Bzk−c+uk/ρ‖22)zk+1:=argminz(g(z)+ρ2‖Aθk+1+Bz−c+uk/ρ‖22)uk+1:=uk+ρ(Aθk+1+Bzk+1−c)

• ρ can be any positive number
• Usually minθf(θ)+ρ2‖Aθ+w‖22 and minzg(z)+ρ2‖Bz+w‖22 are easier to compute than the original problem

• Lasso

βk+1:=(X′X+ρI)−1(X′y+ρzk−uk)zk+1:=Snλ/ρ(βk+1+uk/ρ)uk+1:=uk+ρ(βk+1−zk+1)

• Median Regression

βk+1:=(X′X)−1X′(y+zk−uk/ρ)zk+1:=S1/ρ(Xβk+1−y+uk/ρ)uk+1:=uk+ρ(Xβk+1−zk+1−y)

• where Sκ(a)=max(0,1−κ/|a|)⋅a

• Usually statistical models have the following form: minN∑i=1fi(θ)+g(θ)
• Here fi(θ) are partitions of the neg-log-likelihood function: If sample is cut into N partitions, the log-likelihood can be written as the summation of N parts
• g(θ) can be some penalty terms or prior information
• Models of this form can be easily parallelized

• The previous model can be rewritten as

minN∑i=1fi(θi)+g(z)s.t.θi−z=0, i=1,…,N

• The iteration steps are

θk+1i:=argminθi(fi(θi)+ρ2‖θi−zk+uki/ρ‖22)zk+1:=argminz(g(z)+Nρ2‖z−ˉθk+1−ˉuk/ρ‖22)uk+1i:=uki+ρ(θk+1i−zk+1)

• A diagram of P-ADMM

• https://github.com/yixuan/ADMM-Spark
• Currently used to solve logistic lasso regression
• Can be easily extended to general penalized likelihood model minθl(θ)+λ‖θ‖1

• class PADMML1
• For general L1-penalized likelihood model minθ∑Ni=1l(θ;xi,yi)+λ‖θ‖1
• Implementing the following iteration algorithm θk+1i:=argminθi(l(θi;xi,yi)+ρ2‖θi−zk+uki/ρ‖22)zk+1:=Sλ/(ρN)(ˉθk+1+ˉuk+1/ρ)uk+1i:=uki+ρ(θk+1i−zk+1)
• θ-step is undefined (abstract)

• class PLogisticLasso extends PADMML1
• Solves sparse logistic regression
• Only needs to implement the function for θ-step
• Calls the solver of logistic ridge regression

• class LogisticRidge
• class LogisticRidgeNative
• Solves the problem minθl(θ)+λ‖θ−v‖22
• LogisticRidge: pure Java code
• LogisticRidgeNative: using C++ code to speed up

import statr.stat598bd.PLogisticLasso

// Create RDD objects x and y representing
// the data matrix and response vector
// ...

val model = new PLogisticLasso(x, y, sc)
model.set_lambda(2.0)
model.set_opts(max_iter = 500,
               eps_abs = 1e-3,
               eps_rel = 1e-3,
               logs = true)
model.run()

val beta = model.coef

• Obtained from the Kaggle competition (https://www.kaggle.com/c/tradeshift-text-classification)
• Consists of thousands of documents, in which words are detected and combined to form text blocks

• 1.7 million rows (text blocks) / 145 variables

• Classify each text block to one or more labels (33 in total), such as date, invoice number or monetary amount

• We consider the 33rd label

• Numbers: Continuous/discrete numerical values
• Boolean: The values include YES (true) or NO (false)
• Encrypted text

• Preprocess the data by removing text and missing values, and encoding YES/NO values to 1/0
• Store the cleaned data to HDFS
• Read data from Spark and randomly split into training set (99%) and tuning set (1%)
• Standardize the training set
• Consider 10 λ values ranging from 0.001n to 0.1n, equally spaced in the log scale

• Fit a logistic lasso regression for each λ, and calculate the log-likelihood on tuning set
• Select the “optimal” λ
• Obtain the final regression coefficient based on the selected λ

• With 10 / 20 Executors
• Cleaning data: 41.3 s / 41.3 s
• Reading data in Spark: 23.6 s / 17.6 s
• Data standardization: 7.9 s / 7.7 s
• Setting up model: 9.4 s / 6.6 s
• Model fitting:

• Coding time: FOREVER

“Multithreading: Theory and Practice” (Image from http://wrathematics.github.io/RparallelGuide/)

• Boyd, S., Parikh, N., Chu, E., Peleato, B., & Eckstein, J. (2011). Distributed optimization and statistical learning via the alternating direction method of multipliers. Foundations and Trends® in Machine Learning, 3(1), 1-122.
• Minka, T. P. (2003). Algorithms for maximum-likelihood logistic regression.
• Spark documentation
• Breeze documentation
• Java Native Interface (JNI)