[Haskell-cafe] Automated Differentiation of Matrices (hmatrix)

Dominic Steinitz dominic at steinitz.org
Wed Apr 10 18:39:36 CEST 2013


Hi Edward,

I see now that the issues are deeper than performance.

I took another package that supports matrix operations: repa.

> data MyMatrix a = MyMatrix
>   {
>     myRows :: Int
>   , myCols :: Int
>   , myElts :: [a]
>   } deriving (Show, Functor, Foldable, Traversable)
> 
> f (MyMatrix r c es) = sum es
> 
> g (MyMatrix r c es) = head $ toList $ sumS $ sumS n
>   where
>     n   = fromListUnboxed (Z :. r :. c) es


I can take the grad of f but not of g even though they are the same function:

> *Main> :t grad f
> grad f :: Num a => MyMatrix a -> MyMatrix a
> *Main> :t grad g
> 
> <interactive>:1:6:
>     Could not deduce (repa-3.2.3.1:Data.Array.Repa.Eval.Elt.Elt
>                         (ad-3.4:Numeric.AD.Internal.Types.AD s a))
>       arising from a use of `g'
>     from the context (Num a)
>       bound by the inferred type of
>                it :: Num a => MyMatrix a -> MyMatrix a
>       at Top level
>     or from (Numeric.AD.Internal.Classes.Mode s)
>       bound by a type expected by the context:
>                  Numeric.AD.Internal.Classes.Mode s =>
>                  MyMatrix (ad-3.4:Numeric.AD.Internal.Types.AD s a)
>                  -> ad-3.4:Numeric.AD.Internal.Types.AD s a
>       at <interactive>:1:1-6
>     Possible fix:
>       add an instance declaration for
>       (repa-3.2.3.1:Data.Array.Repa.Eval.Elt.Elt
>          (ad-3.4:Numeric.AD.Internal.Types.AD s a))
>     In the first argument of `grad', namely `g'
>     In the expression: grad g


2. By monomorphic, do you mean that I can do:

> g :: MyMatrix Double -> Double
> g (MyMatrix r c es) = head $ toList $ sumS $ sumS n
>   where
>     n   = fromListUnboxed (Z :. r :. c) es


but not get a type error as I currently do:

> *Main> :t grad g
> 
> <interactive>:1:6:
>     Couldn't match type `Double'
>                   with `ad-3.4:Numeric.AD.Internal.Types.AD s a0'
>     Expected type: MyMatrix (ad-3.4:Numeric.AD.Internal.Types.AD s a0)
>                    -> ad-3.4:Numeric.AD.Internal.Types.AD s a0
>       Actual type: MyMatrix Double -> Double
>     In the first argument of `grad', namely `g'
>     In the expression: grad g


If so that would help as at least I could then convert between e.g. repa and structures that ad can play with. Of course, better would be that I could just apply ad to e.g. repa.

Dominic.

On 9 Apr 2013, at 23:03, Edward Kmett <ekmett at gmail.com> wrote:

> 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/20130410/dd96bede/attachment-0001.htm>


More information about the Haskell-Cafe mailing list