diff --git a/CHANGELOG.md b/CHANGELOG.md
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,8 @@
 # Revision history for downhill
 
+## 0.4.0.0
+* `Manifold` class
+
 ## 0.2.0.0
 * `MetricTensor` is no longer required to be a `VectorSpace`
 * `T2`, `T3` pattern synonyms for `BVar`
diff --git a/downhill.cabal b/downhill.cabal
--- a/downhill.cabal
+++ b/downhill.cabal
@@ -1,7 +1,7 @@
 cabal-version:       2.4
 
 name:                downhill
-version:             0.3.0.0
+version:             0.4.0.0
 synopsis:            Reverse mode automatic differentiation
 homepage:            https://andriusstank.github.io/downhill/
 description:
@@ -39,9 +39,7 @@
   build-depends:       base                  >= 4.12.0.0 && <4.18,
                        containers            >= 0.6.5 && < 0.7,
                        reflection            >= 2.1.6 && < 2.2,
-                       template-haskell      >= 2.16.0 && < 2.20,
                        transformers          >= 0.5.6 && < 0.7,
-                       th-abstraction        >= 0.4.3 && < 0.5,
                        unordered-containers  >= 0.2.14 && < 0.3,
                        vector-space          >= 0.16 && < 0.17,
   hs-source-dirs:      src
