lasso.clj

(ns ml-clj.lasso
  (:require [clojure.algo.generic.math-functions :as math]
            [clojure.core.matrix :as m]
            [clojure.core.matrix.operators :as M]))

(defn max-lambda
  [Y X]
  (second
   (apply max-key second
          (map-indexed vector
                       (map (fn [x]
                              (/ (m/mget (m/dot x (m/transpose Y)) 0)
                                 (m/row-count X)))
                            X)))))

(defn soft-thresholding
  [s lambda]
  (* (math/sgn s)
     (max 0 (- (math/abs s) lambda))))

(defn active-set
  [B]
  (filter (complement nil?)
          (map (fn [b i] (if (not (= 0 b)) i nil))
               B
               (range 0 (count B)))))

(defn convergenced?
  [Bnew Bold]
  (= (active-set Bnew)
     (active-set Bold)))

(defn coordinate-descent
  ([f Bold converged? j]
   (coordinate-descent f 1 Bold converged? j))
  ([f n Bold converged? j]
   (let [Bold-snapshot Bold
         Bnew (f Bold j)]
     (if (converged?
          Bnew
          Bold-snapshot)
       {:betas Bnew :iterations n}
       (recur f (+ n 1) Bnew converged? j)))))

(defn mae
  [Y X B]
  (/ (reduce +
             (flatten
              (m/emap math/abs
                      (M/- Y (map #(reduce + (m/emul % B)) X)))))
     (m/row-count X)))

(defn init-coefs
  [X]
  (vec (repeat (m/column-count X) 0)))

(defn lasso
  [Y X Bold lambda precision]
  (let [inner (fn [Blast j]
                (let [Xj (m/slice X j)
                      Rj (M/+
                          (M/- Y (map #(reduce + (m/emul % Blast)) X))
                          (M/* Xj (get Blast j)))
                      Bj* (/ (m/dot (m/transpose Xj)
                                    (flatten Rj)) ; Rj must be a vector
                             (m/row-count X))
                      Bjnew (soft-thresholding Bj* lambda)
                      Blast (assoc Blast j Bjnew)]
                  (if (zero? j)
                    Blast
                    (recur Blast (dec 1)))))
        coefs (coordinate-descent inner Bold convergenced?
                                  (dec (m/column-count X)))]
    (-> coefs
        (assoc :lambda lambda)
        (assoc :mae (mae Y X
                         (coefs :betas (init-coefs X)))))))

(defn lasso-cv
  [Y X & {:keys [precision step]
          :or {precision 0.0001 step 0.0001}}]
  (let [lambda (max-lambda Y X)]
    (apply min-key :mae
           (pmap #(lasso Y X (init-coefs X) % precision)
                 (range 0 lambda step)))))