Neural Network in Haskell

Question

Following this website I wrote a neural network which uses the MNIST training data to recognize digits. The author writes that it should take a couple of minutes to train the network with 30 epochs of training. My network needs something like 5 minutes alone for 1 epoch of training.

How could I make it process faster?

Furthermore after 1 epoch of training it recognizes about 10 percent of the digits in the test file. The author's network recognizes 90 percent of the digits in the test file after one epoch of training.

Can I make it train better in the first epoch?

{-# LANGUAGE TypeFamilies #-}

module Blueprint where

import Codec.Compression.GZip (decompress)
import qualified Data.ByteString.Lazy as BS

import Prelude
import Numeric.LinearAlgebra
import Control.Monad
import Control.Arrow
import System.Random
import Data.List
import Data.Ord
import Data.VectorSpace
import Data.Array.IO


newtype Network = Network [( Matrix Double, Vector Double)] deriving (Eq,Show)

instance AdditiveGroup Network where
  (Network n1) ^+^ (Network n2) = Network $ zipWith (\(m,v) (n,w) -> (m+n,v+w)) n1 n2
  (Network n1) ^-^ (Network n2) = Network $ zipWith (\(m,v) (n,w) -> (m-n,v-w)) n1 n2
  zeroV = Network [(0,0) | x<-[0..]]
  negateV (Network n) = Network $ (\(m,v) -> (-m,-v)) <$> n

instance VectorSpace Network where
  type Scalar Network = Double
  lambda *^ (Network n) = Network $ (scale lambda Control.Arrow.*** scale lambda) <$> n


part :: Int -> [a] -> [[a]]
part n xs = if length xs >= n then take n xs : part n (drop n xs) else []

randomlist :: Int -> StdGen -> [Int]
randomlist n = take n . unfoldr (Just . random)

shuffle :: [a] -> IO [a]
shuffle xs = do
        ar <- newArray n xs
        forM [1..n] $ \i -> do
            j <- randomRIO (i,n)
            vi <- readArray ar i
            vj <- readArray ar j
            writeArray ar j vi
            return vj
  where
    n = length xs
    newArray :: Int -> [a] -> IO (IOArray Int a)
    newArray n =  newListArray (1,n)

sigmoid :: Double -> Double
sigmoid x = 1 / (1+exp (-x))

sigmoid' :: Double -> Double 
sigmoid' x = sigmoid x / (1 - sigmoid x)

getNetwork :: [Int] -> IO Network
getNetwork as@(_:bs) = 
  do 
    matrices <- zipWithM randn bs as
    seed <- newStdGen
    let rs = randomlist (length bs) seed
    let vectors = map (\(n,m) -> randomVector n Gaussian m) (zip rs bs) 
    return $ Network $ zip matrices vectors    

feed :: Network -> Vector Double -> Vector Double
feed (Network network) input = foldl (\ v (m, w) -> cmap sigmoid (m #> v + w)) input network 

train :: Network ->  [(Vector Double, Vector Double)] -> (Int, Int) -> Double -> IO Network
train network tdata (epochs,batchSize) eta = 
    if epochs == 0 then
        return network
    else 
        do
            shuffledData <- shuffle tdata
            let miniBatches = part batchSize shuffledData
            train (foldl (updateNetwork eta) network miniBatches) tdata (epochs-1, batchSize) eta

updateNetwork :: Double -> Network -> [(Vector Double, Vector Double)] -> Network
updateNetwork eta network miniBatch = 
    let 
      nabla = foldl (^+^) zeroV ((backpropagate network) <$> miniBatch)
      alpha = eta / fromIntegral (length miniBatch)
    in 
      network ^-^ alpha *^ nabla

backpropagate :: Network -> (Vector Double, Vector Double) -> Network
backpropagate (Network network) (input ,output) = 
    let 
      zs = tail $ scanl (\ z (m,v) -> m #> cmap sigmoid z + v) input network
      network1 = zip (tail $ fst <$> network) zs
      as = input : ((cmap sigmoid) <$> zs)
      deltaL = (last as - output) * cmap sigmoid' (last zs)
      deltas = scanr (\ (m,z) delta -> ((tr m) #> delta) * (cmap sigmoid' z)) deltaL network1  
    in 
      Network $ zip (zipWith (*) (asColumn <$> deltas) (tr . asColumn <$> as)) deltas

getInput s n = vector $ (/ 256) . fromIntegral . BS.index s . (n*28^2 + 16 +) <$> [0..28^2-1]
getLabel s n = fromIntegral $ BS.index s (n+8)
getOutput     s n = vector $ fromIntegral . fromEnum . (getLabel s n ==) <$> [0..9]

main = do  
    [inData, outData, inTest, outTest] <- mapM BS.readFile
      [ "train-images-idx3-ubyte"
      , "train-labels-idx1-ubyte"
      ,  "t10k-images-idx3-ubyte"
      ,  "t10k-labels-idx1-ubyte"
      ]
    network <- getNetwork [784, 30, 10]
    let tData = zip (getInput inData <$> [0..49999]) (getOutput outData <$> [0..49999])
    smart <- train network tData (30,10) 3
    let
      bestOf = fst . maximumBy (comparing snd) . zip [0..] . toList
      guesses = bestOf . (\n -> feed smart (getInput inTest n))  <$> [0..9999]
      answers = getLabel outTest <$> [0..9999]
    putStrLn $ show (sum $ fromEnum <$> zipWith (==) guesses answers) ++ " / 10000"

Answer 1

Let's start small and refactor:

type Network = ZipList [(Matrix Double, Vector Double)] and Data.NumInstances.Tuple will allow you to rewrite instance AdditiveGroup Network like so:

(^+^) = liftA2 (+)
(^-^) = liftA2 (-)
zeroV = pure 0
negateV = fmap negate

part is chunksOf from Data.List.Split, except it doesn't discard the remainder.

In getNetwork, letting the zipWithM envelop all removes the need for randomList, length bs and zipping:

getNetwork :: [Int] -> IO Network
getNetwork as@(_:bs) = fmap Network $ zipWithM foo bs as where
  foo b a = do 
    matrix <- randn b a
    n <- randomIO
    let vector = randomVector n Gaussian b
    return (matrix, vector)

shuffle has been done for you, for example in System.Random.Shuffle.

replicateM can take the recursion from train:

train :: Network ->  [(Vector Double, Vector Double)] -> (Int, Int) -> Double -> IO Network
train network tdata (epochs,batchSize) eta = do
  shuffledData <- replicateM epochs $ shuffle tdata
  let miniBatches = concatMap (part batchSize) shuffledData
  return $ foldl (updateNetwork eta) network miniBatches

I would put Network as the last argument when a function turns parameters into a Network transformer. Your mileage may vary.

Not really being proficient with neural networks and just looking at it vaguely, are you sure you want to use last and tail in backpropagate rather than head and tail or init and last?

Edit: To start making sense of backpropagate, I've rewritten it into State form to make the data flow more linear. coerce allows us to ignore Networks newtype constructor. This won't actually work, because my modifys are trying to change the type of the state along the way, but perhaps this'll be instructive?

backpropagate :: (Vector Double, Vector Double) -> Network -> Network
backpropagate (input, output) = coerce $ execState $ do
  modify $ liftA2 (zip `on` tail) (map fst) (scanl (\z (m,v) -> m #> cmap sigmoid z + v) input)
  as <- gets $ (input :) . map (cmap sigmoid) . map snd
  deltaL <- gets $ ((last as - output) *) . cmap sigmoid' . last . map snd
  modify $ scanr (\(m,z) delta -> (tr m #> delta) * cmap sigmoid' z) deltaL
  modify $ zipWith (\a delta -> (asColumn a * tr (asColumn delta), delta)) as