Support Vector Machine implementation translated from course material of "Machine Learning" on Coursera

ex6data1_X.csv

ex6data1_y.csv

1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1

SVM.hs

{-# LANGUAGE BangPatterns, FlexibleContexts, ScopedTypeVariables #-}
{-# OPTIONS_GHC -Wall #-}
-- Translation from svmTrain.m and svmPredict.m from https://www.coursera.org/course/ml
module SVM
  ( SVM (..)
  , TrainOptions (..)
  , Kernel (..)
  , train
  , predict
  ) where

import Control.Monad
import Data.Default.Class -- data-default-class
import Data.IORef
import qualified Data.Vector.Storable as VS
import qualified Data.Vector.Generic as V
import qualified Data.Vector.Generic.Mutable as VM
import Numeric.LinearAlgebra.HMatrix -- hmatrix
import qualified System.Random.MWC as Rand -- mwc-random

import System.Exit
import Text.CSV
import Text.Printf

data Kernel
  = LinearKernel
  | GaussianKernel !Double
  | Kernel (Vector Double -> Vector Double -> Double)

-- We have implemented optimized vectorized version of the Kernels here so
-- that the svm training will run faster.
kernelFunction :: Kernel -> Vector Double -> Vector Double -> Double
kernelFunction LinearKernel x1 x2 = x1 `dot` x2
kernelFunction (GaussianKernel sigma) x1 x2 = exp (- sumElements ((x1 - x2) ** scalar 2) / (2 * sigma**2)) 
kernelFunction (Kernel f) x1 x2 = f x1 x2

-- (m1 >< n) -> (m2 >< n) -> (m1 >< m2)
kernelMatrix :: Kernel -> Matrix Double -> Matrix Double -> Matrix Double
kernelMatrix LinearKernel xs1 xs2 = xs1 <> tr xs2
kernelMatrix (GaussianKernel sigma) xs1 xs2 = scalar (exp (-1 / (2*sigma**2))) ** k
  where
    m1 = rows xs1
    m2 = rows xs2
    xs1_2 = V.fromList $ map V.sum $ toRows $ xs1 ** scalar 2
    xs2_2 = V.fromList $ map V.sum $ toRows $ xs2 ** scalar 2
    k = accum (scale (-2) (xs1 <> tr xs2)) (+) [((i,j), xs1_2!i + xs2_2!j) | i<-[0..m1-1], j<-[0..m2-1]]
kernelMatrix (Kernel f) xs1 xs2
  = assoc (m1,m2) 0 [((i,j), f x1 x2) | (i,x1) <- zip [0..] (toRows xs1), (j,x2) <- zip [0..] (toRows xs2)]
  where
    m1 = rows xs1
    m2 = rows xs2
    
data SVM
  = SVM
  { svmKernel :: !Kernel
  , svmX      :: !(Matrix Double)
  , svmY      :: !(VS.Vector Double)
  , svmAlphas :: !(VS.Vector Double)
  , svmW      :: !(VS.Vector Double)
  , svmB      :: !Double
  }

data TrainOptions
  = TrainOptions
  { trainTol :: !Double
  , trainMaxPasses :: !Double
  , trainC :: !Double
  , trainKernel :: !Kernel
  , trainGen :: !(Maybe Rand.GenIO)
  }
  
instance Default TrainOptions where
  def =
    TrainOptions
    { trainTol = 1e-3
    , trainMaxPasses = 20
    , trainC = 100
    , trainKernel = LinearKernel
    , trainGen = Nothing
    }

{- |
@'train'@ Trains an SVM classifier using a simplified version of the SMO algorithm. 

@'train' opt xs ys@ trains an SVM classifier and returns trained model.
@xs@ is the matrix of training examples.  Each row is a training example,
and the jth column holds the jth feature. @ys@ is a vector containing @True@
for positive examples and @False@ for negative examples.

* @trainC opt@ is the standard SVM regularization  parameter.

* @trainTol opt@ is a tolerance value used for determining equality of
  floating point numbers.
* @trainMaxPasses opt@ controls the number of iterations over the dataset
  (without changes to alpha) before the algorithm quits.

Note: This is a simplified version of the SMO algorithm for training
SVMs. In practice, if you want to train an SVM classifier, we
recommend using an optimized package such as:

* LIBSVM   <http://www.csie.ntu.edu.tw/~cjlin/libsvm/>

* SVMLight <http://svmlight.joachims.org/>
-}
train :: TrainOptions -> Matrix Double -> Vector Bool -> IO SVM
train opt xs ys = do
  let tol = trainTol opt
      maxPasses = trainMaxPasses opt
      c = trainC opt
      -- m : number of samples
      -- n : number of features
      (m,n) = size xs
      ys2 = VS.map (\y -> if y then 1 else -1) ys
  gen <- do
    case trainGen opt of
      Nothing -> Rand.createSystemRandom
      Just gen -> return gen
      
  alphasM <- VM.replicate m 0
  bRef <- newIORef 0  

  --  Pre-compute the Kernel Matrix since our dataset is small
  --  (in practice, optimized SVM packages that handle large datasets
  --   gracefully will _not_ do this)
  let k = kernelMatrix (trainKernel opt) xs xs

  let -- Calculate Ei = f(x(i)) - y(i) using (2).
      computeE :: Int -> IO Double
      computeE i = do
        alphas <- V.unsafeFreeze alphasM
        b <- readIORef bRef
        -- f(x) = wx + b
        --   where w = sum_j alphas_j y_j x_j
        return $! (b + sumElements ((k ? [i]) #> (alphas * ys2))) - ys2!i
                     
  let loop !passes = when (passes < maxPasses) $ do
        numChangedAlphasRef <- newIORef (0::Int)
        forM_ [0..m-1] $ \i -> do
          alpha_i <- VM.read alphasM i
          b <- readIORef bRef

          -- Calculate Ei = f(x(i)) - y(i) using (2). 
          e_i <- computeE i

          when (ys2!i * e_i < - tol && alpha_i < c || ys2!i * e_i > tol && alpha_i > 0) $ do
            -- In practice, there are many heuristics one can use to select
            -- the i and j. In this simplified code, we select them randomly.
            let pickJ = do
                  j <- Rand.uniformR (0, m-1) gen
                  if j /= i then return j else pickJ -- Make sure i \neq j
            j <- pickJ

            alpha_j <- VM.read alphasM j
            -- Calculate Ej = f(x(j)) - y(j) using (2).
            e_j <- computeE j

            let -- Compute L and H by (10) or (11).
                (l,h) =
                  if ys2!i == ys2!j then
                    ( max 0 (alpha_j + alpha_i - c)
                    , min c (alpha_j + alpha_i) )
                  else
                    ( max 0 (alpha_j - alpha_i)
                    , min c (c + alpha_j - alpha_i) )
                -- Compute eta by (14).
                eta = 2 * atIndex k (i,j) - atIndex k (i,i) - atIndex k (j,j)
                -- Compute and clip new value for alpha j using (12) and (15).
                alpha_j_new = clip (l,h) $ alpha_j - (ys2!j * (e_i - e_j)) / eta
                -- Determine value for alpha i using (16). 
                alpha_i_new = alpha_i + ys2!i * ys2!j * (alpha_j  - alpha_j_new)
                -- Compute b1 and b2 using (17) and (18) respectively.
                b1 = b - e_i
                     - ys2!i * (alpha_i_new - alpha_i) * atIndex k (i,j) -- atIndex k (i,i) でなくて良い?
                     - ys2!j * (alpha_j_new - alpha_j) * atIndex k (i,j)
                b2 = b - e_j
                     - ys2!i * (alpha_i_new - alpha_i) * atIndex k (i,j)
                     - ys2!j * (alpha_j_new - alpha_j) * atIndex k (j,j) -- もしくは atIndex k (i,j) でなくて良い?
                -- Compute b by (19).
                b_new
                  | 0 < alpha_i_new && alpha_i_new < c = b1
                  | 0 < alpha_j_new && alpha_j_new < c = b2
                  | otherwise = (b1+b2) / 2

            when (l /= h && eta < 0 && abs (alpha_j_new - alpha_j) >= tol) $ do
              VM.write alphasM i alpha_i_new
              VM.write alphasM j alpha_j_new
              writeIORef bRef b_new
              modifyIORef' numChangedAlphasRef (+1)

              -- alphas <- V.unsafeFreeze alphasM
              -- let obj = V.sum alphas + (1/2) * sum [alphas!i * alphas!j * ys2!i * ys2!j * atIndex k (i,j) | i <- [0..m-1], j <- [0..m-1]]
              -- printf "(dual) obj = %f\n" obj

        numChangedAlphas <- readIORef numChangedAlphasRef
        if numChangedAlphas == 0 then
          loop (passes + 1)
        else
          loop 0

  loop 0

  alphas <- V.unsafeFreeze alphasM
  b <- readIORef bRef
  let idx = V.fromList [i | (i, alpha_i) <- zip [0..] (V.toList alphas), alpha_i > 0]
  return $
    SVM
    { svmKernel = trainKernel opt
    , svmX = xs ? (V.toList idx)
    , svmY = V.map (ys2 !) idx
    , svmAlphas = V.map (alphas !) idx
    , svmW = tr xs #> (alphas * ys2)
    , svmB = b
    }

clip :: Ord a => (a,a) -> a -> a
clip (l,h) x
  | x < l = l
  | h < x = h
  | otherwise = x

{-|
@'predict'@ returns a vector of predictions using a trained SVM model ('train'). 

@predict model xs@ returns a vector of predictions using a 
trained SVM model ('train'). xs is a mxn matrix where there each 
example is a row. model is a svm model returned from 'train'.
The result is a m vector of predictions of @Bool@ values.
-}
predict :: SVM -> Matrix Double -> Vector Bool
predict svm xs = V.map (>=0) p
  where
    (m,n) = size xs
    p = case svmKernel svm of
          LinearKernel -> xs #> svmW svm + scalar (svmB svm)
          GaussianKernel sigma -> kernelMatrix (svmKernel svm) xs (svmX svm) #> (svmAlphas svm * svmY svm) -- + scalar (svmB svm) は不要?
          Kernel f -> kernelMatrix (svmKernel svm) xs (svmX svm) #> (svmAlphas svm * svmY svm) + scalar (svmB svm)

loadData :: FilePath -> FilePath -> IO (Matrix Double, Vector Bool)
loadData xpath ypath = do
  x <- do
    ret <- parseCSVFromFile xpath
    case ret of
      Left err -> error (show err)
      Right rows -> return $ fromRows [V.fromList [read x :: Double | x <- row] | row <- rows, row /= [""]]
  y <- do
    ret <- parseCSVFromFile ypath
    case ret of
      Left err -> error (show err)
      Right rows -> return $ V.fromList [if y' == (1::Int) then True else False | [y] <- rows, y /= "", let y' = read y]
  return (x,y)
                      
ex6_2 :: IO SVM
ex6_2 = do
  (x,y) <- loadData "ex6data1_X.csv" "ex6data1_y.csv"
  svm <- train def x y
  print $ svmX svm
  print $ svmY svm
  print $ svmAlphas svm
  print $ svmW svm
  print $ svmB svm

  let y' = predict svm x
      accuracy :: Double
      accuracy = sum (zipWith (\y1 y2 -> if y1==y2 then 1 else 0) (V.toList y) (V.toList y')) / fromIntegral (V.length y)
  print $ accuracy
  
  return svm

ex6_3 :: Bool         
ex6_3 = abs (actual - expected) <= 1e-4
  where
    x1 = fromList [1,2,1]
    x2 = fromList [0,4,-1]
    sigma = 2
    actual = kernelFunction (GaussianKernel sigma) x1 x2
    expected = 0.324652

ex6_5 :: IO SVM
ex6_5 = do
  (x,y) <- loadData "ex6data2_X.csv" "ex6data2_y.csv"
  let sigma = 0.1
  svm <- train def{ trainC = 1, trainKernel = GaussianKernel sigma }  x y
  print $ svmX svm
  print $ svmY svm
  print $ svmAlphas svm
  print $ svmW svm
  print $ svmB svm

  let y' = predict svm x
      accuracy :: Double
      accuracy = sum (zipWith (\y1 y2 -> if y1==y2 then 1 else 0) (V.toList y) (V.toList y')) / fromIntegral (V.length y)
  print $ accuracy

  return svm