R言語で勾配とヘシアンを出力する関数を作成する

Rには、微分して値を出力したり、ヘシアンを出力する既存関数としてderiv関数というものがあるみたいなのですが、
なんかネットで調べて試してみても、何か使い方がよくわかんなかったので、まどろっこしくなって作りました。

関数形と微分したい変数名、変数に代入したい値を入力すると勾配とヘシアンを出力する関数です。

勾配とヘシアンとは
今回使用するRの既存関数
勾配とヘシアンの表現型を表示する関数の作成
- 勾配
- ヘシアン
勾配とヘシアンの値を出力する関数
- 勾配値を出力する関数
- ヘシアン値を出力する関数

勾配とヘシアンとは

とりあえず、勾配とヘシアンの説明を簡単にしておきます。

ある多変量関数 $f(X)$ があるとします。ここで $X$ は $n$ 次元ベクトルとします。
この時、 $X$ に関する勾配とは次の $G$ を指します。

$G = \frac{\partial f(X)}{\partial X}= ( \frac{\partial f(X)}{\partial X_1},\cdots , \frac{\partial f(X)}{\partial X_n})'$

ただし、 $X_i$ とは $X$ の $i$ 番目の要素とします。

また、 $X$ に関するヘシアン $H$ とは次のような行列を指します。

$H = \begin{pmatrix} \frac{\partial^2 f(X)}{\partial X_1^2} & \frac{\partial^2 f(X)}{\partial X_1 \partial X_2} & \cdots & \frac{\partial^2 f(X)}{\partial X_1 \partial X_n} \\ \frac{\partial^2 f(X)}{\partial X_2 \partial X_1} & \frac{\partial^2 f(X)}{\partial X_2^2} & & \vdots\\ \vdots & &\ddots & \vdots\\ \frac{\partial^2 f(X)}{\partial X_n \partial X_1} &\cdots& & \frac{\partial^2 f(X)}{\partial X_n^2} \\ \end{pmatrix}$

要は1回微分したものを並べたベクトルと、それぞれの組み合わせで二回微分したものを並べた行列です。

これを作成した後適当な値を $X$ に代入することで、数値の入った行列やベクトルを得ます。
この勾配やヘシアンは特に最適化において力を発揮します。

今回使用するRの既存関数

expression関数

Rで、値を代入していない変数を扱う時はexpressionを利用することになります。
expression関数に入れた計算式は計算されることなく(これを評価されないと呼ぶ)、計算式のまま保持されます。
具体例を見てみます。

> expression(X^2)
expression(X^2)
> expression(X^2+1)
expression(X^2 + 1)
> expression(sqrt(X))
expression(sqrt(X))
> siki=expression(sqrt(X)+1)
> siki
expression(sqrt(X) + 1)

eval関数

expression関数に入った値を評価するにはeval関数を利用します。これも具体例を見た方が早いので見てみます。

> siki=expression(sqrt(X)+1)
> siki
expression(sqrt(X) + 1)
> eval(siki)
 eval(siki) でエラー:  オブジェクト 'X' がありません 
> X=4
> eval(siki)
[1] 3

グローバル環境にXの値が無い場合、式は評価不可能なのでエラーを吐きます。
Xの値が入っている場合は式が評価され、値が出力します。

assign関数

文字列(character型)を変数に変換し、そこに値を代入する関数です。これも具体例を見たらすぐわかるかと。

> xx
 エラー:  オブジェクト 'xx' がありません 
> assign("xx",5)
> xx
[1] 5

これも地味に便利です。

勾配とヘシアンの表現型を表示する関数の作成

ある関数が与えられた時に、勾配やヘシアンに関して値を代入する前の関数形を知りたいという需要もあるかもしれないので、まずそういったものを出力する関数を作成します。

勾配

grad=function(expr,par=c()){
	###expr:表現式(expression)
	###par:微分する変数(character)
	G=list()
	g_len=length(par)
	for(i in 1:g_len){
		x=par[i]
		G[[i]]=D(expr,x)
	}
return(G)
}

grad関数は、exprにexpression、parに微分したい変数を突っ込むことで、微分された状態をリストで返します。
具体例としてはこんな感じ

