<html><head><meta http-equiv="Content-Type" content="text/html charset=iso-8859-1"></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; ">Hi Edward,<div><br></div><div>I see now that the issues are deeper than performance.</div><div><br></div><div>I took another package that supports matrix operations: repa.</div><div><br></div><div><blockquote type="cite"><div>data MyMatrix a = MyMatrix</div><div> {</div><div> myRows :: Int</div><div> , myCols :: Int</div><div> , myElts :: [a]</div><div> } deriving (Show, Functor, Foldable, Traversable)</div><div><br></div><div>f (MyMatrix r c es) = sum es</div><div><br></div><div>g (MyMatrix r c es) = head $ toList $ sumS $ sumS n</div><div> where</div><div> n = fromListUnboxed (Z :. r :. c) es</div></blockquote></div><div><div><br></div></div><div>I can take the grad of f but not of g even though they are the same function:</div><div><br></div><div><blockquote type="cite"><div>*Main> :t grad f</div><div>grad f :: Num a => MyMatrix a -> MyMatrix a</div><div>*Main> :t grad g</div><div><br></div><div><interactive>:1:6:</div><div> Could not deduce (repa-3.2.3.1:Data.Array.Repa.Eval.Elt.Elt</div><div> (ad-3.4:Numeric.AD.Internal.Types.AD s a))</div><div> arising from a use of `g'</div><div> from the context (Num a)</div><div> bound by the inferred type of</div><div> it :: Num a => MyMatrix a -> MyMatrix a</div><div> at Top level</div><div> or from (Numeric.AD.Internal.Classes.Mode s)</div><div> bound by a type expected by the context:</div><div> Numeric.AD.Internal.Classes.Mode s =></div><div> MyMatrix (ad-3.4:Numeric.AD.Internal.Types.AD s a)</div><div> -> ad-3.4:Numeric.AD.Internal.Types.AD s a</div><div> at <interactive>:1:1-6</div><div> Possible fix:</div><div> add an instance declaration for</div><div> (repa-3.2.3.1:Data.Array.Repa.Eval.Elt.Elt</div><div> (ad-3.4:Numeric.AD.Internal.Types.AD s a))</div><div> In the first argument of `grad', namely `g'</div><div> In the expression: grad g</div></blockquote></div><div><div><br></div></div><div>2. By monomorphic, do you mean that I can do:</div><div><br></div><div><blockquote type="cite"><div>g :: MyMatrix Double -> Double</div><div>g (MyMatrix r c es) = head $ toList $ sumS $ sumS n</div><div> where</div><div> n = fromListUnboxed (Z :. r :. c) es</div></blockquote></div><div><div><br></div></div><div>but not get a type error as I currently do:</div><div><br></div><div><blockquote type="cite"><div>*Main> :t grad g</div><div><br></div><div><interactive>:1:6:</div><div> Couldn't match type `Double'</div><div> with `ad-3.4:Numeric.AD.Internal.Types.AD s a0'</div><div> Expected type: MyMatrix (ad-3.4:Numeric.AD.Internal.Types.AD s a0)</div><div> -> ad-3.4:Numeric.AD.Internal.Types.AD s a0</div><div> Actual type: MyMatrix Double -> Double</div><div> In the first argument of `grad', namely `g'</div><div> In the expression: grad g</div></blockquote></div><div><br></div><div>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.</div><div><br></div><div>Dominic.</div><div><div><br><div><div>On 9 Apr 2013, at 23:03, Edward Kmett <<a href="mailto:ekmett@gmail.com">ekmett@gmail.com</a>> wrote:</div><br class="Apple-interchange-newline"><blockquote type="cite"><div dir="ltr">hmatrix and ad don't (currently) mix.<div><br></div><div style="">The problem is that hmatrix uses a packed structure that can't hold any of the AD mode variants we have as an Element. =(</div><div style="">
<br></div><div style="">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.<br>
</div><div style=""><br></div><div style="">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. =/</div><div style=""><br></div><div style="">
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.</div>
<div style=""><br></div><div style="">This is a long term goal, but not one you're likely to see support for out of 'ad' in the short term.</div><div style=""><br></div><div style=""><div>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.</div>
</div><div style=""><br></div><div style="">-Edward</div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Tue, Apr 9, 2013 at 10:46 AM, Dominic Steinitz <span dir="ltr"><<a href="mailto:dominic@steinitz.org" target="_blank">dominic@steinitz.org</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hi Cafe,<br>
<br>
Suppose I want to find the grad of a function then it's easy I just<br>
use <a href="http://hackage.haskell.org/package/ad-3.4" target="_blank">http://hackage.haskell.org/package/ad-3.4</a>:<br>
<br>
import Numeric.AD<br>
import Data.Foldable (Foldable)<br>
import Data.Traversable (Traversable)<br>
<br>
data MyMatrix a = MyMatrix (a, a)<br>
deriving (Show, Functor, Foldable, Traversable)<br>
<br>
f :: Floating a => MyMatrix a -> a<br>
f (MyMatrix (x, y)) = exp $ negate $ (x^2 + y^2) / 2.0<br>
<br>
main :: IO ()<br>
main = do<br>
putStrLn $ show $ f $ MyMatrix (0.0, 0.0)<br>
putStrLn $ show $ grad f $ MyMatrix (0.0, 0.0)<br>
<br>
But now suppose I am doing some matrix calculations<br>
<a href="http://hackage.haskell.org/package/hmatrix-0.14.1.0" target="_blank">http://hackage.haskell.org/package/hmatrix-0.14.1.0</a> and I want to find<br>
the grad of a function of a matrix:<br>
<br>
import Numeric.AD<br>
import Numeric.LinearAlgebra<br>
import Data.Foldable (Foldable)<br>
import Data.Traversable (Traversable)<br>
<br>
g :: (Element a, Floating a) => Matrix a -> a<br>
g m = exp $ negate $ (x^2 + y^2) / 2.0<br>
where r = (toLists m)!!0<br>
x = r!!0<br>
y = r!!1<br>
<br>
main :: IO ()<br>
main = do<br>
putStrLn $ show $ g $ (1 >< 2) ([0.0, 0.0] :: [Double])<br>
putStrLn $ show $ grad g $ (1 >< 2) ([0.0, 0.0] :: [Double])<br>
<br>
Then I am in trouble:<br>
<br>
/Users/dom/Dropbox/Private/Whales/MyAD.hs:24:21:<br>
No instance for (Traversable Matrix) arising from a use of `grad'<br>
Possible fix: add an instance declaration for (Traversable Matrix)<br>
In the expression: grad g<br>
In the second argument of `($)', namely<br>
`grad g $ (1 >< 2) ([0.0, 0.0] :: [Double])'<br>
In the second argument of `($)', namely<br>
`show $ grad g $ (1 >< 2) ([0.0, 0.0] :: [Double])'<br>
<br>
/Users/dom/Dropbox/Private/Whales/MyAD.hs:24:26:<br>
Could not deduce (Element<br>
(ad-3.4:<a href="http://numeric.ad.internal.types.ad/" target="_blank">Numeric.AD.Internal.Types.AD</a> s Double))<br>
arising from a use of `g'<br>
from the context (Numeric.AD.Internal.Classes.Mode s)<br>
bound by a type expected by the context:<br>
Numeric.AD.Internal.Classes.Mode s =><br>
Matrix (ad-3.4:<a href="http://numeric.ad.internal.types.ad/" target="_blank">Numeric.AD.Internal.Types.AD</a> s Double)<br>
-> ad-3.4:<a href="http://numeric.ad.internal.types.ad/" target="_blank">Numeric.AD.Internal.Types.AD</a> s Double<br>
at /Users/dom/Dropbox/Private/Whales/MyAD.hs:24:21-26<br>
Possible fix:<br>
add an instance declaration for<br>
(Element (ad-3.4:<a href="http://numeric.ad.internal.types.ad/" target="_blank">Numeric.AD.Internal.Types.AD</a> s Double))<br>
In the first argument of `grad', namely `g'<br>
In the expression: grad g<br>
In the second argument of `($)', namely<br>
`grad g $ (1 >< 2) ([0.0, 0.0] :: [Double])'<br>
<br>
What are my options here? Clearly I can convert my matrix into a list<br>
(which is traversable), find the grad and convert it back into a<br>
matrix but given I am doing numerical calculations and speed is an<br>
important factor, this seems undesirable.<br>
<br>
I think I would have the same problem with:<br>
<br>
<a href="http://hackage.haskell.org/package/repa" target="_blank">http://hackage.haskell.org/package/repa</a><br>
<a href="http://hackage.haskell.org/package/yarr-1.3.1" target="_blank">http://hackage.haskell.org/package/yarr-1.3.1</a><br>
<br>
although I haven'¯t checked.<br>
<br>
Thanks, Dominic.<br>
<br>
<br>
_______________________________________________<br>
Haskell-Cafe mailing list<br>
<a href="mailto:Haskell-Cafe@haskell.org">Haskell-Cafe@haskell.org</a><br>
<a href="http://www.haskell.org/mailman/listinfo/haskell-cafe" target="_blank">http://www.haskell.org/mailman/listinfo/haskell-cafe</a><br>
</blockquote></div><br></div>
</blockquote></div><br></div></div></body></html>