1 はじめに

私たちは最新の注意を払ってこのパッケージを作成してきましたが、 通常のソフトウエアと同様に、結果の完全な正確性を保証する訳ではありません。

出力される記号演算式が、 例えば二つの行列の積のように簡単な場合、数学的知識を活用すれば 目視で結果の正確性を判断することが可能です。 しかし、記号演算式が複雑な場合や、 行列式の演算やスカラー関数の行列による微分の様に、数学的知識を活用しても 簡単に正解が導けない場合には、この方法を用いることが出来ません。

そこで、この文書では、lazy.symbolic パッケージの記号的出力を数値的に確かめる方法 について説明します。

2 パッケージの宣言

options(width=200)
library(lazy.symbolic)
modify_math_operators()

3 結果の数値的検証法

いま、expr を文字列を要素として持つオブジェクトとします。 例えば

 expr <- "A"
 expr <- "Det(A)"
 expr <- matrix( c("a11","a21","a12","a22"), 2,2 )
 expr <- mD( "tr(A%*%X)","X" )

の様なものです。

lazy.symbolic パッケージには文字列を要素として持つオブジェクト expr の中に 現れる行列やベクトルと見なされる文字(以下、行列)やその要素(行列の要素)に、 指定された値を「代入」した結果を評価する関数 Eval が用意されています。

この関数は以下の様な仕様になっており

     Eval( expr, …, values, fullsymb=NULL)

ここで … は

     pat11=pat21, pat12=pat22, pat13=pat23

の様な形をしており これが意味するところは expr の中に現れる pat11 という行列に pat21 の値を「代入」することを意味しています。 したがって pat21 は行列として定義されていなければなりません。 pat12 と pat22、pat13 と pat23 に関しても同様です。

同じことが、

     values=“pat11=pat21; pat12=pat22; pat13=pat23”

という形でも行えますが、この場合は、右辺に現れる pat21 等が Eval が実行される 環境の中で定義されている行列であることが必要です。 また、pat21 等に関数呼び出しを書くことも出来ません。

fullsymb オプションは expr に含まれる行列が 記号行列か完全記号行列であるかを区別するためのもので、 expr が、その要素が文字列として既に定義されている記号行列 からなる式の場合にはデフォールトの NULL を用い、 expr の中に現れる文字が要素は未定義であるが行列と 見なすべきもの、すなわち完全記号行列である場合にには 1 と指定する必要があります。

例えば、以下の例は A という完全記号行列に demomat(2,2) という数値行列を 「代入」し、それを評価するためのものです。

なお、demomat 関数は

     demomat( nrow, ncol, root, shape )

という形をしており、行数が nrow, 列数が ncol の行列を作成するための関数ですが、 root=“文字” を指定すると記号行列が定義されます。また、shape オプションで、 対称行列、上三角行列、対角行列等の指定が出来ます。

suppressWarnings(rm(A))
expr <- "A"
Eval(expr, A=demomat(2,2), fullsymb=1)
##      [,1] [,2]
## [1,]   11   12
## [2,]   21   22

これは

suppressWarnings(rm(A,An))
An=demomat(2,2)
expr <- "A"
Eval(expr, A=An, fullsymb=1)
##      [,1] [,2]
## [1,]   11   12
## [2,]   21   22

または

suppressWarnings(rm(A,An))
An=demomat(2,2)
values="A=An"
expr <- "A"
Eval(expr, values=values, fullsymb=1)
##      [,1] [,2]
## [1,]   11   12
## [2,]   21   22

としても同じことですが、

suppressWarnings(rm(A,An))
values="A=demomat(2,2)"
expr <- "A"
Eval(expr, values=values, fullsymb=1)
## 
##  The input expression cannot be evaluated.
##      A             
## [1,] "demomat(2,2)"

は意図した動作を行いません。

次に、expr が記号的行列からなる式の場合を説明します。 たとえば、行列 A が記号的行列として定義されている場合、その要素に数値的行列 An の対応する要素を代入して評価する場合は以下の様にします。 以下の例の Eval の中の A=An は、 expr の中に現れる A という 行列の要素(例えば a11)を、 対応する An 行列の要素(たとえば 11 )で置き換えて評価することを示しています。 したがって、A ならびに An という行列の双方が Eval 関数が呼ばれる環境に 存在している必要があります。

A <- demomat(2,2)
An <- demomat(2,2,root="a")
expr <- A
Eval( expr, A=An )
##      [,1] [,2]
## [1,]   11   12
## [2,]   21   22

