Compose two functions created using either sourceCppAD or %.% itself

Description

Returns a function with a matrix X as input which computes the value of (f \circ g)(x) = f(g(x)). Note that the order of the composition is such that g is applied to X first. f is then applied to the result of g. The returned function is compatible with both J and H, and if J and H are applied to a function produced by composition, the resulting Jacobian or Hessian matrices are constructed from the Jacobians and Hessians of f and g using a combinatorical form of Faa di Bruno's formula (Hardy 2006). The functions f and g must both be functions of a single argument. Note that a function of multiple arguments can be Curried into a function with a single argument. The R package functional provides the method Curry which is convenient for this purpose.

Usage

f %.% g

Arguments

Value

A function which computes the composition f and g

References

Hardy M (2006). “Combinatorics of Partial Derivatives.” Electr. J. Comb., 13(1). https://dblp.uni-trier.de/db/journals/combinatorics/combinatorics13.html#Hardy06.

Ma T (2009). “Higher Chain Formula proved by Combinatorics.” Electr. J. Comb., 16.

Examples

library(RcppEigenAD)
# define a function to calculate the eigen vectors of a matrix
f<-sourceCppAD('
ADmat f(const ADmat& X)
{
  Eigen::EigenSolver<ADmat> es(X);
  return es.pseudoEigenvectors();

}
')
# define function to calculate the inverse of a matrix
g<-sourceCppAD('
ADmat g(const ADmat& X)
{
  return X.inverse(); 
}
')
# compose f and g to produce a functions to calculate the eigenvectors of the inverse of a matrix
h<-f%.%g  #  h = f o g
x<-matrix(c(1,2,3,4),2,2)
x<-x%*%t(x) # positive definite matrix
J(h)(x) # Jacobian of h = f o g
H(h)(x) # Stacked Hessians of h = f o g
# redefine h as a function to directly calculate the eigenvectors of the inverse of a matrix
h<-sourceCppAD('
ADmat h(const ADmat& X)
{
  Eigen::EigenSolver<ADmat> es(X.inverse());
  return es.pseudoEigenvectors();
}
')
# calculate the Jacobian and Hessian of h to compare with previous result
J(h)(x) # Jacobian of h = f o g
H(h)(x) # Stacked Hessians of h = f o g

Construct a function to calculate the Hessian of a function.

Description

Constructs a function to calculate the Hessian of a function produced either using sourceCppAD or the composition operator %.% . The returned function has the same argument signature as f but returns a matrix representing the blocked Hessians of f evaluated at the functions arguments. The partial derivatives are formed with respect to the argument specified when f was created with sourceCppAD.

In what follows it is assumed that f has a single matrix argument (the one with which f is differentiated with respect to). When this is not the case, the other arguments will be considered constant at the point the Hessian is evaluated at. Consequently, the structure of the output of the function produced by H is unchanged by the additional arguments. The blocked Hessian {\bf H} is organised as follows.

If f:{\bf R}^{n} \rightarrow {\bf R}^{m} where {\bf Y}_{n_{Y} \times m_{Y}} = f({\bf X}_{n_{X} \times m_{X}}) and n=n_{X}m_{X}, m=n_{Y}m_{Y} then by numbering the elements of the matrices row-wise so that,

{\bf Y} = \left[ \begin{array}{ccc} y_{1} & \dots & y_{m_{Y}} \\ y_{m_{Y}+1} & \dots & y_{2m_{Y}} \\ \vdots & \ddots & \vdots \\ y_{(n_{Y}-1)m_{Y}+1} & \dots & y_{n_{Y}m_{Y}} \end{array} \right]

and

{\bf X} = \left[ \begin{array}{ccc} x_{1} & \dots & x_{m_{X}} \\ x_{m_{X}+1} & \dots & x_{2m_{X}} \\ \vdots & \ddots & \vdots \\ x_{(n_{X}-1)m_{X}+1} & \dots & x_{n_{X}m_{X}} \end{array} \right]

then the n \times nm blocked Hessian matrix {\bf H} is structured as

{\bf H} = \left[ \begin{array}{cccc} {\bf H}_{1} & {\bf H}_{2} & \dots & {\bf H}_{m} \end{array} \right]

with {\bf H}_{k} the Hessian matrix for y_{k} i.e.,

\begingroup \renewcommand*{\arraystretch}{1.5} {\bf H}_{k} = \left[ \begin{array}{ccc} \frac{\partial^2 y_{k}}{\partial x_1 \partial x_1} & \dots & \frac{\partial^2 y_{k}}{\partial x_1 \partial x_n} \\ \vdots & \ddots & \vdots \\ \frac{\partial^2 y_{k}}{\partial x_n \partial x_1} & \dots & \frac{\partial^2 y_{k}}{\partial x_n \partial x_n} \end{array} \right] \endgroup

Usage

H(f)

Arguments

Value

A function which computes the Hessian of the function.

Examples

library(RcppEigenAD)
# define f as the eigen vectors of its argument x
# calculated using the Eigen library
f<-sourceCppAD('
ADmat f(const ADmat& X)
{
   Eigen::EigenSolver<ADmat> es(X);
   return es.pseudoEigenvectors();
}
')
Hf<-H(f)
X<-matrix(c(1,2,3,4),2,2)
Hmat<-Hf(X) 
Hmat # the Hessian matrices of second derivatives stacked column wise
Hmat[3,9] # the second derivative of f(X)[2,1] with respect to X[2,1] and X[1,2]

Construct a function to calculate the Jacobian of a function.

Description

Constructs a function to calculate the Jacobian of a function produced either using sourceCppAD or the composition operator %.%. The returned function has the same argument signature as f but returns a matrix representing the Jacobian of f evaluated at the functions arguments. The partial derivatives are formed with respect to the arguments specified when f was created with sourceCppAD.

In what follows it is assumed that f has a single matrix argument (the one with which f is differentiated with respect to). When this is not the case, the other arguments will be considered constant at the point the Jacobian is evaluated at. Consequently, the structure of the output of the function produced by J is unchanged by the additional arguments. The Jacobian matrix {\bf J} is organised as follows.

If f:{\bf R}^{n} \rightarrow {\bf R}^{m} where {\bf Y}_{n_{Y} \times m_{Y}} = f({\bf X}_{n_{X} \times m_{X}}) and n=n_{X}m_{X}, m=n_{Y}m_{Y} then by numbering the elements of the matrices row-wise so that,

{\bf Y} = \left[ \begin{array}{ccc} y_{1} & \dots & y_{m_{Y}} \\ y_{m_{Y}+1} & \dots & y_{2m_{Y}} \\ \vdots & \ddots & \vdots \\ y_{(n_{Y}-1)m_{Y}+1} & \dots & y_{n_{Y}m_{Y}} \end{array} \right]

and

{\bf X} = \left[ \begin{array}{ccc} x_{1} & \dots & x_{m_{X}} \\ x_{m_{X}+1} & \dots & x_{2m_{X}} \\ \vdots & \ddots & \vdots \\ x_{(n_{X}-1)m_{X}+1} & \dots & x_{n_{X}m_{X}} \end{array} \right]

then the m \times n Jacobian matrix is given by

\begingroup \renewcommand*{\arraystretch}{1.5} {\bf J} = \left[ \begin{array}{ccc} \frac{\partial y_{1}}{\partial x_1} & \dots & \frac{\partial y_{1}}{\partial x_n} \\ \vdots & \ddots & \vdots \\ \frac{\partial y_{m}}{\partial x_1} & \dots & \frac{\partial y_{m}}{\partial x_n} \end{array} \right] \endgroup

Usage

J(f)

Arguments

Value

A function which computes the Jacobian of the function

Examples

# define f as the eigen vectors of its argument X
# calculated using the Eigen library
library(RcppEigenAD)
f<-sourceCppAD('
ADmat f(const ADmat& X)
{
   Eigen::EigenSolver<ADmat> es(X);
   return es.pseudoEigenvectors();
}
')
Jf<-J(f)
X<-matrix(c(1,2,3,4),2,2)
Jmat<-Jf(X)
Jmat # the Jacobian matrix of first derivatives
Jmat[2,3] # the derivative of f(X)[1,2] with respect to X[2,1]

Construct a function to calculate the Jacobian of a function.

Description

Constructs an R function that wraps a call to a C++ function and compiles the C++ code used to define the function. The code must define a function having arguments of the type const ADmat& and a return type of ADmat. The ADmat type is a type definition for a dynamically sized matrix from the Eigen linear algebra library. (Guennebaud G, Jacob B and others 2010). The function can employ methods from the Eigen library, such as all of the standard operators from linear algebra, and a wide range of functions for calculating inverses, decompositions and so on. See Guennebaud G, Jacob B and others (2010) for further details.

The returned function can then be used as an argument to J or H which provide functions that apply algorithmic differentiation to calculate the Jacobian or Hessian matrices. The functions produced by J or H evaluate the partial derivatives with respect to the elements of the argument located at the position in the functions argument list that is specified by the wrt argument of sourceCppAD.

Usage

sourceCppAD(code=NULL,file=NULL,wrt=1,output="method")

Arguments

Value

A function which invokes the compiled code or the c++ code that wraps the users function for use with the algorithmic differentiation libraries.

References

Guennebaud G, Jacob B, others (2010). “Eigen v3.” http://eigen.tuxfamily.org.

Examples

library(RcppEigenAD)
# define a function to calculate sin(cos(x)) where x is a matrix
f<-sourceCppAD('
ADmat f(const ADmat& X)
{
  return X.array().cos().sin();
}
')
x<-matrix(c(1,2,3,4),2,2)
# call it
f(x)