scripts/builtin/lasso.dml

#-------------------------------------------------------------
#
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.
#
#-------------------------------------------------------------

# Builtin function for the SpaRSA algorithm to perform lasso regression
# (SpaRSA .. Sparse Reconstruction by Separable Approximation)
#
# INPUTS:
# ---------------------------------------------------------------------------------------------
# NAME            TYPE     DEFAULT     MEANING
# ---------------------------------------------------------------------------------------------
# X               Double   ---         input feature matrix
# y               Double   ---         matrix Y columns of the design matrix
# tol             Double   1e-15       target convergence tolerance
# M               Integer  5           history length
# tau             Double   1           regularization component
# maxi            Integer  100         maximum number of iterations until convergence
# ---------------------------------------------------------------------------------------------
# OUTPUTS
# w               Double   ---        model matrix


m_lasso = function(Matrix[Double] X, Matrix[Double] y, Double tol = 10^(-15), 
  Integer M = 5, Double tau = 1, Integer maxi = 100, Boolean verbose = TRUE)
  return(Matrix[Double] w)
{
  n = nrow(X)
  m = ncol(X)

  #constants
  eta = 2
  sigma = 0.01
  alpha_min = 10^(-30)
  alpha_max = 10^(30)

  #init
  alpha = 1
  w = Rand(rows=m, cols=1, min=0, max=1, pdf="uniform")
  history = -1e30 * matrix(1, rows=M, cols=1)

  r = X %*% w - y
  g = t(X) %*% r
  obj = 0.5 * sum(r*r) + tau*sum(abs(w))

  if( verbose )
    print("Initial OBJ=" + obj)

  history[M,1] = obj

  inactive_set = matrix(1, rows=m, cols=1)
  iter = 0
  continue = TRUE
  while(iter < maxi & continue) {
    dw = matrix(0, rows=m, cols=1)
    dg = matrix(0, rows=m, cols=1)
    relChangeObj = -1.0

    inner_iter = 0
    inner_continue = TRUE
    inner_maxiter = 100
    while(inner_iter < inner_maxiter & inner_continue) {
      u = w - g/alpha
      lambda = tau/alpha

      wnew = sign(u) * (abs(u) - lambda) * ((abs(u) - lambda) > 0)
      dw = wnew - w
      dw2 = sum(dw*dw)

      r = X %*% wnew - y
      gnew = t(X) %*% r
      objnew = 0.5 * sum(r*r) + tau*sum(abs(wnew))
      obj_threshold = max(history) - 0.5*sigma*alpha*dw2

      if(objnew <= obj_threshold) {
        w = wnew
        dg = gnew - g
        g = gnew
        inner_continue = FALSE

        history[1:(M-1),] = history[2:M,]
        history[M,1] = objnew
        relChangeObj = abs(objnew - obj)/obj
        obj = objnew
      }
      else
        alpha = eta*alpha

      inner_iter = inner_iter + 1
    }

    if(inner_continue)
      print("Inner loop did not converge")

    alphanew = sum(dw*dg)/sum(dw*dw)
    alpha = max(alpha_min, min(alpha_max, alphanew))

    old_inactive_set = inactive_set
    inactive_set = w != 0
    diff = sum(abs(old_inactive_set - inactive_set))
    continue = (diff != 0 | relChangeObj >= tol)

    num_inactive = sum(w != 0)
    if( verbose )
      print("ITER=" + iter + " OBJ=" + obj + " relative change=" + relChangeObj + " num_inactive=" + num_inactive)
    iter = iter + 1
  }
}