diff --git a/.ghci b/.ghci
new file mode 100644
--- /dev/null
+++ b/.ghci
@@ -0,0 +1,1 @@
+:set -isrc -idist/build/autogen -optP-include -optPdist/build/autogen/cabal_macros.h
diff --git a/.gitignore b/.gitignore
new file mode 100644
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,13 @@
+dist
+docs
+wiki
+TAGS
+tags
+wip
+.DS_Store
+.*.swp
+.*.swo
+*.o
+*.hi
+*~
+*#
diff --git a/.travis.yml b/.travis.yml
new file mode 100644
--- /dev/null
+++ b/.travis.yml
@@ -0,0 +1,25 @@
+language: haskell
+before_install:
+  # Uncomment whenever hackage is down.
+  # - mkdir -p ~/.cabal && cp travis/config ~/.cabal/config && cabal update
+
+  # Try installing some of the build-deps with apt-get for speed.
+  - travis/cabal-apt-install $mode
+
+install:
+  - cabal configure $mode
+  - cabal build
+
+script:
+  - $script && hlint src --cpp-define HLINT
+
+notifications:
+  irc:
+    channels:
+      - "irc.freenode.org#haskell-lens"
+    skip_join: true
+    template:
+      - "\x0313foo\x03/\x0306%{branch}\x03 \x0314%{commit}\x03 %{build_url} %{message}"
+
+env:
+  - mode="--enable-tests" script="cabal test"
diff --git a/.vim.custom b/.vim.custom
new file mode 100644
--- /dev/null
+++ b/.vim.custom
@@ -0,0 +1,31 @@
+" Add the following to your .vimrc to automatically load this on startup
+
+" if filereadable(".vim.custom")
+"     so .vim.custom
+" endif
+
+function StripTrailingWhitespace()
+  let myline=line(".")
+  let mycolumn = col(".")
+  silent %s/  *$//
+  call cursor(myline, mycolumn)
+endfunction
+
+" enable syntax highlighting
+syntax on
+
+" search for the tags file anywhere between here and /
+set tags=TAGS;/
+
+" highlight tabs and trailing spaces
+set listchars=tab:‗‗,trail:‗
+set list
+
+" f2 runs hasktags
+map <F2> :exec ":!hasktags -x -c --ignore src"<CR><CR>
+
+" strip trailing whitespace before saving
+" au BufWritePre *.hs,*.markdown silent! cal StripTrailingWhitespace()
+
+" rebuild hasktags after saving
+au BufWritePost *.hs silent! :exec ":!hasktags -x -c --ignore src"
diff --git a/CHANGELOG.markdown b/CHANGELOG.markdown
new file mode 100644
--- /dev/null
+++ b/CHANGELOG.markdown
@@ -0,0 +1,3 @@
+0.1
+---
+* Repository initialized
diff --git a/HLint.hs b/HLint.hs
new file mode 100644
--- /dev/null
+++ b/HLint.hs
@@ -0,0 +1,1 @@
+import "hint" HLint.Default
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright 2011 Edward Kmett
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+
+1. Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+
+3. Neither the name of the author nor the names of his contributors
+   may be used to endorse or promote products derived from this software
+   without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR
+IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED.  IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR
+ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
+STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
+ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGE.
diff --git a/README.markdown b/README.markdown
new file mode 100644
--- /dev/null
+++ b/README.markdown
@@ -0,0 +1,27 @@
+optimization
+===
+
+These are a set of implementations of various numerical optimization
+methods in Haskell. Note that these implementations were originally
+written as part of a class project; while at one point they worked
+no attention has been given to numerical stability or the many other
+potential difficulties of implementing robust numerical
+methods. That being said, they should serve to succinctly illustrate
+a number of optimization techniques from the modern optimization
+literature.
+
+Those seeking a high-level overview of some of these methods are
+referred to Stephen Wright's excellent
+[tutorial](http://videolectures.net/nips2010_wright_oaml/) from NIPS
+2010. A deeper introduction can be found in Boyd and Vandenberghe's
+*Complex Optimization* which available freely online.
+
+
+Contact Information
+-------------------
+
+Contributions and bug reports are welcome!
+
+Please feel free to contact me through github or on the #haskell IRC channel on irc.freenode.net.
+
+- Ben Gamari
diff --git a/Setup.lhs b/Setup.lhs
new file mode 100644
--- /dev/null
+++ b/Setup.lhs
@@ -0,0 +1,44 @@
+#!/usr/bin/runhaskell
+\begin{code}
+{-# OPTIONS_GHC -Wall #-}
+module Main (main) where
+
+import Data.List ( nub )
+import Data.Version ( showVersion )
+import Distribution.Package ( PackageName(PackageName), PackageId, InstalledPackageId, packageVersion, packageName )
+import Distribution.PackageDescription ( PackageDescription(), TestSuite(..) )
+import Distribution.Simple ( defaultMainWithHooks, UserHooks(..), simpleUserHooks )
+import Distribution.Simple.Utils ( rewriteFile, createDirectoryIfMissingVerbose )
+import Distribution.Simple.BuildPaths ( autogenModulesDir )
+import Distribution.Simple.Setup ( BuildFlags(buildVerbosity), fromFlag )
+import Distribution.Simple.LocalBuildInfo ( withLibLBI, withTestLBI, LocalBuildInfo(), ComponentLocalBuildInfo(componentPackageDeps) )
+import Distribution.Verbosity ( Verbosity )
+import System.FilePath ( (</>) )
+
+main :: IO ()
+main = defaultMainWithHooks simpleUserHooks
+  { buildHook = \pkg lbi hooks flags -> do
+     generateBuildModule (fromFlag (buildVerbosity flags)) pkg lbi
+     buildHook simpleUserHooks pkg lbi hooks flags
+  }
+
+generateBuildModule :: Verbosity -> PackageDescription -> LocalBuildInfo -> IO ()
+generateBuildModule verbosity pkg lbi = do
+  let dir = autogenModulesDir lbi
+  createDirectoryIfMissingVerbose verbosity True dir
+  withLibLBI pkg lbi $ \_ libcfg -> do
+    withTestLBI pkg lbi $ \suite suitecfg -> do
+      rewriteFile (dir </> "Build_" ++ testName suite ++ ".hs") $ unlines
+        [ "module Build_" ++ testName suite ++ " where"
+        , "deps :: [String]"
+        , "deps = " ++ (show $ formatdeps (testDeps libcfg suitecfg))
+        ]
+  where
+    formatdeps = map (formatone . snd)
+    formatone p = case packageName p of
+      PackageName n -> n ++ "-" ++ showVersion (packageVersion p)
+
+testDeps :: ComponentLocalBuildInfo -> ComponentLocalBuildInfo -> [(InstalledPackageId, PackageId)]
+testDeps xs ys = nub $ componentPackageDeps xs ++ componentPackageDeps ys
+
+\end{code}
diff --git a/optimization.cabal b/optimization.cabal
new file mode 100644
--- /dev/null
+++ b/optimization.cabal
@@ -0,0 +1,85 @@
+name:          optimization
+category:      Math
+version:       0.1
+license:       BSD3
+cabal-version: >= 1.10
+license-file:  LICENSE
+author:        Ben Gamari
+maintainer:    Ben Gamari <bgamari@gmail.com>
+stability:     experimental
+homepage:      http://github.com/bgamari/optimization
+bug-reports:   http://github.com/bgamari/optimization/issues
+copyright:     Copyright (C) 2013 Ben Gamari
+synopsis:      Numerical optimization
+description:
+  These are a set of implementations of various numerical optimization
+  methods in Haskell. Note that these implementations were originally
+  written as part of a class project; while at one point they worked
+  no attention has been given to numerical stability or the many other
+  potential difficulties of implementing robust numerical
+  methods. That being said, they should serve to succinctly illustrate
+  a number of optimization techniques from the modern optimization
+  literature.
+  .
+  Those seeking a high-level overview of some of these methods are
+  referred to Stephen Wright's excellent tutorial from NIPS 2010
+  <http://videolectures.net/nips2010_wright_oaml/>. A deeper
+  introduction can be found in Boyd and Vandenberghe's *Complex
+  Optimization* which available freely online.
+
+build-type:    Custom
+
+extra-source-files:
+  .ghci
+  .gitignore
+  .travis.yml
+  .vim.custom
+  CHANGELOG.markdown
+  HLint.hs
+  README.markdown
+  travis/cabal-apt-install
+  travis/config
+
+source-repository head
+  type: git
+  location: git://github.com/bgamari/optimization.git
+
+library
+  hs-source-dirs: src
+  default-language: Haskell2010
+  ghc-options: -Wall -fno-warn-type-defaults
+  build-depends:
+    base                >= 4.4          && < 5,
+    vector              >= 0.10         && < 1.0,
+    ad                  >= 3.4          && < 4.0,
+    linear              >= 1.0          && < 2.0,
+    semigroupoids       >= 3.0          && < 4.0,
+    distributive        >= 0.3          && < 0.4
+
+  exposed-modules:
+    Optimization.LineSearch
+    Optimization.LineSearch.ConjugateGradient
+    Optimization.LineSearch.BarzilaiBorwein
+    Optimization.LineSearch.SteepestDescent
+    Optimization.LineSearch.MirrorDescent
+    Optimization.LineSearch.BFGS
+    Optimization.TrustRegion.Nesterov1983
+    Optimization.TrustRegion.Fista
+    Optimization.TrustRegion.Newton
+    Optimization.Constrained.Penalty
+    Optimization.Constrained.ProjectedSubgradient
+
+
+test-suite doctests
+  type:    exitcode-stdio-1.0
+  main-is: doctests.hs
+  default-language: Haskell2010
+  build-depends:
+    base,
+    directory >= 1.0,
+    doctest >= 0.9.1,
+    filepath
+  ghc-options: -Wall -threaded
+  if impl(ghc<7.6.1)
+    ghc-options: -Werror
+  hs-source-dirs: tests
diff --git a/src/Optimization/Constrained/Penalty.hs b/src/Optimization/Constrained/Penalty.hs
new file mode 100644
--- /dev/null
+++ b/src/Optimization/Constrained/Penalty.hs
@@ -0,0 +1,88 @@
+{-# LANGUAGE DeriveFunctor, DeriveFoldable, DeriveTraversable, DeriveGeneric,
+    FlexibleInstances, FlexibleContexts, TypeFamilies,
+    KindSignatures, DataKinds, TypeOperators, RankNTypes, ExistentialQuantification #-}
+
+module Optimization.Constrained.Penalty
+  ( -- * Building the problem
+    Opt
+  , FU(..)
+  , optimize
+  , constrainEQ
+  , constrainLT
+  , constrainGT
+    -- * Optimizing the problem
+  , minimize
+  , maximize
+    -- * Finding the Lagrangian
+  , lagrangian
+  ) where
+
+import           Numeric.AD.Types
+
+import qualified Data.Vector as V
+
+newtype FU f a = FU { runFU :: forall s. Mode s => f (AD s a) -> AD s a }
+
+type V = V.Vector
+
+-- | @Opt d f gs hs@ is a Lagrangian optimization problem with objective @f@
+-- equality (@g(x) == 0@) constraints @gs@ and less-than (@h(x) < 0@)
+-- constraints @hs@
+data Opt f a = Opt (FU f a) (V (FU f a)) (V (FU f a))
+
+optimize :: (forall s. Mode s => f (AD s a) -> AD s a) -> Opt f a
+optimize f = Opt (FU f) V.empty V.empty
+
+augment :: a -> V a -> V a
+augment a xs = V.cons a xs
+
+constrainEQ :: (forall s. Mode s => f (AD s a) -> AD s a)
+            -> Opt f a -> Opt f a
+constrainEQ g (Opt f gs hs) = Opt f (augment (FU g) gs) hs
+
+constrainLT :: (forall s. Mode s => f (AD s a) -> AD s a)
+            -> Opt f a -> Opt f a
+constrainLT h (Opt f gs hs) = Opt f gs (augment (FU h) hs)
+
+constrainGT :: (Num a) => (forall s. Mode s => f (AD s a) -> AD s a)
+            -> Opt f a -> Opt f a
+constrainGT h (Opt f gs hs) = Opt f gs (augment (FU $ negate . h) hs)
+
+-- | Minimize the given constrained optimization problem
+-- This is a basic penalty method approach
+minimize :: (Functor f, Num a, Ord a, g ~ V)
+         => (FU f a -> f a -> [f a])   -- ^ Primal minimizer
+         -> Opt f a                    -- ^ The optimization problem of interest
+         -> a                          -- ^ The penalty increase factor
+         -> f a                        -- ^ The primal starting value
+         -> g a                        -- ^ The dual starting value
+         -> [f a]                      -- ^ Optimizing iterates
+minimize minX opt alpha = go
+  where go x0 l0 = let l1 = fmap (*alpha) l0
+                       x1 = head $ drop 100 $ minX (FU $ \x -> augLagrangian opt x (fmap auto l1)) x0
+                   in x1 : go x1 l1
+
+-- | Maximize the given constrained optimization problem
+maximize :: (Functor f, Num a, Ord a, g ~ V)
+         => (FU f a -> f a -> [f a])   -- ^ Primal minimizer
+         -> Opt f a                    -- ^ The optimization problem of interest
+         -> a                          -- ^ The penalty increase factor
+         -> f a                        -- ^ The primal starting value
+         -> g a                        -- ^ The dual starting value
+         -> [f a]                      -- ^ Optimizing iterates
+maximize minX (Opt (FU f) gs hs) alpha =
+    minimize minX (Opt (FU $ negate . f) gs hs) alpha
+
+-- | The Lagrangian for the given constrained optimization
+lagrangian :: (Num a) => Opt f a
+           -> (forall s. Mode s => f (AD s a) -> V (AD s a) -> AD s a)
+lagrangian (Opt (FU f) gs hs) x l =
+    f x - V.sum (V.zipWith (\lamb (FU g)->lamb * g x) l gs)
+
+-- | The augmented Lagrangian for the given constrained optimization
+augLagrangian :: (Num a, Ord a) => Opt f a
+           -> (forall s. Mode s => f (AD s a) -> V (AD s a) -> AD s a)
+augLagrangian (Opt (FU f) gs hs) x l =
+    f x + V.sum (V.zipWith (*) l $ V.concat [gs', hs'])
+  where gs' = V.map (\(FU g) -> (g x)^2) gs
+        hs' = V.map (\(FU h) -> (max 0 $ h x)^2) hs
diff --git a/src/Optimization/Constrained/ProjectedSubgradient.hs b/src/Optimization/Constrained/ProjectedSubgradient.hs
new file mode 100644
--- /dev/null
+++ b/src/Optimization/Constrained/ProjectedSubgradient.hs
@@ -0,0 +1,114 @@
+module Optimization.Constrained.ProjectedSubgradient
+    ( -- * Projected subgradient method
+      projSubgrad
+    , linearProjSubgrad
+      -- * Step schedules
+    , StepSched
+    , optimalStepSched
+    , constStepSched
+    , sqrtKStepSched
+    , invKStepSched
+      -- * Linear constraints
+    , Constraint(..)
+    , linearProjection
+    ) where
+
+import Linear
+import Data.Traversable
+import Data.Function (on)
+import Data.List (maximumBy)
+
+-- | A step size schedule
+-- A list of functions (one per step) which, given a function's
+-- gradient and value, will provide a size for the next step
+type StepSched f a = [f a -> a -> a]
+
+-- | @projSubgrad stepSizes proj a b x0@ minimizes the objective @A
+-- x - b@ potentially projecting iterates into a feasible space with
+-- @proj@ with the given step-size schedule
+projSubgrad :: (Additive f, Traversable f, Metric f, Ord a, Fractional a)
+            => StepSched f a  -- ^ A step size schedule
+            -> (f a -> f a)   -- ^ Function projecting into the feasible space
+            -> (f a -> f a)   -- ^ Gradient of objective function
+            -> (f a -> a)     -- ^ The objective function
+            -> f a            -- ^ Initial solution
+            -> [f a]
+projSubgrad stepSizes proj df f = go stepSizes
+  where go (alpha:rest) x0 =
+            let p = negated $ df x0
+                step = alpha (df x0) (f x0)
+                x1 = proj $ x0 ^+^ step *^ p
+            in x1 : go rest x1
+        go [] _ = []
+
+-- | @linearProjSubgrad stepSizes proj a b x0@ minimizes the objective @A
+-- x - b@ potentially projecting iterates into a feasible space with
+-- @proj@ with the given step-size schedule
+linearProjSubgrad :: (Additive f, Traversable f, Metric f, Ord a, Fractional a)
+                  => StepSched f a  -- ^ A step size schedule
+                  -> (f a -> f a)   -- ^ Function projecting into the feasible space
+                  -> f a            -- ^ Coefficient vector @A@ of objective
+                  -> a              -- ^ Constant @b@ of objective
+                  -> f a            -- ^ Initial solution
+                  -> [f a]
+linearProjSubgrad stepSizes proj a b = go stepSizes
+  where go (alpha:rest) x0 =
+            let p = negated $ df x0
+                step = alpha a (f x0)
+                x1 = proj $ x0 ^+^ step *^ p
+            in x1 : go rest x1
+        go [] _ = []
+        df _ = a
+        f x = a `dot` x - b
+
+-- | The optimal step size schedule when the optimal value of the
+-- objective is known
+optimalStepSched :: (Fractional a, Metric f)
+                 => a    -- ^ The optimal value of the objective
+                 -> StepSched f a
+optimalStepSched fStar =
+    repeat $ \gk fk->(fk - fStar) / quadrance gk
+
+-- | Constant step size schedule
+constStepSched :: a    -- ^ The step size
+               -> StepSched f a
+constStepSched gamma =
+    repeat $ \_ _ -> gamma
+
+-- | A non-summable step size schedule
+sqrtKStepSched :: Floating a
+               => a       -- ^ The size of the first step
+               -> StepSched f a
+sqrtKStepSched gamma =
+    map (\k _ _ -> gamma / sqrt (fromIntegral k)) [0..]
+
+-- | A square-summable step size schedule
+invKStepSched :: Fractional a
+              => a        -- ^ The size of the first step
+              -> StepSched f a
+invKStepSched gamma =
+    map (\k _ _ -> gamma / fromIntegral k) [0..]
+
+-- | A linear constraint. For instance, @Constr LT 2 (V2 1 3)@ defines
+-- the constraint @x_1 + 3 x_2 <= 2@
+data Constraint f a = Constr Ordering a (f a)
+                    deriving (Show)
+
+-- | Project onto a the space of feasible solutions defined by a set
+-- of linear constraints
+linearProjection :: (Fractional a, Ord a, RealFloat a, Metric f)
+                 => [Constraint f a] -- ^ A set of linear constraints
+                 -> f a -> f a
+linearProjection constraints x =
+    case unmet of
+      []   -> x
+      _    -> linearProjection constraints $ fixConstraint x
+              $ maximumBy (flip compare `on` (`ap` x)) unmet
+  where unmet = filter (not . met x) constraints
+        ap (Constr _ b a) c = a `dot` c - b
+        met c (Constr t a constr) = let y = constr `dot` c - a
+                                    in case t of
+                                       EQ -> abs y < 1e-4
+                                       GT -> y >= 0 || abs y < 1e-4
+                                       LT -> y <= 0 || abs y < 1e-4
+        fixConstraint c (Constr _ b a) = c ^-^ (a `dot` c - b) *^ a ^/ quadrance a
diff --git a/src/Optimization/LineSearch.hs b/src/Optimization/LineSearch.hs
new file mode 100644
--- /dev/null
+++ b/src/Optimization/LineSearch.hs
@@ -0,0 +1,100 @@
+-- |
+-- Module      : Optimization.LineSearch
+-- Copyright   : (c) 2012-2013 Ben Gamari
+-- License     : BSD-style (see the file LICENSE)
+-- Maintainer  : Ben Gamari <bgamari@gmail.com>
+-- Stability   : provisional
+-- Portability : portable
+--
+-- Line search algorithms are a class of iterative optimization
+-- methods. These methods are distinguished by the characteristic of,
+-- starting from a point @x0@, choosing a direction @d@ (by some method)
+-- to advance and then finding an optimal distance @a@ (known as the
+-- step-size) to advance in this direction.
+--
+-- Here we provide several methods for determining this optimal
+-- distance. These can be used with any of line-search optimization
+-- algorithms found in this namespace.
+
+module Optimization.LineSearch
+    ( -- * Line search methods
+      LineSearch
+    , backtrackingSearch
+    , armijoSearch
+    , wolfeSearch
+    , newtonSearch
+    , secantSearch
+    , constantSearch
+    ) where
+
+import Prelude hiding (pred)
+import Linear
+
+-- | A 'LineSearch' method 'search df p x' determines a step size
+-- in direction 'p' from point 'x' for function 'f' with gradient 'df'
+type LineSearch f a = (f a -> f a) -> f a -> f a -> a
+
+-- | Armijo condition
+--
+-- The Armijo condition captures the intuition that step should
+-- move far enough from its starting point to change the function enough,
+-- as predicted by its gradient. This often finds its place as a criterion
+-- for line-search
+armijo :: (Num a, Additive f, Ord a, Metric f)
+       => a -> (f a -> a) -> (f a -> f a) -> f a -> f a -> a -> Bool
+armijo c1 f df x p a =
+    f (x ^+^ a *^ p) <= f x + c1 * a * (df x `dot` p)
+
+-- | Curvature condition
+curvature :: (Num a, Ord a, Additive f, Metric f)
+          => a -> (f a -> f a) -> f a -> f a -> a -> Bool
+curvature c2 df x p a =
+    df (x ^+^ a *^ p) `dot` p >= c2 * (df x `dot` p)
+
+-- | Backtracking line search algorithm
+--
+-- @backtrackingSearch gamma alpha pred@ starts with the given step
+-- size @alpha@ and reduces it by a factor of @gamma@ until the given
+-- condition is satisfied.
+backtrackingSearch :: (Num a, Ord a, Metric f)
+                   => a -> a -> (a -> Bool) -> LineSearch f a
+backtrackingSearch gamma alpha pred _ _ _ =
+    head $ dropWhile (not . pred) $ nonzero $ iterate (*gamma) alpha
+  where nonzero (x:xs) | not $ x > 0 = error "Backtracking search failed: alpha=0" -- FIXME
+                       | otherwise   = x : nonzero xs
+        nonzero [] = error "Backtracking search failed: no more iterates"
+
+-- | Armijo backtracking line search algorithm
+--
+-- @armijoSearch gamma alpha c1@ starts with the given step size @alpha@
+-- and reduces it by a factor of @gamma@ until the Armijo condition
+-- is satisfied.
+armijoSearch :: (Num a, Ord a, Metric f)
+             => a -> a -> a -> (f a -> a) -> LineSearch f a
+armijoSearch gamma alpha c1 f df p x =
+    backtrackingSearch gamma alpha (armijo c1 f df x p) df p x
+
+-- | Wolfe backtracking line search algorithm
+--
+-- @wolfeSearch gamma alpha c1@ starts with the given step size @alpha@
+-- and reduces it by a factor of @gamma@ until both the Armijo and
+-- curvature conditions is satisfied.
+wolfeSearch :: (Num a, Ord a, Metric f)
+             => a -> a -> a -> a -> (f a -> a) -> LineSearch f a
+wolfeSearch gamma alpha c1 c2 f df p x =
+    backtrackingSearch gamma alpha wolfe df p x
+  where wolfe a = armijo c1 f df p x a && curvature c2 df x p a
+
+-- | Line search by Newton's method
+newtonSearch :: (Num a) => LineSearch f a
+newtonSearch = undefined
+
+-- | Line search by secant method with given tolerance
+secantSearch :: (Num a, Fractional a) => a -> LineSearch f a
+secantSearch = undefined
+
+-- | Constant line search
+--
+-- @constantSearch c@ always chooses a step-size @c@.
+constantSearch :: a -> LineSearch f a
+constantSearch c _ _ _ = c
diff --git a/src/Optimization/LineSearch/BFGS.hs b/src/Optimization/LineSearch/BFGS.hs
new file mode 100644
--- /dev/null
+++ b/src/Optimization/LineSearch/BFGS.hs
@@ -0,0 +1,31 @@
+{-# LANGUAGE ScopedTypeVariables #-}
+
+module Optimization.LineSearch.BFGS (bfgs) where
+
+import Linear
+import Optimization.LineSearch
+import Control.Applicative
+import Data.Traversable
+import Data.Distributive
+import Data.Foldable
+
+-- | Broyden–Fletcher–Goldfarb–Shanno (BFGS) method
+-- @bfgs search df b0 x0@ where @b0@ is the inverse Hessian (the
+-- identity is often a good initial value).
+bfgs :: ( Metric f, Additive f, Distributive f, Foldable f, Traversable f, Applicative f
+        , Fractional a, Epsilon a)
+     => LineSearch f a -> (f a -> f a) -> f (f a) -> f a -> [f a]
+bfgs search df = go
+    where go b0 x0 = let p1 = negated $ b0 !* df x0
+                         alpha = search df p1 x0
+                         s = alpha *^ p1
+                         x1 = x0 ^+^ s
+                         y = df x1 ^-^ df x0
+                         -- Sherman-Morrison update of inverse Hessian
+                         sy = s `dot` y
+                         rho = if nearZero sy then 1000 else 1 / sy
+                         i = kronecker (pure 1)
+                         u = i !-! rho *!! outer y s
+                         v = i !-! rho *!! outer s y
+                         b1 = u !*! b0 !*! v !+! rho *!! outer s s
+                     in x1 : go b1 x1
diff --git a/src/Optimization/LineSearch/BarzilaiBorwein.hs b/src/Optimization/LineSearch/BarzilaiBorwein.hs
new file mode 100644
--- /dev/null
+++ b/src/Optimization/LineSearch/BarzilaiBorwein.hs
@@ -0,0 +1,17 @@
+module Optimization.LineSearch.BarzilaiBorwein
+    ( barzilaiBorwein
+    ) where
+
+import Linear
+
+-- | Barzilai-Borwein 1988 is a non-monotonic optimization method
+barzilaiBorwein :: (Additive f, Metric f, Functor f, Fractional a, Epsilon a)
+                => (f a -> f a) -> f a -> f a -> [f a]
+barzilaiBorwein df = go
+  where go x0 x1 = let s = x1 ^-^ x0
+                       z = df x1 ^-^ df x0
+                       alpha = (s `dot` z) / (z `dot` z)
+                       x2 = x1 ^-^ alpha *^ df x1
+                   in if nearZero (z `dot` z)
+                        then [x2]
+                        else x2 : go x1 x2
diff --git a/src/Optimization/LineSearch/ConjugateGradient.hs b/src/Optimization/LineSearch/ConjugateGradient.hs
new file mode 100644
--- /dev/null
+++ b/src/Optimization/LineSearch/ConjugateGradient.hs
@@ -0,0 +1,54 @@
+module Optimization.LineSearch.ConjugateGradient
+    ( -- * Conjugate gradient methods
+      conjGrad
+      -- * General line search
+    , module Optimization.LineSearch
+      -- * Beta expressions
+    , Beta
+    , fletcherReeves
+    , polakRibiere
+    , hestenesStiefel
+    ) where
+
+import Optimization.LineSearch
+import Linear
+
+-- | A beta expression 'beta df0 df1 p' is an expression for the
+-- conjugate direction contribution given the derivative 'df0' and
+-- direction 'p' for iteration 'k', 'df1' for iteration 'k+1'
+type Beta f a = f a -> f a -> f a -> a
+
+-- | Conjugate gradient method with given beta and line search method
+--
+-- The conjugate gradient method avoids the trouble encountered by the
+-- steepest descent method on poorly conditioned problems (e.g. those with
+-- a wide range of eigenvalues). It does this by choosing directions which
+-- satisfy a condition of @A@ orthogonality, ensuring that steps in the
+-- "unstretched" search space are orthogonal.
+-- TODO: clarify explanation
+{-# INLINEABLE conjGrad #-}
+conjGrad :: (Num a, RealFloat a, Additive f, Metric f)
+         => LineSearch f a -> Beta f a
+         -> (f a -> f a) -> f a -> [f a]
+conjGrad search beta df x0 = go (negated $ df x0) x0
+  where go p x = let a = search df p x
+                     x' = x ^+^ a *^ p
+                     b = beta (df x) (df x') p
+                     p' = negated (df x') ^+^ b *^ p
+                 in x' : go p' x'
+
+-- | Fletcher-Reeves expression for beta
+{-# INLINEABLE fletcherReeves #-}
+fletcherReeves :: (Num a, RealFloat a, Metric f) => Beta f a
+fletcherReeves df0 df1 _ = norm df1 / norm df0
+
+-- | Polak-Ribiere expression for beta
+{-# INLINEABLE polakRibiere #-}
+polakRibiere :: (Num a, RealFloat a, Metric f) => Beta f a
+polakRibiere df0 df1 _ = df1 `dot` (df1 ^-^ df0) / norm df0
+
+-- | Hestenes-Stiefel expression for beta
+{-# INLINEABLE hestenesStiefel #-}
+hestenesStiefel :: (Num a, RealFloat a, Metric f) => Beta f a
+hestenesStiefel df0 df1 p0 =
+    - (df1 `dot` (df1 ^-^ df0)) / (p0 `dot` (df1 ^-^ df0))
diff --git a/src/Optimization/LineSearch/MirrorDescent.hs b/src/Optimization/LineSearch/MirrorDescent.hs
new file mode 100644
--- /dev/null
+++ b/src/Optimization/LineSearch/MirrorDescent.hs
@@ -0,0 +1,23 @@
+module Optimization.LineSearch.MirrorDescent
+    ( mirrorDescent ) where
+
+import Optimization.LineSearch
+import Linear
+
+-- | Mirror descent method.
+--
+-- Originally described by Nemirovsky and Yudin and later elucidated
+-- by Beck and Teboulle, the mirror descent method is a generalization of
+-- the projected subgradient method for convex optimization.
+-- Mirror descent requires the gradient of a strongly
+-- convex function @psi@ (and its dual) as well as a way to get a
+-- subgradient for each point of the objective function @f@.
+mirrorDescent :: (Num a, Additive f)
+              => LineSearch f a -> (f a -> f a) -> (f a -> f a)
+              -> (f a -> f a) -> f a -> [f a]
+mirrorDescent search dPsi dPsiStar df = go
+  where go y0 = let x0 = dPsiStar y0
+                    t0 = search df (df x0) x0
+                    y1 = dPsi x0 ^-^ t0 *^ df x0
+                    x1 = dPsiStar y1
+                in x1 : go y1
diff --git a/src/Optimization/LineSearch/SteepestDescent.hs b/src/Optimization/LineSearch/SteepestDescent.hs
new file mode 100644
--- /dev/null
+++ b/src/Optimization/LineSearch/SteepestDescent.hs
@@ -0,0 +1,22 @@
+module Optimization.LineSearch.SteepestDescent
+    ( -- * Steepest descent method
+      steepestDescent
+    ) where
+
+import Optimization.LineSearch
+import Linear
+
+-- | Steepest descent method
+--
+-- @steepestDescent search f df x0@ optimizes a function @f@ with gradient @df@
+-- with step size schedule @search@ starting from initial point @x0@
+--
+-- The steepest descent method chooses the negative gradient of the function
+-- as its step direction.
+{-# INLINEABLE steepestDescent #-}
+steepestDescent :: (Num a, Ord a, Additive f, Metric f)
+                => LineSearch f a -> (f a -> f a) -> f a -> [f a]
+steepestDescent search df x0 = iterate go x0
+  where go x = let p = negated (df x)
+                   a = search df p x
+               in x ^+^ a *^ p
diff --git a/src/Optimization/TrustRegion/Fista.hs b/src/Optimization/TrustRegion/Fista.hs
new file mode 100644
--- /dev/null
+++ b/src/Optimization/TrustRegion/Fista.hs
@@ -0,0 +1,17 @@
+module Optimization.TrustRegion.Fista
+    ( -- * Fast Iterative Shrinkage-Thresholding Algorithm
+      fista
+    ) where
+
+import Linear
+
+-- | Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) with
+-- constant stepsize
+{-# INLINEABLE fista #-}
+fista :: (Additive f, Fractional a, Floating a)
+      => a -> (f a -> f a) -> f a -> [f a]
+fista l df x0' = go x0' x0' 1
+  where go x0 y1 t1 = let x1 = y1 ^-^ df y1 ^/ l
+                          t2 = (1 + sqrt (1 + 4 * t1^2)) / 2
+                          y2 = x1 ^+^ (t1-1) / t2 *^ (x1 ^-^ x0)
+                      in x1 : go x1 y2 t2
diff --git a/src/Optimization/TrustRegion/Nesterov1983.hs b/src/Optimization/TrustRegion/Nesterov1983.hs
new file mode 100644
--- /dev/null
+++ b/src/Optimization/TrustRegion/Nesterov1983.hs
@@ -0,0 +1,35 @@
+module Optimization.TrustRegion.Nesterov1983
+    ( -- * Nesterov's Optimal Gradient method
+      optimalGradient
+    ) where
+
+import Linear
+
+-- | Nesterov 1983
+-- @optimalGradient kappa l df alpha0 x0@ is Nesterov's optimal
+-- gradient method, first described in 1983. This method requires
+-- knowledge of the Lipschitz constant @l@ of the gradient, the condition
+-- number @kappa@, as well as an initial step size @alpha0@ in @(0,1)@.
+{-# INLINEABLE optimalGradient #-}
+optimalGradient :: (Additive f, Functor f, Ord a, Floating a, Epsilon a)
+                => a -> a -> (f a -> f a) -> a -> f a -> [f a]
+optimalGradient kappa l df a0' x0' = go a0' x0' x0'
+  where go a0 x0 y0 = let x1 = y0 ^-^ df y0 ^/ l
+                          alphas = quadratic 1 (a0^2 - 1/kappa) (-a0^2)
+                          a1 = case filter (\x->x >= 0 && x <= 1) alphas of
+                                 a:_  -> a
+                                 []   -> error "No solution for alpha_{k+1}"
+                          b1 = a0 * (1 - a0) / (a0^2 + a1)
+                          y1 = x1 ^+^ b1 *^ (x1 ^-^ x0)
+                      in x1 : go a0 x1 y1
+
+-- | 'quadratic a b c' is the real solutions to a quadratic equation
+-- 'a x^2 + b x + c == 0'
+quadratic :: (Ord a, Floating a, Epsilon a)
+          => a -> a -> a -> [a]
+quadratic a b c
+    | discr < 0      = []
+    | nearZero discr = [-b / 2 / a]
+    | otherwise      = [ (-b + sqrt discr) / 2 / a
+                       , (-b - sqrt discr) / 2 / a ]
+  where discr = b^2 - 4*a*c
diff --git a/src/Optimization/TrustRegion/Newton.hs b/src/Optimization/TrustRegion/Newton.hs
new file mode 100644
--- /dev/null
+++ b/src/Optimization/TrustRegion/Newton.hs
@@ -0,0 +1,37 @@
+module Optimization.TrustRegion.Newton
+    ( -- * Newton's method
+      newton
+      -- * Matrix inversion methods
+    , bicInv
+    , bicInv'
+    ) where
+
+import Control.Applicative
+import Data.Distributive (Distributive)
+import Data.Functor.Bind (Apply)
+import Data.Foldable (Foldable)
+import Linear
+
+-- | Newton's method
+{-# INLINEABLE newton #-}
+newton :: (Num a, Ord a, Additive f, Metric f, Foldable f)
+       => (f a -> f a) -> (f a -> f (f a)) -> f a -> [f a]
+newton df ddfInv x0 = iterate go x0
+  where go x = x ^-^ ddfInv x !* df x
+
+-- | Inverse by iterative method of Ben-Israel and Cohen
+-- with given starting condition
+bicInv' :: (Functor m, Distributive m, Additive m,
+            Applicative m, Apply m, Foldable m, Conjugate a)
+        => m (m a) -> m (m a) -> [m (m a)]
+bicInv' a0 a = iterate go a0
+  where go ak = 2 *!! ak !-! ak !*! a !*! ak
+
+-- | Inverse by iterative method of Ben-Israel and Cohen
+-- starting from 'alpha A^T'. Alpha should be set such that
+-- 0 < alpha < 2/sigma^2 where sigma denotes the largest singular
+-- value of A
+bicInv :: (Functor m, Distributive m, Additive m,
+           Applicative m, Apply m, Foldable m, Conjugate a)
+       => a -> m (m a) -> [m (m a)]
+bicInv alpha a = bicInv' (alpha *!! adjoint a) a
diff --git a/tests/doctests.hsc b/tests/doctests.hsc
new file mode 100644
--- /dev/null
+++ b/tests/doctests.hsc
@@ -0,0 +1,73 @@
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE ForeignFunctionInterface #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Main (doctests)
+-- Copyright   :  (C) 2012-13 Edward Kmett
+-- License     :  BSD-style (see the file LICENSE)
+-- Maintainer  :  Edward Kmett <ekmett@gmail.com>
+-- Stability   :  provisional
+-- Portability :  portable
+--
+-- This module provides doctests for a project based on the actual versions
+-- of the packages it was built with. It requires a corresponding Setup.lhs
+-- to be added to the project
+-----------------------------------------------------------------------------
+module Main where
+
+import Build_doctests (deps)
+import Control.Applicative
+import Control.Monad
+import Data.List
+import System.Directory
+import System.FilePath
+import Test.DocTest
+
+##ifdef mingw32_HOST_ARCH
+##ifdef i386_HOST_ARCH
+##define USE_CP
+import Control.Applicative
+import Control.Exception
+import Foreign.C.Types
+foreign import stdcall "windows.h SetConsoleCP" c_SetConsoleCP :: CUInt -> IO Bool
+foreign import stdcall "windows.h GetConsoleCP" c_GetConsoleCP :: IO CUInt
+##elif defined(x86_64_HOST_ARCH)
+##define USE_CP
+import Control.Applicative
+import Control.Exception
+import Foreign.C.Types
+foreign import ccall "windows.h SetConsoleCP" c_SetConsoleCP :: CUInt -> IO Bool
+foreign import ccall "windows.h GetConsoleCP" c_GetConsoleCP :: IO CUInt
+##endif
+##endif
+
+-- | Run in a modified codepage where we can print UTF-8 values on Windows.
+withUnicode :: IO a -> IO a
+##ifdef USE_CP
+withUnicode m = do
+  cp <- c_GetConsoleCP
+  (c_SetConsoleCP 65001 >> m) `finally` c_SetConsoleCP cp
+##else
+withUnicode m = m
+##endif
+
+main :: IO ()
+main = withUnicode $ getSources >>= \sources -> doctest $
+    "-isrc"
+  : "-idist/build/autogen"
+  : "-optP-include"
+  : "-optPdist/build/autogen/cabal_macros.h"
+  : "-hide-all-packages"
+  : map ("-package="++) deps ++ sources
+
+getSources :: IO [FilePath]
+getSources = filter (isSuffixOf ".hs") <$> go "src"
+  where
+    go dir = do
+      (dirs, files) <- getFilesAndDirectories dir
+      (files ++) . concat <$> mapM go dirs
+
+getFilesAndDirectories :: FilePath -> IO ([FilePath], [FilePath])
+getFilesAndDirectories dir = do
+  c <- map (dir </>) . filter (`notElem` ["..", "."]) <$> getDirectoryContents dir
+  (,) <$> filterM doesDirectoryExist c <*> filterM doesFileExist c
diff --git a/travis/cabal-apt-install b/travis/cabal-apt-install
new file mode 100644
--- /dev/null
+++ b/travis/cabal-apt-install
@@ -0,0 +1,27 @@
+#! /bin/bash
+set -eu
+
+APT="sudo apt-get -q -y"
+CABAL_INSTALL_DEPS="cabal install --only-dependencies --force-reinstall"
+
+$APT update
+$APT install dctrl-tools
+
+# Find potential system packages to satisfy cabal dependencies
+deps()
+{
+	local M='^\([^ ]\+\)-[0-9.]\+ (.*$'
+	local G=' -o ( -FPackage -X libghc-\L\1\E-dev )'
+	local E="$($CABAL_INSTALL_DEPS "$@" --dry-run -v 2> /dev/null \
+		| sed -ne "s/$M/$G/p" | sort -u)"
+	grep-aptavail -n -sPackage \( -FNone -X None \) $E | sort -u
+}
+
+$APT install $(deps "$@") libghc-quickcheck2-dev # QuickCheck is special
+$CABAL_INSTALL_DEPS "$@" # Install the rest via Hackage
+
+if ! $APT install hlint ; then
+	$APT install $(deps hlint)
+	cabal install hlint
+fi
+
diff --git a/travis/config b/travis/config
new file mode 100644
--- /dev/null
+++ b/travis/config
@@ -0,0 +1,16 @@
+-- This provides a custom ~/.cabal/config file for use when hackage is down that should work on unix
+--
+-- This is particularly useful for travis-ci to get it to stop complaining
+-- about a broken build when everything is still correct on our end.
+--
+-- This uses Luite Stegeman's mirror of hackage provided by his 'hdiff' site instead
+--
+-- To enable this, uncomment the before_script in .travis.yml
+
+remote-repo: hdiff.luite.com:http://hdiff.luite.com/packages/archive
+remote-repo-cache: ~/.cabal/packages
+world-file: ~/.cabal/world
+build-summary: ~/.cabal/logs/build.log
+remote-build-reporting: anonymous
+install-dirs user
+install-dirs global