これは、実は、以下の形の省略形となっていますが、このように指定する場合には A 行列の全ての値を指定する必要があります。

expr <- A
Eval( expr, a11=11, a12=12, a21=21, a22=22, fullsymb=1 )
##      [,1] [,2]
## [1,]   11   12
## [2,]   21   22

なお、expr は記号行列からなる式ですから以下の様なものでも構いません。

A <- demomat(3,3,root="a")
An <- demomat(3,3)
expr <- A+t(A)
resn <- Eval( expr, A=An )
printm( A, An , expr, resn )
##  A                 An             expr                          resn                
##      c1  c2  c3       c1 c2 c3            c1      c2      c3         [,1] [,2] [,3] 
##  r1 a11 a12 a13    r1 11 12 13    r1   2*a11 a12+a21 a13+a31    [1,]   22   33   44 
##  r2 a21 a22 a23    r2 21 22 23    r2 a21+a12   2*a22 a23+a32    [2,]   33   44   55 
##  r3 a31 a32 a33    r3 31 32 33    r3 a31+a13 a32+a23   2*a33    [3,]   44   55   66

もう少し複雑な場合は以下の様になります。 この例では A だけではなく expr の中に現れる B 行列の要素も対応する Bn 行列の 要素で置き換えられた後に評価されています。

A <- demomat(2,2,root="a")
B <- demomat(2,2,root="b")
An <- demomat(2,2)
Bn <- demomat(2,2)
expr <- A%*%B
resn <- Eval( expr, A=An, B=Bn )
printm( A, An, B, Bn, "/", expr, resn )
##  A             An          B             Bn       
##      c1  c2       c1 c2        c1  c2       c1 c2 
##  r1 a11 a12    r1 11 12    r1 b11 b12    r1 11 12 
##  r2 a21 a22    r2 21 22    r2 b21 b22    r2 21 22 
## 
##  expr                                  resn           
##                  c1              c2         [,1] [,2] 
##  r1 a11*b11+a12*b21 a11*b12+a12*b22    [1,]  373  396 
##  r2 a21*b11+a22*b21 a21*b12+a22*b22    [2,]  693  736

またこれは以下と同じです。

A <- demomat(2,2,root="a")
B <- demomat(2,2,root="b")
An <- demomat(2,2)
Bn <- demomat(2,2)
vals <- "A=An; B=Bn"
expr <- A%*%B
resn <- Eval( expr, values=vals )
printm( A, An, B, Bn, "/", expr, resn )
##  A             An          B             Bn       
##      c1  c2       c1 c2        c1  c2       c1 c2 
##  r1 a11 a12    r1 11 12    r1 b11 b12    r1 11 12 
##  r2 a21 a22    r2 21 22    r2 b21 b22    r2 21 22 
## 
##  expr                                  resn           
##                  c1              c2         [,1] [,2] 
##  r1 a11*b11+a12*b21 a11*b12+a12*b22    [1,]  373  396 
##  r2 a21*b11+a22*b21 a21*b12+a22*b22    [2,]  693  736

これまで見てきたように expr が記号行列の要素からなる式の場合には fullsymb オプションを指定する必要はありません。

4 数値的検証の具体例

4.1 四則演算の検証

以下に、簡単な記号的な四則演算が正しく処理されているかを調べる例を示します。

A <- demomat(2,2,root="a"); An <- demomat(2,2)
B <- demomat(2,2,root="b"); Bn <- demomat(2,2)
C <- demomat(2,2,root="c"); Cn <- demomat(2,2)
vals <- "A=An; B=Bn; C=Cn"
( expr <- A%*%(B+C) )
##    c1                            c2                           
## r1 "a11*(b11+c11)+a12*(b21+c21)" "a11*(b12+c12)+a12*(b22+c22)"
## r2 "a21*(b11+c11)+a22*(b21+c21)" "a21*(b12+c12)+a22*(b22+c22)"

この例では行列 A, B, C が全て 2 x 2 の行列の場合の A%*%(B+C) という 行列演算の結果が示されていますが、これが本当に正しい結果かどうかを知りたい場合を 考えます。数学的知識があればこれは正しい結果だということは自明ですが、 そうでない場合を想定しています。 その場合、R の数値演算の結果が正しいと仮定すれば、上記の expr に出てくる 行列の要素を、数値的行列の要素で置き換えて評価したものが、最初から数値的 行列を用いて R で計算した結果と一致すれば、 (この数値例に関してだけではありますが)結果は正しいということになります。

そのためには以下の様にします。

A <- demomat(2,2,root="a"); An <- demomat(2,2)
B <- demomat(2,2,root="b"); Bn <- demomat(2,2)
C <- demomat(2,2,root="c"); Cn <- demomat(2,2)
vals <- "A=An; B=Bn; C=Cn"
expr <- A%*%(B+C)
resc <- Eval(expr, values=vals )
resn <- An%*%(Bn+Cn)
printm( resc, resn )
##  resc              resn         
##       [,1] [,2]         c1   c2 
##  [1,]  746  792    r1  746  792 
##  [2,] 1386 1472    r2 1386 1472

この操作は、一般性を持たせると、以下の様に書くことが出来ます。 まず、完全記号行列 A, B, C を用いて expr を定義し、 それに数値行列 An, Bn, Cn を「代入」したものを resn とします。 これが R が出してくる正解です。 次に、expr に記号行列 Ac, Bc, Cc をを「代入」したものを resc としますが、 これは記号行列となり、その結果の正当性が問われているものです。 そして、その記号行列 resc に現れる行列の要素に対応する数値行列の 要素を「代入」して評価したものが rescn です。 したがって、resn と rescn が等しければ、lazy.symbolic パッケージで得られた resc の表現は正しい結果であると言えることになります。 (実際は、一つの数値例だけの結果ですから、 単に間違いであるとは言えないというだけですが。)

suppressWarnings(rm(A,B,C))
Ac <- demomat(2,2,root="a"); An <- demomat(2,2)
Bc <- demomat(2,2,root="b"); Bn <- demomat(2,2)
Cc <- demomat(2,2,root="c"); Cn <- demomat(2,2)
valsc <- "A=Ac; B=Bc; C=Cc"
valsn <- "A=An; B=Bn; C=Cn"
valscn <- "Ac=An; Bc=Bn; Cc=Cn"
expr <- "A%*%(B+C)"
resn <- Eval( expr, values=valsn, fullsymb=1 )
resc <- Eval( expr, values=valsc, fullsymb=1 )
## 
##  The input expression cannot be evaluated.
rescn <- Eval( resc, values=valscn )
printm(expr, resc, "/", resn,rescn)
##  expr          resc                                                       
##                                            c1                          c2 
##   A%*%(B+C)    r1 a11*(b11+c11)+a12*(b21+c21) a11*(b12+c12)+a12*(b22+c22) 
##                r2 a21*(b11+c11)+a22*(b21+c21) a21*(b12+c12)+a22*(b22+c22) 
## 
##  resn              rescn          
##       [,1] [,2]         [,1] [,2] 
##  [1,]  746  792    [1,]  746  792 
##  [2,] 1386 1472    [2,] 1386 1472

なお、上記の例では、expr の定義が完全記号行列であることを強調するために 敢えて A, B, C という行列を用いませんでしたが、expr は単なる文字列として 定義されていますから、その中に現れる行列は完全記号行列となります。 以下の様に書いた方が簡潔になります。

suppressWarnings(rm(A,B,C))
A <- demomat(2,2,root="a"); An <- demomat(2,2)
B <- demomat(2,2,root="b"); Bn <- demomat(2,2)
C <- demomat(2,2,root="c"); Cn <- demomat(2,2)
valsc <- "A=A; B=B; C=C"
valsn <- "A=An; B=Bn; C=Cn"
expr <- "A%*%(B+C)"
resn <- Eval( expr, values=valsn, fullsymb=1 )
resc <- Eval( expr, values=valsc, fullsymb=1 )
## 
##  The input expression cannot be evaluated.
rescn <- Eval( resc, values=valsn )
printm(expr, resc, "/", resn,rescn)
##  expr          resc                                                       
##                                            c1                          c2 
##   A%*%(B+C)    r1 a11*(b11+c11)+a12*(b21+c21) a11*(b12+c12)+a12*(b22+c22) 
##                r2 a21*(b11+c11)+a22*(b21+c21) a21*(b12+c12)+a22*(b22+c22) 
## 
##  resn              rescn          
##       [,1] [,2]         [,1] [,2] 
##  [1,]  746  792    [1,]  746  792 
##  [2,] 1386 1472    [2,] 1386 1472

以下ではこの表記を用います。

4.2 逆行列の検証

行列 A の逆行列を記号的に求め、それを数値的に求めた逆行列と比較します。 Inv は lazy.symbolic パッケージの逆行列を求める関数です。

A <- demomat(2,2,root="a"); An <- demomat(2,2)
valsc <- "A=A"
valsn <- "A=An"
expr <- "Inv(A)"
resn <- Eval( expr, values=valsn, fullsymb=1 )
resc <- Eval( expr, values=valsc, fullsymb=1 )
## 
##  The input expression cannot be evaluated.
rescn <- Eval( resc, values=valsn )
printm(expr, resc, "/", resn,rescn)
##  expr       resc                                                                         
##                                                          c1                           c2 
##   Inv(A)    r1 1/a11-(a12/a11)*(-a21/a11)/(a22-a21*a12/a11) -(a12/a11)/(a22-a21*a12/a11) 
##             r2                 (-a21/a11)/(a22-a21*a12/a11)          1/(a22-a21*a12/a11) 
## 
##  resn              rescn          
##       [,1] [,2]         [,1] [,2] 
##  [1,] -2.2  1.2    [1,] -2.2  1.2 
##  [2,]  2.1 -1.1    [2,]  2.1 -1.1

逆行列の検証 その2

行列 A の逆行列 Inv(A) を記号的に求め、A%Inv(A) が単位行列になるかを調べます。 demomat(3,3) はランク落ちをするので使えないため、乱数を用いています。

A <- demomat(3,3,root="a")
set.seed(1701)
An <- matrix(sample(0:9,3*3,replace=TRUE),3,3)
expr <- "A%*%Inv(A)"
resn <- Eval( expr, A=An, fullsymb=1 )
## 
##  The input expression cannot be evaluated.
resc <- Eval( expr, A=A, fullsymb=1 )
## 
##  The input expression cannot be evaluated.
rescn <- Eval( resc, A=An )
printm(expr, resc, "/", resn,rescn, fmt="12.10")
##  expr           
##   A%*%Inv(A)    
## 
##  resc                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                c1 
##  r1 a11*(1/a11-(a12/a11)*(-a21/a11)/(a22-a21*a12/a11)-(a13/a11-(a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11))*(-a31/a11-(a32-a31*a12/a11)*(-a21/a11)/(a22-a21*a12/a11))/(a33-a31*a13/a11-(a32-a31*a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11)))+a12*((-a21/a11)/(a22-a21*a12/a11)-((a23-a21*a13/a11)/(a22-a21*a12/a11))*(-a31/a11-(a32-a31*a12/a11)*(-a21/a11)/(a22-a21*a12/a11))/(a33-a31*a13/a11-(a32-a31*a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11)))+a13*(-a31/a11-(a32-a31*a12/a11)*(-a21/a11)/(a22-a21*a12/a11))/(a33-a31*a13/a11-(a32-a31*a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11)) 
##  r2 a21*(1/a11-(a12/a11)*(-a21/a11)/(a22-a21*a12/a11)-(a13/a11-(a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11))*(-a31/a11-(a32-a31*a12/a11)*(-a21/a11)/(a22-a21*a12/a11))/(a33-a31*a13/a11-(a32-a31*a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11)))+a22*((-a21/a11)/(a22-a21*a12/a11)-((a23-a21*a13/a11)/(a22-a21*a12/a11))*(-a31/a11-(a32-a31*a12/a11)*(-a21/a11)/(a22-a21*a12/a11))/(a33-a31*a13/a11-(a32-a31*a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11)))+a23*(-a31/a11-(a32-a31*a12/a11)*(-a21/a11)/(a22-a21*a12/a11))/(a33-a31*a13/a11-(a32-a31*a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11)) 
##  r3 a31*(1/a11-(a12/a11)*(-a21/a11)/(a22-a21*a12/a11)-(a13/a11-(a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11))*(-a31/a11-(a32-a31*a12/a11)*(-a21/a11)/(a22-a21*a12/a11))/(a33-a31*a13/a11-(a32-a31*a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11)))+a32*((-a21/a11)/(a22-a21*a12/a11)-((a23-a21*a13/a11)/(a22-a21*a12/a11))*(-a31/a11-(a32-a31*a12/a11)*(-a21/a11)/(a22-a21*a12/a11))/(a33-a31*a13/a11-(a32-a31*a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11)))+a33*(-a31/a11-(a32-a31*a12/a11)*(-a21/a11)/(a22-a21*a12/a11))/(a33-a31*a13/a11-(a32-a31*a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11)) 
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                              c2                                                                                   
##  r1 a11*(-(a12/a11)/(a22-a21*a12/a11)-(a13/a11-(a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11))*(-(a32-a31*a12/a11)/(a22-a21*a12/a11))/(a33-a31*a13/a11-(a32-a31*a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11)))+a12*(1/(a22-a21*a12/a11)-((a23-a21*a13/a11)/(a22-a21*a12/a11))*(-(a32-a31*a12/a11)/(a22-a21*a12/a11))/(a33-a31*a13/a11-(a32-a31*a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11)))+a13*(-(a32-a31*a12/a11)/(a22-a21*a12/a11))/(a33-a31*a13/a11-(a32-a31*a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11))                                                                                   
##  r2 a21*(-(a12/a11)/(a22-a21*a12/a11)-(a13/a11-(a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11))*(-(a32-a31*a12/a11)/(a22-a21*a12/a11))/(a33-a31*a13/a11-(a32-a31*a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11)))+a22*(1/(a22-a21*a12/a11)-((a23-a21*a13/a11)/(a22-a21*a12/a11))*(-(a32-a31*a12/a11)/(a22-a21*a12/a11))/(a33-a31*a13/a11-(a32-a31*a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11)))+a23*(-(a32-a31*a12/a11)/(a22-a21*a12/a11))/(a33-a31*a13/a11-(a32-a31*a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11))                                                                                   
##  r3 a31*(-(a12/a11)/(a22-a21*a12/a11)-(a13/a11-(a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11))*(-(a32-a31*a12/a11)/(a22-a21*a12/a11))/(a33-a31*a13/a11-(a32-a31*a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11)))+a32*(1/(a22-a21*a12/a11)-((a23-a21*a13/a11)/(a22-a21*a12/a11))*(-(a32-a31*a12/a11)/(a22-a21*a12/a11))/(a33-a31*a13/a11-(a32-a31*a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11)))+a33*(-(a32-a31*a12/a11)/(a22-a21*a12/a11))/(a33-a31*a13/a11-(a32-a31*a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11))                                                                                   
##                                                                                                                                                                                                                                                                                                                                     c3                                                                                                                                                                                                                                                            
##  r1 -a11*(a13/a11-(a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11))/(a33-a31*a13/a11-(a32-a31*a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11))-a12*((a23-a21*a13/a11)/(a22-a21*a12/a11))/(a33-a31*a13/a11-(a32-a31*a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11))+a13/(a33-a31*a13/a11-(a32-a31*a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11))                                                                                                                                                                                                                                                            
##  r2 -a21*(a13/a11-(a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11))/(a33-a31*a13/a11-(a32-a31*a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11))-a22*((a23-a21*a13/a11)/(a22-a21*a12/a11))/(a33-a31*a13/a11-(a32-a31*a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11))+a23/(a33-a31*a13/a11-(a32-a31*a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11))                                                                                                                                                                                                                                                            
##  r3 -a31*(a13/a11-(a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11))/(a33-a31*a13/a11-(a32-a31*a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11))-a32*((a23-a21*a13/a11)/(a22-a21*a12/a11))/(a33-a31*a13/a11-(a32-a31*a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11))+a33/(a33-a31*a13/a11-(a32-a31*a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11))                                                                                                                                                                                                                                                            
## 
##  resn                                              rescn                                          
##       [,1]          [,2]          [,3]                  [,1]          [,2]          [,3]          
##  [1,]  1.0000000000  0.0000000000  0.0000000000    [1,]  1.0000000000 -0.0000000000  0.0000000000 
##  [2,] -0.0000000000  1.0000000000  0.0000000000    [2,] -0.0000000000  1.0000000000  0.0000000000 
##  [3,] -0.0000000000  0.0000000000  1.0000000000    [3,] -0.0000000000  0.0000000000  1.0000000000

二つの結果は誤差範囲で一致していることが分かります。

4.3 行列式の検証

行列 A の逆行列を記号的に求め、それを数値的に求めた逆行列と比較します。 Det は lazy.symbolic パッケージの行列式を求める関数です。

A <- demomat(3,3,root="a")
set.seed(1701)
An <- matrix(sample(0:9,3*3,replace=TRUE),3,3)
valsc <- "A=A"
valsn <- "A=An"
expr <- "Det(A)"
resn <- Eval( expr, values=valsn, fullsymb=1 )
resc <- Eval( expr, values=valsc, fullsymb=1 )
## 
##  The input expression cannot be evaluated.
rescn <- Eval( resc, values=valsn )
printm(expr, resc, "/", resn,rescn)
##  expr       resc                                                                                               
##                  [,1]                                                                                          
##   Det(A)    [1,] (a33-a31*a13/a11-(a32-a31*a12/a11)*(a23-a21*a13/a11)/(a22-a21*a12/a11))*(a22-a21*a12/a11)*a11 
## 
##  resn         rescn     
##       [,1]         [,1] 
##  [1,]  -82    [1,]  -82

4.4 総和記号演算の検証

lazy.symbolic パッケージの mat2sum 関数は、完全記号行列で定義された行列の 式の i-j 要素が、総和記号で書くとどのように書けるかを示してくれる関数です。

mat2sum 関数は以下の様な形をしていますが、

     mat2sum( mexpr, simple=0,s sublist=NA )

mexpr は完全記号行列を使った行列の式、 simple は総和記号の和の範囲を簡易表示するか否かの指定です。 数値的検証が必要な場合は simple=0 を指定しなければなりません。

mat2sum の結果は以下の様な総和記号式で表されますが、

     s( expr(i), {i,1,n} )

これは expr(i) を添字 i を含む行列の要素の式として \[\begin{align} \sum_{i=1}^{n} \text{expr}(i) \end{align}\] を表しています。 すなわち、中括弧の中が添字の名前とその範囲となりますが、デフォールトの添字は s1, s2,… が使われます。これを別のものに変えるには、 sublist=c(“k”,“m”) の用に指定します。

以下の例では、行列 A, B, C がそれぞれ 2 x 2 の行列であった場合に A%%B%%C の i-j 要素がどのようになるかを示し、それを数値的に検証しています。 mexpr を数値的に評価したものが resm でこれが正解になります。 mexpr の i-j 要素を総和記号式で表現したものが sexpr で、それを数値的に 評価したのが、ressn になります。 ただし、総和記号式の数値的評価には Eval 関数ではなく sumEval 関数を利用します。

# summation
nvar1 <- 2; nvar2 <- 2
set.seed(1701)
An <- matrix(rnorm(nvar2*nvar1),nvar2,nvar1)
Bn <- matrix(rnorm(nvar2*nvar1),nvar2,nvar1)
Cn <- matrix(rnorm(nvar2*nvar1),nvar2,nvar1)
mexpr <- "A%*%B%*%C"
sexpr <- mat2sum(mexpr, simple=0, sublist=c("k","m"))
resm <- Eval(mexpr, A=An,B=Bn,C=Cn, fullsymb=1)
ressn <- sumEval(sexpr, A=An,B=Bn,C=Cn)
printm(mexpr,sexpr, "/", resm,ressn, fmt="10.8")
##  mexpr         sexpr                                                   
##   A%*%B%*%C     s(A[i,k]*s(B[k,m]*C[m,j],{m,1,ncol(B)}),{k,1,ncol(A)}) 
## 
##  resm                            ressn                      
##       [,1]        [,2]                       j1          j2 
##  [1,]  2.03854669 -1.24886685    i1  2.03854669 -1.24886685 
##  [2,]  2.15141735 -0.45954943    i2  2.15141735 -0.45954943
printm(max(abs(resm-ressn))/mean(abs(resm)), fmt="12.10")
##  max(abs(resm - ressn))/mean(abs(resm))  
##    0.0000000000

結果は誤差の範囲で一致しています。

もう一つの総和記号の例です。

mexpr <- "t(A)%*%Diag(B)%*%C"
sexpr <- mat2sum(mexpr, simple=0)
resm <- Eval(mexpr, A=An,B=Bn,C=Cn, fullsymb=1)
resc <- sumEval(sexpr, A=An,B=Bn,C=Cn)
printm(mexpr,sexpr, "/", resm,resc, fmt="10.8")
##  mexpr                  sexpr                                       
##   t(A)%*%Diag(B)%*%C     s(A[s1,i]*B[s1,s1]*C[s1,j],{s1,1,nrow(A)}) 
## 
##  resm                            resc                       
##       [,1]        [,2]                       j1          j2 
##  [1,]  1.47196175  0.94328262    i1  1.47196175  0.94328262 
##  [2,] -1.90385820 -0.06876642    i2 -1.90385820 -0.06876642
printm(max(abs(resm-resc))/mean(abs(resm)), fmt="12.10")
##  max(abs(resm - resc))/mean(abs(resm))  
##    0.0000000000

また、lazy.symbolic パッケージには sum2mat という総和記号式から、 それを i-j 要素として持つ行列の式を表示する関数が含まれますが、 それの評価は直接に行えないので、 sexpr として総和記号式を作成し、それを sum2mat に渡した結果 の行列の式から添字を除いた形で mexpr を作り、 上記の方法で評価することになります。

4.5 スカラー関数の行列による微分

lazy.symbolic パッケージの mD 関数は、 tr や det の様にスカラー値を返す行列の関数の、その関数に含まれる一つの行列に 関する一次偏微分を完全記号行列の式として返す関数です。

mD 関数は以下の様な形をしていますが、

     mD( expr, X_, sym=0 )

expr は完全記号行列を使った tr や det のような行列の式、 X_ は式の中に現れる行列の名称で、その行列に関する偏微分が計算されます。

なお、スカラー値を返す行列の関数 \(f\) を行列 \(X\) で偏微分した結果とは、 \(\frac{\displaystyle \partial \; f }{ \displaystyle \partial \; x_{ij} }\) を行列 \(X\) の形に並べたもののことです。

この場合、通常の R の関数では正解を計算できませんので、\(f\) を行列 \(X\) の要素で 数値的に偏微分したものを行列 \(X\) の形に並べたものを正解とする必要があります。 そのためには、以下の gradmn 関数を用います。

gradmn は、

     gradmn( expr, X_, …, values = NULL, sym = 0 )

という形をしていますが、 その第1引数 expr として与えられた完全記号行列の関数 \(f\) を その第2引数の X_ として与えられる変数の 与えられた数値(以下の例では Xn)付近で数値的に微分した値 を返す関数です。必要となる定数の値は、適宜、第3変数以降に

    変数名=変数の値

という形で与えます。もしくは Eval 関数同様に values オプションを利用します。

結果は第2引数の名前と同じサイズの行列になります。

数値微分には lazy.symbolic::JacobianMat 関数を使っています。

以下に tr(t(X)%%A%%X) を X で偏微分した結果の数値的検証例を示します。 dexpr が mD が出してくる完全記号行列としての偏微分の結果で、 それを数値的に評価したものが Ga です。 また、Gn は gradmn 関数で計算した数値微分の値です。 以下の例では Ga と Gn が誤差の範囲で一致していることが分かります。

expr <- "tr(t(X)%*%A%*%X)"
dexpr <- mD(expr,"X")
## P1: using product rule.....
## S2: tr(A%*%X)
## S2: tr(A%*%X)
printm(dexpr)
##  dexpr           
##   t(A)%*%X+A%*%X
nvar1 <- 2; nvar2 <- 2
set.seed(1701)
Xn <- matrix(rnorm(nvar1*nvar1),nvar1,nvar1)
An <- matrix(rnorm(nvar2*nvar1),nvar2,nvar1)
values <- "X=Xn; A=An"
Ga <- 0; Gn <- 9999
Ga <- Eval( dexpr, dexpr=dexpr, X=Xn, A=An, fullsymb=1 )
Gn <- gradmn( expr, X=Xn, A=An)
printm(Ga,Gn, fmt="12.10")
##  Ga                                  Gn                               
##       [,1]          [,2]                  [,1]          [,2]          
##  [1,]  4.6682279593 -8.2196764252    [1,]  4.6682279606 -8.2196764257 
##  [2,] -3.0594730085  2.1828383095    [2,] -3.0594730092  2.1828383101
printm(max(abs(Ga-Gn))/mean(abs(Gn)), fmt="12.10")
##  max(abs(Ga - Gn))/mean(abs(Gn))  
##    0.0000000003

微分の例をもう一つ。 tr(A%%inv(t(X)%%B%%X)%%C) を X で偏微分した結果を X が対称行列で無い場合と ある場合に分けて検証したものです。 微分する変数が対称行列の場合、mD と gradm に sym=1 を指定し Xn として 対称行列を用ます。

まず、X が対称では無い場合、

# inverse: tr(A%*%inv(t(X)%*%B%*%X)%*%C)
expr <- "tr(A%*%inv(t(X)%*%B%*%X)%*%C)"
dexpr <- mD(expr,"X")
## C1: using chain rule...
## S3.1 or 3.2: tr(A%*%inv(X))
## P1: using product rule.....
## S2: tr(A%*%X)
## S2: tr(A%*%X)
printm(dexpr)
##  dexpr                                                                                                               
##   -t(B)%*%X%*%t(Inv(t(X)%*%B%*%X)%*%C%*%A%*%Inv(t(X)%*%B%*%X))-B%*%X%*%Inv(t(X)%*%B%*%X)%*%C%*%A%*%Inv(t(X)%*%B%*%X)
nvar1 <- 2; nvar2 <- 3
set.seed(1701)
Xn <- matrix(rnorm(nvar1*nvar1),nvar1,nvar1)
Xn <- Xn+t(Xn)
An <- matrix(rnorm(nvar2*nvar1),nvar2,nvar1)
Bn <- matrix(rnorm(nvar1*nvar1),nvar1,nvar1)
Cn <- matrix(rnorm(nvar1*nvar2),nvar1,nvar2)
values <- "X=Xn; A=An; B=Bn; C=Cn"
Ga <- 0; Gn <- 9999
Ga <- Eval( dexpr, values=values, fullsymb=1 )
Gn <- gradmn( expr, values=values )
printm(Ga,Gn, fmt="12.10")
##  Ga                                  Gn                               
##       [,1]          [,2]                  [,1]          [,2]          
##  [1,] -0.1374688946 -0.2058239818    [1,] -0.1374688947 -0.2058239817 
##  [2,] -0.2464479337  0.8749811608    [2,] -0.2464479336  0.8749811604
printm(max(abs(Ga-Gn))/mean(abs(Gn)), fmt="12.10")
##  max(abs(Ga - Gn))/mean(abs(Gn))  
##    0.0000000011

次に対称な場合の例です。

# sym inverse: tr(A%*%inv(t(X)%*%B%*%X)%*%C)
expr <- "tr(A%*%inv(t(X)%*%B%*%X)%*%C)"
dexpr <- mD(expr,"X", sym=1)
## C1: using chain rule...
## S3.1 or 3.2: tr(A%*%inv(X))
## P1: using product rule.....
## S2: tr(A%*%X)
## S2: tr(A%*%X)
printm(dexpr)
##  dexpr                                                                                                                                                                                                                                                                                                                                                                          
##   -t(B)%*%X%*%t(Inv(t(X)%*%B%*%X)%*%C%*%A%*%Inv(t(X)%*%B%*%X))-B%*%X%*%Inv(t(X)%*%B%*%X)%*%C%*%A%*%Inv(t(X)%*%B%*%X)+-Inv(t(X)%*%B%*%X)%*%C%*%A%*%Inv(t(X)%*%B%*%X)%*%t(X)%*%B-t(Inv(t(X)%*%B%*%X))%*%t(A)%*%t(C)%*%t(Inv(t(X)%*%B%*%X))%*%t(X)%*%t(B)-Diag(-t(B)%*%X%*%t(Inv(t(X)%*%B%*%X)%*%C%*%A%*%Inv(t(X)%*%B%*%X))-B%*%X%*%Inv(t(X)%*%B%*%X)%*%C%*%A%*%Inv(t(X)%*%B%*%X))
values <- "X=Xn; A=An; B=Bn; C=Cn"
Ga <- 0; Gn <- 9999
Ga <- Eval( dexpr, values=values, fullsymb=1 )
Gn <- gradmn( expr, values=values, sym=1)
printm(Ga,Gn, fmt="12.10")
##  Ga                                  Gn                               
##       [,1]          [,2]                  [,1]          [,2]          
##  [1,] -0.1374688946 -0.4522719156    [1,] -0.1374688947 -0.4522719156 
##  [2,] -0.4522719156  0.8749811608    [2,] -0.4522719156  0.8749811604
printm(max(abs(Ga-Gn))/mean(abs(Gn)), fmt="12.10")
##  max(abs(Ga - Gn))/mean(abs(Gn))  
##    0.0000000008

両者とも誤差の範囲内で等しいと見ませます。

5 まとめ

lazy.symbolic パッケージは その正確性に関して細心の注意を払って作成されていますが、 結果に疑問がある場合には、この文書で示した方法により、 その結果を数値的に検証することが可能です。

また、上記の例では一つの数値例しか用いていませんが、乱数を用いて 多数の数値例を用いて検証することにより、より正確性への確信が高まると思われます。 ただ、乱数を用いる場合には、例えば逆行列や行列式の検証に際しては 作成された行列がフルランクであることを確かめてから利用するようにしてください。

########################################################################
# Stop Knitting Here !!!!!!
knitr::knit_exit()
########################################################################