[Haskell-cafe] Automated Differentiation of Matrices (hmatrix)
Edward Kmett
ekmett at gmail.com
Wed Apr 10 00:03:22 CEST 2013
hmatrix and ad don't (currently) mix.
The problem is that hmatrix uses a packed structure that can't hold any of
the AD mode variants we have as an Element. =(
I've been working with Alex Lang to explore in ad 4.0 ways that we can
support monomorphic AD modes and still present largely the same API. We've
seen a number of performance gains off of this refactoring already, but it
doesn't go far enough to address what you need.
A goal a bit farther out is to support AD on vector/matrix operations, but
it is a much bigger refactoring than the one currently in the pipeline. =/
To support automatic differentiation on vector-based operations in a form
that works with standard BLAS-like storage like the packed matrix rep used
in hmatrix we need to convert from a 'matrix of AD variables' to an 'AD
mode over of matrices'. This is similar to the difference between a matrix
of complex numbers and a real matrix plus an imaginary matrix.
This is a long term goal, but not one you're likely to see support for out
of 'ad' in the short term.
I can't build AD on hmatrix itself due in part to licensing restrictions
and differing underlying storage requirements, so there are a lot of little
issues in making that latter vision a reality.
-Edward
On Tue, Apr 9, 2013 at 10:46 AM, Dominic Steinitz <dominic at steinitz.org>wrote:
> Hi Cafe,
>
> Suppose I want to find the grad of a function then it's easy I just
> use http://hackage.haskell.org/package/ad-3.4:
>
> import Numeric.AD
> import Data.Foldable (Foldable)
> import Data.Traversable (Traversable)
>
> data MyMatrix a = MyMatrix (a, a)
> deriving (Show, Functor, Foldable, Traversable)
>
> f :: Floating a => MyMatrix a -> a
> f (MyMatrix (x, y)) = exp $ negate $ (x^2 + y^2) / 2.0
>
> main :: IO ()
> main = do
> putStrLn $ show $ f $ MyMatrix (0.0, 0.0)
> putStrLn $ show $ grad f $ MyMatrix (0.0, 0.0)
>
> But now suppose I am doing some matrix calculations
> http://hackage.haskell.org/package/hmatrix-0.14.1.0 and I want to find
> the grad of a function of a matrix:
>
> import Numeric.AD
> import Numeric.LinearAlgebra
> import Data.Foldable (Foldable)
> import Data.Traversable (Traversable)
>
> g :: (Element a, Floating a) => Matrix a -> a
> g m = exp $ negate $ (x^2 + y^2) / 2.0
> where r = (toLists m)!!0
> x = r!!0
> y = r!!1
>
> main :: IO ()
> main = do
> putStrLn $ show $ g $ (1 >< 2) ([0.0, 0.0] :: [Double])
> putStrLn $ show $ grad g $ (1 >< 2) ([0.0, 0.0] :: [Double])
>
> Then I am in trouble:
>
> /Users/dom/Dropbox/Private/Whales/MyAD.hs:24:21:
> No instance for (Traversable Matrix) arising from a use of `grad'
> Possible fix: add an instance declaration for (Traversable Matrix)
> In the expression: grad g
> In the second argument of `($)', namely
> `grad g $ (1 >< 2) ([0.0, 0.0] :: [Double])'
> In the second argument of `($)', namely
> `show $ grad g $ (1 >< 2) ([0.0, 0.0] :: [Double])'
>
> /Users/dom/Dropbox/Private/Whales/MyAD.hs:24:26:
> Could not deduce (Element
> (ad-3.4:Numeric.AD.Internal.Types.AD s Double))
> arising from a use of `g'
> from the context (Numeric.AD.Internal.Classes.Mode s)
> bound by a type expected by the context:
> Numeric.AD.Internal.Classes.Mode s =>
> Matrix (ad-3.4:Numeric.AD.Internal.Types.AD s Double)
> -> ad-3.4:Numeric.AD.Internal.Types.AD s Double
> at /Users/dom/Dropbox/Private/Whales/MyAD.hs:24:21-26
> Possible fix:
> add an instance declaration for
> (Element (ad-3.4:Numeric.AD.Internal.Types.AD s Double))
> In the first argument of `grad', namely `g'
> In the expression: grad g
> In the second argument of `($)', namely
> `grad g $ (1 >< 2) ([0.0, 0.0] :: [Double])'
>
> What are my options here? Clearly I can convert my matrix into a list
> (which is traversable), find the grad and convert it back into a
> matrix but given I am doing numerical calculations and speed is an
> important factor, this seems undesirable.
>
> I think I would have the same problem with:
>
> http://hackage.haskell.org/package/repa
> http://hackage.haskell.org/package/yarr-1.3.1
>
> although I haven'¯t checked.
>
> Thanks, Dominic.
>
>
> _______________________________________________
> Haskell-Cafe mailing list
> Haskell-Cafe at haskell.org
> http://www.haskell.org/mailman/listinfo/haskell-cafe
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.haskell.org/pipermail/haskell-cafe/attachments/20130409/fcb1be92/attachment.htm>
More information about the Haskell-Cafe
mailing list