diff --git a/Build.hs b/Build.hs --- a/Build.hs +++ b/Build.hs @@ -14,22 +14,28 @@ data Doc = Lab main :: IO () -main = getDirectoryFilesIO "samples" ["/*.lhs"] >>= \allSamps -> - getDirectoryFilesIO "src" ["//*.hs"] >>= \allSrc -> - getDirectoryFilesIO "app" ["//*.hs"] >>= \allApp -> +main = getDirectoryFilesIO "samples" ["/*.lhs", "/*.hs"] >>= \allSamps -> + getDirectoryFilesIO "src" ["//*.hs"] >>= \allSrc -> shakeArgs opts $ do want ["all"] "all" ~> - need ["pdf", "md", "haddocks", "gentags", "install"] + need ["pdf", "md", "haddocks", "gentags", "install", "exe"] "pdf" ~> - need (map (\f -> "renders" takeFileName f -<.> "pdf") allSamps) + need [ "renders" takeFileName f -<.> ".pdf" + | f <- allSamps, takeExtension f == ".lhs" + ] "md" ~> - need (map (\f -> "renders" takeFileName f -<.> "md") allSamps) + need [ "renders" takeFileName f -<.> ".md" + | f <- allSamps, takeExtension f == ".lhs" + ] + "exe" ~> + need (map (\f -> "samples-exe" dropExtension f) allSamps) + "haddocks" ~> do need (("src" ) <$> allSrc) cmd "jle-git-haddocks" @@ -37,7 +43,6 @@ "install" ~> do need . concat $ [ ("src" ) <$> allSrc , ("samples" ) <$> allSamps - , ("app" ) <$> allApp ] cmd "stack install" @@ -58,6 +63,23 @@ "-o" f src + "samples-exe/*" %> \f -> do + need ["install"] + [src] <- getDirectoryFiles "samples" $ (takeFileName f <.>) <$> ["hs","lhs"] + liftIO $ do + createDirectoryIfMissing True "samples-exe" + createDirectoryIfMissing True ".build" + removeFilesAfter "samples" ["/*.o"] + cmd "stack ghc --" ("samples" src) + "-o" f + "-hidir" ".build" + "-threaded" + "-rtsopts" + "-with-rtsopts=-N" + "-Wall" + "-O2" + "-package backprop" + ["tags","TAGS"] &%> \_ -> do need (("src" ) <$> allSrc) cmd "hasktags" "src/" @@ -65,4 +87,5 @@ "clean" ~> do unit $ cmd "stack clean" removeFilesAfter ".shake" ["//*"] + removeFilesAfter ".build" ["//*"] diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,10 +1,24 @@ Changelog ========= +Version 0.0.3.0 +--------------- + + + +* Removed samples as registered executables in the cabal file, to reduce + dependences to a bare minimum. For convenience, build script now also + compiles the samples into the local directory if *stack* is installed. + +* Added experimental (unsafe) combinators for working with GADTs with + existential types, `withGADT`, to *Numeric.Backprop* module. + +* Fixed broken links in Changelog. + Version 0.0.2.0 --------------- - + * Added optimized numeric `Op`s, and re-write `Num`/`Fractional`/`Floating` instances in terms of them. @@ -19,7 +33,7 @@ Version 0.0.1.0 --------------- - + * Initial pre-release, as a request for comments. API is in a usable form and everything is fully documented, but there are definitely some things diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -39,10 +39,21 @@ literate haskell file][mnist-lhs], or [rendered here as a PDF][mnist-pdf]! **Read this first!** -[mnist-lhs]: https://github.com/mstksg/backprop/blob/master/samples/MNIST.lhs -[mnist-pdf]: https://github.com/mstksg/backprop/blob/master/renders/MNIST.pdf +[mnist-lhs]: https://github.com/mstksg/backprop/blob/master/samples/backprop-mnist.lhs +[mnist-pdf]: https://github.com/mstksg/backprop/blob/master/renders/backprop-mnist.pdf +The [literate haskell file][mnist-lhs] is a standalone haskell file that you +can compile (preferably with `-O2`) on its own with stack or some other +dependency manager. It can also be compiled with the build script in the +project directory (if [stack][] is installed, and appropriate dependencies are +installed), using +[stack]: http://haskellstack.org/ + +~~~bash +$ ./Build.hs exe +~~~ + Brief example ------------- @@ -192,16 +203,22 @@ 5. Some open questions: - a. Is it possible to offer pattern matching on sum types/with different - constructors for implicit-graph backprop? It's possible for - explicit-graph versions already, with `choicesVar`, but not yet with - the implicit-graph interface. Could be similar to an "Applicative vs. - Monad" issue where you can only have pre-determined fixed computation - paths when using `Applicative`, but I'm not sure. Still, it would be - nice, because if this was possible, we could possibly do away with - explicit-graph mode completely. + a. Is it possible to offer pattern matching on sum types/with different + constructors for implicit-graph backprop? It's possible for + explicit-graph versions already, with `choicesVar`, but not yet with + the implicit-graph interface. Could be similar to an "Applicative vs. + Monad" issue where you can only have pre-determined fixed computation + paths when using `Applicative`, but I'm not sure. Still, it would be + nice, because if this was possible, we could possibly do away with + explicit-graph mode completely. - b. Though we already have sum type support with explicit-graph mode, we - can't support GADTs yet. It'd be nice to see if this is possible, - because a lot of dependently typed neural network stuff is made much - simpler with GADTs. + b. Though we already have safe sum type support with explicit-graph mode, + we can't support GADTs yet safely. It'd be nice to see if this is + possible, because a lot of dependently typed neural network stuff is + made much simpler with GADTs. + + As of v0.0.3.0, we have a way of dealing with GADTs in explicit-graph + mode (using `withGADT`) that is *unsafe*, and requires some ugly manual + plumbing by the user that could potentially be confusing. But it would + still be nice to have a way that is safe and doesn't require the manual + plumbing and isn't as easy to mess up. diff --git a/backprop.cabal b/backprop.cabal --- a/backprop.cabal +++ b/backprop.cabal @@ -1,5 +1,5 @@ name: backprop -version: 0.0.2.0 +version: 0.0.3.0 synopsis: Heterogeneous, type-safe automatic backpropagation in Haskell description: See . @@ -21,10 +21,13 @@ extra-source-files: README.md CHANGELOG.md Build.hs - renders/MNIST.md - renders/MNIST.pdf - renders/NeuralTest.md - renders/NeuralTest.pdf + renders/backprop-mnist.md + renders/backprop-mnist.pdf + renders/backprop-neural-test.md + renders/backprop-neural-test.pdf + samples/backprop-mnist.lhs + samples/backprop-monotest.hs + samples/backprop-neural-test.lhs cabal-version: >=1.10 library @@ -53,47 +56,6 @@ default-language: Haskell2010 ghc-options: -Wall -executable backprop-monotest - hs-source-dirs: samples - main-is: MonoTest.hs - ghc-options: -threaded -rtsopts -with-rtsopts=-N -Wall -O2 - build-depends: base - , backprop - default-language: Haskell2010 - -executable backprop-neuraltest - hs-source-dirs: samples - main-is: NeuralTest.lhs - ghc-options: -threaded -rtsopts -with-rtsopts=-N -Wall -O2 - build-depends: base - , ad - , backprop - , generics-sop - , hmatrix >= 0.18 - , mwc-random - , primitive - , singletons - , type-combinators - default-language: Haskell2010 - -executable backprop-mnist - hs-source-dirs: samples - main-is: MNIST.lhs - ghc-options: -threaded -rtsopts -with-rtsopts=-N -Wall -O2 - build-depends: base - , backprop - , bifunctors - , deepseq - , generics-sop - , hmatrix >= 0.18 - , mnist-idx - , mwc-random - , split - , time - , transformers - , vector - default-language: Haskell2010 - benchmark backprop-mnist-bench type: exitcode-stdio-1.0 hs-source-dirs: bench @@ -114,6 +76,17 @@ , type-combinators , vector default-language: Haskell2010 + +-- test-suite backprop-doctest +-- type: exitcode-stdio-1.0 +-- hs-source-dirs: doctest +-- main-is: doctest.hs +-- build-depends: base +-- , backprop +-- , doctest +-- , Glob +-- ghc-options: -threaded -rtsopts -with-rtsopts=-N +-- default-language: Haskell2010 -- test-suite backprop-test -- type: exitcode-stdio-1.0 diff --git a/renders/MNIST.md b/renders/MNIST.md deleted file mode 100644 --- a/renders/MNIST.md +++ /dev/null @@ -1,554 +0,0 @@ ---- -author: -- Justin Le -fontfamily: 'palatino,cmtt' -geometry: margin=1in -links-as-notes: true -title: Learning MNIST with Neural Networks with backprop library ---- - -The *backprop* library performs back-propagation over a *hetereogeneous* -system of relationships. It offers both an implicit (*[ad]*-like) and -explicit graph building usage style. Let’s use it to build neural -networks and learn mnist! - - [ad]: http://hackage.haskell.org/package/ad - -Repository source is [on github], and docs are [on hackage]. - - [on github]: https://github.com/mstksg/backprop - [on hackage]: http://hackage.haskell.org/package/backprop - -If you’re reading this as a literate haskell file, you should know that -a [rendered pdf version is available on github.]. If you are reading -this as a pdf file, you should know that a [literate haskell version -that you can run] is also available on github! - - [rendered pdf version is available on github.]: https://github.com/mstksg/backprop/blob/master/renders/MNIST.pdf - [literate haskell version that you can run]: https://github.com/mstksg/backprop/blob/master/samples/MNIST.lhs - -``` {.sourceCode .literate .haskell} -{-# LANGUAGE BangPatterns #-} -{-# LANGUAGE DataKinds #-} -{-# LANGUAGE DeriveGeneric #-} -{-# LANGUAGE GADTs #-} -{-# LANGUAGE LambdaCase #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE TupleSections #-} -{-# LANGUAGE TypeApplications #-} -{-# LANGUAGE ViewPatterns #-} -{-# OPTIONS_GHC -fno-warn-orphans #-} -{-# OPTIONS_GHC -fno-warn-incomplete-patterns #-} -{-# OPTIONS_GHC -fno-warn-unused-top-binds #-} - -import Control.DeepSeq -import Control.Exception -import Control.Monad -import Control.Monad.IO.Class -import Control.Monad.Trans.Maybe -import Control.Monad.Trans.State -import Data.Bitraversable -import Data.Foldable -import Data.IDX -import Data.List.Split -import Data.Maybe -import Data.Time.Clock -import Data.Traversable -import Data.Tuple -import GHC.Generics (Generic) -import GHC.TypeLits -import Numeric.Backprop -import Numeric.LinearAlgebra.Static hiding (dot) -import Text.Printf -import qualified Data.Vector as V -import qualified Data.Vector.Generic as VG -import qualified Data.Vector.Unboxed as VU -import qualified Generics.SOP as SOP -import qualified Numeric.LinearAlgebra as HM -import qualified System.Random.MWC as MWC -import qualified System.Random.MWC.Distributions as MWC -``` - -Types -===== - -For the most part, we’re going to be using the great *[hmatrix]* library -and its vector and matrix types. It offers a type `L m n` for -$m \times n$ matrices, and a type `R n` for an $n$ vector. - - [hmatrix]: http://hackage.haskell.org/package/hmatrix - -First things first: let’s define our neural networks as simple -containers of parameters (weight matrices and bias vectors). - -First, a type for layers: - -``` {.sourceCode .literate .haskell} -data Layer i o = - Layer { _lWeights :: !(L o i) - , _lBiases :: !(R o) - } - deriving (Show, Generic) - -instance SOP.Generic (Layer i o) -instance NFData (Layer i o) -``` - -And a type for a simple feed-forward network with two hidden layers: - -``` {.sourceCode .literate .haskell} -data Network i h1 h2 o = - Net { _nLayer1 :: !(Layer i h1) - , _nLayer2 :: !(Layer h1 h2) - , _nLayer3 :: !(Layer h2 o) - } - deriving (Show, Generic) - -instance SOP.Generic (Network i h1 h2 o) -instance NFData (Network i h1 h2 o) -``` - -These are pretty straightforward container types…pretty much exactly the -type you’d make to represent these networks! Note that, following true -Haskell form, we separate out logic from data. This should be all we -need. - -We derive an instance of `SOP.Generic` from the *[generics-sop]* -package, which *backprop* uses to propagate derivatives on values inside -product types. - - [generics-sop]: http://hackage.haskell.org/package/generics-sop - -Instances ---------- - -Things are much simplier if we had `Num` and `Fractional` instances for -everything, so let’s just go ahead and define that now, as well. Just a -little bit of boilerplate. - -``` {.sourceCode .literate .haskell} -instance (KnownNat i, KnownNat o) => Num (Layer i o) where - Layer w1 b1 + Layer w2 b2 = Layer (w1 + w2) (b1 + b2) - Layer w1 b1 - Layer w2 b2 = Layer (w1 - w2) (b1 - b2) - Layer w1 b1 * Layer w2 b2 = Layer (w1 * w2) (b1 * b2) - abs (Layer w b) = Layer (abs w) (abs b) - signum (Layer w b) = Layer (signum w) (signum b) - negate (Layer w b) = Layer (negate w) (negate b) - fromInteger x = Layer (fromInteger x) (fromInteger x) - -instance (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) => Num (Network i h1 h2 o) where - Net a b c + Net d e f = Net (a + d) (b + e) (c + f) - Net a b c - Net d e f = Net (a - d) (b - e) (c - f) - Net a b c * Net d e f = Net (a * d) (b * e) (c * f) - abs (Net a b c) = Net (abs a) (abs b) (abs c) - signum (Net a b c) = Net (signum a) (signum b) (signum c) - negate (Net a b c) = Net (negate a) (negate b) (negate c) - fromInteger x = Net (fromInteger x) (fromInteger x) (fromInteger x) - -instance (KnownNat i, KnownNat o) => Fractional (Layer i o) where - Layer w1 b1 / Layer w2 b2 = Layer (w1 / w2) (b1 / b2) - recip (Layer w b) = Layer (recip w) (recip b) - fromRational x = Layer (fromRational x) (fromRational x) - -instance (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) => Fractional (Network i h1 h2 o) where - Net a b c / Net d e f = Net (a / d) (b / e) (c / f) - recip (Net a b c) = Net (recip a) (recip b) (recip c) - fromRational x = Net (fromRational x) (fromRational x) (fromRational x) -``` - -`KnownNat` comes from *base*; it’s a typeclass that *hmatrix* uses to -refer to the numbers in its type and use it to go about its normal -hmatrixy business. - -Ops -=== - -Now, *backprop* does require *primitive* differentiable operations on -our relevant types to be defined. *backprop* uses these primitive `Op`s -to tie everything together. Ideally we’d import these from a library -that implements these for you, and the end-user never has to make `Op` -primitives. - -But in this case, I’m going to put the definitions here to show that -there isn’t any magic going on. If you’re curious, refer to -[documentation for `Op`] for more details on how `Op` is implemented and -how this works. - - [documentation for `Op`]: http://hackage.haskell.org/package/backprop/docs/Numeric-Backprop-Op.html - -First, matrix-vector multiplication primitive, giving an explicit -gradient function. - -``` {.sourceCode .literate .haskell} -matVec - :: (KnownNat m, KnownNat n) - => Op '[ L m n, R n ] (R m) -matVec = op2' $ \m v -> - ( m #> v, \(fromMaybe 1 -> g) -> - (g `outer` v, tr m #> g) - ) -``` - -Dot products would be nice too. - -``` {.sourceCode .literate .haskell} -dot :: KnownNat n - => Op '[ R n, R n ] Double -dot = op2' $ \x y -> - ( x <.> y, \case Nothing -> (y, x) - Just g -> (konst g * y, x * konst g) - ) -``` - -Also a “scaling” function, scales a vector by a given factor. - -``` {.sourceCode .literate .haskell} -scale - :: KnownNat n - => Op '[ Double, R n ] (R n) -scale = op2' $ \a x -> - ( konst a * x - , \case Nothing -> (HM.sumElements (extract x ), konst a ) - Just g -> (HM.sumElements (extract (x * g)), konst a * g) - ) -``` - -Finally, an operation to sum all of the items in the vector. - -``` {.sourceCode .literate .haskell} -vsum - :: KnownNat n - => Op '[ R n ] Double -vsum = op1' $ \x -> (HM.sumElements (extract x), maybe 1 konst) -``` - -And why not, here’s the [logistic function], which we’ll use as an -activation function for internal layers. We don’t need to define this as -an `Op` up-front right now, because the library can automatically -promote any numeric polymorphic function (an `a -> a` or `a -> a -> a`, -etc.) to an `Op` anyways. - - [logistic function]: https://en.wikipedia.org/wiki/Logistic_function - -``` {.sourceCode .literate .haskell} -logistic :: Floating a => a -> a -logistic x = 1 / (1 + exp (-x)) -``` - -Running our Network -=================== - -Now that we have our primitives in place, let’s actually write a -function to run our network! - -``` {.sourceCode .literate .haskell} -runLayer - :: (KnownNat i, KnownNat o) - => BPOp s '[ R i, Layer i o ] (R o) -runLayer = withInps $ \(x :< l :< Ø) -> do - w :< b :< Ø <- gTuple #<~ l - y <- matVec ~$ (w :< x :< Ø) - return $ y + b -``` - -A `BPOp s '[ R i, Layer i o ] (R o)` is a backpropagatable function that -produces an `R o` (a vector with `o` elements, from the *[hmatrix]* -library) given an input environment of an `R i` (the “input” of the -layer) and a layer. - - [hmatrix]: http://hackage.haskell.org/package/hmatrix - -We use `withInps` to bring the environment into scope as a bunch of -`BVar`s. `x` is a `BVar` containing the input vector, and `l` is a -`BVar` containing the layer. - -The first thing we do is split out the parts of the layer so we can work -with the internal matrices. We can use `#<~` to “split out” the -components of a `BVar`, splitting on `gTuple` (which uses `GHC.Generics` -to automatically figure out how to split up a product type). - -Then we apply `matVec` (our primitive `Op` that does matrix-vector -multiplication) to `w` and `x`, and then the result is that added to the -bias vector `b`. - -We can write the `runNetwork` function pretty much the same way. - -``` {.sourceCode .literate .haskell} -runNetwork - :: (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) - => BPOp s '[ R i, Network i h1 h2 o ] (R o) -runNetwork = withInps $ \(x :< n :< Ø) -> do - l1 :< l2 :< l3 :< Ø <- gTuple #<~ n - y <- runLayer -$ (x :< l1 :< Ø) - z <- runLayer -$ (logistic y :< l2 :< Ø) - r <- runLayer -$ (logistic z :< l3 :< Ø) - softmax -$ (r :< Ø) - where - softmax :: KnownNat n => BPOp s '[ R n ] (R n) - softmax = withInps $ \(x :< Ø) -> do - expX <- bindVar (exp x) - totX <- vsum ~$ (expX :< Ø) - scale ~$ (1/totX :< expX :< Ø) -``` - -After splitting out the layers in the input `Network`, we run each layer -successively using our previously defined `runLayer`, giving inputs -using `-$`. We can directly apply `logistic` to `BVar`s. At the end, we -run a [softmax function] because MNIST is a classification challenge. -The softmax is done by applying $e^x$ for every item in the input -vector, and dividing each element by the total. - - [softmax function]: https://en.wikipedia.org/wiki/Softmax_function - -The Magic ---------- - -What did we just define? Well, with a `BPOp s rs a`, we can *run* it and -get the output: - -``` {.sourceCode .literate .haskell} -runNetOnInp - :: (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) - => Network i h1 h2 o - -> R i - -> R o -runNetOnInp n x = evalBPOp runNetwork (x ::< n ::< Ø) -``` - -But, the magic part is that we can also get the gradient! - -``` {.sourceCode .literate .haskell} -gradNet - :: (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) - => Network i h1 h2 o - -> R i - -> Network i h1 h2 o -gradNet n x = case gradBPOp runNetwork (x ::< n ::< Ø) of - _gradX ::< gradN ::< Ø -> gradN -``` - -This gives the gradient of all of the parameters in the matrices and -vectors inside the `Network`, which we can use to “train”! - -Training -======== - -Now for the real work. To train a network, we can do gradient descent -based on the gradient of some type of *error function* with respect to -the network parameters. Let’s use the [cross entropy], which is popular -for classification problems. - - [cross entropy]: https://en.wikipedia.org/wiki/Cross_entropy - -``` {.sourceCode .literate .haskell} -crossEntropy - :: KnownNat n - => R n - -> BPOpI s '[ R n ] Double -crossEntropy targ (r :< Ø) = negate (dot .$ (log r :< t :< Ø)) - where - t = constVar targ -``` - -Given a target vector and a `BVar` referring to the result of the -network, we can directly apply: - -$$ -H(\mathbf{r}, \mathbf{t}) = - (log(\mathbf{r}) \cdot \mathbf{t}) -$$ - -Just for fun, I implemented `crossEntropy` in “implicit-graph” mode, so -you don’t see any binds or returns. - -Now, a function to make one gradient descent step based on an input -vector and a target, using `gradBPOp`: - -``` {.sourceCode .literate .haskell} -trainStep - :: forall i h1 h2 o. (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) - => Double - -> R i - -> R o - -> Network i h1 h2 o - -> Network i h1 h2 o -trainStep r !x !t !n = case gradBPOp o (x ::< n ::< Ø) of - _ ::< gN ::< Ø -> - n - (realToFrac r * gN) - where - o :: BPOp s '[ R i, Network i h1 h2 o ] Double - o = do - y <- runNetwork - implicitly (crossEntropy t) -$ (y :< Ø) -``` - -A convenient wrapper for training over all of the observations in a -list: - -``` {.sourceCode .literate .haskell} -trainList - :: (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) - => Double - -> [(R i, R o)] - -> Network i h1 h2 o - -> Network i h1 h2 o -trainList r = flip $ foldl' (\n (x,y) -> trainStep r x y n) -``` - -Pulling it all together -======================= - -`testNet` will be a quick way to test our net by computing the -percentage of correct guesses: (mostly using *hmatrix* stuff) - -``` {.sourceCode .literate .haskell} -testNet - :: forall i h1 h2 o. (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) - => [(R i, R o)] - -> Network i h1 h2 o - -> Double -testNet xs n = sum (map (uncurry test) xs) / fromIntegral (length xs) - where - test :: R i -> R o -> Double - test x (extract->t) - | HM.maxIndex t == HM.maxIndex (extract r) = 1 - | otherwise = 0 - where - r :: R o - r = evalBPOp runNetwork (x ::< n ::< Ø) -``` - -And now, a main loop! - -If you are following along at home, download the [mnist data set files] -and uncompress them into the folder `data`, and everything should work -fine. - - [mnist data set files]: http://yann.lecun.com/exdb/mnist/ - -``` {.sourceCode .literate .haskell} -main :: IO () -main = MWC.withSystemRandom $ \g -> do - Just train <- loadMNIST "data/train-images-idx3-ubyte" "data/train-labels-idx1-ubyte" - Just test <- loadMNIST "data/t10k-images-idx3-ubyte" "data/t10k-labels-idx1-ubyte" - putStrLn "Loaded data." - net0 <- MWC.uniformR @(Network 784 300 100 9) (-0.5, 0.5) g - flip evalStateT net0 . forM_ [1..] $ \e -> do - train' <- liftIO . fmap V.toList $ MWC.uniformShuffle (V.fromList train) g - liftIO $ printf "[Epoch %d]\n" (e :: Int) - - forM_ ([1..] `zip` chunksOf batch train') $ \(b, chnk) -> StateT $ \n0 -> do - printf "(Batch %d)\n" (b :: Int) - - t0 <- getCurrentTime - n' <- evaluate . force $ trainList rate chnk n0 - t1 <- getCurrentTime - printf "Trained on %d points in %s.\n" batch (show (t1 `diffUTCTime` t0)) - - let trainScore = testNet chnk n' - testScore = testNet test n' - printf "Training error: %.2f%%\n" ((1 - trainScore) * 100) - printf "Validation error: %.2f%%\n" ((1 - testScore ) * 100) - - return ((), n') - where - rate = 0.02 - batch = 5000 -``` - -Each iteration of the loop: - -1. Shuffles the training set -2. Splits it into chunks of `batch` size -3. Uses `trainList` to train over the batch -4. Computes the score based on `testNet` based on the training set and - the test set -5. Prints out the results - -And, that’s really it! - -Result ------- - -I haven’t put much into optimizing the library yet, but the network -(with hidden layer sizes 300 and 100) seems to take 25s on my computer -to finish a batch of 5000 training points. It’s slow (five minutes per -60000 point epooch), but it’s a first unoptimized run and a proof of -concept! It’s my goal to get this down to a point where the result has -the same performance characteristics as the actual backend (*hmatrix*), -and so overhead is 0. - -Main takeaways -============== - -Most of the actual heavy lifting/logic actually came from the *hmatrix* -library itself. We just created simple types to wrap up our bare -matrices. - -Basically, all that *backprop* did was give you an API to define *how to -run* a neural net — how to *run* a net based on a `Network` and `R i` -input you were given. The goal of the library is to let you write down -how to run things in as natural way as possible. - -And then, after things are run, we can just get the gradient and roll -from there! - -Because the heavy lifting is done by the data types themselves, we can -presumably plug in *any* type and any tensor/numerical backend, and reap -the benefits of those libraries’ optimizations and parallelizations. -*Any* type can be backpropagated! :D - -What now? ---------- - -Check out the docs for the [Numeric.Backprop] module for a more detailed -picture of what’s going on, or find more examples at the [github repo]! - - [Numeric.Backprop]: http://hackage.haskell.org/package/backprop/docs/Numeric-Backprop.html - [github repo]: https://github.com/mstksg/backprop - -Boring stuff -============ - -Here is a small wrapper function over the [mnist-idx] library loading -the contents of the idx files into *hmatrix* vectors: - - [mnist-idx]: http://hackage.haskell.org/package/mnist-idx - -``` {.sourceCode .literate .haskell} -loadMNIST - :: FilePath - -> FilePath - -> IO (Maybe [(R 784, R 9)]) -loadMNIST fpI fpL = runMaybeT $ do - i <- MaybeT $ decodeIDXFile fpI - l <- MaybeT $ decodeIDXLabelsFile fpL - d <- MaybeT . return $ labeledIntData l i - r <- MaybeT . return $ for d (bitraverse mkImage mkLabel . swap) - liftIO . evaluate $ force r - where - mkImage :: VU.Vector Int -> Maybe (R 784) - mkImage = create . VG.convert . VG.map (\i -> fromIntegral i / 255) - mkLabel :: Int -> Maybe (R 9) - mkLabel n = create $ HM.build 9 (\i -> if round i == n then 1 else 0) -``` - -And here are instances to generating random -vectors/matrices/layers/networks, used for the initialization step. - -``` {.sourceCode .literate .haskell} -instance KnownNat n => MWC.Variate (R n) where - uniform g = randomVector <$> MWC.uniform g <*> pure Uniform - uniformR (l, h) g = (\x -> x * (h - l) + l) <$> MWC.uniform g - -instance (KnownNat m, KnownNat n) => MWC.Variate (L m n) where - uniform g = uniformSample <$> MWC.uniform g <*> pure 0 <*> pure 1 - uniformR (l, h) g = (\x -> x * (h - l) + l) <$> MWC.uniform g - -instance (KnownNat i, KnownNat o) => MWC.Variate (Layer i o) where - uniform g = Layer <$> MWC.uniform g <*> MWC.uniform g - uniformR (l, h) g = (\x -> x * (h - l) + l) <$> MWC.uniform g - -instance (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) => MWC.Variate (Network i h1 h2 o) where - uniform g = Net <$> MWC.uniform g <*> MWC.uniform g <*> MWC.uniform g - uniformR (l, h) g = (\x -> x * (h - l) + l) <$> MWC.uniform g -``` diff --git a/renders/MNIST.pdf b/renders/MNIST.pdf deleted file mode 100644 Binary files a/renders/MNIST.pdf and /dev/null differ diff --git a/renders/NeuralTest.md b/renders/NeuralTest.md deleted file mode 100644 --- a/renders/NeuralTest.md +++ /dev/null @@ -1,447 +0,0 @@ ---- -author: -- Justin Le -fontfamily: 'palatino,cmtt' -geometry: margin=1in -links-as-notes: true -title: Neural networks with backprop library ---- - -The *backprop* library performs back-propagation over a *hetereogeneous* -system of relationships. It offers both an implicit ([ad]-like) and -explicit graph building usage style. Let’s use it to build neural -networks! - - [ad]: http://hackage.haskell.org/package/ad - -Repository source is [on github], and so are the [rendered unstable -docs]. - - [on github]: https://github.com/mstksg/backprop - [rendered unstable docs]: https://mstksg.github.io/backprop - -``` {.sourceCode .literate .haskell} -{-# LANGUAGE DeriveGeneric #-} -{-# LANGUAGE GADTs #-} -{-# LANGUAGE LambdaCase #-} -{-# LANGUAGE RankNTypes #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE StandaloneDeriving #-} -{-# LANGUAGE TypeApplications #-} -{-# LANGUAGE TypeInType #-} -{-# LANGUAGE TypeOperators #-} -{-# LANGUAGE ViewPatterns #-} -{-# OPTIONS_GHC -fno-warn-orphans #-} -{-# OPTIONS_GHC -fno-warn-unused-top-binds #-} - -import Data.Functor -import Data.Kind -import Data.Maybe -import Data.Singletons -import Data.Singletons.Prelude -import Data.Singletons.TypeLits -import Data.Type.Combinator -import Data.Type.Product -import GHC.Generics (Generic) -import Numeric.Backprop -import Numeric.Backprop.Iso -import Numeric.LinearAlgebra.Static hiding (dot) -import System.Random.MWC -import qualified Generics.SOP as SOP -``` - -Ops -=== - -First, we define values of `Op` for the operations we want to do. `Op`s -are bundles of functions packaged with their hetereogeneous gradients. -For simple numeric functions, *backprop* can derive `Op`s automatically. -But for matrix operations, we have to derive them ourselves. - -The types help us with matching up the dimensions, but we still need to -be careful that our gradients are calculated correctly. - -`L` and `R` are matrix and vector types from the great *hmatrix* -library. - -First, matrix-vector multiplication: - -``` {.sourceCode .literate .haskell} -matVec - :: (KnownNat m, KnownNat n) - => Op '[ L m n, R n ] (R m) -matVec = op2' $ \m v -> ( m #> v - , \(fromMaybe 1 -> g) -> - (g `outer` v, tr m #> g) - ) -``` - -Now, dot products: - -``` {.sourceCode .literate .haskell} -dot :: KnownNat n - => Op '[ R n, R n ] Double -dot = op2' $ \x y -> ( x <.> y - , \case Nothing -> (y, x) - Just g -> (konst g * y, x * konst g) - ) -``` - -Polymorphic functions can be easily turned into `Op`s with `op1`/`op2` -etc., but they can also be run directly on graph nodes. - -``` {.sourceCode .literate .haskell} -logistic :: Floating a => a -> a -logistic x = 1 / (1 + exp (-x)) -``` - -A Simple Complete Example -========================= - -At this point, we already have enough to train a simple -single-hidden-layer neural network: - -``` {.sourceCode .literate .haskell} -simpleOp - :: (KnownNat m, KnownNat n, KnownNat o) - => R m - -> BPOpI s '[ L n m, R n, L o n, R o ] (R o) -simpleOp inp = \(w1 :< b1 :< w2 :< b2 :< Ø) -> - let z = logistic $ liftB2 matVec w1 x + b1 - in logistic $ liftB2 matVec w2 z + b2 - where - x = constVar inp -``` - -Here, `simpleOp` is defined in implicit (non-monadic) style, given a -tuple of inputs and returning outputs. Now `simpleOp` can be “run” with -the input vectors and parameters (a `L n m`, `R n`, `L o n`, and `R o`) -and calculate the output of the neural net. - -``` {.sourceCode .literate .haskell} -runSimple - :: (KnownNat m, KnownNat n, KnownNat o) - => R m - -> Tuple '[ L n m, R n, L o n, R o ] - -> R o -runSimple inp = evalBPOp (implicitly $ simpleOp inp) -``` - -Alternatively, we can define `simpleOp` in explicit monadic style, were -we specify our graph nodes explicitly. The results should be the same. - -``` {.sourceCode .literate .haskell} -simpleOpExplicit - :: (KnownNat m, KnownNat n, KnownNat o) - => R m - -> BPOp s '[ L n m, R n, L o n, R o ] (R o) -simpleOpExplicit inp = withInps $ \(w1 :< b1 :< w2 :< b2 :< Ø) -> do - -- First layer - y1 <- matVec ~$ (w1 :< x1 :< Ø) - let x2 = logistic (y1 + b1) - -- Second layer - y2 <- matVec ~$ (w2 :< x2 :< Ø) - return $ logistic (y2 + b2) - where - x1 = constVar inp -``` - -Now, for the magic of *backprop*: the library can now take advantage of -the implicit (or explicit) graph and use it to do back-propagation, too! - -``` {.sourceCode .literate .haskell} -simpleGrad - :: forall m n o. (KnownNat m, KnownNat n, KnownNat o) - => R m - -> R o - -> Tuple '[ L n m, R n, L o n, R o ] - -> Tuple '[ L n m, R n, L o n, R o ] -simpleGrad inp targ params = gradBPOp opError params - where - opError :: BPOp s '[ L n m, R n, L o n, R o ] Double - opError = do - res <- implicitly $ simpleOp inp - -- we explicitly bind err to prevent recomputation - err <- bindVar $ res - t - dot ~$ (err :< err :< Ø) - where - t = constVar targ -``` - -The result is the gradient of the input tuple’s components, with respect -to the `Double` result of `opError` (the squared error). We can then use -this gradient to do gradient descent. - -With Parameter Containers -========================= - -This method doesn’t quite scale, because we might want to make networks -with multiple layers and parameterize networks by layers. Let’s make -some basic container data types to help us organize our types, including -a recursive `Network` type that lets us chain multiple layers. - -``` {.sourceCode .literate .haskell} -data Layer :: Nat -> Nat -> Type where - Layer :: { _lWeights :: L m n - , _lBiases :: R m - } - -> Layer n m - deriving (Show, Generic) - - -data Network :: Nat -> [Nat] -> Nat -> Type where - NØ :: !(Layer a b) -> Network a '[] b - (:&) :: !(Layer a b) -> Network b bs c -> Network a (b ': bs) c -``` - -A `Layer n m` is a layer taking an n-vector and returning an m-vector. A -`Network a '[b, c, d] e` would be a Network that takes in an a-vector -and outputs an e-vector, with hidden layers of sizes b, c, and d. - -Isomorphisms ------------- - -The *backprop* library lets you apply operations on “parts” of data -types (like on the weights and biases of a `Layer`) by using `Iso`’s -(isomorphisms), like the ones from the *lens* library. The library -doesn’t depend on lens, but it can use the `Iso`s from the library and -also custom-defined ones. - -First, we can auto-generate isomorphisms using the *generics-sop* -library: - -``` {.sourceCode .literate .haskell} -instance SOP.Generic (Layer n m) -``` - -And then can create isomorphisms by hand for the two `Network` -constructors: - -``` {.sourceCode .literate .haskell} -netExternal :: Iso' (Network a '[] b) (Tuple '[Layer a b]) -netExternal = iso (\case NØ x -> x ::< Ø) - (\case I x :< Ø -> NØ x ) - -netInternal :: Iso' (Network a (b ': bs) c) (Tuple '[Layer a b, Network b bs c]) -netInternal = iso (\case x :& xs -> x ::< xs ::< Ø) - (\case I x :< I xs :< Ø -> x :& xs ) -``` - -An `Iso' a (Tuple as)` means that an `a` can really just be seen as a -tuple of `as`. - -Running a network -================= - -Now, we can write the `BPOp` that reprenents running the network and -getting a result. We pass in a `Sing bs` (a singleton list of the hidden -layer sizes) so that we can “pattern match” on the list and handle the -different network constructors differently. - -``` {.sourceCode .literate .haskell} -netOp - :: forall s a bs c. (KnownNat a, KnownNat c) - => Sing bs - -> BPOp s '[ R a, Network a bs c ] (R c) -netOp sbs = go sbs - where - go :: forall d es. KnownNat d - => Sing es - -> BPOp s '[ R d, Network d es c ] (R c) - go = \case - SNil -> withInps $ \(x :< n :< Ø) -> do - -- peek into the NØ using netExternal iso - l :< Ø <- netExternal #<~ n - -- run the 'layerOp' BP, with x and l as inputs - bpOp layerOp ~$ (x :< l :< Ø) - SNat `SCons` ses -> withInps $ \(x :< n :< Ø) -> withSingI ses $ do - -- peek into the (:&) using the netInternal iso - l :< n' :< Ø <- netInternal #<~ n - -- run the 'layerOp' BP, with x and l as inputs - z <- bpOp layerOp ~$ (x :< l :< Ø) - -- run the 'go ses' BP, with z and n as inputs - bpOp (go ses) ~$ (z :< n' :< Ø) - layerOp - :: forall d e. (KnownNat d, KnownNat e) - => BPOp s '[ R d, Layer d e ] (R e) - layerOp = withInps $ \(x :< l :< Ø) -> do - -- peek into the layer using the gTuple iso, auto-generated with SOP.Generic - w :< b :< Ø <- gTuple #<~ l - y <- matVec ~$ (w :< x :< Ø) - return $ logistic (y + b) -``` - -There’s some singletons work going on here, but it’s fairly standard -singletons stuff. Most of the complexity here is from the static typing -in our neural network type, and *not* from *backprop*. - -From *backprop* specifically, the only elements are `#<~` lets you -“split” an input ref with the given iso, and `bpOp`, which converts a -`BPOp` into an `Op` that you can bind with `~$`. - -Note that this library doesn’t support truly pattern matching on GADTs, -and that we had to pass in `Sing bs` as a reference to the structure of -our networks. - -Gradient Descent ----------------- - -Now we can do simple gradient descent. Defining an error function: - -``` {.sourceCode .literate .haskell} -errOp - :: KnownNat m - => R m - -> BVar s rs (R m) - -> BPOp s rs Double -errOp targ r = do - err <- bindVar $ r - t - dot ~$ (err :< err :< Ø) - where - t = constVar targ -``` - -And now, we can use `backprop` to generate the gradient, and shift the -`Network`! Things are made a bit cleaner from the fact that -`Network a bs c` has a `Num` instance, so we can use `(-)` and `(*)` -etc. - -``` {.sourceCode .literate .haskell} -train - :: (KnownNat a, SingI bs, KnownNat c) - => Double - -> R a - -> R c - -> Network a bs c - -> Network a bs c -train r x t n = case backprop (errOp t =<< netOp sing) (x ::< n ::< Ø) of - (_, _ :< I g :< Ø) -> n - (realToFrac r * g) -``` - -(`(::<)` is cons and `Ø` is nil for tuples.) - -Main -==== - -`main`, which will train on sample data sets, is still in progress! -Right now it just generates a random network using the *mwc-random* -library and prints each internal layer. - -``` {.sourceCode .literate .haskell} -main :: IO () -main = withSystemRandom $ \g -> do - n <- uniform @(Network 4 '[3,2] 1) g - void $ traverseNetwork sing (\l -> l <$ print l) n -``` - -Appendix: Boilerplate -===================== - -And now for some typeclass instances and boilerplates unrelated to the -*backprop* library that makes our custom types easier to use. - -``` {.sourceCode .literate .haskell} -instance KnownNat n => Variate (R n) where - uniform g = randomVector <$> uniform g <*> pure Uniform - uniformR (l, h) g = (\x -> x * (h - l) + l) <$> uniform g - -instance (KnownNat m, KnownNat n) => Variate (L m n) where - uniform g = uniformSample <$> uniform g <*> pure 0 <*> pure 1 - uniformR (l, h) g = (\x -> x * (h - l) + l) <$> uniform g - -instance (KnownNat n, KnownNat m) => Variate (Layer n m) where - uniform g = subtract 1 . (* 2) <$> (Layer <$> uniform g <*> uniform g) - uniformR (l, h) g = (\x -> x * (h - l) + l) <$> uniform g - -instance (KnownNat m, KnownNat n) => Num (Layer n m) where - Layer w1 b1 + Layer w2 b2 = Layer (w1 + w2) (b1 + b2) - Layer w1 b1 - Layer w2 b2 = Layer (w1 - w2) (b1 - b2) - Layer w1 b1 * Layer w2 b2 = Layer (w1 * w2) (b1 * b2) - abs (Layer w b) = Layer (abs w) (abs b) - signum (Layer w b) = Layer (signum w) (signum b) - negate (Layer w b) = Layer (negate w) (negate b) - fromInteger x = Layer (fromInteger x) (fromInteger x) - -instance (KnownNat m, KnownNat n) => Fractional (Layer n m) where - Layer w1 b1 / Layer w2 b2 = Layer (w1 / w2) (b1 / b2) - recip (Layer w b) = Layer (recip w) (recip b) - fromRational x = Layer (fromRational x) (fromRational x) - -instance (KnownNat a, SingI bs, KnownNat c) => Variate (Network a bs c) where - uniform g = genNet sing (uniform g) - uniformR (l, h) g = (\x -> x * (h - l) + l) <$> uniform g - -genNet - :: forall f a bs c. (Applicative f, KnownNat a, KnownNat c) - => Sing bs - -> (forall d e. (KnownNat d, KnownNat e) => f (Layer d e)) - -> f (Network a bs c) -genNet sbs f = go sbs - where - go :: forall d es. KnownNat d => Sing es -> f (Network d es c) - go = \case - SNil -> NØ <$> f - SNat `SCons` ses -> (:&) <$> f <*> go ses - -mapNetwork0 - :: forall a bs c. (KnownNat a, KnownNat c) - => Sing bs - -> (forall d e. (KnownNat d, KnownNat e) => Layer d e) - -> Network a bs c -mapNetwork0 sbs f = getI $ genNet sbs (I f) - -traverseNetwork - :: forall a bs c f. (KnownNat a, KnownNat c, Applicative f) - => Sing bs - -> (forall d e. (KnownNat d, KnownNat e) => Layer d e -> f (Layer d e)) - -> Network a bs c - -> f (Network a bs c) -traverseNetwork sbs f = go sbs - where - go :: forall d es. KnownNat d => Sing es -> Network d es c -> f (Network d es c) - go = \case - SNil -> \case - NØ x -> NØ <$> f x - SNat `SCons` ses -> \case - x :& xs -> (:&) <$> f x <*> go ses xs - -mapNetwork1 - :: forall a bs c. (KnownNat a, KnownNat c) - => Sing bs - -> (forall d e. (KnownNat d, KnownNat e) => Layer d e -> Layer d e) - -> Network a bs c - -> Network a bs c -mapNetwork1 sbs f = getI . traverseNetwork sbs (I . f) - -mapNetwork2 - :: forall a bs c. (KnownNat a, KnownNat c) - => Sing bs - -> (forall d e. (KnownNat d, KnownNat e) => Layer d e -> Layer d e -> Layer d e) - -> Network a bs c - -> Network a bs c - -> Network a bs c -mapNetwork2 sbs f = go sbs - where - go :: forall d es. KnownNat d => Sing es -> Network d es c -> Network d es c -> Network d es c - go = \case - SNil -> \case - NØ x -> \case - NØ y -> NØ (f x y) - SNat `SCons` ses -> \case - x :& xs -> \case - y :& ys -> f x y :& go ses xs ys - -instance (KnownNat a, SingI bs, KnownNat c) => Num (Network a bs c) where - (+) = mapNetwork2 sing (+) - (-) = mapNetwork2 sing (-) - (*) = mapNetwork2 sing (*) - negate = mapNetwork1 sing negate - abs = mapNetwork1 sing abs - signum = mapNetwork1 sing signum - fromInteger x = mapNetwork0 sing (fromInteger x) - -instance (KnownNat a, SingI bs, KnownNat c) => Fractional (Network a bs c) where - (/) = mapNetwork2 sing (/) - recip = mapNetwork1 sing recip - fromRational x = mapNetwork0 sing (fromRational x) -``` diff --git a/renders/NeuralTest.pdf b/renders/NeuralTest.pdf deleted file mode 100644 Binary files a/renders/NeuralTest.pdf and /dev/null differ diff --git a/renders/backprop-mnist.md b/renders/backprop-mnist.md new file mode 100644 --- /dev/null +++ b/renders/backprop-mnist.md @@ -0,0 +1,554 @@ +--- +author: +- Justin Le +fontfamily: 'palatino,cmtt' +geometry: margin=1in +links-as-notes: true +title: Learning MNIST with Neural Networks with backprop library +--- + +The *backprop* library performs back-propagation over a *hetereogeneous* +system of relationships. It offers both an implicit (*[ad]*-like) and +explicit graph building usage style. Let’s use it to build neural +networks and learn mnist! + + [ad]: http://hackage.haskell.org/package/ad + +Repository source is [on github], and docs are [on hackage]. + + [on github]: https://github.com/mstksg/backprop + [on hackage]: http://hackage.haskell.org/package/backprop + +If you’re reading this as a literate haskell file, you should know that +a [rendered pdf version is available on github.]. If you are reading +this as a pdf file, you should know that a [literate haskell version +that you can run] is also available on github! + + [rendered pdf version is available on github.]: https://github.com/mstksg/backprop/blob/master/renders/backprop-mnist.pdf + [literate haskell version that you can run]: https://github.com/mstksg/backprop/blob/master/samples/backprop-mnist.lhs + +``` {.sourceCode .literate .haskell} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE DeriveGeneric #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE LambdaCase #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TupleSections #-} +{-# LANGUAGE TypeApplications #-} +{-# LANGUAGE ViewPatterns #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} +{-# OPTIONS_GHC -fno-warn-incomplete-patterns #-} +{-# OPTIONS_GHC -fno-warn-unused-top-binds #-} + +import Control.DeepSeq +import Control.Exception +import Control.Monad +import Control.Monad.IO.Class +import Control.Monad.Trans.Maybe +import Control.Monad.Trans.State +import Data.Bitraversable +import Data.Foldable +import Data.IDX +import Data.List.Split +import Data.Maybe +import Data.Time.Clock +import Data.Traversable +import Data.Tuple +import GHC.Generics (Generic) +import GHC.TypeLits +import Numeric.Backprop +import Numeric.LinearAlgebra.Static hiding (dot) +import Text.Printf +import qualified Data.Vector as V +import qualified Data.Vector.Generic as VG +import qualified Data.Vector.Unboxed as VU +import qualified Generics.SOP as SOP +import qualified Numeric.LinearAlgebra as HM +import qualified System.Random.MWC as MWC +import qualified System.Random.MWC.Distributions as MWC +``` + +Types +===== + +For the most part, we’re going to be using the great *[hmatrix]* library +and its vector and matrix types. It offers a type `L m n` for +$m \times n$ matrices, and a type `R n` for an $n$ vector. + + [hmatrix]: http://hackage.haskell.org/package/hmatrix + +First things first: let’s define our neural networks as simple +containers of parameters (weight matrices and bias vectors). + +First, a type for layers: + +``` {.sourceCode .literate .haskell} +data Layer i o = + Layer { _lWeights :: !(L o i) + , _lBiases :: !(R o) + } + deriving (Show, Generic) + +instance SOP.Generic (Layer i o) +instance NFData (Layer i o) +``` + +And a type for a simple feed-forward network with two hidden layers: + +``` {.sourceCode .literate .haskell} +data Network i h1 h2 o = + Net { _nLayer1 :: !(Layer i h1) + , _nLayer2 :: !(Layer h1 h2) + , _nLayer3 :: !(Layer h2 o) + } + deriving (Show, Generic) + +instance SOP.Generic (Network i h1 h2 o) +instance NFData (Network i h1 h2 o) +``` + +These are pretty straightforward container types…pretty much exactly the +type you’d make to represent these networks! Note that, following true +Haskell form, we separate out logic from data. This should be all we +need. + +We derive an instance of `SOP.Generic` from the *[generics-sop]* +package, which *backprop* uses to propagate derivatives on values inside +product types. + + [generics-sop]: http://hackage.haskell.org/package/generics-sop + +Instances +--------- + +Things are much simplier if we had `Num` and `Fractional` instances for +everything, so let’s just go ahead and define that now, as well. Just a +little bit of boilerplate. + +``` {.sourceCode .literate .haskell} +instance (KnownNat i, KnownNat o) => Num (Layer i o) where + Layer w1 b1 + Layer w2 b2 = Layer (w1 + w2) (b1 + b2) + Layer w1 b1 - Layer w2 b2 = Layer (w1 - w2) (b1 - b2) + Layer w1 b1 * Layer w2 b2 = Layer (w1 * w2) (b1 * b2) + abs (Layer w b) = Layer (abs w) (abs b) + signum (Layer w b) = Layer (signum w) (signum b) + negate (Layer w b) = Layer (negate w) (negate b) + fromInteger x = Layer (fromInteger x) (fromInteger x) + +instance (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) => Num (Network i h1 h2 o) where + Net a b c + Net d e f = Net (a + d) (b + e) (c + f) + Net a b c - Net d e f = Net (a - d) (b - e) (c - f) + Net a b c * Net d e f = Net (a * d) (b * e) (c * f) + abs (Net a b c) = Net (abs a) (abs b) (abs c) + signum (Net a b c) = Net (signum a) (signum b) (signum c) + negate (Net a b c) = Net (negate a) (negate b) (negate c) + fromInteger x = Net (fromInteger x) (fromInteger x) (fromInteger x) + +instance (KnownNat i, KnownNat o) => Fractional (Layer i o) where + Layer w1 b1 / Layer w2 b2 = Layer (w1 / w2) (b1 / b2) + recip (Layer w b) = Layer (recip w) (recip b) + fromRational x = Layer (fromRational x) (fromRational x) + +instance (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) => Fractional (Network i h1 h2 o) where + Net a b c / Net d e f = Net (a / d) (b / e) (c / f) + recip (Net a b c) = Net (recip a) (recip b) (recip c) + fromRational x = Net (fromRational x) (fromRational x) (fromRational x) +``` + +`KnownNat` comes from *base*; it’s a typeclass that *hmatrix* uses to +refer to the numbers in its type and use it to go about its normal +hmatrixy business. + +Ops +=== + +Now, *backprop* does require *primitive* differentiable operations on +our relevant types to be defined. *backprop* uses these primitive `Op`s +to tie everything together. Ideally we’d import these from a library +that implements these for you, and the end-user never has to make `Op` +primitives. + +But in this case, I’m going to put the definitions here to show that +there isn’t any magic going on. If you’re curious, refer to +[documentation for `Op`] for more details on how `Op` is implemented and +how this works. + + [documentation for `Op`]: http://hackage.haskell.org/package/backprop/docs/Numeric-Backprop-Op.html + +First, matrix-vector multiplication primitive, giving an explicit +gradient function. + +``` {.sourceCode .literate .haskell} +matVec + :: (KnownNat m, KnownNat n) + => Op '[ L m n, R n ] (R m) +matVec = op2' $ \m v -> + ( m #> v, \(fromMaybe 1 -> g) -> + (g `outer` v, tr m #> g) + ) +``` + +Dot products would be nice too. + +``` {.sourceCode .literate .haskell} +dot :: KnownNat n + => Op '[ R n, R n ] Double +dot = op2' $ \x y -> + ( x <.> y, \case Nothing -> (y, x) + Just g -> (konst g * y, x * konst g) + ) +``` + +Also a “scaling” function, scales a vector by a given factor. + +``` {.sourceCode .literate .haskell} +scale + :: KnownNat n + => Op '[ Double, R n ] (R n) +scale = op2' $ \a x -> + ( konst a * x + , \case Nothing -> (HM.sumElements (extract x ), konst a ) + Just g -> (HM.sumElements (extract (x * g)), konst a * g) + ) +``` + +Finally, an operation to sum all of the items in the vector. + +``` {.sourceCode .literate .haskell} +vsum + :: KnownNat n + => Op '[ R n ] Double +vsum = op1' $ \x -> (HM.sumElements (extract x), maybe 1 konst) +``` + +And why not, here’s the [logistic function], which we’ll use as an +activation function for internal layers. We don’t need to define this as +an `Op` up-front right now, because the library can automatically +promote any numeric polymorphic function (an `a -> a` or `a -> a -> a`, +etc.) to an `Op` anyways. + + [logistic function]: https://en.wikipedia.org/wiki/Logistic_function + +``` {.sourceCode .literate .haskell} +logistic :: Floating a => a -> a +logistic x = 1 / (1 + exp (-x)) +``` + +Running our Network +=================== + +Now that we have our primitives in place, let’s actually write a +function to run our network! + +``` {.sourceCode .literate .haskell} +runLayer + :: (KnownNat i, KnownNat o) + => BPOp s '[ R i, Layer i o ] (R o) +runLayer = withInps $ \(x :< l :< Ø) -> do + w :< b :< Ø <- gTuple #<~ l + y <- matVec ~$ (w :< x :< Ø) + return $ y + b +``` + +A `BPOp s '[ R i, Layer i o ] (R o)` is a backpropagatable function that +produces an `R o` (a vector with `o` elements, from the *[hmatrix]* +library) given an input environment of an `R i` (the “input” of the +layer) and a layer. + + [hmatrix]: http://hackage.haskell.org/package/hmatrix + +We use `withInps` to bring the environment into scope as a bunch of +`BVar`s. `x` is a `BVar` containing the input vector, and `l` is a +`BVar` containing the layer. + +The first thing we do is split out the parts of the layer so we can work +with the internal matrices. We can use `#<~` to “split out” the +components of a `BVar`, splitting on `gTuple` (which uses `GHC.Generics` +to automatically figure out how to split up a product type). + +Then we apply `matVec` (our primitive `Op` that does matrix-vector +multiplication) to `w` and `x`, and then the result is that added to the +bias vector `b`. + +We can write the `runNetwork` function pretty much the same way. + +``` {.sourceCode .literate .haskell} +runNetwork + :: (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) + => BPOp s '[ R i, Network i h1 h2 o ] (R o) +runNetwork = withInps $ \(x :< n :< Ø) -> do + l1 :< l2 :< l3 :< Ø <- gTuple #<~ n + y <- runLayer -$ (x :< l1 :< Ø) + z <- runLayer -$ (logistic y :< l2 :< Ø) + r <- runLayer -$ (logistic z :< l3 :< Ø) + softmax -$ (r :< Ø) + where + softmax :: KnownNat n => BPOp s '[ R n ] (R n) + softmax = withInps $ \(x :< Ø) -> do + expX <- bindVar (exp x) + totX <- vsum ~$ (expX :< Ø) + scale ~$ (1/totX :< expX :< Ø) +``` + +After splitting out the layers in the input `Network`, we run each layer +successively using our previously defined `runLayer`, giving inputs +using `-$`. We can directly apply `logistic` to `BVar`s. At the end, we +run a [softmax function] because MNIST is a classification challenge. +The softmax is done by applying $e^x$ for every item in the input +vector, and dividing each element by the total. + + [softmax function]: https://en.wikipedia.org/wiki/Softmax_function + +The Magic +--------- + +What did we just define? Well, with a `BPOp s rs a`, we can *run* it and +get the output: + +``` {.sourceCode .literate .haskell} +runNetOnInp + :: (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) + => Network i h1 h2 o + -> R i + -> R o +runNetOnInp n x = evalBPOp runNetwork (x ::< n ::< Ø) +``` + +But, the magic part is that we can also get the gradient! + +``` {.sourceCode .literate .haskell} +gradNet + :: (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) + => Network i h1 h2 o + -> R i + -> Network i h1 h2 o +gradNet n x = case gradBPOp runNetwork (x ::< n ::< Ø) of + _gradX ::< gradN ::< Ø -> gradN +``` + +This gives the gradient of all of the parameters in the matrices and +vectors inside the `Network`, which we can use to “train”! + +Training +======== + +Now for the real work. To train a network, we can do gradient descent +based on the gradient of some type of *error function* with respect to +the network parameters. Let’s use the [cross entropy], which is popular +for classification problems. + + [cross entropy]: https://en.wikipedia.org/wiki/Cross_entropy + +``` {.sourceCode .literate .haskell} +crossEntropy + :: KnownNat n + => R n + -> BPOpI s '[ R n ] Double +crossEntropy targ (r :< Ø) = negate (dot .$ (log r :< t :< Ø)) + where + t = constVar targ +``` + +Given a target vector and a `BVar` referring to the result of the +network, we can directly apply: + +$$ +H(\mathbf{r}, \mathbf{t}) = - (log(\mathbf{r}) \cdot \mathbf{t}) +$$ + +Just for fun, I implemented `crossEntropy` in “implicit-graph” mode, so +you don’t see any binds or returns. + +Now, a function to make one gradient descent step based on an input +vector and a target, using `gradBPOp`: + +``` {.sourceCode .literate .haskell} +trainStep + :: forall i h1 h2 o. (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) + => Double + -> R i + -> R o + -> Network i h1 h2 o + -> Network i h1 h2 o +trainStep r !x !t !n = case gradBPOp o (x ::< n ::< Ø) of + _ ::< gN ::< Ø -> + n - (realToFrac r * gN) + where + o :: BPOp s '[ R i, Network i h1 h2 o ] Double + o = do + y <- runNetwork + implicitly (crossEntropy t) -$ (y :< Ø) +``` + +A convenient wrapper for training over all of the observations in a +list: + +``` {.sourceCode .literate .haskell} +trainList + :: (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) + => Double + -> [(R i, R o)] + -> Network i h1 h2 o + -> Network i h1 h2 o +trainList r = flip $ foldl' (\n (x,y) -> trainStep r x y n) +``` + +Pulling it all together +======================= + +`testNet` will be a quick way to test our net by computing the +percentage of correct guesses: (mostly using *hmatrix* stuff) + +``` {.sourceCode .literate .haskell} +testNet + :: forall i h1 h2 o. (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) + => [(R i, R o)] + -> Network i h1 h2 o + -> Double +testNet xs n = sum (map (uncurry test) xs) / fromIntegral (length xs) + where + test :: R i -> R o -> Double + test x (extract->t) + | HM.maxIndex t == HM.maxIndex (extract r) = 1 + | otherwise = 0 + where + r :: R o + r = evalBPOp runNetwork (x ::< n ::< Ø) +``` + +And now, a main loop! + +If you are following along at home, download the [mnist data set files] +and uncompress them into the folder `data`, and everything should work +fine. + + [mnist data set files]: http://yann.lecun.com/exdb/mnist/ + +``` {.sourceCode .literate .haskell} +main :: IO () +main = MWC.withSystemRandom $ \g -> do + Just train <- loadMNIST "data/train-images-idx3-ubyte" "data/train-labels-idx1-ubyte" + Just test <- loadMNIST "data/t10k-images-idx3-ubyte" "data/t10k-labels-idx1-ubyte" + putStrLn "Loaded data." + net0 <- MWC.uniformR @(Network 784 300 100 9) (-0.5, 0.5) g + flip evalStateT net0 . forM_ [1..] $ \e -> do + train' <- liftIO . fmap V.toList $ MWC.uniformShuffle (V.fromList train) g + liftIO $ printf "[Epoch %d]\n" (e :: Int) + + forM_ ([1..] `zip` chunksOf batch train') $ \(b, chnk) -> StateT $ \n0 -> do + printf "(Batch %d)\n" (b :: Int) + + t0 <- getCurrentTime + n' <- evaluate . force $ trainList rate chnk n0 + t1 <- getCurrentTime + printf "Trained on %d points in %s.\n" batch (show (t1 `diffUTCTime` t0)) + + let trainScore = testNet chnk n' + testScore = testNet test n' + printf "Training error: %.2f%%\n" ((1 - trainScore) * 100) + printf "Validation error: %.2f%%\n" ((1 - testScore ) * 100) + + return ((), n') + where + rate = 0.02 + batch = 5000 +``` + +Each iteration of the loop: + +1. Shuffles the training set +2. Splits it into chunks of `batch` size +3. Uses `trainList` to train over the batch +4. Computes the score based on `testNet` based on the training set and + the test set +5. Prints out the results + +And, that’s really it! + +Result +------ + +I haven’t put much into optimizing the library yet, but the network +(with hidden layer sizes 300 and 100) seems to take 25s on my computer +to finish a batch of 5000 training points. It’s slow (five minutes per +60000 point epooch), but it’s a first unoptimized run and a proof of +concept! It’s my goal to get this down to a point where the result has +the same performance characteristics as the actual backend (*hmatrix*), +and so overhead is 0. + +Main takeaways +============== + +Most of the actual heavy lifting/logic actually came from the *hmatrix* +library itself. We just created simple types to wrap up our bare +matrices. + +Basically, all that *backprop* did was give you an API to define *how to +run* a neural net — how to *run* a net based on a `Network` and `R i` +input you were given. The goal of the library is to let you write down +how to run things in as natural way as possible. + +And then, after things are run, we can just get the gradient and roll +from there! + +Because the heavy lifting is done by the data types themselves, we can +presumably plug in *any* type and any tensor/numerical backend, and reap +the benefits of those libraries’ optimizations and parallelizations. +*Any* type can be backpropagated! :D + +What now? +--------- + +Check out the docs for the [Numeric.Backprop] module for a more detailed +picture of what’s going on, or find more examples at the [github repo]! + + [Numeric.Backprop]: http://hackage.haskell.org/package/backprop/docs/Numeric-Backprop.html + [github repo]: https://github.com/mstksg/backprop + +Boring stuff +============ + +Here is a small wrapper function over the [mnist-idx] library loading +the contents of the idx files into *hmatrix* vectors: + + [mnist-idx]: http://hackage.haskell.org/package/mnist-idx + +``` {.sourceCode .literate .haskell} +loadMNIST + :: FilePath + -> FilePath + -> IO (Maybe [(R 784, R 9)]) +loadMNIST fpI fpL = runMaybeT $ do + i <- MaybeT $ decodeIDXFile fpI + l <- MaybeT $ decodeIDXLabelsFile fpL + d <- MaybeT . return $ labeledIntData l i + r <- MaybeT . return $ for d (bitraverse mkImage mkLabel . swap) + liftIO . evaluate $ force r + where + mkImage :: VU.Vector Int -> Maybe (R 784) + mkImage = create . VG.convert . VG.map (\i -> fromIntegral i / 255) + mkLabel :: Int -> Maybe (R 9) + mkLabel n = create $ HM.build 9 (\i -> if round i == n then 1 else 0) +``` + +And here are instances to generating random +vectors/matrices/layers/networks, used for the initialization step. + +``` {.sourceCode .literate .haskell} +instance KnownNat n => MWC.Variate (R n) where + uniform g = randomVector <$> MWC.uniform g <*> pure Uniform + uniformR (l, h) g = (\x -> x * (h - l) + l) <$> MWC.uniform g + +instance (KnownNat m, KnownNat n) => MWC.Variate (L m n) where + uniform g = uniformSample <$> MWC.uniform g <*> pure 0 <*> pure 1 + uniformR (l, h) g = (\x -> x * (h - l) + l) <$> MWC.uniform g + +instance (KnownNat i, KnownNat o) => MWC.Variate (Layer i o) where + uniform g = Layer <$> MWC.uniform g <*> MWC.uniform g + uniformR (l, h) g = (\x -> x * (h - l) + l) <$> MWC.uniform g + +instance (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) => MWC.Variate (Network i h1 h2 o) where + uniform g = Net <$> MWC.uniform g <*> MWC.uniform g <*> MWC.uniform g + uniformR (l, h) g = (\x -> x * (h - l) + l) <$> MWC.uniform g +``` diff --git a/renders/backprop-mnist.pdf b/renders/backprop-mnist.pdf new file mode 100644 Binary files /dev/null and b/renders/backprop-mnist.pdf differ diff --git a/renders/backprop-neural-test.md b/renders/backprop-neural-test.md new file mode 100644 --- /dev/null +++ b/renders/backprop-neural-test.md @@ -0,0 +1,447 @@ +--- +author: +- Justin Le +fontfamily: 'palatino,cmtt' +geometry: margin=1in +links-as-notes: true +title: Neural networks with backprop library +--- + +The *backprop* library performs back-propagation over a *hetereogeneous* +system of relationships. It offers both an implicit ([ad]-like) and +explicit graph building usage style. Let’s use it to build neural +networks! + + [ad]: http://hackage.haskell.org/package/ad + +Repository source is [on github], and so are the [rendered unstable +docs]. + + [on github]: https://github.com/mstksg/backprop + [rendered unstable docs]: https://mstksg.github.io/backprop + +``` {.sourceCode .literate .haskell} +{-# LANGUAGE DeriveGeneric #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE LambdaCase #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE StandaloneDeriving #-} +{-# LANGUAGE TypeApplications #-} +{-# LANGUAGE TypeInType #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE ViewPatterns #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} +{-# OPTIONS_GHC -fno-warn-unused-top-binds #-} + +import Data.Functor +import Data.Kind +import Data.Maybe +import Data.Singletons +import Data.Singletons.Prelude +import Data.Singletons.TypeLits +import Data.Type.Combinator +import Data.Type.Product +import GHC.Generics (Generic) +import Numeric.Backprop +import Numeric.Backprop.Iso +import Numeric.LinearAlgebra.Static hiding (dot) +import System.Random.MWC +import qualified Generics.SOP as SOP +``` + +Ops +=== + +First, we define values of `Op` for the operations we want to do. `Op`s +are bundles of functions packaged with their hetereogeneous gradients. +For simple numeric functions, *backprop* can derive `Op`s automatically. +But for matrix operations, we have to derive them ourselves. + +The types help us with matching up the dimensions, but we still need to +be careful that our gradients are calculated correctly. + +`L` and `R` are matrix and vector types from the great *hmatrix* +library. + +First, matrix-vector multiplication: + +``` {.sourceCode .literate .haskell} +matVec + :: (KnownNat m, KnownNat n) + => Op '[ L m n, R n ] (R m) +matVec = op2' $ \m v -> ( m #> v + , \(fromMaybe 1 -> g) -> + (g `outer` v, tr m #> g) + ) +``` + +Now, dot products: + +``` {.sourceCode .literate .haskell} +dot :: KnownNat n + => Op '[ R n, R n ] Double +dot = op2' $ \x y -> ( x <.> y + , \case Nothing -> (y, x) + Just g -> (konst g * y, x * konst g) + ) +``` + +Polymorphic functions can be easily turned into `Op`s with `op1`/`op2` +etc., but they can also be run directly on graph nodes. + +``` {.sourceCode .literate .haskell} +logistic :: Floating a => a -> a +logistic x = 1 / (1 + exp (-x)) +``` + +A Simple Complete Example +========================= + +At this point, we already have enough to train a simple +single-hidden-layer neural network: + +``` {.sourceCode .literate .haskell} +simpleOp + :: (KnownNat m, KnownNat n, KnownNat o) + => R m + -> BPOpI s '[ L n m, R n, L o n, R o ] (R o) +simpleOp inp = \(w1 :< b1 :< w2 :< b2 :< Ø) -> + let z = logistic $ liftB2 matVec w1 x + b1 + in logistic $ liftB2 matVec w2 z + b2 + where + x = constVar inp +``` + +Here, `simpleOp` is defined in implicit (non-monadic) style, given a +tuple of inputs and returning outputs. Now `simpleOp` can be “run” with +the input vectors and parameters (a `L n m`, `R n`, `L o n`, and `R o`) +and calculate the output of the neural net. + +``` {.sourceCode .literate .haskell} +runSimple + :: (KnownNat m, KnownNat n, KnownNat o) + => R m + -> Tuple '[ L n m, R n, L o n, R o ] + -> R o +runSimple inp = evalBPOp (implicitly $ simpleOp inp) +``` + +Alternatively, we can define `simpleOp` in explicit monadic style, were +we specify our graph nodes explicitly. The results should be the same. + +``` {.sourceCode .literate .haskell} +simpleOpExplicit + :: (KnownNat m, KnownNat n, KnownNat o) + => R m + -> BPOp s '[ L n m, R n, L o n, R o ] (R o) +simpleOpExplicit inp = withInps $ \(w1 :< b1 :< w2 :< b2 :< Ø) -> do + -- First layer + y1 <- matVec ~$ (w1 :< x1 :< Ø) + let x2 = logistic (y1 + b1) + -- Second layer + y2 <- matVec ~$ (w2 :< x2 :< Ø) + return $ logistic (y2 + b2) + where + x1 = constVar inp +``` + +Now, for the magic of *backprop*: the library can now take advantage of +the implicit (or explicit) graph and use it to do back-propagation, too! + +``` {.sourceCode .literate .haskell} +simpleGrad + :: forall m n o. (KnownNat m, KnownNat n, KnownNat o) + => R m + -> R o + -> Tuple '[ L n m, R n, L o n, R o ] + -> Tuple '[ L n m, R n, L o n, R o ] +simpleGrad inp targ params = gradBPOp opError params + where + opError :: BPOp s '[ L n m, R n, L o n, R o ] Double + opError = do + res <- implicitly $ simpleOp inp + -- we explicitly bind err to prevent recomputation + err <- bindVar $ res - t + dot ~$ (err :< err :< Ø) + where + t = constVar targ +``` + +The result is the gradient of the input tuple’s components, with respect +to the `Double` result of `opError` (the squared error). We can then use +this gradient to do gradient descent. + +With Parameter Containers +========================= + +This method doesn’t quite scale, because we might want to make networks +with multiple layers and parameterize networks by layers. Let’s make +some basic container data types to help us organize our types, including +a recursive `Network` type that lets us chain multiple layers. + +``` {.sourceCode .literate .haskell} +data Layer :: Nat -> Nat -> Type where + Layer :: { _lWeights :: L m n + , _lBiases :: R m + } + -> Layer n m + deriving (Show, Generic) + + +data Network :: Nat -> [Nat] -> Nat -> Type where + NØ :: !(Layer a b) -> Network a '[] b + (:&) :: !(Layer a b) -> Network b bs c -> Network a (b ': bs) c +``` + +A `Layer n m` is a layer taking an n-vector and returning an m-vector. A +`Network a '[b, c, d] e` would be a Network that takes in an a-vector +and outputs an e-vector, with hidden layers of sizes b, c, and d. + +Isomorphisms +------------ + +The *backprop* library lets you apply operations on “parts” of data +types (like on the weights and biases of a `Layer`) by using `Iso`’s +(isomorphisms), like the ones from the *lens* library. The library +doesn’t depend on lens, but it can use the `Iso`s from the library and +also custom-defined ones. + +First, we can auto-generate isomorphisms using the *generics-sop* +library: + +``` {.sourceCode .literate .haskell} +instance SOP.Generic (Layer n m) +``` + +And then can create isomorphisms by hand for the two `Network` +constructors: + +``` {.sourceCode .literate .haskell} +netExternal :: Iso' (Network a '[] b) (Tuple '[Layer a b]) +netExternal = iso (\case NØ x -> x ::< Ø) + (\case I x :< Ø -> NØ x ) + +netInternal :: Iso' (Network a (b ': bs) c) (Tuple '[Layer a b, Network b bs c]) +netInternal = iso (\case x :& xs -> x ::< xs ::< Ø) + (\case I x :< I xs :< Ø -> x :& xs ) +``` + +An `Iso' a (Tuple as)` means that an `a` can really just be seen as a +tuple of `as`. + +Running a network +================= + +Now, we can write the `BPOp` that reprenents running the network and +getting a result. We pass in a `Sing bs` (a singleton list of the hidden +layer sizes) so that we can “pattern match” on the list and handle the +different network constructors differently. + +``` {.sourceCode .literate .haskell} +netOp + :: forall s a bs c. (KnownNat a, KnownNat c) + => Sing bs + -> BPOp s '[ R a, Network a bs c ] (R c) +netOp sbs = go sbs + where + go :: forall d es. KnownNat d + => Sing es + -> BPOp s '[ R d, Network d es c ] (R c) + go = \case + SNil -> withInps $ \(x :< n :< Ø) -> do + -- peek into the NØ using netExternal iso + l :< Ø <- netExternal #<~ n + -- run the 'layerOp' BP, with x and l as inputs + bpOp layerOp ~$ (x :< l :< Ø) + SNat `SCons` ses -> withInps $ \(x :< n :< Ø) -> withSingI ses $ do + -- peek into the (:&) using the netInternal iso + l :< n' :< Ø <- netInternal #<~ n + -- run the 'layerOp' BP, with x and l as inputs + z <- bpOp layerOp ~$ (x :< l :< Ø) + -- run the 'go ses' BP, with z and n as inputs + bpOp (go ses) ~$ (z :< n' :< Ø) + layerOp + :: forall d e. (KnownNat d, KnownNat e) + => BPOp s '[ R d, Layer d e ] (R e) + layerOp = withInps $ \(x :< l :< Ø) -> do + -- peek into the layer using the gTuple iso, auto-generated with SOP.Generic + w :< b :< Ø <- gTuple #<~ l + y <- matVec ~$ (w :< x :< Ø) + return $ logistic (y + b) +``` + +There’s some singletons work going on here, but it’s fairly standard +singletons stuff. Most of the complexity here is from the static typing +in our neural network type, and *not* from *backprop*. + +From *backprop* specifically, the only elements are `#<~` lets you +“split” an input ref with the given iso, and `bpOp`, which converts a +`BPOp` into an `Op` that you can bind with `~$`. + +Note that this library doesn’t support truly pattern matching on GADTs, +and that we had to pass in `Sing bs` as a reference to the structure of +our networks. + +Gradient Descent +---------------- + +Now we can do simple gradient descent. Defining an error function: + +``` {.sourceCode .literate .haskell} +errOp + :: KnownNat m + => R m + -> BVar s rs (R m) + -> BPOp s rs Double +errOp targ r = do + err <- bindVar $ r - t + dot ~$ (err :< err :< Ø) + where + t = constVar targ +``` + +And now, we can use `backprop` to generate the gradient, and shift the +`Network`! Things are made a bit cleaner from the fact that +`Network a bs c` has a `Num` instance, so we can use `(-)` and `(*)` +etc. + +``` {.sourceCode .literate .haskell} +train + :: (KnownNat a, SingI bs, KnownNat c) + => Double + -> R a + -> R c + -> Network a bs c + -> Network a bs c +train r x t n = case backprop (errOp t =<< netOp sing) (x ::< n ::< Ø) of + (_, _ :< I g :< Ø) -> n - (realToFrac r * g) +``` + +(`(::<)` is cons and `Ø` is nil for tuples.) + +Main +==== + +`main`, which will train on sample data sets, is still in progress! +Right now it just generates a random network using the *mwc-random* +library and prints each internal layer. + +``` {.sourceCode .literate .haskell} +main :: IO () +main = withSystemRandom $ \g -> do + n <- uniform @(Network 4 '[3,2] 1) g + void $ traverseNetwork sing (\l -> l <$ print l) n +``` + +Appendix: Boilerplate +===================== + +And now for some typeclass instances and boilerplates unrelated to the +*backprop* library that makes our custom types easier to use. + +``` {.sourceCode .literate .haskell} +instance KnownNat n => Variate (R n) where + uniform g = randomVector <$> uniform g <*> pure Uniform + uniformR (l, h) g = (\x -> x * (h - l) + l) <$> uniform g + +instance (KnownNat m, KnownNat n) => Variate (L m n) where + uniform g = uniformSample <$> uniform g <*> pure 0 <*> pure 1 + uniformR (l, h) g = (\x -> x * (h - l) + l) <$> uniform g + +instance (KnownNat n, KnownNat m) => Variate (Layer n m) where + uniform g = subtract 1 . (* 2) <$> (Layer <$> uniform g <*> uniform g) + uniformR (l, h) g = (\x -> x * (h - l) + l) <$> uniform g + +instance (KnownNat m, KnownNat n) => Num (Layer n m) where + Layer w1 b1 + Layer w2 b2 = Layer (w1 + w2) (b1 + b2) + Layer w1 b1 - Layer w2 b2 = Layer (w1 - w2) (b1 - b2) + Layer w1 b1 * Layer w2 b2 = Layer (w1 * w2) (b1 * b2) + abs (Layer w b) = Layer (abs w) (abs b) + signum (Layer w b) = Layer (signum w) (signum b) + negate (Layer w b) = Layer (negate w) (negate b) + fromInteger x = Layer (fromInteger x) (fromInteger x) + +instance (KnownNat m, KnownNat n) => Fractional (Layer n m) where + Layer w1 b1 / Layer w2 b2 = Layer (w1 / w2) (b1 / b2) + recip (Layer w b) = Layer (recip w) (recip b) + fromRational x = Layer (fromRational x) (fromRational x) + +instance (KnownNat a, SingI bs, KnownNat c) => Variate (Network a bs c) where + uniform g = genNet sing (uniform g) + uniformR (l, h) g = (\x -> x * (h - l) + l) <$> uniform g + +genNet + :: forall f a bs c. (Applicative f, KnownNat a, KnownNat c) + => Sing bs + -> (forall d e. (KnownNat d, KnownNat e) => f (Layer d e)) + -> f (Network a bs c) +genNet sbs f = go sbs + where + go :: forall d es. KnownNat d => Sing es -> f (Network d es c) + go = \case + SNil -> NØ <$> f + SNat `SCons` ses -> (:&) <$> f <*> go ses + +mapNetwork0 + :: forall a bs c. (KnownNat a, KnownNat c) + => Sing bs + -> (forall d e. (KnownNat d, KnownNat e) => Layer d e) + -> Network a bs c +mapNetwork0 sbs f = getI $ genNet sbs (I f) + +traverseNetwork + :: forall a bs c f. (KnownNat a, KnownNat c, Applicative f) + => Sing bs + -> (forall d e. (KnownNat d, KnownNat e) => Layer d e -> f (Layer d e)) + -> Network a bs c + -> f (Network a bs c) +traverseNetwork sbs f = go sbs + where + go :: forall d es. KnownNat d => Sing es -> Network d es c -> f (Network d es c) + go = \case + SNil -> \case + NØ x -> NØ <$> f x + SNat `SCons` ses -> \case + x :& xs -> (:&) <$> f x <*> go ses xs + +mapNetwork1 + :: forall a bs c. (KnownNat a, KnownNat c) + => Sing bs + -> (forall d e. (KnownNat d, KnownNat e) => Layer d e -> Layer d e) + -> Network a bs c + -> Network a bs c +mapNetwork1 sbs f = getI . traverseNetwork sbs (I . f) + +mapNetwork2 + :: forall a bs c. (KnownNat a, KnownNat c) + => Sing bs + -> (forall d e. (KnownNat d, KnownNat e) => Layer d e -> Layer d e -> Layer d e) + -> Network a bs c + -> Network a bs c + -> Network a bs c +mapNetwork2 sbs f = go sbs + where + go :: forall d es. KnownNat d => Sing es -> Network d es c -> Network d es c -> Network d es c + go = \case + SNil -> \case + NØ x -> \case + NØ y -> NØ (f x y) + SNat `SCons` ses -> \case + x :& xs -> \case + y :& ys -> f x y :& go ses xs ys + +instance (KnownNat a, SingI bs, KnownNat c) => Num (Network a bs c) where + (+) = mapNetwork2 sing (+) + (-) = mapNetwork2 sing (-) + (*) = mapNetwork2 sing (*) + negate = mapNetwork1 sing negate + abs = mapNetwork1 sing abs + signum = mapNetwork1 sing signum + fromInteger x = mapNetwork0 sing (fromInteger x) + +instance (KnownNat a, SingI bs, KnownNat c) => Fractional (Network a bs c) where + (/) = mapNetwork2 sing (/) + recip = mapNetwork1 sing recip + fromRational x = mapNetwork0 sing (fromRational x) +``` diff --git a/renders/backprop-neural-test.pdf b/renders/backprop-neural-test.pdf new file mode 100644 Binary files /dev/null and b/renders/backprop-neural-test.pdf differ diff --git a/samples/MNIST.lhs b/samples/MNIST.lhs deleted file mode 100644 --- a/samples/MNIST.lhs +++ /dev/null @@ -1,504 +0,0 @@ -% Learning MNIST with Neural Networks with backprop library -% Justin Le - -The *backprop* library performs back-propagation over a *hetereogeneous* -system of relationships. It offers both an implicit (*[ad][]*-like) and explicit graph -building usage style. Let's use it to build neural networks and learn -mnist! - -[ad]: http://hackage.haskell.org/package/ad - -Repository source is [on github][repo], and docs are [on hackage][hackage]. - -[repo]: https://github.com/mstksg/backprop -[hackage]: http://hackage.haskell.org/package/backprop - -If you're reading this as a literate haskell file, you should know that a -[rendered pdf version is available on github.][rendered]. If you are reading -this as a pdf file, you should know that a [literate haskell version that -you can run][lhs] is also available on github! - -[rendered]: https://github.com/mstksg/backprop/blob/master/renders/MNIST.pdf -[lhs]: https://github.com/mstksg/backprop/blob/master/samples/MNIST.lhs - - -> {-# LANGUAGE BangPatterns #-} -> {-# LANGUAGE DataKinds #-} -> {-# LANGUAGE DeriveGeneric #-} -> {-# LANGUAGE GADTs #-} -> {-# LANGUAGE LambdaCase #-} -> {-# LANGUAGE ScopedTypeVariables #-} -> {-# LANGUAGE TupleSections #-} -> {-# LANGUAGE TypeApplications #-} -> {-# LANGUAGE ViewPatterns #-} -> {-# OPTIONS_GHC -fno-warn-orphans #-} -> {-# OPTIONS_GHC -fno-warn-incomplete-patterns #-} -> {-# OPTIONS_GHC -fno-warn-unused-top-binds #-} -> -> import Control.DeepSeq -> import Control.Exception -> import Control.Monad -> import Control.Monad.IO.Class -> import Control.Monad.Trans.Maybe -> import Control.Monad.Trans.State -> import Data.Bitraversable -> import Data.Foldable -> import Data.IDX -> import Data.List.Split -> import Data.Maybe -> import Data.Time.Clock -> import Data.Traversable -> import Data.Tuple -> import GHC.Generics (Generic) -> import GHC.TypeLits -> import Numeric.Backprop -> import Numeric.LinearAlgebra.Static hiding (dot) -> import Text.Printf -> import qualified Data.Vector as V -> import qualified Data.Vector.Generic as VG -> import qualified Data.Vector.Unboxed as VU -> import qualified Generics.SOP as SOP -> import qualified Numeric.LinearAlgebra as HM -> import qualified System.Random.MWC as MWC -> import qualified System.Random.MWC.Distributions as MWC - -Types -===== - -For the most part, we're going to be using the great *[hmatrix][]* library -and its vector and matrix types. It offers a type `L m n` for $m \times n$ -matrices, and a type `R n` for an $n$ vector. - -[hmatrix]: http://hackage.haskell.org/package/hmatrix - -First things first: let's define our neural networks as simple containers -of parameters (weight matrices and bias vectors). - -First, a type for layers: - -> data Layer i o = -> Layer { _lWeights :: !(L o i) -> , _lBiases :: !(R o) -> } -> deriving (Show, Generic) -> -> instance SOP.Generic (Layer i o) -> instance NFData (Layer i o) - -And a type for a simple feed-forward network with two hidden layers: - -> data Network i h1 h2 o = -> Net { _nLayer1 :: !(Layer i h1) -> , _nLayer2 :: !(Layer h1 h2) -> , _nLayer3 :: !(Layer h2 o) -> } -> deriving (Show, Generic) -> -> instance SOP.Generic (Network i h1 h2 o) -> instance NFData (Network i h1 h2 o) - -These are pretty straightforward container types...pretty much exactly the -type you'd make to represent these networks! Note that, following true -Haskell form, we separate out logic from data. This should be all we need. - -We derive an instance of `SOP.Generic` from the *[generics-sop][]* package, -which *backprop* uses to propagate derivatives on values inside product -types. - -[generics-sop]: http://hackage.haskell.org/package/generics-sop - -Instances ---------- - -Things are much simplier if we had `Num` and `Fractional` instances for -everything, so let's just go ahead and define that now, as well. Just a -little bit of boilerplate. - -> instance (KnownNat i, KnownNat o) => Num (Layer i o) where -> Layer w1 b1 + Layer w2 b2 = Layer (w1 + w2) (b1 + b2) -> Layer w1 b1 - Layer w2 b2 = Layer (w1 - w2) (b1 - b2) -> Layer w1 b1 * Layer w2 b2 = Layer (w1 * w2) (b1 * b2) -> abs (Layer w b) = Layer (abs w) (abs b) -> signum (Layer w b) = Layer (signum w) (signum b) -> negate (Layer w b) = Layer (negate w) (negate b) -> fromInteger x = Layer (fromInteger x) (fromInteger x) -> -> instance (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) => Num (Network i h1 h2 o) where -> Net a b c + Net d e f = Net (a + d) (b + e) (c + f) -> Net a b c - Net d e f = Net (a - d) (b - e) (c - f) -> Net a b c * Net d e f = Net (a * d) (b * e) (c * f) -> abs (Net a b c) = Net (abs a) (abs b) (abs c) -> signum (Net a b c) = Net (signum a) (signum b) (signum c) -> negate (Net a b c) = Net (negate a) (negate b) (negate c) -> fromInteger x = Net (fromInteger x) (fromInteger x) (fromInteger x) -> -> instance (KnownNat i, KnownNat o) => Fractional (Layer i o) where -> Layer w1 b1 / Layer w2 b2 = Layer (w1 / w2) (b1 / b2) -> recip (Layer w b) = Layer (recip w) (recip b) -> fromRational x = Layer (fromRational x) (fromRational x) -> -> instance (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) => Fractional (Network i h1 h2 o) where -> Net a b c / Net d e f = Net (a / d) (b / e) (c / f) -> recip (Net a b c) = Net (recip a) (recip b) (recip c) -> fromRational x = Net (fromRational x) (fromRational x) (fromRational x) - -`KnownNat` comes from *base*; it's a typeclass that *hmatrix* uses to refer -to the numbers in its type and use it to go about its normal hmatrixy -business. - -Ops -=== - -Now, *backprop* does require *primitive* differentiable operations on our -relevant types to be defined. *backprop* uses these primitive `Op`s to tie -everything together. Ideally we'd import these from a library that -implements these for you, and the end-user never has to make `Op` -primitives. - -But in this case, I'm going to put the definitions here to show that there -isn't any magic going on. If you're curious, refer to [documentation for -`Op`][opdoc] for more details on how `Op` is implemented and how this -works. - -[opdoc]: http://hackage.haskell.org/package/backprop/docs/Numeric-Backprop-Op.html - -First, matrix-vector multiplication primitive, giving an explicit gradient -function. - -> matVec -> :: (KnownNat m, KnownNat n) -> => Op '[ L m n, R n ] (R m) -> matVec = op2' $ \m v -> -> ( m #> v, \(fromMaybe 1 -> g) -> -> (g `outer` v, tr m #> g) -> ) - -Dot products would be nice too. - -> dot :: KnownNat n -> => Op '[ R n, R n ] Double -> dot = op2' $ \x y -> -> ( x <.> y, \case Nothing -> (y, x) -> Just g -> (konst g * y, x * konst g) -> ) - -Also a "scaling" function, scales a vector by a given factor. - -> scale -> :: KnownNat n -> => Op '[ Double, R n ] (R n) -> scale = op2' $ \a x -> -> ( konst a * x -> , \case Nothing -> (HM.sumElements (extract x ), konst a ) -> Just g -> (HM.sumElements (extract (x * g)), konst a * g) -> ) - -Finally, an operation to sum all of the items in the vector. - -> vsum -> :: KnownNat n -> => Op '[ R n ] Double -> vsum = op1' $ \x -> (HM.sumElements (extract x), maybe 1 konst) - -And why not, here's the [logistic function][], which we'll use as an -activation function for internal layers. We don't need to define this as -an `Op` up-front right now, because the library can automatically promote -any numeric polymorphic function (an `a -> a` or `a -> a -> a`, etc.) to an -`Op` anyways. - -[logistic function]: https://en.wikipedia.org/wiki/Logistic_function - -> logistic :: Floating a => a -> a -> logistic x = 1 / (1 + exp (-x)) - -Running our Network -=================== - -Now that we have our primitives in place, let's actually write a function -to run our network! - -> runLayer -> :: (KnownNat i, KnownNat o) -> => BPOp s '[ R i, Layer i o ] (R o) -> runLayer = withInps $ \(x :< l :< Ø) -> do -> w :< b :< Ø <- gTuple #<~ l -> y <- matVec ~$ (w :< x :< Ø) -> return $ y + b - -A `BPOp s '[ R i, Layer i o ] (R o)` is a backpropagatable function that -produces an `R o` (a vector with `o` elements, from the *[hmatrix][]* -library) given an input environment of an `R i` (the "input" of the layer) -and a layer. - -We use `withInps` to bring the environment into scope as a bunch of -`BVar`s. `x` is a `BVar` containing the input vector, and `l` is a `BVar` -containing the layer. - -The first thing we do is split out the parts of the layer so we can work -with the internal matrices. We can use `#<~` to "split out" the components -of a `BVar`, splitting on `gTuple` (which uses `GHC.Generics` to -automatically figure out how to split up a product type). - -Then we apply `matVec` (our primitive `Op` that does matrix-vector -multiplication) to `w` and `x`, and then the result is that added to the -bias vector `b`. - -We can write the `runNetwork` function pretty much the same way. - -> runNetwork -> :: (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) -> => BPOp s '[ R i, Network i h1 h2 o ] (R o) -> runNetwork = withInps $ \(x :< n :< Ø) -> do -> l1 :< l2 :< l3 :< Ø <- gTuple #<~ n -> y <- runLayer -$ (x :< l1 :< Ø) -> z <- runLayer -$ (logistic y :< l2 :< Ø) -> r <- runLayer -$ (logistic z :< l3 :< Ø) -> softmax -$ (r :< Ø) -> where -> softmax :: KnownNat n => BPOp s '[ R n ] (R n) -> softmax = withInps $ \(x :< Ø) -> do -> expX <- bindVar (exp x) -> totX <- vsum ~$ (expX :< Ø) -> scale ~$ (1/totX :< expX :< Ø) - - -After splitting out the layers in the input `Network`, we run each layer -successively using our previously defined `runLayer`, giving inputs using -`-$`. We can directly apply `logistic` to `BVar`s. At the end, we run a -[softmax function][] because MNIST is a classification challenge. The softmax -is done by applying $e^x$ for every item in the input vector, and dividing -each element by the total. - -[softmax function]: https://en.wikipedia.org/wiki/Softmax_function - - -The Magic ---------- - -What did we just define? Well, with a `BPOp s rs a`, we can *run* it and -get the output: - -> runNetOnInp -> :: (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) -> => Network i h1 h2 o -> -> R i -> -> R o -> runNetOnInp n x = evalBPOp runNetwork (x ::< n ::< Ø) - -But, the magic part is that we can also get the gradient! - -> gradNet -> :: (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) -> => Network i h1 h2 o -> -> R i -> -> Network i h1 h2 o -> gradNet n x = case gradBPOp runNetwork (x ::< n ::< Ø) of -> _gradX ::< gradN ::< Ø -> gradN - -This gives the gradient of all of the parameters in the matrices and -vectors inside the `Network`, which we can use to "train"! - -Training -======== - -Now for the real work. To train a network, we can do gradient descent -based on the gradient of some type of *error function* with respect to the -network parameters. Let's use the [cross entropy][], which is popular for -classification problems. - -[cross entropy]: https://en.wikipedia.org/wiki/Cross_entropy - -> crossEntropy -> :: KnownNat n -> => R n -> -> BPOpI s '[ R n ] Double -> crossEntropy targ (r :< Ø) = negate (dot .$ (log r :< t :< Ø)) -> where -> t = constVar targ - -Given a target vector and a `BVar` referring to the result of the network, -we can directly apply: - -$$ -H(\mathbf{r}, \mathbf{t}) = - (log(\mathbf{r}) \cdot \mathbf{t}) -$$ - -Just for fun, I implemented `crossEntropy` in "implicit-graph" mode, so you -don't see any binds or returns. - -Now, a function to make one gradient descent step based on an input vector -and a target, using `gradBPOp`: - -> trainStep -> :: forall i h1 h2 o. (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) -> => Double -> -> R i -> -> R o -> -> Network i h1 h2 o -> -> Network i h1 h2 o -> trainStep r !x !t !n = case gradBPOp o (x ::< n ::< Ø) of -> _ ::< gN ::< Ø -> -> n - (realToFrac r * gN) -> where -> o :: BPOp s '[ R i, Network i h1 h2 o ] Double -> o = do -> y <- runNetwork -> implicitly (crossEntropy t) -$ (y :< Ø) - -A convenient wrapper for training over all of the observations in a list: - -> trainList -> :: (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) -> => Double -> -> [(R i, R o)] -> -> Network i h1 h2 o -> -> Network i h1 h2 o -> trainList r = flip $ foldl' (\n (x,y) -> trainStep r x y n) - -Pulling it all together -======================= - -`testNet` will be a quick way to test our net by computing the percentage -of correct guesses: (mostly using *hmatrix* stuff) - -> testNet -> :: forall i h1 h2 o. (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) -> => [(R i, R o)] -> -> Network i h1 h2 o -> -> Double -> testNet xs n = sum (map (uncurry test) xs) / fromIntegral (length xs) -> where -> test :: R i -> R o -> Double -> test x (extract->t) -> | HM.maxIndex t == HM.maxIndex (extract r) = 1 -> | otherwise = 0 -> where -> r :: R o -> r = evalBPOp runNetwork (x ::< n ::< Ø) - -And now, a main loop! - -If you are following along at home, download the [mnist data set -files][mnist] and uncompress them into the folder `data`, and everything -should work fine. - -[mnist]: http://yann.lecun.com/exdb/mnist/ - -> main :: IO () -> main = MWC.withSystemRandom $ \g -> do -> Just train <- loadMNIST "data/train-images-idx3-ubyte" "data/train-labels-idx1-ubyte" -> Just test <- loadMNIST "data/t10k-images-idx3-ubyte" "data/t10k-labels-idx1-ubyte" -> putStrLn "Loaded data." -> net0 <- MWC.uniformR @(Network 784 300 100 9) (-0.5, 0.5) g -> flip evalStateT net0 . forM_ [1..] $ \e -> do -> train' <- liftIO . fmap V.toList $ MWC.uniformShuffle (V.fromList train) g -> liftIO $ printf "[Epoch %d]\n" (e :: Int) -> -> forM_ ([1..] `zip` chunksOf batch train') $ \(b, chnk) -> StateT $ \n0 -> do -> printf "(Batch %d)\n" (b :: Int) -> -> t0 <- getCurrentTime -> n' <- evaluate . force $ trainList rate chnk n0 -> t1 <- getCurrentTime -> printf "Trained on %d points in %s.\n" batch (show (t1 `diffUTCTime` t0)) -> -> let trainScore = testNet chnk n' -> testScore = testNet test n' -> printf "Training error: %.2f%%\n" ((1 - trainScore) * 100) -> printf "Validation error: %.2f%%\n" ((1 - testScore ) * 100) -> -> return ((), n') -> where -> rate = 0.02 -> batch = 5000 - -Each iteration of the loop: - -1. Shuffles the training set -2. Splits it into chunks of `batch` size -3. Uses `trainList` to train over the batch -4. Computes the score based on `testNet` based on the training set and the - test set -5. Prints out the results - -And, that's really it! - -Result ------- - -I haven't put much into optimizing the library yet, but the network (with -hidden layer sizes 300 and 100) seems to take 25s on my computer to finish -a batch of 5000 training points. It's slow (five minutes per 60000 point -epooch), but it's a first unoptimized run and a proof of concept! It's my -goal to get this down to a point where the result has the same performance -characteristics as the actual backend (*hmatrix*), and so overhead is 0. - -Main takeaways -============== - -Most of the actual heavy lifting/logic actually came from the *hmatrix* -library itself. We just created simple types to wrap up our bare matrices. - -Basically, all that *backprop* did was give you an API to define *how to -run* a neural net --- how to *run* a net based on a `Network` and `R i` input -you were given. The goal of the library is to let you write down how to -run things in as natural way as possible. - -And then, after things are run, we can just get the gradient and roll from -there! - -Because the heavy lifting is done by the data types themselves, we can -presumably plug in *any* type and any tensor/numerical backend, and reap -the benefits of those libraries' optimizations and parallelizations. *Any* -type can be backpropagated! :D - -What now? ---------- - -Check out the docs for the [Numeric.Backprop][] module for a more detailed -picture of what's going on, or find more examples at the [github repo][repo]! - -[Numeric.Backprop]: http://hackage.haskell.org/package/backprop/docs/Numeric-Backprop.html - -Boring stuff -============ - -Here is a small wrapper function over the [mnist-idx][] library loading the -contents of the idx files into *hmatrix* vectors: - -[mnist-idx]: http://hackage.haskell.org/package/mnist-idx - -> loadMNIST -> :: FilePath -> -> FilePath -> -> IO (Maybe [(R 784, R 9)]) -> loadMNIST fpI fpL = runMaybeT $ do -> i <- MaybeT $ decodeIDXFile fpI -> l <- MaybeT $ decodeIDXLabelsFile fpL -> d <- MaybeT . return $ labeledIntData l i -> r <- MaybeT . return $ for d (bitraverse mkImage mkLabel . swap) -> liftIO . evaluate $ force r -> where -> mkImage :: VU.Vector Int -> Maybe (R 784) -> mkImage = create . VG.convert . VG.map (\i -> fromIntegral i / 255) -> mkLabel :: Int -> Maybe (R 9) -> mkLabel n = create $ HM.build 9 (\i -> if round i == n then 1 else 0) - -And here are instances to generating random -vectors/matrices/layers/networks, used for the initialization step. - -> instance KnownNat n => MWC.Variate (R n) where -> uniform g = randomVector <$> MWC.uniform g <*> pure Uniform -> uniformR (l, h) g = (\x -> x * (h - l) + l) <$> MWC.uniform g -> -> instance (KnownNat m, KnownNat n) => MWC.Variate (L m n) where -> uniform g = uniformSample <$> MWC.uniform g <*> pure 0 <*> pure 1 -> uniformR (l, h) g = (\x -> x * (h - l) + l) <$> MWC.uniform g -> -> instance (KnownNat i, KnownNat o) => MWC.Variate (Layer i o) where -> uniform g = Layer <$> MWC.uniform g <*> MWC.uniform g -> uniformR (l, h) g = (\x -> x * (h - l) + l) <$> MWC.uniform g -> -> instance (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) => MWC.Variate (Network i h1 h2 o) where -> uniform g = Net <$> MWC.uniform g <*> MWC.uniform g <*> MWC.uniform g -> uniformR (l, h) g = (\x -> x * (h - l) + l) <$> MWC.uniform g diff --git a/samples/MonoTest.hs b/samples/MonoTest.hs deleted file mode 100644 --- a/samples/MonoTest.hs +++ /dev/null @@ -1,18 +0,0 @@ -{-# LANGUAGE GADTs #-} - -import Numeric.Backprop.Mono - -testImplicit :: BPOp s N3 Double Double -testImplicit = implicitly $ \(x :* y :* z :* ØV) -> - ((x * y) + y) * z - -testExplicit :: BPOp s N3 Double Double -testExplicit = withInps $ \(x :* y :* z :* ØV) -> do - xy <- op2 (*) ~$ (x :* y :* ØV) - xyy <- op2 (+) ~$ (xy :* y :* ØV) - op2 (*) ~$ (xyy :* z :* ØV) - -main :: IO () -main = do - print $ backprop testImplicit (2 :+ 3 :+ 4 :+ ØV) - print $ backprop testExplicit (2 :+ 3 :+ 4 :+ ØV) diff --git a/samples/NeuralTest.lhs b/samples/NeuralTest.lhs deleted file mode 100644 --- a/samples/NeuralTest.lhs +++ /dev/null @@ -1,405 +0,0 @@ -% Neural networks with backprop library -% Justin Le - -The *backprop* library performs back-propagation over a *hetereogeneous* -system of relationships. It offers both an implicit ([ad][]-like) and explicit graph -building usage style. Let's use it to build neural networks! - -[ad]: http://hackage.haskell.org/package/ad - -Repository source is [on github][repo], and so are the [rendered unstable -docs][docs]. - -[repo]: https://github.com/mstksg/backprop -[docs]: https://mstksg.github.io/backprop - -> {-# LANGUAGE DeriveGeneric #-} -> {-# LANGUAGE GADTs #-} -> {-# LANGUAGE LambdaCase #-} -> {-# LANGUAGE RankNTypes #-} -> {-# LANGUAGE ScopedTypeVariables #-} -> {-# LANGUAGE StandaloneDeriving #-} -> {-# LANGUAGE TypeApplications #-} -> {-# LANGUAGE TypeInType #-} -> {-# LANGUAGE TypeOperators #-} -> {-# LANGUAGE ViewPatterns #-} -> {-# OPTIONS_GHC -fno-warn-orphans #-} -> {-# OPTIONS_GHC -fno-warn-unused-top-binds #-} -> -> import Data.Functor -> import Data.Kind -> import Data.Maybe -> import Data.Singletons -> import Data.Singletons.Prelude -> import Data.Singletons.TypeLits -> import Data.Type.Combinator -> import Data.Type.Product -> import GHC.Generics (Generic) -> import Numeric.Backprop -> import Numeric.Backprop.Iso -> import Numeric.LinearAlgebra.Static hiding (dot) -> import System.Random.MWC -> import qualified Generics.SOP as SOP - -Ops -=== - -First, we define values of `Op` for the operations we want to do. `Op`s -are bundles of functions packaged with their hetereogeneous gradients. For -simple numeric functions, *backprop* can derive `Op`s automatically. But -for matrix operations, we have to derive them ourselves. - -The types help us with matching up the dimensions, but we still need to be -careful that our gradients are calculated correctly. - -`L` and `R` are matrix and vector types from the great *hmatrix* library. - -First, matrix-vector multiplication: - -> matVec -> :: (KnownNat m, KnownNat n) -> => Op '[ L m n, R n ] (R m) -> matVec = op2' $ \m v -> ( m #> v -> , \(fromMaybe 1 -> g) -> -> (g `outer` v, tr m #> g) -> ) - -Now, dot products: - -> dot :: KnownNat n -> => Op '[ R n, R n ] Double -> dot = op2' $ \x y -> ( x <.> y -> , \case Nothing -> (y, x) -> Just g -> (konst g * y, x * konst g) -> ) - -Polymorphic functions can be easily turned into `Op`s with `op1`/`op2` -etc., but they can also be run directly on graph nodes. - -> logistic :: Floating a => a -> a -> logistic x = 1 / (1 + exp (-x)) - -A Simple Complete Example -========================= - -At this point, we already have enough to train a simple single-hidden-layer -neural network: - -> simpleOp -> :: (KnownNat m, KnownNat n, KnownNat o) -> => R m -> -> BPOpI s '[ L n m, R n, L o n, R o ] (R o) -> simpleOp inp = \(w1 :< b1 :< w2 :< b2 :< Ø) -> -> let z = logistic $ liftB2 matVec w1 x + b1 -> in logistic $ liftB2 matVec w2 z + b2 -> where -> x = constVar inp - -Here, `simpleOp` is defined in implicit (non-monadic) style, given a tuple -of inputs and returning outputs. Now `simpleOp` can be "run" with the -input vectors and parameters (a `L n m`, `R n`, `L o n`, and `R o`) and -calculate the output of the neural net. - -> runSimple -> :: (KnownNat m, KnownNat n, KnownNat o) -> => R m -> -> Tuple '[ L n m, R n, L o n, R o ] -> -> R o -> runSimple inp = evalBPOp (implicitly $ simpleOp inp) - -Alternatively, we can define `simpleOp` in explicit monadic style, were we -specify our graph nodes explicitly. The results should be the same. - -> simpleOpExplicit -> :: (KnownNat m, KnownNat n, KnownNat o) -> => R m -> -> BPOp s '[ L n m, R n, L o n, R o ] (R o) -> simpleOpExplicit inp = withInps $ \(w1 :< b1 :< w2 :< b2 :< Ø) -> do -> -- First layer -> y1 <- matVec ~$ (w1 :< x1 :< Ø) -> let x2 = logistic (y1 + b1) -> -- Second layer -> y2 <- matVec ~$ (w2 :< x2 :< Ø) -> return $ logistic (y2 + b2) -> where -> x1 = constVar inp - -Now, for the magic of *backprop*: the library can now take advantage of -the implicit (or explicit) graph and use it to do back-propagation, too! - -> simpleGrad -> :: forall m n o. (KnownNat m, KnownNat n, KnownNat o) -> => R m -> -> R o -> -> Tuple '[ L n m, R n, L o n, R o ] -> -> Tuple '[ L n m, R n, L o n, R o ] -> simpleGrad inp targ params = gradBPOp opError params -> where -> opError :: BPOp s '[ L n m, R n, L o n, R o ] Double -> opError = do -> res <- implicitly $ simpleOp inp -> -- we explicitly bind err to prevent recomputation -> err <- bindVar $ res - t -> dot ~$ (err :< err :< Ø) -> where -> t = constVar targ - -The result is the gradient of the input tuple's components, with respect -to the `Double` result of `opError` (the squared error). We can then use -this gradient to do gradient descent. - -With Parameter Containers -========================= - -This method doesn't quite scale, because we might want to make networks -with multiple layers and parameterize networks by layers. Let's make some -basic container data types to help us organize our types, including a -recursive `Network` type that lets us chain multiple layers. - -> data Layer :: Nat -> Nat -> Type where -> Layer :: { _lWeights :: L m n -> , _lBiases :: R m -> } -> -> Layer n m -> deriving (Show, Generic) -> -> -> data Network :: Nat -> [Nat] -> Nat -> Type where -> NØ :: !(Layer a b) -> Network a '[] b -> (:&) :: !(Layer a b) -> Network b bs c -> Network a (b ': bs) c - -A `Layer n m` is a layer taking an n-vector and returning an m-vector. A -`Network a '[b, c, d] e` would be a Network that takes in an a-vector and -outputs an e-vector, with hidden layers of sizes b, c, and d. - -Isomorphisms ------------- - -The *backprop* library lets you apply operations on "parts" of data types -(like on the weights and biases of a `Layer`) by using `Iso`'s -(isomorphisms), like the ones from the *lens* library. The library doesn't -depend on lens, but it can use the `Iso`s from the library and also -custom-defined ones. - -First, we can auto-generate isomorphisms using the *generics-sop* library: - -> instance SOP.Generic (Layer n m) - -And then can create isomorphisms by hand for the two `Network` -constructors: - -> netExternal :: Iso' (Network a '[] b) (Tuple '[Layer a b]) -> netExternal = iso (\case NØ x -> x ::< Ø) -> (\case I x :< Ø -> NØ x ) -> -> netInternal :: Iso' (Network a (b ': bs) c) (Tuple '[Layer a b, Network b bs c]) -> netInternal = iso (\case x :& xs -> x ::< xs ::< Ø) -> (\case I x :< I xs :< Ø -> x :& xs ) - -An `Iso' a (Tuple as)` means that an `a` can really just be seen as a tuple -of `as`. - -Running a network -================= - -Now, we can write the `BPOp` that reprenents running the network and -getting a result. We pass in a `Sing bs` (a singleton list of the hidden -layer sizes) so that we can "pattern match" on the list and handle the -different network constructors differently. - -> netOp -> :: forall s a bs c. (KnownNat a, KnownNat c) -> => Sing bs -> -> BPOp s '[ R a, Network a bs c ] (R c) -> netOp sbs = go sbs -> where -> go :: forall d es. KnownNat d -> => Sing es -> -> BPOp s '[ R d, Network d es c ] (R c) -> go = \case -> SNil -> withInps $ \(x :< n :< Ø) -> do -> -- peek into the NØ using netExternal iso -> l :< Ø <- netExternal #<~ n -> -- run the 'layerOp' BP, with x and l as inputs -> bpOp layerOp ~$ (x :< l :< Ø) -> SNat `SCons` ses -> withInps $ \(x :< n :< Ø) -> withSingI ses $ do -> -- peek into the (:&) using the netInternal iso -> l :< n' :< Ø <- netInternal #<~ n -> -- run the 'layerOp' BP, with x and l as inputs -> z <- bpOp layerOp ~$ (x :< l :< Ø) -> -- run the 'go ses' BP, with z and n as inputs -> bpOp (go ses) ~$ (z :< n' :< Ø) -> layerOp -> :: forall d e. (KnownNat d, KnownNat e) -> => BPOp s '[ R d, Layer d e ] (R e) -> layerOp = withInps $ \(x :< l :< Ø) -> do -> -- peek into the layer using the gTuple iso, auto-generated with SOP.Generic -> w :< b :< Ø <- gTuple #<~ l -> y <- matVec ~$ (w :< x :< Ø) -> return $ logistic (y + b) - -There's some singletons work going on here, but it's fairly standard -singletons stuff. Most of the complexity here is from the static typing in -our neural network type, and *not* from *backprop*. - -From *backprop* specifically, the only elements are `#<~` lets you "split" an -input ref with the given iso, and `bpOp`, which converts a `BPOp` into an `Op` -that you can bind with `~$`. - -Note that this library doesn't support truly pattern matching on GADTs, and -that we had to pass in `Sing bs` as a reference to the structure of our -networks. - -Gradient Descent ----------------- - -Now we can do simple gradient descent. Defining an error function: - -> errOp -> :: KnownNat m -> => R m -> -> BVar s rs (R m) -> -> BPOp s rs Double -> errOp targ r = do -> err <- bindVar $ r - t -> dot ~$ (err :< err :< Ø) -> where -> t = constVar targ - -And now, we can use `backprop` to generate the gradient, and shift the -`Network`! Things are made a bit cleaner from the fact that `Network a bs c` -has a `Num` instance, so we can use `(-)` and `(*)` etc. - -> train -> :: (KnownNat a, SingI bs, KnownNat c) -> => Double -> -> R a -> -> R c -> -> Network a bs c -> -> Network a bs c -> train r x t n = case backprop (errOp t =<< netOp sing) (x ::< n ::< Ø) of -> (_, _ :< I g :< Ø) -> n - (realToFrac r * g) - -(`(::<)` is cons and `Ø` is nil for tuples.) - -Main -==== - -`main`, which will train on sample data sets, is still in progress! Right -now it just generates a random network using the *mwc-random* library and -prints each internal layer. - -> main :: IO () -> main = withSystemRandom $ \g -> do -> n <- uniform @(Network 4 '[3,2] 1) g -> void $ traverseNetwork sing (\l -> l <$ print l) n - -Appendix: Boilerplate -===================== - -And now for some typeclass instances and boilerplates unrelated to the -*backprop* library that makes our custom types easier to use. - -> instance KnownNat n => Variate (R n) where -> uniform g = randomVector <$> uniform g <*> pure Uniform -> uniformR (l, h) g = (\x -> x * (h - l) + l) <$> uniform g -> -> instance (KnownNat m, KnownNat n) => Variate (L m n) where -> uniform g = uniformSample <$> uniform g <*> pure 0 <*> pure 1 -> uniformR (l, h) g = (\x -> x * (h - l) + l) <$> uniform g -> -> instance (KnownNat n, KnownNat m) => Variate (Layer n m) where -> uniform g = subtract 1 . (* 2) <$> (Layer <$> uniform g <*> uniform g) -> uniformR (l, h) g = (\x -> x * (h - l) + l) <$> uniform g -> -> instance (KnownNat m, KnownNat n) => Num (Layer n m) where -> Layer w1 b1 + Layer w2 b2 = Layer (w1 + w2) (b1 + b2) -> Layer w1 b1 - Layer w2 b2 = Layer (w1 - w2) (b1 - b2) -> Layer w1 b1 * Layer w2 b2 = Layer (w1 * w2) (b1 * b2) -> abs (Layer w b) = Layer (abs w) (abs b) -> signum (Layer w b) = Layer (signum w) (signum b) -> negate (Layer w b) = Layer (negate w) (negate b) -> fromInteger x = Layer (fromInteger x) (fromInteger x) -> -> instance (KnownNat m, KnownNat n) => Fractional (Layer n m) where -> Layer w1 b1 / Layer w2 b2 = Layer (w1 / w2) (b1 / b2) -> recip (Layer w b) = Layer (recip w) (recip b) -> fromRational x = Layer (fromRational x) (fromRational x) -> -> instance (KnownNat a, SingI bs, KnownNat c) => Variate (Network a bs c) where -> uniform g = genNet sing (uniform g) -> uniformR (l, h) g = (\x -> x * (h - l) + l) <$> uniform g -> -> genNet -> :: forall f a bs c. (Applicative f, KnownNat a, KnownNat c) -> => Sing bs -> -> (forall d e. (KnownNat d, KnownNat e) => f (Layer d e)) -> -> f (Network a bs c) -> genNet sbs f = go sbs -> where -> go :: forall d es. KnownNat d => Sing es -> f (Network d es c) -> go = \case -> SNil -> NØ <$> f -> SNat `SCons` ses -> (:&) <$> f <*> go ses -> -> mapNetwork0 -> :: forall a bs c. (KnownNat a, KnownNat c) -> => Sing bs -> -> (forall d e. (KnownNat d, KnownNat e) => Layer d e) -> -> Network a bs c -> mapNetwork0 sbs f = getI $ genNet sbs (I f) -> -> traverseNetwork -> :: forall a bs c f. (KnownNat a, KnownNat c, Applicative f) -> => Sing bs -> -> (forall d e. (KnownNat d, KnownNat e) => Layer d e -> f (Layer d e)) -> -> Network a bs c -> -> f (Network a bs c) -> traverseNetwork sbs f = go sbs -> where -> go :: forall d es. KnownNat d => Sing es -> Network d es c -> f (Network d es c) -> go = \case -> SNil -> \case -> NØ x -> NØ <$> f x -> SNat `SCons` ses -> \case -> x :& xs -> (:&) <$> f x <*> go ses xs -> -> mapNetwork1 -> :: forall a bs c. (KnownNat a, KnownNat c) -> => Sing bs -> -> (forall d e. (KnownNat d, KnownNat e) => Layer d e -> Layer d e) -> -> Network a bs c -> -> Network a bs c -> mapNetwork1 sbs f = getI . traverseNetwork sbs (I . f) -> -> mapNetwork2 -> :: forall a bs c. (KnownNat a, KnownNat c) -> => Sing bs -> -> (forall d e. (KnownNat d, KnownNat e) => Layer d e -> Layer d e -> Layer d e) -> -> Network a bs c -> -> Network a bs c -> -> Network a bs c -> mapNetwork2 sbs f = go sbs -> where -> go :: forall d es. KnownNat d => Sing es -> Network d es c -> Network d es c -> Network d es c -> go = \case -> SNil -> \case -> NØ x -> \case -> NØ y -> NØ (f x y) -> SNat `SCons` ses -> \case -> x :& xs -> \case -> y :& ys -> f x y :& go ses xs ys -> -> instance (KnownNat a, SingI bs, KnownNat c) => Num (Network a bs c) where -> (+) = mapNetwork2 sing (+) -> (-) = mapNetwork2 sing (-) -> (*) = mapNetwork2 sing (*) -> negate = mapNetwork1 sing negate -> abs = mapNetwork1 sing abs -> signum = mapNetwork1 sing signum -> fromInteger x = mapNetwork0 sing (fromInteger x) -> -> instance (KnownNat a, SingI bs, KnownNat c) => Fractional (Network a bs c) where -> (/) = mapNetwork2 sing (/) -> recip = mapNetwork1 sing recip -> fromRational x = mapNetwork0 sing (fromRational x) diff --git a/samples/backprop-mnist.lhs b/samples/backprop-mnist.lhs new file mode 100644 --- /dev/null +++ b/samples/backprop-mnist.lhs @@ -0,0 +1,504 @@ +% Learning MNIST with Neural Networks with backprop library +% Justin Le + +The *backprop* library performs back-propagation over a *hetereogeneous* +system of relationships. It offers both an implicit (*[ad][]*-like) and explicit graph +building usage style. Let's use it to build neural networks and learn +mnist! + +[ad]: http://hackage.haskell.org/package/ad + +Repository source is [on github][repo], and docs are [on hackage][hackage]. + +[repo]: https://github.com/mstksg/backprop +[hackage]: http://hackage.haskell.org/package/backprop + +If you're reading this as a literate haskell file, you should know that a +[rendered pdf version is available on github.][rendered]. If you are reading +this as a pdf file, you should know that a [literate haskell version that +you can run][lhs] is also available on github! + +[rendered]: https://github.com/mstksg/backprop/blob/master/renders/backprop-mnist.pdf +[lhs]: https://github.com/mstksg/backprop/blob/master/samples/backprop-mnist.lhs + + +> {-# LANGUAGE BangPatterns #-} +> {-# LANGUAGE DataKinds #-} +> {-# LANGUAGE DeriveGeneric #-} +> {-# LANGUAGE GADTs #-} +> {-# LANGUAGE LambdaCase #-} +> {-# LANGUAGE ScopedTypeVariables #-} +> {-# LANGUAGE TupleSections #-} +> {-# LANGUAGE TypeApplications #-} +> {-# LANGUAGE ViewPatterns #-} +> {-# OPTIONS_GHC -fno-warn-orphans #-} +> {-# OPTIONS_GHC -fno-warn-incomplete-patterns #-} +> {-# OPTIONS_GHC -fno-warn-unused-top-binds #-} +> +> import Control.DeepSeq +> import Control.Exception +> import Control.Monad +> import Control.Monad.IO.Class +> import Control.Monad.Trans.Maybe +> import Control.Monad.Trans.State +> import Data.Bitraversable +> import Data.Foldable +> import Data.IDX +> import Data.List.Split +> import Data.Maybe +> import Data.Time.Clock +> import Data.Traversable +> import Data.Tuple +> import GHC.Generics (Generic) +> import GHC.TypeLits +> import Numeric.Backprop +> import Numeric.LinearAlgebra.Static hiding (dot) +> import Text.Printf +> import qualified Data.Vector as V +> import qualified Data.Vector.Generic as VG +> import qualified Data.Vector.Unboxed as VU +> import qualified Generics.SOP as SOP +> import qualified Numeric.LinearAlgebra as HM +> import qualified System.Random.MWC as MWC +> import qualified System.Random.MWC.Distributions as MWC + +Types +===== + +For the most part, we're going to be using the great *[hmatrix][]* library +and its vector and matrix types. It offers a type `L m n` for $m \times n$ +matrices, and a type `R n` for an $n$ vector. + +[hmatrix]: http://hackage.haskell.org/package/hmatrix + +First things first: let's define our neural networks as simple containers +of parameters (weight matrices and bias vectors). + +First, a type for layers: + +> data Layer i o = +> Layer { _lWeights :: !(L o i) +> , _lBiases :: !(R o) +> } +> deriving (Show, Generic) +> +> instance SOP.Generic (Layer i o) +> instance NFData (Layer i o) + +And a type for a simple feed-forward network with two hidden layers: + +> data Network i h1 h2 o = +> Net { _nLayer1 :: !(Layer i h1) +> , _nLayer2 :: !(Layer h1 h2) +> , _nLayer3 :: !(Layer h2 o) +> } +> deriving (Show, Generic) +> +> instance SOP.Generic (Network i h1 h2 o) +> instance NFData (Network i h1 h2 o) + +These are pretty straightforward container types...pretty much exactly the +type you'd make to represent these networks! Note that, following true +Haskell form, we separate out logic from data. This should be all we need. + +We derive an instance of `SOP.Generic` from the *[generics-sop][]* package, +which *backprop* uses to propagate derivatives on values inside product +types. + +[generics-sop]: http://hackage.haskell.org/package/generics-sop + +Instances +--------- + +Things are much simplier if we had `Num` and `Fractional` instances for +everything, so let's just go ahead and define that now, as well. Just a +little bit of boilerplate. + +> instance (KnownNat i, KnownNat o) => Num (Layer i o) where +> Layer w1 b1 + Layer w2 b2 = Layer (w1 + w2) (b1 + b2) +> Layer w1 b1 - Layer w2 b2 = Layer (w1 - w2) (b1 - b2) +> Layer w1 b1 * Layer w2 b2 = Layer (w1 * w2) (b1 * b2) +> abs (Layer w b) = Layer (abs w) (abs b) +> signum (Layer w b) = Layer (signum w) (signum b) +> negate (Layer w b) = Layer (negate w) (negate b) +> fromInteger x = Layer (fromInteger x) (fromInteger x) +> +> instance (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) => Num (Network i h1 h2 o) where +> Net a b c + Net d e f = Net (a + d) (b + e) (c + f) +> Net a b c - Net d e f = Net (a - d) (b - e) (c - f) +> Net a b c * Net d e f = Net (a * d) (b * e) (c * f) +> abs (Net a b c) = Net (abs a) (abs b) (abs c) +> signum (Net a b c) = Net (signum a) (signum b) (signum c) +> negate (Net a b c) = Net (negate a) (negate b) (negate c) +> fromInteger x = Net (fromInteger x) (fromInteger x) (fromInteger x) +> +> instance (KnownNat i, KnownNat o) => Fractional (Layer i o) where +> Layer w1 b1 / Layer w2 b2 = Layer (w1 / w2) (b1 / b2) +> recip (Layer w b) = Layer (recip w) (recip b) +> fromRational x = Layer (fromRational x) (fromRational x) +> +> instance (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) => Fractional (Network i h1 h2 o) where +> Net a b c / Net d e f = Net (a / d) (b / e) (c / f) +> recip (Net a b c) = Net (recip a) (recip b) (recip c) +> fromRational x = Net (fromRational x) (fromRational x) (fromRational x) + +`KnownNat` comes from *base*; it's a typeclass that *hmatrix* uses to refer +to the numbers in its type and use it to go about its normal hmatrixy +business. + +Ops +=== + +Now, *backprop* does require *primitive* differentiable operations on our +relevant types to be defined. *backprop* uses these primitive `Op`s to tie +everything together. Ideally we'd import these from a library that +implements these for you, and the end-user never has to make `Op` +primitives. + +But in this case, I'm going to put the definitions here to show that there +isn't any magic going on. If you're curious, refer to [documentation for +`Op`][opdoc] for more details on how `Op` is implemented and how this +works. + +[opdoc]: http://hackage.haskell.org/package/backprop/docs/Numeric-Backprop-Op.html + +First, matrix-vector multiplication primitive, giving an explicit gradient +function. + +> matVec +> :: (KnownNat m, KnownNat n) +> => Op '[ L m n, R n ] (R m) +> matVec = op2' $ \m v -> +> ( m #> v, \(fromMaybe 1 -> g) -> +> (g `outer` v, tr m #> g) +> ) + +Dot products would be nice too. + +> dot :: KnownNat n +> => Op '[ R n, R n ] Double +> dot = op2' $ \x y -> +> ( x <.> y, \case Nothing -> (y, x) +> Just g -> (konst g * y, x * konst g) +> ) + +Also a "scaling" function, scales a vector by a given factor. + +> scale +> :: KnownNat n +> => Op '[ Double, R n ] (R n) +> scale = op2' $ \a x -> +> ( konst a * x +> , \case Nothing -> (HM.sumElements (extract x ), konst a ) +> Just g -> (HM.sumElements (extract (x * g)), konst a * g) +> ) + +Finally, an operation to sum all of the items in the vector. + +> vsum +> :: KnownNat n +> => Op '[ R n ] Double +> vsum = op1' $ \x -> (HM.sumElements (extract x), maybe 1 konst) + +And why not, here's the [logistic function][], which we'll use as an +activation function for internal layers. We don't need to define this as +an `Op` up-front right now, because the library can automatically promote +any numeric polymorphic function (an `a -> a` or `a -> a -> a`, etc.) to an +`Op` anyways. + +[logistic function]: https://en.wikipedia.org/wiki/Logistic_function + +> logistic :: Floating a => a -> a +> logistic x = 1 / (1 + exp (-x)) + +Running our Network +=================== + +Now that we have our primitives in place, let's actually write a function +to run our network! + +> runLayer +> :: (KnownNat i, KnownNat o) +> => BPOp s '[ R i, Layer i o ] (R o) +> runLayer = withInps $ \(x :< l :< Ø) -> do +> w :< b :< Ø <- gTuple #<~ l +> y <- matVec ~$ (w :< x :< Ø) +> return $ y + b + +A `BPOp s '[ R i, Layer i o ] (R o)` is a backpropagatable function that +produces an `R o` (a vector with `o` elements, from the *[hmatrix][]* +library) given an input environment of an `R i` (the "input" of the layer) +and a layer. + +We use `withInps` to bring the environment into scope as a bunch of +`BVar`s. `x` is a `BVar` containing the input vector, and `l` is a `BVar` +containing the layer. + +The first thing we do is split out the parts of the layer so we can work +with the internal matrices. We can use `#<~` to "split out" the components +of a `BVar`, splitting on `gTuple` (which uses `GHC.Generics` to +automatically figure out how to split up a product type). + +Then we apply `matVec` (our primitive `Op` that does matrix-vector +multiplication) to `w` and `x`, and then the result is that added to the +bias vector `b`. + +We can write the `runNetwork` function pretty much the same way. + +> runNetwork +> :: (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) +> => BPOp s '[ R i, Network i h1 h2 o ] (R o) +> runNetwork = withInps $ \(x :< n :< Ø) -> do +> l1 :< l2 :< l3 :< Ø <- gTuple #<~ n +> y <- runLayer -$ (x :< l1 :< Ø) +> z <- runLayer -$ (logistic y :< l2 :< Ø) +> r <- runLayer -$ (logistic z :< l3 :< Ø) +> softmax -$ (r :< Ø) +> where +> softmax :: KnownNat n => BPOp s '[ R n ] (R n) +> softmax = withInps $ \(x :< Ø) -> do +> expX <- bindVar (exp x) +> totX <- vsum ~$ (expX :< Ø) +> scale ~$ (1/totX :< expX :< Ø) + + +After splitting out the layers in the input `Network`, we run each layer +successively using our previously defined `runLayer`, giving inputs using +`-$`. We can directly apply `logistic` to `BVar`s. At the end, we run a +[softmax function][] because MNIST is a classification challenge. The softmax +is done by applying $e^x$ for every item in the input vector, and dividing +each element by the total. + +[softmax function]: https://en.wikipedia.org/wiki/Softmax_function + + +The Magic +--------- + +What did we just define? Well, with a `BPOp s rs a`, we can *run* it and +get the output: + +> runNetOnInp +> :: (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) +> => Network i h1 h2 o +> -> R i +> -> R o +> runNetOnInp n x = evalBPOp runNetwork (x ::< n ::< Ø) + +But, the magic part is that we can also get the gradient! + +> gradNet +> :: (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) +> => Network i h1 h2 o +> -> R i +> -> Network i h1 h2 o +> gradNet n x = case gradBPOp runNetwork (x ::< n ::< Ø) of +> _gradX ::< gradN ::< Ø -> gradN + +This gives the gradient of all of the parameters in the matrices and +vectors inside the `Network`, which we can use to "train"! + +Training +======== + +Now for the real work. To train a network, we can do gradient descent +based on the gradient of some type of *error function* with respect to the +network parameters. Let's use the [cross entropy][], which is popular for +classification problems. + +[cross entropy]: https://en.wikipedia.org/wiki/Cross_entropy + +> crossEntropy +> :: KnownNat n +> => R n +> -> BPOpI s '[ R n ] Double +> crossEntropy targ (r :< Ø) = negate (dot .$ (log r :< t :< Ø)) +> where +> t = constVar targ + +Given a target vector and a `BVar` referring to the result of the network, +we can directly apply: + +$$ +H(\mathbf{r}, \mathbf{t}) = - (log(\mathbf{r}) \cdot \mathbf{t}) +$$ + +Just for fun, I implemented `crossEntropy` in "implicit-graph" mode, so you +don't see any binds or returns. + +Now, a function to make one gradient descent step based on an input vector +and a target, using `gradBPOp`: + +> trainStep +> :: forall i h1 h2 o. (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) +> => Double +> -> R i +> -> R o +> -> Network i h1 h2 o +> -> Network i h1 h2 o +> trainStep r !x !t !n = case gradBPOp o (x ::< n ::< Ø) of +> _ ::< gN ::< Ø -> +> n - (realToFrac r * gN) +> where +> o :: BPOp s '[ R i, Network i h1 h2 o ] Double +> o = do +> y <- runNetwork +> implicitly (crossEntropy t) -$ (y :< Ø) + +A convenient wrapper for training over all of the observations in a list: + +> trainList +> :: (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) +> => Double +> -> [(R i, R o)] +> -> Network i h1 h2 o +> -> Network i h1 h2 o +> trainList r = flip $ foldl' (\n (x,y) -> trainStep r x y n) + +Pulling it all together +======================= + +`testNet` will be a quick way to test our net by computing the percentage +of correct guesses: (mostly using *hmatrix* stuff) + +> testNet +> :: forall i h1 h2 o. (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) +> => [(R i, R o)] +> -> Network i h1 h2 o +> -> Double +> testNet xs n = sum (map (uncurry test) xs) / fromIntegral (length xs) +> where +> test :: R i -> R o -> Double +> test x (extract->t) +> | HM.maxIndex t == HM.maxIndex (extract r) = 1 +> | otherwise = 0 +> where +> r :: R o +> r = evalBPOp runNetwork (x ::< n ::< Ø) + +And now, a main loop! + +If you are following along at home, download the [mnist data set +files][mnist] and uncompress them into the folder `data`, and everything +should work fine. + +[mnist]: http://yann.lecun.com/exdb/mnist/ + +> main :: IO () +> main = MWC.withSystemRandom $ \g -> do +> Just train <- loadMNIST "data/train-images-idx3-ubyte" "data/train-labels-idx1-ubyte" +> Just test <- loadMNIST "data/t10k-images-idx3-ubyte" "data/t10k-labels-idx1-ubyte" +> putStrLn "Loaded data." +> net0 <- MWC.uniformR @(Network 784 300 100 9) (-0.5, 0.5) g +> flip evalStateT net0 . forM_ [1..] $ \e -> do +> train' <- liftIO . fmap V.toList $ MWC.uniformShuffle (V.fromList train) g +> liftIO $ printf "[Epoch %d]\n" (e :: Int) +> +> forM_ ([1..] `zip` chunksOf batch train') $ \(b, chnk) -> StateT $ \n0 -> do +> printf "(Batch %d)\n" (b :: Int) +> +> t0 <- getCurrentTime +> n' <- evaluate . force $ trainList rate chnk n0 +> t1 <- getCurrentTime +> printf "Trained on %d points in %s.\n" batch (show (t1 `diffUTCTime` t0)) +> +> let trainScore = testNet chnk n' +> testScore = testNet test n' +> printf "Training error: %.2f%%\n" ((1 - trainScore) * 100) +> printf "Validation error: %.2f%%\n" ((1 - testScore ) * 100) +> +> return ((), n') +> where +> rate = 0.02 +> batch = 5000 + +Each iteration of the loop: + +1. Shuffles the training set +2. Splits it into chunks of `batch` size +3. Uses `trainList` to train over the batch +4. Computes the score based on `testNet` based on the training set and the + test set +5. Prints out the results + +And, that's really it! + +Result +------ + +I haven't put much into optimizing the library yet, but the network (with +hidden layer sizes 300 and 100) seems to take 25s on my computer to finish +a batch of 5000 training points. It's slow (five minutes per 60000 point +epooch), but it's a first unoptimized run and a proof of concept! It's my +goal to get this down to a point where the result has the same performance +characteristics as the actual backend (*hmatrix*), and so overhead is 0. + +Main takeaways +============== + +Most of the actual heavy lifting/logic actually came from the *hmatrix* +library itself. We just created simple types to wrap up our bare matrices. + +Basically, all that *backprop* did was give you an API to define *how to +run* a neural net --- how to *run* a net based on a `Network` and `R i` input +you were given. The goal of the library is to let you write down how to +run things in as natural way as possible. + +And then, after things are run, we can just get the gradient and roll from +there! + +Because the heavy lifting is done by the data types themselves, we can +presumably plug in *any* type and any tensor/numerical backend, and reap +the benefits of those libraries' optimizations and parallelizations. *Any* +type can be backpropagated! :D + +What now? +--------- + +Check out the docs for the [Numeric.Backprop][] module for a more detailed +picture of what's going on, or find more examples at the [github repo][repo]! + +[Numeric.Backprop]: http://hackage.haskell.org/package/backprop/docs/Numeric-Backprop.html + +Boring stuff +============ + +Here is a small wrapper function over the [mnist-idx][] library loading the +contents of the idx files into *hmatrix* vectors: + +[mnist-idx]: http://hackage.haskell.org/package/mnist-idx + +> loadMNIST +> :: FilePath +> -> FilePath +> -> IO (Maybe [(R 784, R 9)]) +> loadMNIST fpI fpL = runMaybeT $ do +> i <- MaybeT $ decodeIDXFile fpI +> l <- MaybeT $ decodeIDXLabelsFile fpL +> d <- MaybeT . return $ labeledIntData l i +> r <- MaybeT . return $ for d (bitraverse mkImage mkLabel . swap) +> liftIO . evaluate $ force r +> where +> mkImage :: VU.Vector Int -> Maybe (R 784) +> mkImage = create . VG.convert . VG.map (\i -> fromIntegral i / 255) +> mkLabel :: Int -> Maybe (R 9) +> mkLabel n = create $ HM.build 9 (\i -> if round i == n then 1 else 0) + +And here are instances to generating random +vectors/matrices/layers/networks, used for the initialization step. + +> instance KnownNat n => MWC.Variate (R n) where +> uniform g = randomVector <$> MWC.uniform g <*> pure Uniform +> uniformR (l, h) g = (\x -> x * (h - l) + l) <$> MWC.uniform g +> +> instance (KnownNat m, KnownNat n) => MWC.Variate (L m n) where +> uniform g = uniformSample <$> MWC.uniform g <*> pure 0 <*> pure 1 +> uniformR (l, h) g = (\x -> x * (h - l) + l) <$> MWC.uniform g +> +> instance (KnownNat i, KnownNat o) => MWC.Variate (Layer i o) where +> uniform g = Layer <$> MWC.uniform g <*> MWC.uniform g +> uniformR (l, h) g = (\x -> x * (h - l) + l) <$> MWC.uniform g +> +> instance (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o) => MWC.Variate (Network i h1 h2 o) where +> uniform g = Net <$> MWC.uniform g <*> MWC.uniform g <*> MWC.uniform g +> uniformR (l, h) g = (\x -> x * (h - l) + l) <$> MWC.uniform g diff --git a/samples/backprop-monotest.hs b/samples/backprop-monotest.hs new file mode 100644 --- /dev/null +++ b/samples/backprop-monotest.hs @@ -0,0 +1,18 @@ +{-# LANGUAGE GADTs #-} + +import Numeric.Backprop.Mono + +testImplicit :: BPOp s N3 Double Double +testImplicit = implicitly $ \(x :* y :* z :* ØV) -> + ((x * y) + y) * z + +testExplicit :: BPOp s N3 Double Double +testExplicit = withInps $ \(x :* y :* z :* ØV) -> do + xy <- op2 (*) ~$ (x :* y :* ØV) + xyy <- op2 (+) ~$ (xy :* y :* ØV) + op2 (*) ~$ (xyy :* z :* ØV) + +main :: IO () +main = do + print $ backprop testImplicit (2 :+ 3 :+ 4 :+ ØV) + print $ backprop testExplicit (2 :+ 3 :+ 4 :+ ØV) diff --git a/samples/backprop-neural-test.lhs b/samples/backprop-neural-test.lhs new file mode 100644 --- /dev/null +++ b/samples/backprop-neural-test.lhs @@ -0,0 +1,405 @@ +% Neural networks with backprop library +% Justin Le + +The *backprop* library performs back-propagation over a *hetereogeneous* +system of relationships. It offers both an implicit ([ad][]-like) and explicit graph +building usage style. Let's use it to build neural networks! + +[ad]: http://hackage.haskell.org/package/ad + +Repository source is [on github][repo], and so are the [rendered unstable +docs][docs]. + +[repo]: https://github.com/mstksg/backprop +[docs]: https://mstksg.github.io/backprop + +> {-# LANGUAGE DeriveGeneric #-} +> {-# LANGUAGE GADTs #-} +> {-# LANGUAGE LambdaCase #-} +> {-# LANGUAGE RankNTypes #-} +> {-# LANGUAGE ScopedTypeVariables #-} +> {-# LANGUAGE StandaloneDeriving #-} +> {-# LANGUAGE TypeApplications #-} +> {-# LANGUAGE TypeInType #-} +> {-# LANGUAGE TypeOperators #-} +> {-# LANGUAGE ViewPatterns #-} +> {-# OPTIONS_GHC -fno-warn-orphans #-} +> {-# OPTIONS_GHC -fno-warn-unused-top-binds #-} +> +> import Data.Functor +> import Data.Kind +> import Data.Maybe +> import Data.Singletons +> import Data.Singletons.Prelude +> import Data.Singletons.TypeLits +> import Data.Type.Combinator +> import Data.Type.Product +> import GHC.Generics (Generic) +> import Numeric.Backprop +> import Numeric.Backprop.Iso +> import Numeric.LinearAlgebra.Static hiding (dot) +> import System.Random.MWC +> import qualified Generics.SOP as SOP + +Ops +=== + +First, we define values of `Op` for the operations we want to do. `Op`s +are bundles of functions packaged with their hetereogeneous gradients. For +simple numeric functions, *backprop* can derive `Op`s automatically. But +for matrix operations, we have to derive them ourselves. + +The types help us with matching up the dimensions, but we still need to be +careful that our gradients are calculated correctly. + +`L` and `R` are matrix and vector types from the great *hmatrix* library. + +First, matrix-vector multiplication: + +> matVec +> :: (KnownNat m, KnownNat n) +> => Op '[ L m n, R n ] (R m) +> matVec = op2' $ \m v -> ( m #> v +> , \(fromMaybe 1 -> g) -> +> (g `outer` v, tr m #> g) +> ) + +Now, dot products: + +> dot :: KnownNat n +> => Op '[ R n, R n ] Double +> dot = op2' $ \x y -> ( x <.> y +> , \case Nothing -> (y, x) +> Just g -> (konst g * y, x * konst g) +> ) + +Polymorphic functions can be easily turned into `Op`s with `op1`/`op2` +etc., but they can also be run directly on graph nodes. + +> logistic :: Floating a => a -> a +> logistic x = 1 / (1 + exp (-x)) + +A Simple Complete Example +========================= + +At this point, we already have enough to train a simple single-hidden-layer +neural network: + +> simpleOp +> :: (KnownNat m, KnownNat n, KnownNat o) +> => R m +> -> BPOpI s '[ L n m, R n, L o n, R o ] (R o) +> simpleOp inp = \(w1 :< b1 :< w2 :< b2 :< Ø) -> +> let z = logistic $ liftB2 matVec w1 x + b1 +> in logistic $ liftB2 matVec w2 z + b2 +> where +> x = constVar inp + +Here, `simpleOp` is defined in implicit (non-monadic) style, given a tuple +of inputs and returning outputs. Now `simpleOp` can be "run" with the +input vectors and parameters (a `L n m`, `R n`, `L o n`, and `R o`) and +calculate the output of the neural net. + +> runSimple +> :: (KnownNat m, KnownNat n, KnownNat o) +> => R m +> -> Tuple '[ L n m, R n, L o n, R o ] +> -> R o +> runSimple inp = evalBPOp (implicitly $ simpleOp inp) + +Alternatively, we can define `simpleOp` in explicit monadic style, were we +specify our graph nodes explicitly. The results should be the same. + +> simpleOpExplicit +> :: (KnownNat m, KnownNat n, KnownNat o) +> => R m +> -> BPOp s '[ L n m, R n, L o n, R o ] (R o) +> simpleOpExplicit inp = withInps $ \(w1 :< b1 :< w2 :< b2 :< Ø) -> do +> -- First layer +> y1 <- matVec ~$ (w1 :< x1 :< Ø) +> let x2 = logistic (y1 + b1) +> -- Second layer +> y2 <- matVec ~$ (w2 :< x2 :< Ø) +> return $ logistic (y2 + b2) +> where +> x1 = constVar inp + +Now, for the magic of *backprop*: the library can now take advantage of +the implicit (or explicit) graph and use it to do back-propagation, too! + +> simpleGrad +> :: forall m n o. (KnownNat m, KnownNat n, KnownNat o) +> => R m +> -> R o +> -> Tuple '[ L n m, R n, L o n, R o ] +> -> Tuple '[ L n m, R n, L o n, R o ] +> simpleGrad inp targ params = gradBPOp opError params +> where +> opError :: BPOp s '[ L n m, R n, L o n, R o ] Double +> opError = do +> res <- implicitly $ simpleOp inp +> -- we explicitly bind err to prevent recomputation +> err <- bindVar $ res - t +> dot ~$ (err :< err :< Ø) +> where +> t = constVar targ + +The result is the gradient of the input tuple's components, with respect +to the `Double` result of `opError` (the squared error). We can then use +this gradient to do gradient descent. + +With Parameter Containers +========================= + +This method doesn't quite scale, because we might want to make networks +with multiple layers and parameterize networks by layers. Let's make some +basic container data types to help us organize our types, including a +recursive `Network` type that lets us chain multiple layers. + +> data Layer :: Nat -> Nat -> Type where +> Layer :: { _lWeights :: L m n +> , _lBiases :: R m +> } +> -> Layer n m +> deriving (Show, Generic) +> +> +> data Network :: Nat -> [Nat] -> Nat -> Type where +> NØ :: !(Layer a b) -> Network a '[] b +> (:&) :: !(Layer a b) -> Network b bs c -> Network a (b ': bs) c + +A `Layer n m` is a layer taking an n-vector and returning an m-vector. A +`Network a '[b, c, d] e` would be a Network that takes in an a-vector and +outputs an e-vector, with hidden layers of sizes b, c, and d. + +Isomorphisms +------------ + +The *backprop* library lets you apply operations on "parts" of data types +(like on the weights and biases of a `Layer`) by using `Iso`'s +(isomorphisms), like the ones from the *lens* library. The library doesn't +depend on lens, but it can use the `Iso`s from the library and also +custom-defined ones. + +First, we can auto-generate isomorphisms using the *generics-sop* library: + +> instance SOP.Generic (Layer n m) + +And then can create isomorphisms by hand for the two `Network` +constructors: + +> netExternal :: Iso' (Network a '[] b) (Tuple '[Layer a b]) +> netExternal = iso (\case NØ x -> x ::< Ø) +> (\case I x :< Ø -> NØ x ) +> +> netInternal :: Iso' (Network a (b ': bs) c) (Tuple '[Layer a b, Network b bs c]) +> netInternal = iso (\case x :& xs -> x ::< xs ::< Ø) +> (\case I x :< I xs :< Ø -> x :& xs ) + +An `Iso' a (Tuple as)` means that an `a` can really just be seen as a tuple +of `as`. + +Running a network +================= + +Now, we can write the `BPOp` that reprenents running the network and +getting a result. We pass in a `Sing bs` (a singleton list of the hidden +layer sizes) so that we can "pattern match" on the list and handle the +different network constructors differently. + +> netOp +> :: forall s a bs c. (KnownNat a, KnownNat c) +> => Sing bs +> -> BPOp s '[ R a, Network a bs c ] (R c) +> netOp sbs = go sbs +> where +> go :: forall d es. KnownNat d +> => Sing es +> -> BPOp s '[ R d, Network d es c ] (R c) +> go = \case +> SNil -> withInps $ \(x :< n :< Ø) -> do +> -- peek into the NØ using netExternal iso +> l :< Ø <- netExternal #<~ n +> -- run the 'layerOp' BP, with x and l as inputs +> bpOp layerOp ~$ (x :< l :< Ø) +> SNat `SCons` ses -> withInps $ \(x :< n :< Ø) -> withSingI ses $ do +> -- peek into the (:&) using the netInternal iso +> l :< n' :< Ø <- netInternal #<~ n +> -- run the 'layerOp' BP, with x and l as inputs +> z <- bpOp layerOp ~$ (x :< l :< Ø) +> -- run the 'go ses' BP, with z and n as inputs +> bpOp (go ses) ~$ (z :< n' :< Ø) +> layerOp +> :: forall d e. (KnownNat d, KnownNat e) +> => BPOp s '[ R d, Layer d e ] (R e) +> layerOp = withInps $ \(x :< l :< Ø) -> do +> -- peek into the layer using the gTuple iso, auto-generated with SOP.Generic +> w :< b :< Ø <- gTuple #<~ l +> y <- matVec ~$ (w :< x :< Ø) +> return $ logistic (y + b) + +There's some singletons work going on here, but it's fairly standard +singletons stuff. Most of the complexity here is from the static typing in +our neural network type, and *not* from *backprop*. + +From *backprop* specifically, the only elements are `#<~` lets you "split" an +input ref with the given iso, and `bpOp`, which converts a `BPOp` into an `Op` +that you can bind with `~$`. + +Note that this library doesn't support truly pattern matching on GADTs, and +that we had to pass in `Sing bs` as a reference to the structure of our +networks. + +Gradient Descent +---------------- + +Now we can do simple gradient descent. Defining an error function: + +> errOp +> :: KnownNat m +> => R m +> -> BVar s rs (R m) +> -> BPOp s rs Double +> errOp targ r = do +> err <- bindVar $ r - t +> dot ~$ (err :< err :< Ø) +> where +> t = constVar targ + +And now, we can use `backprop` to generate the gradient, and shift the +`Network`! Things are made a bit cleaner from the fact that `Network a bs c` +has a `Num` instance, so we can use `(-)` and `(*)` etc. + +> train +> :: (KnownNat a, SingI bs, KnownNat c) +> => Double +> -> R a +> -> R c +> -> Network a bs c +> -> Network a bs c +> train r x t n = case backprop (errOp t =<< netOp sing) (x ::< n ::< Ø) of +> (_, _ :< I g :< Ø) -> n - (realToFrac r * g) + +(`(::<)` is cons and `Ø` is nil for tuples.) + +Main +==== + +`main`, which will train on sample data sets, is still in progress! Right +now it just generates a random network using the *mwc-random* library and +prints each internal layer. + +> main :: IO () +> main = withSystemRandom $ \g -> do +> n <- uniform @(Network 4 '[3,2] 1) g +> void $ traverseNetwork sing (\l -> l <$ print l) n + +Appendix: Boilerplate +===================== + +And now for some typeclass instances and boilerplates unrelated to the +*backprop* library that makes our custom types easier to use. + +> instance KnownNat n => Variate (R n) where +> uniform g = randomVector <$> uniform g <*> pure Uniform +> uniformR (l, h) g = (\x -> x * (h - l) + l) <$> uniform g +> +> instance (KnownNat m, KnownNat n) => Variate (L m n) where +> uniform g = uniformSample <$> uniform g <*> pure 0 <*> pure 1 +> uniformR (l, h) g = (\x -> x * (h - l) + l) <$> uniform g +> +> instance (KnownNat n, KnownNat m) => Variate (Layer n m) where +> uniform g = subtract 1 . (* 2) <$> (Layer <$> uniform g <*> uniform g) +> uniformR (l, h) g = (\x -> x * (h - l) + l) <$> uniform g +> +> instance (KnownNat m, KnownNat n) => Num (Layer n m) where +> Layer w1 b1 + Layer w2 b2 = Layer (w1 + w2) (b1 + b2) +> Layer w1 b1 - Layer w2 b2 = Layer (w1 - w2) (b1 - b2) +> Layer w1 b1 * Layer w2 b2 = Layer (w1 * w2) (b1 * b2) +> abs (Layer w b) = Layer (abs w) (abs b) +> signum (Layer w b) = Layer (signum w) (signum b) +> negate (Layer w b) = Layer (negate w) (negate b) +> fromInteger x = Layer (fromInteger x) (fromInteger x) +> +> instance (KnownNat m, KnownNat n) => Fractional (Layer n m) where +> Layer w1 b1 / Layer w2 b2 = Layer (w1 / w2) (b1 / b2) +> recip (Layer w b) = Layer (recip w) (recip b) +> fromRational x = Layer (fromRational x) (fromRational x) +> +> instance (KnownNat a, SingI bs, KnownNat c) => Variate (Network a bs c) where +> uniform g = genNet sing (uniform g) +> uniformR (l, h) g = (\x -> x * (h - l) + l) <$> uniform g +> +> genNet +> :: forall f a bs c. (Applicative f, KnownNat a, KnownNat c) +> => Sing bs +> -> (forall d e. (KnownNat d, KnownNat e) => f (Layer d e)) +> -> f (Network a bs c) +> genNet sbs f = go sbs +> where +> go :: forall d es. KnownNat d => Sing es -> f (Network d es c) +> go = \case +> SNil -> NØ <$> f +> SNat `SCons` ses -> (:&) <$> f <*> go ses +> +> mapNetwork0 +> :: forall a bs c. (KnownNat a, KnownNat c) +> => Sing bs +> -> (forall d e. (KnownNat d, KnownNat e) => Layer d e) +> -> Network a bs c +> mapNetwork0 sbs f = getI $ genNet sbs (I f) +> +> traverseNetwork +> :: forall a bs c f. (KnownNat a, KnownNat c, Applicative f) +> => Sing bs +> -> (forall d e. (KnownNat d, KnownNat e) => Layer d e -> f (Layer d e)) +> -> Network a bs c +> -> f (Network a bs c) +> traverseNetwork sbs f = go sbs +> where +> go :: forall d es. KnownNat d => Sing es -> Network d es c -> f (Network d es c) +> go = \case +> SNil -> \case +> NØ x -> NØ <$> f x +> SNat `SCons` ses -> \case +> x :& xs -> (:&) <$> f x <*> go ses xs +> +> mapNetwork1 +> :: forall a bs c. (KnownNat a, KnownNat c) +> => Sing bs +> -> (forall d e. (KnownNat d, KnownNat e) => Layer d e -> Layer d e) +> -> Network a bs c +> -> Network a bs c +> mapNetwork1 sbs f = getI . traverseNetwork sbs (I . f) +> +> mapNetwork2 +> :: forall a bs c. (KnownNat a, KnownNat c) +> => Sing bs +> -> (forall d e. (KnownNat d, KnownNat e) => Layer d e -> Layer d e -> Layer d e) +> -> Network a bs c +> -> Network a bs c +> -> Network a bs c +> mapNetwork2 sbs f = go sbs +> where +> go :: forall d es. KnownNat d => Sing es -> Network d es c -> Network d es c -> Network d es c +> go = \case +> SNil -> \case +> NØ x -> \case +> NØ y -> NØ (f x y) +> SNat `SCons` ses -> \case +> x :& xs -> \case +> y :& ys -> f x y :& go ses xs ys +> +> instance (KnownNat a, SingI bs, KnownNat c) => Num (Network a bs c) where +> (+) = mapNetwork2 sing (+) +> (-) = mapNetwork2 sing (-) +> (*) = mapNetwork2 sing (*) +> negate = mapNetwork1 sing negate +> abs = mapNetwork1 sing abs +> signum = mapNetwork1 sing signum +> fromInteger x = mapNetwork0 sing (fromInteger x) +> +> instance (KnownNat a, SingI bs, KnownNat c) => Fractional (Network a bs c) where +> (/) = mapNetwork2 sing (/) +> recip = mapNetwork1 sing recip +> fromRational x = mapNetwork0 sing (fromRational x) diff --git a/src/Data/Type/Util.hs b/src/Data/Type/Util.hs --- a/src/Data/Type/Util.hs +++ b/src/Data/Type/Util.hs @@ -18,6 +18,7 @@ Replicate , unzipP , zipP + , tagSum , indexP , vecToProd , prodToVec' @@ -37,7 +38,7 @@ import Control.Applicative import Data.Bifunctor -import Data.Kind +-- import Data.Kind import Data.Monoid hiding (Sum) import Data.Type.Conjunction import Data.Type.Fin @@ -49,9 +50,9 @@ import Data.Type.Vector import Lens.Micro import Type.Class.Higher -import Type.Class.Known +-- import Type.Class.Known import Type.Class.Witness -import Type.Family.List +-- import Type.Family.List import Type.Family.Nat -- | @'Replicate' n a@ is a list of @a@s repeated @n@ times. @@ -185,6 +186,7 @@ ØV -> Z_ _ :* xs -> S_ (vecLength xs) +-- | Currently not used tagSum :: Prod f as -> Sum g as diff --git a/src/Numeric/Backprop.hs b/src/Numeric/Backprop.hs --- a/src/Numeric/Backprop.hs +++ b/src/Numeric/Backprop.hs @@ -83,6 +83,8 @@ , Sum(..) -- *** As sums of products , sopVar, gSplits, gSOP + -- *** As GADTs + , withGADT, BPCont(..) -- ** Combining , liftB, (.$), liftB1, liftB2, liftB3 -- * Op @@ -105,6 +107,7 @@ import Control.Monad.Reader import Control.Monad.ST import Control.Monad.State +import Data.Kind import Data.Maybe import Data.Monoid ((<>)) import Data.STRef @@ -1082,6 +1085,7 @@ fr | isTerminal = FRTerminal gOut | otherwise = FRInternal (IRConst <$> maybeToList gOut) +-- | WARNING: the gradient continuation must only be run ONCE! backpropWith :: Every Num rs => BPOp s rs a @@ -1386,15 +1390,73 @@ -> BVar s rs d liftB3 o x y z = liftB o (x :< y :< z :< Ø) - - - - - - - - +-- | For usage with 'withGADT', to handle constructors of a GADT. See +-- documentation for 'withGADT' for more information. +data BPCont :: Type -> [Type] -> Type -> Type -> Type where + BPC :: Every Num as + => Tuple as + -> (Tuple as -> a) + -> (Prod (BVar s rs) as -> BP s rs b) + -> BPCont s rs a b +-- | Special __unsafe__ combinator that lets you pattern match and work on +-- GADTs. +-- +-- @ +-- data MyGADT :: Bool -> Type where +-- A :: String -> Int -> MyGADT 'True +-- B :: Bool -> Double -> MyGADT 'False +-- +-- +-- foo :: BP s '[ MyGADT b ] a +-- foo = 'withInps' $ \\( gVar :< Ø ) -\> +-- withGADT gVar $ \\case +-- A s i -\> BPC (s ::< i ::< Ø) (\\(s' ::< i' ::< Ø) -\> A s i) $ +-- \\(sVar :< iVar) -> do +-- -- .. in this 'BP' action, sVar and iVar are 'BPVar's that +-- -- refer to the String and Int inside the A constructor in +-- -- gVar +-- B b d -\> BPC (b ::< d ::< Ø) (\\(b' ::< d' ::< Ø) -\> B b d) $ +-- \\(bVar :< dVar) -> do +-- -- .. in this 'BP' action, bVar and dVar are 'BPVar's that +-- -- refer to the Bool and DOuble inside the B constructor in +-- -- gVar +-- @ +-- +-- 'withGADT' lets to directly pattern match on the GADT, but as soon as +-- you pattern match, you must handle the results with a 'BPCont' +-- containing: +-- +-- 1. /All/ of the items inside the GADT constructor, in a 'Tuple' +-- 2. A function from a 'Tuple' of items inside the GADT constructor that +-- assembles them back into the original /same/ constructor. +-- 3. A function from a 'Prod' of 'BVar's (that contain the items inside +-- the constructor) and doing whatever you wanted to do with it, +-- inside 'BP'. +-- +-- If you don't provide all of the items inside the GADT into the 'BPC', or +-- if your "re-assembling" function doesn't properly reassemble things +-- correctly or changes some of the values, this will not work. +-- +withGADT + :: forall s rs a b. () + => BVar s rs a + -> (a -> BPCont s rs a b) + -> BP s rs b +withGADT v f = do + x <- BP (resolveVar v) + case f x of + BPC (xs :: Tuple as) g h -> do + let bp :: BPNode s rs '[a] as + bp = BPN { _bpnOut = map1 (const (FRInternal [])) xs + , _bpnRes = xs + , _bpnGradFunc = return . only_ . g + . imap1 (\ix -> every @_ @Num ix // maybe (I 1) I) + , _bpnGradCache = Nothing + } + r <- BP . liftBase $ newSTRef bp + registerVar (IRNode IZ r) v + h $ imap1 (\ix _ -> BVNode ix r) xs -- | Apply a function to the contents of an STRef, and cache the results -- using the given lens. If already calculated, simply returned the cached diff --git a/src/Numeric/Backprop/Internal.hs b/src/Numeric/Backprop/Internal.hs --- a/src/Numeric/Backprop/Internal.hs +++ b/src/Numeric/Backprop/Internal.hs @@ -141,11 +141,9 @@ -- 'Op' (or more precisely, an 'Numeric.Backprop.OpB', which is a subtype of -- 'OpM'). So, once you create your fancy 'BP' computation, you can -- transform it into an 'OpM' using 'Numeric.Backprop.bpOp'. -newtype BP s rs a = BP { bpST :: ReaderT (Tuple rs) (StateT (BPState s rs) (ST s)) a } - deriving ( Functor - , Applicative - , Monad - ) +newtype BP s rs a + = BP { bpST :: ReaderT (Tuple rs) (StateT (BPState s rs) (ST s)) a } + deriving (Functor, Applicative, Monad) -- | The basic unit of manipulation inside 'BP' (or inside an -- implicit-graph backprop function). Instead of directly working with diff --git a/src/Numeric/Backprop/Mono/Implicit.hs b/src/Numeric/Backprop/Mono/Implicit.hs --- a/src/Numeric/Backprop/Mono/Implicit.hs +++ b/src/Numeric/Backprop/Mono/Implicit.hs @@ -106,7 +106,7 @@ -- in z + x ** y -- @ -- --- >>> 'backprop' foo (2 :+ 3 :+ ØV) +-- >>> backprop foo (2 :+ 3 :+ ØV) -- (11.46, 13.73 :+ 6.12 :+ ØV) backprop :: forall n a b. (Num a, Known Nat n) @@ -125,7 +125,7 @@ -- in z + x ** y -- @ -- --- >>> 'grad' foo (2 :+ 3 :+ ØV) +-- >>> grad foo (2 :+ 3 :+ ØV) -- 13.73 :+ 6.12 :+ ØV grad :: forall n a b. (Num a, Known Nat n) @@ -144,7 +144,7 @@ -- in z + x ** y -- @ -- --- >>> 'eval' foo (2 :+ 3 :+ ØV) +-- >>> eval foo (2 :+ 3 :+ ØV) -- 11.46 eval :: forall n a b. (Num a, Known Nat n) diff --git a/src/Numeric/Backprop/Op.hs b/src/Numeric/Backprop/Op.hs --- a/src/Numeric/Backprop/Op.hs +++ b/src/Numeric/Backprop/Op.hs @@ -782,27 +782,27 @@ -- 'Numeric.Backprop.liftB2' ('.+') v1 v2 -- @ --- | Optimized version of @'op1' ('+.')@. +-- | Optimized version of @'op1' ('+')@. (+.) :: Num a => Op '[a, a] a (+.) = op2' $ \x y -> (x + y, maybe (1, 1) (\g -> (g, g))) {-# INLINE (+.) #-} --- | Optimized version of @'op1' ('-.')@. +-- | Optimized version of @'op1' ('-')@. (-.) :: Num a => Op '[a, a] a (-.) = op2' $ \x y -> (x - y, maybe (1, -1) (\g -> (g, -g))) {-# INLINE (-.) #-} --- | Optimized version of @'op1' ('*.')@. +-- | Optimized version of @'op1' ('*')@. (*.) :: Num a => Op '[a, a] a (*.) = op2' $ \x y -> (x * y, maybe (y, x) (\g -> (y*g, x*g))) {-# INLINE (*.) #-} --- | Optimized version of @'op1' ('/.')@. +-- | Optimized version of @'op1' ('/')@. (/.) :: Fractional a => Op '[a, a] a (/.) = op2' $ \x y -> (x / y, maybe (1/y, -x/(y*y)) (\g -> (g/y, -g*x/(y*y)))) {-# INLINE (/.) #-} --- | Optimized version of @'op1' ('**.')@. +-- | Optimized version of @'op1' ('**')@. (**.) :: Floating a => Op '[a, a] a (**.) = op2' $ \x y -> (x ** y, let dx = y*x**(y-1) dy = x**y*log(x) @@ -810,42 +810,42 @@ ) {-# INLINE (**.) #-} --- | Optimized version of @'op1' 'negateOp'@. +-- | Optimized version of @'op1' 'negate'@. negateOp :: Num a => Op '[a] a negateOp = op1' $ \x -> (negate x, maybe (-1) negate) {-# INLINE negateOp #-} --- | Optimized version of @'op1' 'signumOp'@. +-- | Optimized version of @'op1' 'signum'@. signumOp :: Num a => Op '[a] a signumOp = op1' $ \x -> (signum x, const 0) {-# INLINE signumOp #-} --- | Optimized version of @'op1' 'absOp'@. +-- | Optimized version of @'op1' 'abs'@. absOp :: Num a => Op '[a] a absOp = op1' $ \x -> (abs x, maybe (signum x) (* signum x)) {-# INLINE absOp #-} --- | Optimized version of @'op1' 'recipOp'@. +-- | Optimized version of @'op1' 'recip'@. recipOp :: Fractional a => Op '[a] a recipOp = op1' $ \x -> (recip x, maybe (-1/(x*x)) ((/(x*x)) . negate)) {-# INLINE recipOp #-} --- | Optimized version of @'op1' 'expOp'@. +-- | Optimized version of @'op1' 'exp'@. expOp :: Floating a => Op '[a] a expOp = op1' $ \x -> (exp x, maybe (exp x) (exp x *)) {-# INLINE expOp #-} --- | Optimized version of @'op1' 'logOp'@. +-- | Optimized version of @'op1' 'log'@. logOp :: Floating a => Op '[a] a logOp = op1' $ \x -> (log x, (/x) . fromMaybe 1) {-# INLINE logOp #-} --- | Optimized version of @'op1' 'sqrtOp'@. +-- | Optimized version of @'op1' 'sqrt'@. sqrtOp :: Floating a => Op '[a] a sqrtOp = op1' $ \x -> (sqrt x, maybe (0.5 * sqrt x) (/ (2 * sqrt x))) {-# INLINE sqrtOp #-} --- | Optimized version of @'op2' 'logBaseOp'@. +-- | Optimized version of @'op2' 'logBase'@. logBaseOp :: Floating a => Op '[a, a] a logBaseOp = op2' $ \x y -> (logBase x y, let dx = - logBase x y / (log x * x) in maybe (dx, 1/(y * log x)) @@ -853,62 +853,62 @@ ) {-# INLINE logBaseOp #-} --- | Optimized version of @'op1' 'sinOp'@. +-- | Optimized version of @'op1' 'sin'@. sinOp :: Floating a => Op '[a] a sinOp = op1' $ \x -> (sin x, maybe (cos x) (* cos x)) {-# INLINE sinOp #-} --- | Optimized version of @'op1' 'cosOp'@. +-- | Optimized version of @'op1' 'cos'@. cosOp :: Floating a => Op '[a] a cosOp = op1' $ \x -> (cos x, maybe (-sin x) (* (-sin x))) {-# INLINE cosOp #-} --- | Optimized version of @'op1' 'tanOp'@. +-- | Optimized version of @'op1' 'tan'@. tanOp :: Floating a => Op '[a] a tanOp = op1' $ \x -> (tan x, (/ cos x^(2::Int)) . fromMaybe 1) {-# INLINE tanOp #-} --- | Optimized version of @'op1' 'asinOp'@. +-- | Optimized version of @'op1' 'asin'@. asinOp :: Floating a => Op '[a] a asinOp = op1' $ \x -> (asin x, (/ sqrt(1 - x*x)) . fromMaybe 1) {-# INLINE asinOp #-} --- | Optimized version of @'op1' 'acosOp'@. +-- | Optimized version of @'op1' 'acos'@. acosOp :: Floating a => Op '[a] a acosOp = op1' $ \x -> (acos x, (/ sqrt (1 - x*x)) . maybe (-1) negate) {-# INLINE acosOp #-} --- | Optimized version of @'op1' 'atanOp'@. +-- | Optimized version of @'op1' 'atan'@. atanOp :: Floating a => Op '[a] a atanOp = op1' $ \x -> (atan x, (/ (x*x + 1)) . fromMaybe 1) {-# INLINE atanOp #-} --- | Optimized version of @'op1' 'sinhOp'@. +-- | Optimized version of @'op1' 'sinh'@. sinhOp :: Floating a => Op '[a] a sinhOp = op1' $ \x -> (sinh x, maybe (cosh x) (* cosh x)) {-# INLINE sinhOp #-} --- | Optimized version of @'op1' 'coshOp'@. +-- | Optimized version of @'op1' 'cosh'@. coshOp :: Floating a => Op '[a] a coshOp = op1' $ \x -> (cosh x, maybe (sinh x) (* sinh x)) {-# INLINE coshOp #-} --- | Optimized version of @'op1' 'tanhOp'@. +-- | Optimized version of @'op1' 'tanh'@. tanhOp :: Floating a => Op '[a] a tanhOp = op1' $ \x -> (tanh x, (/ cosh x^(2::Int)) . fromMaybe 1) {-# INLINE tanhOp #-} --- | Optimized version of @'op1' 'asinhOp'@. +-- | Optimized version of @'op1' 'asinh'@. asinhOp :: Floating a => Op '[a] a asinhOp = op1' $ \x -> (asinh x, (/ sqrt (x*x + 1)) . fromMaybe 1) {-# INLINE asinhOp #-} --- | Optimized version of @'op1' 'acoshOp'@. +-- | Optimized version of @'op1' 'acosh'@. acoshOp :: Floating a => Op '[a] a acoshOp = op1' $ \x -> (acosh x, (/ sqrt (x*x - 1)) . fromMaybe 1) {-# INLINE acoshOp #-} --- | Optimized version of @'op1' 'atanhOp'@. +-- | Optimized version of @'op1' 'atanh'@. atanhOp :: Floating a => Op '[a] a atanhOp = op1' $ \x -> (atanh x, (/ (1 - x*x)) . fromMaybe 1) {-# INLINE atanhOp #-} diff --git a/src/Numeric/Backprop/Op/Mono.hs b/src/Numeric/Backprop/Op/Mono.hs --- a/src/Numeric/Backprop/Op/Mono.hs +++ b/src/Numeric/Backprop/Op/Mono.hs @@ -576,78 +576,78 @@ negateOp :: Num a => Op N1 a a negateOp = BP.negateOp --- | Optimized version of @'op1' 'signumOp'@. +-- | Optimized version of @'op1' 'signum'@. signumOp :: Num a => Op N1 a a signumOp = BP.signumOp --- | Optimized version of @'op1' 'absOp'@. +-- | Optimized version of @'op1' 'abs'@. absOp :: Num a => Op N1 a a absOp = BP.absOp --- | Optimized version of @'op1' 'recipOp'@. +-- | Optimized version of @'op1' 'recip'@. recipOp :: Fractional a => Op N1 a a recipOp = BP.recipOp --- | Optimized version of @'op1' 'expOp'@. +-- | Optimized version of @'op1' 'exp'@. expOp :: Floating a => Op N1 a a expOp = BP.expOp --- | Optimized version of @'op1' 'logOp'@. +-- | Optimized version of @'op1' 'log'@. logOp :: Floating a => Op N1 a a logOp = BP.logOp --- | Optimized version of @'op1' 'sqrtOp'@. +-- | Optimized version of @'op1' 'sqrt'@. sqrtOp :: Floating a => Op N1 a a sqrtOp = BP.sqrtOp --- | Optimized version of @'op2' 'logBaseOp'@. +-- | Optimized version of @'op2' 'logBase'@. logBaseOp :: Floating a => Op N2 a a logBaseOp = BP.logBaseOp --- | Optimized version of @'op1' 'sinOp'@. +-- | Optimized version of @'op1' 'sin'@. sinOp :: Floating a => Op N1 a a sinOp = BP.sinOp --- | Optimized version of @'op1' 'cosOp'@. +-- | Optimized version of @'op1' 'cos'@. cosOp :: Floating a => Op N1 a a cosOp = BP.cosOp --- | Optimized version of @'op1' 'tanOp'@. +-- | Optimized version of @'op1' 'tan'@. tanOp :: Floating a => Op N1 a a tanOp = BP.tanOp --- | Optimized version of @'op1' 'asinOp'@. +-- | Optimized version of @'op1' 'asin'@. asinOp :: Floating a => Op N1 a a asinOp = BP.asinOp --- | Optimized version of @'op1' 'acosOp'@. +-- | Optimized version of @'op1' 'acos'@. acosOp :: Floating a => Op N1 a a acosOp = BP.acosOp --- | Optimized version of @'op1' 'atanOp'@. +-- | Optimized version of @'op1' 'atan'@. atanOp :: Floating a => Op N1 a a atanOp = BP.atanOp --- | Optimized version of @'op1' 'sinhOp'@. +-- | Optimized version of @'op1' 'sinh'@. sinhOp :: Floating a => Op N1 a a sinhOp = BP.sinhOp --- | Optimized version of @'op1' 'coshOp'@. +-- | Optimized version of @'op1' 'cosh'@. coshOp :: Floating a => Op N1 a a coshOp = BP.coshOp --- | Optimized version of @'op1' 'tanhOp'@. +-- | Optimized version of @'op1' 'tanh'@. tanhOp :: Floating a => Op N1 a a tanhOp = BP.tanhOp --- | Optimized version of @'op1' 'asinhOp'@. +-- | Optimized version of @'op1' 'asinh'@. asinhOp :: Floating a => Op N1 a a asinhOp = BP.asinhOp --- | Optimized version of @'op1' 'acoshOp'@. +-- | Optimized version of @'op1' 'acosh'@. acoshOp :: Floating a => Op N1 a a acoshOp = BP.acoshOp --- | Optimized version of @'op1' 'atanhOp'@. +-- | Optimized version of @'op1' 'atanh'@. atanhOp :: Floating a => Op N1 a a atanhOp = BP.atanhOp