Derivatives of a multivariative function and corresponding Jacobian with vector-space package

Question

I am having a problem with the vector-space package again. I received a very helpful answer from @mnish in a recent post, but there I only dealt with a function which depends on only 1 variable. What happens when I have, for instance, a function which maps from polar coordinates to cartesians

f:(0,oo) x [0,2pi] -> R²
(r,phi) -> (r*cos(phi),r*sin(phi))

which depends on 2 variables.

I have tried this out, with quite a naive approach:

polar :: Double -> Double -> ((Double,Double) :~> (Double,Double))
polar r phi = \(r,phi) ->  (((idD) r)*cos( idD phi),((idD) r)*sin( idD phi))

I get the following error:

Couldn't match expected type `(Double, Double) :> (Double, Double)'
            with actual type `(t0, t1)'
In the expression:
  (((idD) r) * cos (idD phi), ((idD) r) * sin (idD phi))
In the expression:
  \ (r, phi)
    -> (((idD) r) * cos (idD phi), ((idD) r) * sin (idD phi))
In an equation for `polar':
    polar r phi
      = \ (r, phi)
          -> (((idD) r) * cos (idD phi), ((idD) r) * sin (idD phi))

For one component

polarx :: Double -> Double -> ((Double,Double) :~> Double)
polarx r phi = \(r,phi) ->  ((idD) r)*cos( idD phi)

I get

Couldn't match expected type `Double'
            with actual type `(Double, Double)'
Expected type: (Double, Double) :> Double
  Actual type: (Double, Double) :> (Double, Double)
In the return type of a call of `idD'
In the first argument of `(*)', namely `((idD) r)'

Apparently there is some type disorder, but I can't figure out what is wrong.

Another question arises, when I want to calculate the Jacobian of such a mapping. As the name suggests, it has something to do with linear maps, which is, of course, covered by the package, actually it is based on those maps. But again, my Haskell knowledge is insufficient, to derive a solution on my own.

Answer 1

I finally found a solution to my problem, it was not that hard, but still it took me a while to figure it out. In case anyone else is interested I present the details.

First here is my code for the polar case:

polarCoordD :: ((Double,Double) :~> (Double,Double))
polarCoordD = \(r,phi) ->  pairD (polarx (r,phi), polary (r,phi))
where polarx :: (Double,Double) :~> Double
      polarx = \(r,phi) -> (fst . unpairD $ (idD) (r,phi))*cos( snd . unpairD $ idD (r, phi))
      polary :: (Double,Double) :~> Double
      polary = \(r,phi) -> (fst . unpairD $ (idD) (r,phi))*sin( snd . unpairD $ idD (r, phi))

The key was to make the "derivation variable" (idD) aware of the tuple (r, phi) which holds the two variables I want to differentiate. Then I have to unpack the tuple via unpairD and chose the first and the second part of the resulting pair (in polarx and polary). Both are packed again into a pair. Maybe there is a more elegant way to do this, but that's how I understood it finally.

From here it is not hard to go further to cylindrical coordinates or, in fact, to any other curved orthogonal coordinate system. For cylindrical coordinates I obtain:

cylCoordD :: (Vec3 Double :~> Vec3 Double)
cylCoordD = \(r,phi,z) ->  tripleD (cylx (r,phi,z), cyly (r,phi,z),cylz (0,0,z))
where cylx :: (Double,Double,Double) :~> Double
      cylx = \(r,phi,z) -> (fst' . untripleD $ (idD) (r,phi,z))*cos( snd' . untripleD $ idD (r, phi,z))
      cyly :: (Double,Double,Double) :~> Double
      cyly = \(r,phi,z) -> (fst' . untripleD $ (idD) (r,phi,z))*sin( snd' . untripleD $ idD (r, phi,z))
      cylz :: (Double,Double,Double) :~> Double
      cylz = \(_,_,z) -> third . untripleD $ idD (0,0,z)
      fst' :: (a,b,c) -> a
      fst' (x,_,_) = x
      snd' :: (a,b,c) -> b
      snd' (_,y,_) = y
      third :: (a,b,c) -> c
      third (_,_,z) = z

where Vec3 Double belongs to type Vec3 a = (a, a, a). Now we can even build a transformation matrix:

let transmat = \(r,phi,z) -> powVal $ liftD3 (,,) (normalized $ derivAtBasis (cylCoordD (r,phi,z)) (Left ())) (normalized $ derivAtBasis (cylCoordD (r,phi,z)) (Right (Left ()))) (normalized $ derivAtBasis (cylCoordD (r,phi,z)) (Right (Right ())))

*Main> transmat (2, rad 0, 0)
((1.0,0.0,0.0),(0.0,1.0,0.0),(0.0,0.0,1.0))

*Main> transmat (2, rad 90, 0)
((6.123233995736766e-17,1.0,0.0),(-1.0,6.123233995736766e-17,0.0),(0.0,0.0,1.0))

rad is a convenience function

rad :: Double -> Double
rad = (pi*) . (*recip 180)

Now it would be interesting to convert this "matrix" to the matrix type of Numeric Prelude and/or hmatrix, but I am not sure if this would be even useful. But still, it would be a nice example for the use of the vector-space -package.

I still have to figure out the use and especially the application of linear maps.

Answer 2

Just saw this followup question. I'm not sure what you want:

• the Jacobian matrix
• a Jacobian-vector product
• a Jacobian-transpose-vector product

In such a low-dimensional system, I'll assume the first. (The others come in handy mainly when the system is high-dimensional enough that you don't want to store or compute the Jacobian per-se, but instead treat it as a generalized sparse matrix.) In any case:

Prelude> :m + Numeric.AD
Prelude Numeric.AD> let f [r,phi] = map (r*) [cos phi, sin phi]
Prelude Numeric.AD> jacobian f [2,3::Float]
[[-0.9899925,-0.28224],[0.14112,-1.979985]]