rderivativeautodiff

Compute partial derivatives with `madness`


The madness packages, as mentioned here, is nice for autodiff in R.

I would like to compute now a derivative wrt x of a derivative wrt y.

$\frac{\partial}{\partial x}\frac{\partial}{\partial y}xy$

how can this be done using madness?

update: actually here I guess it factors.. maybe this will be ok by just multiplying the two derivatives? Maybe this will only be difficult if x is a function of y.


Solution

  • Here's a way using the numderiv function in madness:

    library(madness)
    
    dxdy <- function(x, y, f) {
      dy <- function(x, y) {
        dvdx(f(x, y))
      }
      
      numderiv(dy, x, y = y)
    }
    
    x <- matrix(1:3, nrow = 1)
    y <- matrix(3:1, ncol = 1)
    # identity matrix, as expected
    dvdx(dxdy(madness(x), madness(y), function(x, y) x%*%y))
    #>      [,1] [,2] [,3]
    #> [1,]    1    0    0
    #> [2,]    0    1    0
    #> [3,]    0    0    1
    x <- matrix(2, ncol = 1)
    y <- matrix(3, ncol = 1)
    dvdx(dxdy(madness(x), madness(y), function(x, y) y^x))
    #>          [,1]
    #> [1,] 9.591674
    # compare to analytical solution
    y^(x-1)*(x*log(y) + 1)
    #>          [,1]
    #> [1,] 9.591674
    x <- matrix(1:3, ncol = 1)
    y <- matrix(3:1, ncol = 1)
    dvdx(dxdy(madness(x), madness(y), function(x, y) sum(y^x)))
    #>          [,1]     [,2] [,3]
    #> [1,] 2.098612 0.000000    0
    #> [2,] 0.000000 4.772589    0
    #> [3,] 0.000000 0.000000    1
    # compare to analytical solution
    y^(x-1)*(x*log(y) + 1)
    #>          [,1]
    #> [1,] 2.098612
    #> [2,] 4.772589
    #> [3,] 1.000000
    x <- matrix(1:3, ncol = 1)
    y <- matrix(3:1, ncol = 1)
    dvdx(dxdy(madness(x), madness(y), function(x, y) sum(sin(x*y))))
    #>           [,1]     [,2]      [,3]
    #> [1,] -1.413352 0.000000  0.000000
    #> [2,]  0.000000 2.373566  0.000000
    #> [3,]  0.000000 0.000000 -1.413353
    # compare to analytical solution
    cos(x*y) - x*y*sin(x*y)
    #>           [,1]
    #> [1,] -1.413353
    #> [2,]  2.373566
    #> [3,] -1.413353
    

    Note that the madness derivatives are the partials of each x (columns) with respect to the y partials of each xy (rows).