R中计算jacobian和hessian矩阵有两个包，一个是rootSolve，一个是numDeriv，前者好像只能计算jacobian，我们这里只讲numDeriv的计算过程。

Jacobian

例1. single function. $f(x,y)=3x^2y$
我们手动计算结果为$J_{f}(x,y)=\left[\frac{\partial f1}{\partial x},\frac{\partial f1}{\partial y}\right]=[6xy,3x^2]$. R语言代码如下：

f = function(x){
3*x[1]^2*x[2]
}
library(numDeriv)
jm = jacobian(func=f,x=c(1,2))

结果如下：

例2. scalar function. $f(x,y)=\begin{bmatrix}x^2y\\5x+sin(y)\end{bmatrix}$
我们手动计算结果如下：

$f1(x,y) = x^2y$

$f2(x,y) = 5x+sin(y)$

$J_{f}(x,y)=\begin{bmatrix}\frac{\partial f1}{\partial x} & \frac{\partial f1}{\partial y}\\\frac{\partial f2}{\partial x} & \frac{\partial f2}{\partial y}\end{bmatrix}=\begin{bmatrix}2xy & x^2\\5 & cos(y) \end{bmatrix}$
R语言代码如下：

f = function(x,d=1){
y1 = d*x[1]^2*x[2]
y2 = 5*x[1] + sin(x[2])
c(y1,y2)
}
library(numDeriv)
jm = jacobian(func=f,x=c(1,2),d=1)

结果如下

Hessian

例1. single function. $f(x,y)=3x^2y$
我们手动计算结果为:

$J_{f}(x,y)=\left[\frac{\partial f1}{\partial x},\frac{\partial f1}{\partial y}\right]=[6xy,3x^2]$

$H_{f}(x,y)=\begin{bmatrix}\frac{\partial J_{f1}}{\partial x} & \frac{\partial J_{f1}}{\partial y}\\\frac{\partial J_{f2}}{\partial x} & \frac{\partial J_{f2}}{\partial y} \end{bmatrix}=\begin{bmatrix}6y & 6x\\6x & 0 \end{bmatrix}$
R语言代码如下：

f = function(x){
3*x[1]^2*x[2]
}
library(numDeriv)
hm=hessian(func=f,x=c(1,2))

注意：以上为什么没有scalar function的hessian计算，因为hessian是两次求导的过程，当f(x,y)为single function时候，第一次求导，得到的jacobian matrix为scalar function，而scalar function可以继续进行jacobian计算。当f(x,y)为scalar function时候，第一次求导，得到的jacobian matrix就不是scalar function，是一个matrix function，而matrix function不能利用jacobian计算。