diff --git a/src/Downhill/BVar.hs b/src/Downhill/BVar.hs
--- a/src/Downhill/BVar.hs
+++ b/src/Downhill/BVar.hs
@@ -6,13 +6,13 @@
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE GADTs #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE PatternSynonyms #-}
 {-# LANGUAGE PolyKinds #-}
 {-# LANGUAGE RankNTypes #-}
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE TypeOperators #-}
 {-# LANGUAGE UndecidableInstances #-}
-{-# LANGUAGE PatternSynonyms #-}
 {-# LANGUAGE ViewPatterns #-}
 
 module Downhill.BVar
@@ -20,9 +20,10 @@
     var,
     constant,
     backprop,
+
     -- * Pattern synonyms
     pattern T2,
-    pattern T3
+    pattern T3,
   )
 where
 
@@ -32,23 +33,25 @@
 import Data.VectorSpace
   ( AdditiveGroup (..),
     InnerSpace ((<.>)),
-    VectorSpace ((*^)),
+    VectorSpace (Scalar, (*^)),
   )
-import qualified Data.VectorSpace as VectorSpace
 import Downhill.Grad
   ( Dual (evalGrad),
-    HasGrad (Grad, Tang),
-    HasGradAffine, MScalar
+    HasGrad,
+    HasGradAffine,
+    HilbertSpace (coriesz, riesz),
+    MScalar,
+    Manifold (Grad, Tang),
   )
 import Downhill.Linear.BackGrad
   ( BackGrad (..),
     realNode,
   )
 import qualified Downhill.Linear.Backprop as BP
-import Downhill.Linear.Expr (BasicVector, Expr (ExprVar))
-import Downhill.Linear.Lift (lift2_dense)
-import Prelude hiding (id, (.))
+import Downhill.Linear.Expr (BasicVector (..), Expr (ExprVar))
+import Downhill.Linear.Lift (lift1_dense, lift2_dense)
 import qualified Downhill.Linear.Prelude as Linear
+import Prelude hiding (id, (.))
 
 -- | Variable is a value paired with derivative.
 data BVar r a = BVar
@@ -103,8 +106,8 @@
   ( VectorSpace v,
     HasGrad v,
     Tang v ~ v,
-    BasicVector (MScalar v),
-    Grad (MScalar v) ~ MScalar v
+    HasGrad (MScalar v),
+    Grad (Scalar v) ~ Scalar v
   ) =>
   VectorSpace (BVar r v)
   where
@@ -121,11 +124,60 @@
   BVar y0 dy .+^ BVar z0 dz = BVar (y0 .+^ z0) (dy ^+^ dz)
   BVar y0 dy .-. BVar z0 dz = BVar (y0 .-. z0) (dy ^-^ dz)
 
+-- maybe move all those equality constraints to Dual class?
 instance
+  ( HasGrad (Scalar v),
+    HasGrad v,
+    HasGrad dv,
+    Dual v dv,
+    Grad dv ~ v,
+    Grad v ~ dv,
+    Tang v ~ v,
+    Tang dv ~ dv,
+    Grad (Scalar dv) ~ Scalar dv
+  ) =>
+  Dual (BVar r v) (BVar r dv)
+  where
+  evalGrad (BVar dv d_dv) (BVar v d_v) = BVar (evalGrad dv v) (lift2_dense (*^ v) (*^ dv) d_dv d_v)
+
+instance
+  ( HasGrad (MScalar p),
+    HasGrad (Tang p),
+    HasGrad (Grad p),
+    Grad (Grad p) ~ Tang p,
+    Tang (Grad p) ~ Grad p,
+    Tang (Tang p) ~ Tang p,
+    Grad (Tang p) ~ Grad p,
+    Grad (MScalar p) ~ MScalar p,
+    Scalar (Grad p) ~ Scalar (Tang p),
+    Manifold p
+  ) =>
+  Manifold (BVar r p)
+  where
+  type Tang (BVar r p) = BVar r (Tang p)
+  type Grad (BVar r p) = BVar r (Grad p)
+
+instance
+  ( HilbertSpace v dv,
+    HasGrad (Scalar v),
+    HasGrad v,
+    HasGrad dv,
+    Grad dv ~ v,
+    Grad v ~ dv,
+    Tang v ~ v,
+    Tang dv ~ dv,
+    Grad (Scalar dv) ~ Scalar dv
+  ) =>
+  HilbertSpace (BVar r v) (BVar r dv)
+  where
+  riesz (BVar v dv) = BVar (riesz v) (lift1_dense riesz dv)
+  coriesz (BVar v dv) = BVar (coriesz v) (lift1_dense coriesz dv)
+
+instance
   ( VectorSpace v,
     HasGrad v,
-    Grad v ~ v,
     Tang v ~ v,
+    HilbertSpace (Tang v) (Grad v),
     BasicVector (MScalar v),
     Grad (MScalar v) ~ MScalar v,
     InnerSpace v,
@@ -136,9 +188,9 @@
   BVar u du <.> BVar v dv = BVar (u <.> v) (lift2_dense bpU bpV du dv)
     where
       bpU :: MScalar v -> Grad v
-      bpU dz = dz *^ v
+      bpU dz = dz *^ riesz v
       bpV :: MScalar v -> Grad v
-      bpV dz = dz *^ u
+      bpV dz = dz *^ riesz u
 
 -- | A variable with derivative of zero.
 constant :: forall r a. (BasicVector (Grad a), AdditiveGroup (Grad a)) => a -> BVar r a
@@ -148,27 +200,30 @@
 var :: a -> BVar (Grad a) a
 var x = BVar x (realNode ExprVar)
 
---backprop :: forall a p. (HasGrad p, BasicVector a) => BVar a p -> GradBuilder p -> a
---backprop (BVar _y0 x) = BP.backprop x
-
 -- | Reverse mode differentiation.
---
--- 
 backprop :: forall r a. (HasGrad a, BasicVector r) => BVar r a -> Grad a -> r
 backprop (BVar _y0 x) = BP.backprop x
 
-
 splitPair :: (BasicVector (Grad a), BasicVector (Grad b)) => BVar r (a, b) -> (BVar r a, BVar r b)
 splitPair (BVar (a, b) (Linear.T2 da db)) = (BVar a da, BVar b db)
 
 pattern T2 :: forall r a b. (BasicVector (Grad a), BasicVector (Grad b)) => BVar r a -> BVar r b -> BVar r (a, b)
-pattern T2 a b <- (splitPair -> (a, b))
-  where T2 (BVar a da) (BVar b db) = BVar (a, b) (Linear.T2 da db)
+pattern T2 a b <-
+  (splitPair -> (a, b))
+  where
+    T2 (BVar a da) (BVar b db) = BVar (a, b) (Linear.T2 da db)
 
 splitTriple :: (BasicVector (Grad a), BasicVector (Grad b), BasicVector (Grad c)) => BVar r (a, b, c) -> (BVar r a, BVar r b, BVar r c)
 splitTriple (BVar (a, b, c) (Linear.T3 da db dc)) = (BVar a da, BVar b db, BVar c dc)
 
-pattern T3 :: forall r a b c. (BasicVector (Grad a), BasicVector (Grad b), BasicVector (Grad c))
- => BVar r a -> BVar r b -> BVar r c -> BVar r (a, b, c)
-pattern T3 a b c <- (splitTriple -> (a, b, c))
-  where T3 (BVar a da) (BVar b db) (BVar c dc) = BVar (a, b, c) (Linear.T3 da db dc)
+pattern T3 ::
+  forall r a b c.
+  (BasicVector (Grad a), BasicVector (Grad b), BasicVector (Grad c)) =>
+  BVar r a ->
+  BVar r b ->
+  BVar r c ->
+  BVar r (a, b, c)
+pattern T3 a b c <-
+  (splitTriple -> (a, b, c))
+  where
+    T3 (BVar a da) (BVar b db) (BVar c dc) = BVar (a, b, c) (Linear.T3 da db dc)
diff --git a/src/Downhill/BVar/Num.hs b/src/Downhill/BVar/Num.hs
--- a/src/Downhill/BVar/Num.hs
+++ b/src/Downhill/BVar/Num.hs
@@ -34,7 +34,7 @@
 import qualified Downhill.BVar as BVar
 import Downhill.Grad
   ( Dual (evalGrad),
-    HasGrad (Grad, Tang)
+    Manifold (Grad, Tang)
   )
 import Downhill.Linear.Expr (BasicVector (..))
 import Downhill.Metric (MetricTensor (evalMetric))
@@ -49,7 +49,7 @@
 instance Num a => Dual (AsNum a) (AsNum a) where
   evalGrad = (*)
 
-instance Num a => HasGrad (AsNum a) where
+instance Num a => Manifold (AsNum a) where
   type Grad (AsNum a) = AsNum a
   type Tang (AsNum a) = AsNum a
 
diff --git a/src/Downhill/BVar/Traversable.hs b/src/Downhill/BVar/Traversable.hs
--- a/src/Downhill/BVar/Traversable.hs
+++ b/src/Downhill/BVar/Traversable.hs
@@ -63,7 +63,7 @@
 import Downhill.BVar (BVar (BVar, bvarGrad, bvarValue), backprop, var)
 import Downhill.Grad
   ( Dual (evalGrad),
-    HasGrad (Grad, Tang)
+    Manifold (Grad, Tang), HasGrad
   )
 import Downhill.Linear.BackGrad (BackGrad (BackGrad), castBackGrad, realNode)
 import Downhill.Linear.Expr
@@ -92,13 +92,17 @@
   evalMetric (TraversableMetric m) (IntmapVector da) =
     IntmapVector (IntMap.map (evalMetric @p @g m) da)
 
-instance HasGrad a => HasGrad (TraversableVar f a) where
+instance Manifold a => Manifold (TraversableVar f a) where
   type Tang (TraversableVar f a) = IntmapVector f (Tang a)
   type Grad (TraversableVar f a) = IntmapVector f (Grad a)
 
 -- | @IntmapVector@ serves as a gradient of 'TraversableVar'.
 newtype IntmapVector (f :: Type -> Type) v = IntmapVector {unIntmapVector :: IntMap v}
   deriving (Show)
+
+instance Manifold v => Manifold (IntmapVector f v) where
+  type Tang (IntmapVector f v) = IntmapVector f (Tang v)
+  type Grad (IntmapVector f v) = IntmapVector f (Grad v)
 
 instance AdditiveGroup a => AdditiveGroup (IntmapVector f a) where
   zeroV = IntmapVector IntMap.empty
diff --git a/src/Downhill/Grad.hs b/src/Downhill/Grad.hs
--- a/src/Downhill/Grad.hs
+++ b/src/Downhill/Grad.hs
@@ -13,8 +13,8 @@
 {-# LANGUAGE UndecidableInstances #-}
 
 module Downhill.Grad
-  ( Dual (..),
-    HasGrad (..), MScalar,
+  ( Dual (..), HilbertSpace(..),
+    Manifold(..), HasGrad, MScalar,
     GradBuilder,
     HasGradAffine,
   )
@@ -41,16 +41,20 @@
   default evalGrad :: (GDual (Scalar v) (Rep v) (Rep dv), Generic dv, Generic v) => dv -> v -> Scalar v
   evalGrad dv v = gevalGrad (from dv) (from v)
 
+-- | u <.> v = evalDual (riesz u) v
+-- | du <.> dv = evalDual du (coriesz dv)
+class Dual v dv => HilbertSpace v dv where
+  riesz :: v -> dv
+  coriesz :: dv -> v
+
+
 type MScalar p = Scalar (Tang p)
 
--- | Differentiable functions don't need to be constrained to vector spaces, they
--- can be defined on other smooth manifolds, too.
 class
   ( Dual (Tang p) (Grad p),
-    BasicVector (Grad p),
     Scalar (Tang p) ~ Scalar (Grad p)
   ) =>
-  HasGrad p
+  Manifold p
   where
   -- | Tangent space.
   type Tang p :: Type
@@ -58,6 +62,10 @@
   -- | Cotangent space.
   type Grad p :: Type
 
+-- | Differentiable functions don't need to be constrained to vector spaces, they
+-- can be defined on other smooth manifolds, too.
+type HasGrad p = (Manifold p, BasicVector (Grad p))
+
 type GradBuilder v = VecBuilder (Grad v)
 
 type HasGradAffine p =
@@ -72,7 +80,7 @@
 instance Dual Integer Integer where
   evalGrad = (*)
 
-instance HasGrad Integer where
+instance Manifold Integer where
   type Tang Integer = Integer
   type Grad Integer = Integer
 
@@ -87,7 +95,7 @@
     HasGrad b,
     MScalar b ~ MScalar a
   ) =>
-  HasGrad (a, b)
+  Manifold (a, b)
   where
   type Grad (a, b) = (Grad a, Grad b)
   type Tang (a, b) = (Tang a, Tang b)
@@ -99,7 +107,7 @@
     MScalar b ~ MScalar a,
     MScalar c ~ MScalar a
   ) =>
-  HasGrad (a, b, c)
+  Manifold (a, b, c)
   where
   type Grad (a, b, c) = (Grad a, Grad b, Grad c)
   type Tang (a, b, c) = (Tang a, Tang b, Tang c)
@@ -107,14 +115,14 @@
 instance Dual Float Float where
   evalGrad = (*)
 
-instance HasGrad Float where
+instance Manifold Float where
   type Grad Float = Float
   type Tang Float = Float
 
 instance Dual Double Double where
   evalGrad = (*)
 
-instance HasGrad Double where
+instance Manifold Double where
   type Grad Double = Double
   type Tang Double = Double
 
diff --git a/src/Downhill/Metric.hs b/src/Downhill/Metric.hs
--- a/src/Downhill/Metric.hs
+++ b/src/Downhill/Metric.hs
@@ -14,7 +14,7 @@
 where
 
 import Data.VectorSpace ((^+^))
-import Downhill.Grad (Dual (evalGrad), HasGrad (Grad, Tang), MScalar)
+import Downhill.Grad (Dual (evalGrad), Manifold(..), MScalar)
 
 -- | @MetricTensor@ converts gradients to vectors.
 --
diff --git a/test/DownhillTest/Bilinear.hs b/test/DownhillTest/Bilinear.hs
--- a/test/DownhillTest/Bilinear.hs
+++ b/test/DownhillTest/Bilinear.hs
@@ -1,40 +1,73 @@
-{-# LANGUAGE FlexibleContexts #-}
-{-# LANGUAGE RankNTypes #-}
-{-# LANGUAGE TypeApplications #-}
-{-# LANGUAGE DerivingVia #-}
 {-# LANGUAGE DeriveGeneric #-}
-{-# LANGUAGE UndecidableInstances #-}
-{-# LANGUAGE TemplateHaskell #-}
-{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE DerivingVia #-}
+{-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE OverloadedStrings #-}
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE TypeApplications #-}
 {-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE UndecidableInstances #-}
 
 module DownhillTest.Bilinear where
 
 import Data.AffineSpace ((.+^))
-import Data.VectorSpace (AdditiveGroup, VectorSpace ((*^), Scalar), (^+^))
+import Data.VectorSpace (AdditiveGroup, VectorSpace (Scalar, (*^)), (^+^))
 import Downhill.BVar (BVar (bvarValue))
 import qualified Downhill.BVar as BVar
-import Downhill.Grad (Dual (evalGrad), HasGrad (Grad), MScalar)
+import Downhill.Grad (Dual (evalGrad), HasGrad, MScalar, Manifold (Grad))
+import Downhill.Linear.Expr (BasicVector, DenseVector (DenseVector))
+import GHC.Base (VecElem (DoubleElemRep))
+import GHC.Generics (Generic)
 import Hedgehog
   ( Gen,
     Property,
+    PropertyT,
     forAll,
     property,
     (===),
   )
-import Test.Tasty (TestTree, testGroup)
-import Test.Tasty.Hedgehog (testProperty, testPropertyNamed)
-import GHC.Base (VecElem(DoubleElemRep))
-import Hedgehog.Internal.Show (Value(Integer))
-import qualified Hedgehog.Internal.Show as Gen
 import qualified Hedgehog.Gen as Gen
+import Hedgehog.Internal.Show (Value (Integer))
+import qualified Hedgehog.Internal.Show as Gen
 import qualified Hedgehog.Range as Range
-import GHC.Generics (Generic)
-import Downhill.Linear.Expr (BasicVector, DenseVector (DenseVector))
+import Test.Tasty (TestTree, testGroup)
+import Test.Tasty.Hedgehog (testProperty, testPropertyNamed)
 
+testLinear ::
+  forall m r u z.
+  ( Show u,
+    HasGrad u,
+    Show (Grad z),
+    HasGrad z,
+    Eq z,
+    AdditiveGroup u,
+    Show z,
+    AdditiveGroup z,
+    Dual (Grad u) u,
+    Eq (Scalar u),
+    Show (Scalar u),
+    Scalar u ~ Scalar z,
+    Dual (Grad z) z,
+    Show (MScalar z),
+    Monad m
+  ) =>
+  (u -> z) ->
+  (forall r. BVar r u -> BVar r z) ->
+  Gen u ->
+  Gen (Grad z) ->
+  PropertyT m ()
+testLinear f bf genU genDZ = do
+  u <- forAll genU
+  dz <- forAll genDZ
+  let z = f u
+      bu = BVar.var u
+      bz = bf bu
+      du = BVar.backprop bz dz
+  bvarValue bz === z -- check that `f` and `bf` is the same function
+  evalGrad u du === evalGrad z dz
+
 testBilinear ::
   ( Show u,
     Show v,
@@ -51,7 +84,10 @@
     Show (Scalar u),
     Scalar u ~ Scalar z,
     Scalar v ~ Scalar z,
-    Dual (Grad z) z, Show (MScalar z), Dual (Grad v) v) =>
+    Dual (Grad z) z,
+    Show (MScalar z),
+    Dual (Grad v) v
+  ) =>
   (u -> v -> z) ->
   (forall r. BVar r u -> BVar r v -> BVar r z) ->
   Gen u ->
@@ -70,10 +106,52 @@
   evalGrad u du === evalGrad z dz
   evalGrad v dv === evalGrad z dz
 
+testBilinear' ::
+  ( Show u,
+    Show v,
+    HasGrad u,
+    HasGrad v,
+    Show (Grad z),
+    HasGrad z,
+    Eq z,
+    AdditiveGroup u,
+    Show z,
+    AdditiveGroup z,
+    Dual (Grad u) u,
+    Eq (Scalar u),
+    Show (Scalar u),
+    Scalar u ~ Scalar z,
+    Scalar v ~ Scalar z,
+    Dual (Grad z) z,
+    Show (MScalar z),
+    Dual (Grad v) v
+  ) =>
+  (u -> v -> z) ->
+  (forall r. BVar r u -> BVar r v -> BVar r z) ->
+  Gen u ->
+  Gen v ->
+  Gen (Grad z) ->
+  Property
+testBilinear' f bf genU genV genDZ = property $ do
+  u <- forAll genU
+  v <- forAll genV
+  dz <- forAll genDZ
+  let z = f u v
+      BVar.T2 bu bv = BVar.var (u, v)
+      bz = bf bu bv
+      (du, dv) = BVar.backprop bz dz
+  bvarValue bz === z -- check that `f` and `bf` is the same function
+  evalGrad u du === evalGrad z dz
+  evalGrad v dv === evalGrad z dz
+
+  testLinear (\x -> f x v) (\bx -> bf bx (BVar.constant v)) genU genDZ
+  testLinear (\x -> f u x) (\bx -> bf (BVar.constant u) bx) genV genDZ
+
 data Vector = Vector Integer Integer
-  deriving Generic
+  deriving (Generic)
 
 instance AdditiveGroup Vector
+
 instance VectorSpace Vector
 
 bilinearIntMulProperty :: Property
@@ -84,10 +162,10 @@
     genInt :: Gen Integer
     genInt = Gen.integral (Range.linear (-100) 100)
 
-
 bilinearTests :: TestTree
 bilinearTests =
-   testGroup "Bilinear operations"
-     [ testPropertyNamed "Scalar multiplication" "bilinearIntMulProperty" bilinearIntMulProperty
-       -- TODO: scalar-vector product, inner product
-     ]
+  testGroup
+    "Bilinear operations"
+    [ testPropertyNamed "Scalar multiplication" "bilinearIntMulProperty" bilinearIntMulProperty
+    -- TODO: scalar-vector product, inner product
+    ]
diff --git a/test/DownhillTest/Traversable.hs b/test/DownhillTest/Traversable.hs
--- a/test/DownhillTest/Traversable.hs
+++ b/test/DownhillTest/Traversable.hs
@@ -8,7 +8,7 @@
 
 import Downhill.BVar.Traversable (TraversableVar (TraversableVar), backpropTraversable, splitTraversable)
 import Downhill.BVar (BVar (BVar), backprop, var)
-import Downhill.Grad (HasGrad (Grad))
+import Downhill.Grad (Manifold (Grad), HasGrad)
 import Test.Tasty (TestTree)
 import Test.Tasty.HUnit (testCase, (@?=))
 
@@ -18,7 +18,7 @@
   }
   deriving (Eq, Functor, Foldable, Traversable, Show)
 
-deriving via (TraversableVar MyRecord a) instance HasGrad a => HasGrad (MyRecord a)
+deriving via (TraversableVar MyRecord a) instance Manifold a => Manifold (MyRecord a)
 
 test_r :: MyRecord Integer
 test_r = MyRecord (10, 11) [12, 13, 14]