> siki=expression(x1^2+x2^3)
> grad(expr=siki,par=c("x1","x2"))
[[1]]
2 * x1

[[2]]
3 * x2^2

ヘシアン

Hes=function(expr,par=c()){
	###expr:表現式(expression)
	###par:微分する変数(character)
		res=matrix(,length(par),length(par))
		G=list()
		GG=list()
		g_len=length(par)
		
		##1階微分後の表現型のリストGを作成
		for(i in 1:g_len){
			x=par[i]
			G[[i]]=D(expr,x)
		}
		##2階微分後の表現型のリストGGを作成
		for(i in 1:g_len){
			GG_inner=list()
			for(j in 1:g_len){
				GG_inner[[j]]=D(G[[i]],par[j])
			}
			GG[[i]]=GG_inner
		}
	return(GG)
}

Hes関数も先ほどのgradと同様の引数をもってヘシアンの関数形を返します。matrix型で表示することは出来ないので、リストでi行j列目を表現しています。

> siki=expression(x1^2+x2^3)
> Hes(expr=siki,par=c("x1","x2"))
[[1]]
[[1]][[1]]
[1] 2

[[1]][[2]]
[1] 0


[[2]]
[[2]][[1]]
[1] 0

[[2]][[2]]
3 * (2 * x2)

勾配とヘシアンの値を出力する関数

先ほどの関数も応用しながら、変数に値を代入して、値が出力されるとこまでの関数を作ってみようと思います。

勾配値を出力する関数

cul_grad=function(expr,par=c(),sub=c(),var=c()){
	###expr:表現式(expression)
	###par:微分する変数(character)
	###sub:変数に代入する値(real)
	###var:表現式内の全ての変数(character)
	if(length(var)==length(sub)){
		res=c()
		G=list()
		g_len=length(par)
		
		##微分後の表現型のリストGを作成
		for(i in 1:g_len){
			x=par[i]
			G[[i]]=D(expr,x)
		}
		##各変数に値を代入していく
		for(ii in sub){
			num=which(sub==ii)
			assign(var[num],ii)
		}
		##リストGを評価したものをベクトルとして表示
		for(i in 1:g_len){
			res[i]=eval(G[[i]])
		}
		return(res)
	}
	else{
	print("Error,\"length(var) == length(sub)\" must be TRUE.")
	}
}

exprにはexpression関数で書いた関数の表現形を入れます。
parには微分したい変数をcharacter型で入力します。
subには変数に代入する値を入れます。次に説明するvarと順番を対応させる必要があります。
varには、exprで使用している全ての変数を識別するために入力します。

具体例

> siki=expression(x1^2+x2^3)
> cul_grad(expr=siki,par=c("x1","x2"),sub=c(2,3),var=c("x1","x2"))
[1]  4 27

ヘシアン値を出力する関数

cul_Hes=function(expr,par=c(),sub=c(),var=c()){
	###expr:表現式(expression)
	###par:微分する変数(character)
	###sub:変数に代入する値(real)
	###var:表現式内の全ての変数(character)
	if(length(var)==length(sub)){
		res=matrix(,length(par),length(par))
		G=list()
		GG=list()
		g_len=length(par)
		
		##1階微分後の表現型のリストGを作成
		for(i in 1:g_len){
			x=par[i]
			G[[i]]=D(expr,x)
		}
		##2階微分後の表現型のリストGGを作成
		for(i in 1:g_len){
			GG_inner=list()
			for(j in 1:g_len){
				GG_inner[[j]]=D(G[[i]],par[j])
			}
			GG[[i]]=GG_inner
		}
		##各変数に値を代入していく
		for(ii in sub){
			num=which(sub==ii)
			assign(var[num],ii)
		}
		##リストGを評価したものをベクトルとして表示
		for(i in 1:g_len){
			for(j in 1:g_len){
				res[i,j]=eval(GG[[i]][[j]])
			}
		}
		return(res)
	}
	else{
	print("Error,\"length(var) == length(sub)\" must be TRUE.")
	}
}

先ほどのcul_grad関数と同じように利用します

> siki=expression(x1^3+x2^4)
> cul_Hes(expr=siki,par=c("x1","x2"),sub=c(2,3),var=c("x1","x2"))
     [,1] [,2]
[1,]   12    0
[2,]    0  108