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

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 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 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@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe