Salad is a package for Automatic Differentiation. Its name could stand for Simple And Light Automatic Differentiation, but it stands well by itself (eat salad, salad is good).
Lots of efforts have been done to allow re-using with salad functions written without automatic differentiation in mind.
We are going to illustrate some of the things salad can do by a series of short examples.
Consider for example the following function, \(f(x) = \sin(x^2)\) :
function(x) sin(x**2) f1 <-
The value of its derivative for a given value of \(x\). We just need to apply the function to a dual number created by salad::dual
:
dual(pi)
x <- f1(x)
y <- y
## [1] -0.4303012
## [has derivatives in x1]
d(y)
## $x1
## [1] -5.671739
And it works with vectors too:
dual(c(0, 1, sqrt(pi)))
x <- f1(x)
y <- y
## [1] 0.000000e+00 8.414710e-01 5.665539e-16
## [has derivatives in x1 x2 x3]
# get value and derivative
value(y)
## [1] 0.000000e+00 8.414710e-01 5.665539e-16
d(y)
## $x1
## [1] 0 0 0
##
## $x2
## [1] 0.000000 1.080605 0.000000
##
## $x3
## [1] 0.000000 0.000000 -3.544908
A second example will be given by matrix arithmetic. First define a dual object from a matrix.
dual( matrix( c(1, 2, 4, 7), 2, 2))
x <- x
## [,1] [,2]
## [1,] 1 4
## [2,] 2 7
## [has derivatives in x1.1 x2.1 x1.2 x2.2]
The default behavior of dual
is to name the variables x1.1
, x1.2
, etc.
varnames(x)
## [1] "x1.1" "x2.1" "x1.2" "x2.2"
# derivative along x1.1
d(x, "x1.1")
## $x1.1
## [,1] [,2]
## [1,] 1 0
## [2,] 0 0
Methods have been defined in salad to handle matrix product:
x %*% x
y <- y
## [,1] [,2]
## [1,] 9 32
## [2,] 16 57
## [has derivatives in x1.1 x2.1 x1.2 x2.2]
d(y, "x1.1")
## $x1.1
## [,1] [,2]
## [1,] 2 4
## [2,] 2 0
The determinant can be computed as well:
det(x)
## [1] -1
## [has derivatives in x1.1 x2.1 x1.2 x2.2]
d(det(x))
## $x1.1
## [1] 7
##
## $x2.1
## [1] -4
##
## $x1.2
## [1] -2
##
## $x2.2
## [1] 1
And the inverse:
solve(x)
z <- z
## [,1] [,2]
## [1,] -7 4
## [2,] 2 -1
## [has derivatives in x1.1 x2.1 x1.2 x2.2]
d(z, "x1.1")
## $x1.1
## [,1] [,2]
## [1,] -49 28
## [2,] 14 -8
ifelse
, apply
etc.As a last example, consider this function, which does nothing very interesting, except using ifelse
, cbind
, apply
, and crossprod
:
function(x) {
f2 <- x**(1:2)
a <- ifelse(a > 1, sin(a), 1 - a)
b <- crossprod( cbind(a,b) )
C <-apply(C, 2, function(x) sum(x^2))
}
It works ok:
# creating a dual number for x = 0.2
dual(0.2)
x <- f2(x)
y <- y
## a b
## 0.04109312 2.47795712
## [has derivatives in x1]
# get value and the derivative
value(y)
## a b
## 0.04109312 2.47795712
d(y)
## $x1
## a b
## 0.4200448 -7.0116352
You need to be aware of the following limitations of salad.
Checking the variable along which the derivatives are defined would have slowed the computations an awful lot. The consequence is that if you don’t take care of that yourselves, it may give inconsistent results.
Let’s define to dual numbers with derivatives along x
and y
.
dual(c(1,2), dx = list("x" = c(1,1)))
a <- dual(c(2,1), dx = list("y" = c(2,1))) b <-
It would be neat if a + b
had a derivative along x
and one along y
. It doesn’t.
+ b a
## [1] 3 3
## [has derivatives in x]
d(a + b)
## $x
## [1] 3 2
A simple way to deal with this is to define a
and b
with the appropriate list of derivatives.
dual(c(1,2), dx = list("x" = c(1,1), "y" = c(0,0)))
a <- dual(c(2,1), dx = list("x" = c(0,0), "y" = c(2,1))) b <-
It now works as intended:
+ b a
## [1] 3 3
## [has derivatives in x y]
d(a + b)
## $x
## [1] 1 1
##
## $y
## [1] 2 1
My prefered solution is to first define a dual vector with the two variables x
and y
, this can be done like this:
dual( c(1,1), varnames = c("x", "y"))
v <- v
## [1] 1 1
## [has derivatives in x y]
d(v)
## $x
## [1] 1 0
##
## $y
## [1] 0 1
Equivalently, one could define v
as dual(c(x = 1, y = 1))
. Once this is done, you can proceed as follows. First isolate the variables x
and y
:
v[1]
x <- v[2] y <-
Then create a
with the wanted derivatives:
c(x, x+1)
a <-d(a)
## $x
## [1] 1 1
##
## $y
## [1] 0 0
Same thing for b
:
c(2*y, y)
b <-d(b)
## $x
## [1] 0 0
##
## $y
## [1] 2 1
And now eveything works ok.
+ b a
## [1] 3 3
## [has derivatives in x y]
d(a + b)
## $x
## [1] 1 1
##
## $y
## [1] 2 1
As a general advice, a computation should begin with a single dual()
call, creating all the variables along which one needs to derive. This should avoid all problems related to this limitation of salad.
as.vector
and as.matrix
The functions as.vector
and as.matrix
return base (constant) vector and matrix objects.
dual( matrix( c(1, 2, 4, 7), 2, 2))
x <-as.vector(x)
## Warning in as.vector(x, mode): Dropping derivatives in as.vector. See ?salad to change this behaviour
## [1] 1 2 4 7
This behavior can be changed using salad(drop.derivatives = FALSE)
, but this may break some things. The prefered way to changed the shape of an object is to use `dim<-’ :
dim(x) <- NULL
x
## [1] 1 2 4 7
## [has derivatives in x1.1 x2.1 x1.2 x2.2]
You may need to rewrite partially some functions due to this behavior.
abs
, max
, min
The derivative of abs
has been defined to sign
. It might not be a good idea:
dual(0) + c(-1,0,1)
x <- x
## [1] -1 0 1
## [has derivatives in x1]
d(x)
## $x1
## [1] 1 1 1
abs(x)
## [1] 1 0 1
## [has derivatives in x1]
d(abs(x))
## $x1
## [1] -1 0 1
Also, the derivative of max
relies on which.max
: it works well when there are no ties, that is, when it is well defined. In the presence of ties, it is false.
When there are no ties, the result is correct:
max( dual(c(1, 2)) )
y <- y
## [1] 2
## [has derivatives in x1 x2]
d(y)
## $x1
## [1] 0
##
## $x2
## [1] 1
But in presence of ties, the derivatives should be undefined.
max( dual(c(2, 2)) )
y <- y
## [1] 2
## [has derivatives in x1 x2]
d(y)
## $x1
## [1] 1
##
## $x2
## [1] 0
It must be clear from the previous examples that salad can handle both vector and matrices. In addition to the simple arithmetic operations, most mathematical functions have been defined in salad
: trigonometic functions, hyperbolic trigonometric functions, etc (see the manual for an exhaustive list of functions of method). Many functions such as ifelse
, apply
, outer
, etc, have been defined.
In addition to the simple matrix arithmetic, salad can also handle det
and solve
. Currently, matrix decomposition with eigen
and qr
are not implemented (but this may change in the near future).
Assume you’re using salad to compute the derivative of a quadratic function:
function(x) x**2 + x + 1
f <- dual(4)
x <-f(x)
## [1] 21
## [has derivatives in x1]
d(f(x))
## $x1
## [1] 9
This works, but in the other hand, you know that this derivative is 2*x + 1
. You can tell salad about it with dualFun1
.
dualFun1(f, \(x) 2*x + 1)
f1 <-f1(x)
## [1] 21
## [has derivatives in x1]
d(f1(x))
## $x1
## [1] 9
This allows you to use special functions that salad can’t handle; moreover, even for simple functions like this one, it saves some time:
system.time( for(i in 1:500) f(x) )
## user system elapsed
## 0.018 0.000 0.019
system.time( for(i in 1:500) f1(x) )
## user system elapsed
## 0.004 0.000 0.004
It can thus be useful to define the derivatives of the some of the functions you are using in this way.
You may e-mail the author if for bug reports, feature requests, or contributions. The source of the package is on github.