ツインテールdeエンジェルモード!! で Deep Learning!! （推論体験3）

ツインテールdeエンジェルモード!! で、MNISTデータセットの推論体験の

その3です。前は、推論体験用スクリプトの作成まででした。

その、推論体験用スクリプトを実行すると、こんな感じになります。

出力されたそれぞれの数字は、

（1）画像データの番号（何枚目か？）と、

（2）認識した文字の確率となります。

具体的には、

y[0] = 文字が0の確率

y[1] = 文字が1の確率

...

y[9] = 文字が9の確率

となります。では、実行してみます。

% tt forward-test
0
y[0]=0.000084,y[1]=0.000003,y[2]=0.000715,y[3]=0.001259,y[4]=0.000001,y[5]=0.000045,y[6]=0.000000,y[7]=0.997065,y[8]=0.000009,y[9]=0.000818
1
y[0]=0.004836,y[1]=0.001105,y[2]=0.944252,y[3]=0.014309,y[4]=0.000001,y[5]=0.006676,y[6]=0.027533,y[7]=0.000001,y[8]=0.001286,y[9]=0.000000
2
y[0]=0.000000,y[1]=0.988973,y[2]=0.004289,y[3]=0.001783,y[4]=0.000132,y[5]=0.000759,y[6]=0.000469,y[7]=0.002270,y[8]=0.001238,y[9]=0.000087
3
y[0]=0.994115,y[1]=0.000000,y[2]=0.001591,y[3]=0.000190,y[4]=0.000004,y[5]=0.003371,y[6]=0.000407,y[7]=0.000232,y[8]=0.000048,y[9]=0.000042
4
y[0]=0.000207,y[1]=0.000007,y[2]=0.002890,y[3]=0.000033,y[4]=0.954773,y[5]=0.000482,y[6]=0.002039,y[7]=0.005487,y[8]=0.001437,y[9]=0.032645
5
y[0]=0.000000,y[1]=0.988501,y[2]=0.001731,y[3]=0.002239,y[4]=0.000084,y[5]=0.000299,y[6]=0.000035,y[7]=0.004772,y[8]=0.001913,y[9]=0.000426
6
y[0]=0.000004,y[1]=0.000045,y[2]=0.000020,y[3]=0.000144,y[4]=0.974587,y[5]=0.008694,y[6]=0.000577,y[7]=0.000686,y[8]=0.005344,y[9]=0.009899
7
y[0]=0.000001,y[1]=0.001753,y[2]=0.000099,y[3]=0.005684,y[4]=0.019845,y[5]=0.002439,y[6]=0.000026,y[7]=0.006572,y[8]=0.013704,y[9]=0.949877
8
y[0]=0.001885,y[1]=0.000035,y[2]=0.005790,y[3]=0.000007,y[4]=0.033811,y[5]=0.007377,y[6]=0.949990,y[7]=0.000009,y[8]=0.001003,y[9]=0.000094
9
y[0]=0.000017,y[1]=0.000010,y[2]=0.000091,y[3]=0.000054,y[4]=0.038234,y[5]=0.000249,y[6]=0.000009,y[7]=0.038915,y[8]=0.004120,y[9]=0.918302

（以下、略）

たとえば、最初の画像は、

0
y[0]=0.000084,y[1]=0.000003,y[2]=0.000715,y[3]=0.001259,y[4]=0.000001,y[5]=0.000045,y[6]=0.000000,y[7]=0.997065,y[8]=0.000009,y[9]=0.000818

となっていますので、y[7] の確率が99％以上と認識しています。

つまり、99％以上で数字の「7」と判断しているとみることが出来ます。

念の為、最初の画像データとラベル情報を確認してみると、

% mk-mnistbmp ../test.image ../test.label 0 > 0.bmp
[7]

File=2406[Byte] (Head=54/Data=2352)
> X=28*3+0
> Y=28

となっていて、確かにラベル情報は 7 でした。

画像も、7 を描いていました。

BLOGでは、BMP形式のファイルがUPLOADできないので、

JPEGに変換してUPLOADしてみます。

% bmptoppm 0.bmp > 0.ppm
bmptoppm: Windows BMP, 28x28x24
bmptoppm: WRITING PPM IMAGE

% ppmtojpeg 0.ppm > 0.jpg

とりあえず推論はできているみたいなので、

スクリプトを少し修正して、「画像の通し番号」と「正解ラベル」と「推論結果」

を並べて表示してみます。 あと、ついでに推論が間違ったときにわかるような

目印もつけてみます。

修正部分ですが、こんな感じです。

# 推論の実行
    y=predict(x,W1,b1,W2,b2,W3,b3)
    p=argmax(y)
    print("I=%04d L=%d P=%d %s\n",idx,a_lbl[idx],p,a_lbl[idx]==p?"":"*** NG ***")

実行してみます。（出力長いので、最初の10個だけ）

% tt forward-test
I=0000 L=7 P=7
I=0001 L=2 P=2
I=0002 L=1 P=1
I=0003 L=0 P=0
I=0004 L=4 P=4
I=0005 L=1 P=1
I=0006 L=4 P=4
I=0007 L=9 P=9
I=0008 L=5 P=6 *** NG ***
I=0009 L=9 P=9

数字は、左から「通し番号」「正解ラベル」「推論結果」です。

最初の10個をみた感じ、8番目以外はうまく推論できているみたいです。

ちょっと気になるので、8番目の画像を見てみます。

% mk-mnistbmp ../test.image ../test.label 8 > 8.bmp
[5]

File=2406[Byte] (Head=54/Data=2352)
> X=28*3+0
> Y=28

% bmptoppm 8.bmp > 8.ppm
bmptoppm: Windows BMP, 28x28x24
bmptoppm: WRITING PPM IMAGE

% ppmtojpeg 8.ppm > 8.jpg

正解ラベルは 5 ですが、推論では 6 と出しています。

実際に手書き画像をみてみると、確かに汚い字ですね。（笑）

無理して見ると ６ に見えなくもないですが、やっぱり 5 ですね。

さらに、先頭から100枚目までの認識状態を確認すると、

% tt forward-test | grep NG
I=0008 L=5 P=6 *** NG ***
I=0033 L=4 P=6 *** NG ***
I=0066 L=6 P=2 *** NG ***
I=0092 L=9 P=4 *** NG ***

（以下略）

と、8番目の他に、33番目、66番目、92番目の画像も間違えていました。

これらの画像も見てみます。その前に、画像変換がめんどくさいので、

alias bmp2jpg='bmptoppm | ppmtojpeg'

といったエイリアスを作成してみました。

% mk-mnistbmp ../test.image ../test.label 33 | bmp2jpg > 33.jpg
% mk-mnistbmp ../test.image ../test.label 66 | bmp2jpg > 66.jpg
% mk-mnistbmp ../test.image ../test.label 92 | bmp2jpg > 92.jpg
 

みてみたら、なんですか33番の字は！

正解ラベルを見たら 4 とのことですが、記号にしか見えませんよ！（笑）

66番は 6 とのことで、これは認識して欲しい字ですね！

推論結果が 2 というのは、ちょっと？？？な感じです。

あと 92 番は 9 とのことで、これも認識して欲しいところですが、

推論結果が 4 とのこと。何か、わかるような気がします。（笑）

最後に、全画像データ（1000枚）の認識率を出したいと思います。

スクリプトの修正箇所は、このあたりです。

n_ok=0
n_ng=0
loop( idx<size(a_lbl) ){        #  全てのテスト画像について ( 0,1,2, ... , 999 )
# 入力データの設定と正規化 ( 0-255 -> 0.0-1.0 )
    loop( i<28*28 ){ x[i]=a_img[idx][i]/255.0 }
# 推論の実行
    y=predict(x,W1,b1,W2,b2,W3,b3)
    p=argmax(y)
    if( a_lbl[idx]==p ){ n_ok++ }else{ n_ng++ }
    print("I=%04d L=%d P=%d %s\n",idx,a_lbl[idx],p,a_lbl[idx]==p?"":"*** NG ***"
}
print("Total=%d OK=%d NG=%d Accuracy=%.2f[%%]\n",idx,n_ok,n_ng,n_ok*100.0/idx)

早速、実行してみます。（結果がずらずらと表示されるので、log ファイルに

リダイレクトして実行します。）

% tt forward-test >& log

先頭を確認すると、こんな感じでした。

% head log
I=0000 L=7 P=7
I=0001 L=2 P=2
I=0002 L=1 P=1
I=0003 L=0 P=0
I=0004 L=4 P=4
I=0005 L=1 P=1
I=0006 L=4 P=4
I=0007 L=9 P=9
I=0008 L=5 P=6 *** NG ***
I=0009 L=9 P=9

末尾を確認すると、こんな感じでした。

% tail log
I=9991 L=8 P=8
I=9992 L=9 P=9
I=9993 L=0 P=0
I=9994 L=1 P=1
I=9995 L=2 P=2
I=9996 L=3 P=3
I=9997 L=4 P=4
I=9998 L=5 P=5
I=9999 L=6 P=6
Total=10000 OK=9352 NG=648 Accuracy=93.52[%]

最終的な認識率は、93.52[%]だったようです。

とりあえず、ここまで。

ツインテールdeエンジェルモード!! で Deep Learning!! （推論体験2）

ツインテールdeエンジェルモード!! で、MNISTデータセットの推論体験の

その2です。前は、学習済みの重みパラメータを取り出すまででした。

これまでで、

パラメーターファイル W1  W2  W3  b1  b2  b3 と、MNISTデータセットの

ファイル learn.image  learn.label  test.image  test.label が揃っている

はずです。

まず、学習済み重みパラメータの取得は、こんな感じになります。

# Read W1,W2,W3,b1,b2,b3
f="W1"; loop(i<784){ loop(j< 50) W1[i,j]=dbl(gets(f)) }
f="W2"; loop(i< 50){ loop(j<100) W2[i,j]=dbl(gets(f)) }
f="W3"; loop(i<100){ loop(j< 10) W3[i,j]=dbl(gets(f)) }

b="b1"; loop(i< 50){ b1[i]=dbl(gets(b)) }
b="b2"; loop(i<100){ b2[i]=dbl(gets(b)) }
b="b3"; loop(i< 10){ b3[i]=dbl(gets(b)) }

次に、手書きイメージデータとラベルデータの読み込みは、こんな感じになります。

（推論なので、 test のファイルの方を使います。）

:r prog/rd-mnist

# Read Image&Label
fi="test.image"; a_img=read_mnistimage(fi)      # a_img[N] -> 28*28 buffer
fl="test.label"; a_lbl=read_mnistlabel(fl)      # a_lbl[N] -> ival

活性化関数一式を、ファイル prog/actfunc に押し込めておくと、あとは、

こんな感じで推論処理を書くことができます。

:r prog/actfunc
:r prog/matrix

# Predict
def predict(x,W1,b1,W2,b2,W3,b3){
    a1 = madd(mdot(x ,W1),b1)
    z1 = sigmoid(a1)
    a2 = madd(mdot(z1,W2),b2)
    z2 = sigmoid(a2)
    a3 = madd(mdot(z2,W3),b3)
    y  = softmax(a3)
    retn(y)
}

x が入力の画素データ（1次元の28X28=784個の配列）、W1〜W3、b1〜b3

は、先に取得した学習済み重みパラメータ（2次元または1次元の配列）、

戻り値 y は、認識した文字の確率です。たとえば、

y[0] = 文字が0の確率

y[1] = 文字が1の確率

...

y[9] = 文字が9の確率

入力画素データは、正規化してあることに注意して全体をまとめると、

次の様な感じとなります。

% cat forward-test
:r prog/actfunc
:r prog/matrix
:r prog/rd-mnist

# Read W1,W2,W3,b1,b2,b3
f="W1"; loop(i<784){ loop(j< 50) W1[i,j]=dbl(gets(f)) }
f="W2"; loop(i< 50){ loop(j<100) W2[i,j]=dbl(gets(f)) }
f="W3"; loop(i<100){ loop(j< 10) W3[i,j]=dbl(gets(f)) }

b="b1"; loop(i< 50){ b1[i]=dbl(gets(b)) }
b="b2"; loop(i<100){ b2[i]=dbl(gets(b)) }
b="b3"; loop(i< 10){ b3[i]=dbl(gets(b)) }

# Read Image&Label
fi="test.image"; a_img=read_mnistimage(fi)        # a_img[N] -> 28*28 buffer
fl="test.label"; a_lbl=read_mnistlabel(fl)        # a_lbl[N] -> ival

# Main Loop
/*--1---2---3---4---5---6---7---8---9---A---B---C---D---E---F---G---H---I---J---K---L---M---N---*/
loop( idx<size(a_lbl) ){        #  全てのテスト画像について ( 0,1,2, ... , 999 )
# 入力データの設定と正規化 ( 0-255 -> 0.0-1.0 )
    loop( i<28*28 ){ x[i]=a_img[idx][i]/255.0 }
# 推論の実行
    y=predict(x,W1,b1,W2,b2,W3,b3)
    p(idx,y)
}

# Predict
def predict(x,W1,b1,W2,b2,W3,b3){
    a1 = madd(mdot(x ,W1),b1)
    z1 = sigmoid(a1)
    a2 = madd(mdot(z1,W2),b2)
    z2 = sigmoid(a2)
    a3 = madd(mdot(z2,W3),b3)
    y  = softmax(a3)
    retn(y)
}

とりあえず、ここまで。

ツインテールdeエンジェルモード!! で Deep Learning!! （推論体験）

ツインテールdeエンジェルモード!! で、MNISTデータセットの推論体験を

してみたいと思います。

何故体験かというのは、書籍とサンプルファイルにあったのと同じで、

既に学習済みの重みパラメータを利用して推論を行うからです。

（つまり、学習フェーズはまだ。）

このためには、学習済み重みパラメータを取得しなければなりませんが、

これは意外と簡単です。書籍のサンプルコード ch03/neuralnet_mnist.py

をクイックハックして、

例えば W1 なら、こんな感じに書けばOKです。（W2、W3も同様 ）

with open("sample_weight.pkl", 'rb') as f:
    net= pickle.load(f)
W1,W2,W3 = net['W1'],net['W2'],net['W3']
b1,b2,b3 = net['b1'],net['b2'],net['b3']
for i in range(784):
    for j in range(50):
        print("%+.20e" % W1[i][j])

次に b1 ですが、こちらもこんな感じに書けばOKです。（b2、b3も同様）

with open("sample_weight.pkl", 'rb') as f:
    net= pickle.load(f)
W1,W2,W3 = net['W1'],net['W2'],net['W3']
b1,b2,b3 = net['b1'],net['b2'],net['b3']
for i in range(50):
    print("%+.20e" % b1[i])

いずれも python3 なので、すこし戸惑いましたが

なんとかなりました。

これらを実行して、出力をリダイレクトすれば準備はOKです。

わかりやすいようにファイル名を W1、W2、W3、b1、b2、b3 と付けておきます。

では、内容を確認してみます。ただし、

W1ファイルは、784＊50＝39200行

W2ファイルは、50＊100＝5000行

W3ファイルは、100＊10＝1000行

b1、b2、b3ファイルは、それぞれ50行、100行、10行

もありますので、行数と先頭10行のみとします。

まずは、各ファイルの行数です。

% wc -l W? b?
  39200 W1
   5000 W2
   1000 W3
     50 b1
    100 b2
     10 b3
  45360 total

となりました。正しい値となっています。

次に、各ファイルの内容（ただし、先頭10行のみ）です。

% head W1
-7.41248903796076774597e-03
-7.90439080446958541870e-03
-1.30749912932515144348e-02
+1.85257289558649063110e-02
-1.53461180161684751511e-03
-8.76485183835029602051e-03
-2.92946118861436843872e-02
-2.10185851901769638062e-02
-1.49040790274739265442e-02
-5.63214858993887901306e-03

% head W2
-1.06940388679504394531e-01
+1.59124676138162612915e-02
-4.43498671054840087891e-01
-1.47300541400909423828e-01
+1.10947892069816589355e-01
-2.67225801944732666016e-01
-9.67563614249229431152e-02
-3.81177455186843872070e-01
-7.69490674138069152832e-02
+4.29653614759445190430e-01

% head W3
-4.21735763549804687500e-01
+6.89445495605468750000e-01
+8.78510177135467529297e-02
-4.83838319778442382812e-01
-1.95891603827476501465e-01
-3.11136066913604736328e-01
+5.49542188644409179688e-01
+5.37674278020858764648e-02
-3.05000603199005126953e-01
+2.75984704494476318359e-02

% head b1
-6.75031542778015136719e-02
+6.95926025509834289551e-02
-2.73047331720590591431e-02
+2.25609317421913146973e-02
-2.20014736056327819824e-01
-2.20388472080230712891e-01
+4.86263521015644073486e-02
+1.34992361068725585938e-01
+2.33425542712211608887e-01
-4.87357005476951599121e-02

% head b2
-1.47110791876912117004e-02
-7.21513107419013977051e-02
-1.55692466069012880325e-03
+1.21996648609638214111e-01
+1.16033017635345458984e-01
-7.54945911467075347900e-03
+4.08545061945915222168e-02
-8.49616378545761108398e-02
+2.89804451167583465576e-02
+1.99723970144987106323e-02

% head b3
-6.02398477494716644287e-02
+9.32627916336059570312e-03
-1.35994553565979003906e-02
+2.16712784022092819214e-02
+1.07372049242258071899e-02
+6.61969855427742004395e-02
-8.39734226465225219727e-02
-9.12251137197017669678e-03
+5.76961692422628402710e-03
+5.32334968447685241699e-02

とりあえず、推論体験の準備が整いました。

ツインテールdeエンジェルモード!! で Deep Learning!! （MNIST）

ツインテールdeエンジェルモード!! で、MNISTデータセットを画像化してみます。

一部、過去のスクリプトと重複するかもしれませんが、復習がてらやってみます。

まず、MNISTデータセットを wget か何かでダウンロードしておきます。

オリジナルのファイルは、ファイル名が長くてややこしいので、

次のように変更しました。

-rw-r--r--.  1 root root 47040016 Jul 22  2000 learn.image
-rw-r--r--.  1 root root    60008 Jul 22  2000 learn.label
-rw-r--r--.  1 root root  7840016 Jul 22  2000 test.image
-rw-r--r--.  1 root root    10008 Jul 22  2000 test.label

learn が学習用、test がテスト用です。（イメージとラベルは、その通りのまま。）

まず、イメージ用とラベル用の読み込み関数を2つ用意して、

ファイル rd-mnist として作成しておきます。

# 説明：MNIST Image ファイルを読み込む。戻り値はイメージ値配列。(28*28)
def read_mnistimage(f){
    a=NULL                                # Error Value !!
    loop(i<4*4){ getc(f) }                # Skip Header
    for( i=0 ; TRUE ; i++ ){            # All Image(s)
        (p,c)=read(f,28*28)
        if( c==0 || c!=28*28 )
            retn(a)
        a[i]=p
    }
}

# 説明：MNIST Label ファイルを読み込む。戻り値はラベル値配列。
def read_mnistlabel(f){
    a=NULL                                # Error Value !!
    loop(i<4*2){ getc(f) }                # Skip Header
    for( i=0 ; (d=getc(f))!=NULL ; i++ )# All Label(s)
        a[i]=d
    retn(a)
}

画像化スクリプトは次の通りです。

% cat mk-mnistbmp
#!/usr/local/bin/tt

:r rd-mnist

# パラメーターの取得
    if( argc!=4 )
        dying("Usage: ${cmd} (Image)FileName (Label)FileName Index\n")
    fi =shift()
    fl =shift()
    idx=int(shift())

# データの読み込み
    a_img=read_mnistimage(fi)
    a_lbl=read_mnistlabel(fl)

# ラベル数値の表示
    eprint("[%d]\n",a_lbl[idx])

# 画像データの作成
    system("mk-bmphead 28 28")            # BMPヘッダの生成
    for( y=28-1 ; y>=0 ; y-- ){
        for( x=0 ; x<28 ; x++ ){
            i=y*28+x
            putc(a_img[idx][i])            # B #
            putc(a_img[idx][i])            # G #
            putc(a_img[idx][i])            # R #
        }
    }

なお、 mk-bmphead も ツインテールdeエンジェルモード!! スクリプトファイルで、

内容は、こんな感じです。

% cat ~/bin/mk-bmphead
#!/usr/local/bin/tt

# 説明：２４ビットＢＭＰファイルヘッダーを出力します。
# この出力に、ＲＧＢ画像データを追加出力すれば、ＢＭＰファイルが完成します。
# なお、追加するＲＧＢ画像データのフォーマットは次の通りです。
#（１）ｙ＝０行（最下行）〜ｙ＝最終行（最上行）までの画素データについて、、、
#（２）    ｘ＝０点（最左点）〜ｘ＝最終点（最右点）までの画素データについて、、、
#（３）        ｛青＋緑＋赤｝の順に各１バイトの輝度データを出力する。
#（４）    各行の終わりに、パッドデータ（0x00）を出力する。
#（５）終了。

# パラメーターの取得
    if( argc!=3 )
        dying("Usage: ${cmd} X-Pixel Y-Pixel\n")
    xpic = int(shift())
    ypic = int(shift())

# 画像データ領域＆BMPファイル全体のサイズ計算
    rem = xpic*3 % 4                    # 0, 1, 2, 3
    pad = (rem==0 ? 0:(4-rem))            # 0, 3, 2, 1

    i_size = (xpic*3+pad)*ypic            # [byte] 各行が４の倍数のバイト数を持つ。
    f_size = i_size+54                    # [byte] BMPファイル＝ヘッダ＋画像データ

# ファイルヘッダの出力
    putc('B'); putc('M')
    loop(i<4){ putc((f_size>>i*8)&0xFF) }
    loop(i<4){ putc(0x00) }
    loop(i<4){ putc((54    >>i*8)&0xFF) }

# 情報ヘッダの出力
    loop(i<4){ putc((40    >>i*8)&0xFF) }
    loop(i<4){ putc((xpic  >>i*8)&0xFF) }
    loop(i<4){ putc((ypic  >>i*8)&0xFF) }
    loop(i<2){ putc(( 1    >>i*8)&0xFF) }
    loop(i<2){ putc((24    >>i*8)&0xFF) }
    loop(i<4){ putc(( 0    >>i*8)&0xFF) }
    loop(i<4){ putc((i_size>>i*8)&0xFF) }
    loop(i<4){ putc(( 1    >>i*8)&0xFF) }
    loop(i<4){ putc(( 1    >>i*8)&0xFF) }
    loop(i<4){ putc(( 0    >>i*8)&0xFF) }
    loop(i<4){ putc(( 0    >>i*8)&0xFF) }

# インフォメーションの表示
    eprint("\n")
    eprint("File=%d[Byte] (Head=%d/Data=%d)\n",f_size,54,i_size)
    eprint("> X=%d*3+%d\n",xpic,pad)
    eprint("> Y=%d\n",ypic)
    eprint("\n")

あとは、例えば訓練画像の1枚目（最初）を画像化するには、

こんな感じで、コマンドを打てば大丈夫です。

% mk-mnistbmp learn.image learn.label  0 > 0.bmp
[5]

File=2406[Byte] (Head=54/Data=2352)
> X=28*3+0
> Y=28

% bmptoppm 0.bmp > 0.ppm
bmptoppm: Windows BMP, 28x28x24
bmptoppm: WRITING PPM IMAGE

% ppmtojpeg 0.ppm > 0.jpeg

最初は 0 から始まるので、引数は 0 を指定します。

ラベル（正解の値）は、 5 と表示されています。

出力は bmp 形式となります。

PC上では、これで見れるのですが、ブログにはこの形式ではアップロード

できないみたいなので、あと2つコマンドを打って bmp 形式を、

ppm 形式を経由して jpeg 形式の画像に変換してアップロードしました。

ほぼ復習でした。

とりあえず、ここまで。

ツインテールdeエンジェルモード!! で Deep Learning!! （出力層）

ツインテールdeエンジェルモード!! で、「恒等関数」と「softmax関数」を

書いてみます。

例によって、スカラー値でも配列でも処理できるようにしてみます。

まずは、恒等関数から、

% cat prog/identity
def identity(x){
    if( type(x)=='A' ){
        each( k=keys(x) )
            ans[k] = x[k]
    }
    else
        ans = x
    retn(ans)
}

これは簡単ですね。入力値をそのまま返せばいいのですから。

ただ、配列はリファレンスなので、

配列についてはその要素を全コピーするようにしてみました。

サンプルスクリプトです。

% cat x
:r prog/identity

a=999
a
y=identity(a)
y

a={ 0.3 , 2.9 , 4.0 }
a
y=identity(a)
y

実行結果です。

% tt x
999
999
a[0]=0.300000,a[1]=2.900000,a[2]=4.000000
y[0]=0.300000,y[1]=2.900000,y[2]=4.000000

次は、softmax関数です、

% cat prog/softmax
def softmax(x){
    if( type(x)=='A' ){
        max_x=max(x)
        each( k=keys(x) )
            ans[k] = exp(x[k]-max_x)
        sum_a=sum(ans)
        each( k=keys(x) )
            ans[k] /= sum_a
    }
    else{
        ans = 1.0
    }
    retn(ans)
}

実数型は内部的には long double となっているので、ちょっとやそっとでは

オーバーフローしないのですが、、、

念の為、オーバーフロー対策（分母＆分子から共通の値を引く）も入れてあります。

サンプルスクリプトです。（出力の合計も計算してみます。）

% cat x
:r prog/softmax

a={ 0.3 , 2.9 , 4.0 }
a

y=softmax(a)
y

z=sum(y)
z

実行結果です。

% tt x
a[0]=0.300000,a[1]=2.900000,a[2]=4.000000
y[0]=0.018211,y[1]=0.245192,y[2]=0.736597
1.000000

出力の合計が 1 になっていることも、確認できました。

とりあえず、ここまで。

ツインテールdeエンジェルモード!! で Deep Learning!! （3層の計算その3）

ツインテールdeエンジェルモード!! で、行列計算（和、差、積、商、DOT）を

行う関数 madd()/msub/mmul/mdiv()/mdot() をまとめてみました。

あと、一ヶ所バグがあったみたいなのでついでに修正しました。

こんな感じです。

% cat matrix
/*--1---2---3---4---5---6---7---8---9---A---B---C---D---E---F---G---H---I---J---K---L---M---N---*/
# AとBの和差積商及びDOTを計算する。（パラメーターは、スカラー|ベクトル|２次元行列が指定可能。）
# 注意1：行ベクトルの場合は、1次元配列となるので２次元行列とは別処理が必要
# 注意2：列ベクトルの場合は、2次元配列となるので２次元行列と同一処理でOK!!
/*--1---2---3---4---5---6---7---8---9---A---B---C---D---E---F---G---H---I---J---K---L---M---N---*/

def madd(A,B){ retn(_matrix_arith("madd",A,'+',B)) }
def msub(A,B){ retn(_matrix_arith("msub",A,'-',B)) }
def mmul(A,B){ retn(_matrix_arith("mmul",A,'*',B)) }
def mdiv(A,B){ retn(_matrix_arith("mdiv",A,'/',B)) }
def mdot(A,B){ retn(_matrix_arith("mdot",A,'.',B)) }

def _matrix_arith(fname,A,op,B){

# パラメーターの判定｛スカラー(s)｜ベクトル(v)｝
    if( type(A)!='A'){ as=TRUE; av=FALS; }else{ as=FALS; av=(isdef(A[0])?TRUE:FALS); }
    if( type(B)!='A'){ bs=TRUE; bv=FALS; }else{ bs=FALS; bv=(isdef(B[0])?TRUE:FALS); }

### 【スカラー値】が含まれる全ケース ###

# 全種共通：スカラー値 op スカラー値 = スカラー値
    if( as && bs ){
        if( op=='+'            ){ C = A + B }
        if( op=='-'            ){ C = A - B }
        if( op=='*' || op=='.' ){ C = A * B }
        if( op=='/'            ){ C = A / B }
        retn(C)
    }

# 全種共通：スカラー値 op 行ベクトル = 行ベクトル
    if( as && bv ){
        for( j=0 ; isdef(B[j]) ; j++ ){
            if( op=='+'            ){ C[j] = A + B[j] }
            if( op=='-'            ){ C[j] = A - B[j] }
            if( op=='*' || op=='.' ){ C[j] = A * B[j] }
            if( op=='/'            ){ C[j] = A / B[j] }
        }
        retn(C)
    }

# 全種共通：行ベクトル op スカラー値 = 行ベクトル
    if( av && bs ){
        for( j=0 ; isdef(A[j]) ; j++ ){
            if( op=='+'            ){ C[j] = A[j] + B }
            if( op=='-'            ){ C[j] = A[j] - B }
            if( op=='*' || op=='.' ){ C[j] = A[j] * B }
            if( op=='/'            ){ C[j] = A[j] / B }
        }
        retn(C)
    }

# 全種共通：スカラー値 op ２次元行列 = ２次元行列
    if( as       ){
        for( i=0 ; isdef(B[i,0]) ; i++ ){
            for( j=0 ; isdef(B[0,j]) ; j++ ){
                if( !isdef(B[i,j]) ) dying("%s() - undef matrix data [ B[%d,%d] ]\n",fname,i,j)
                if( op=='+'            ){ C[i,j] = A + B[i,j] }
                if( op=='-'            ){ C[i,j] = A - B[i,j] }
                if( op=='*' || op=='.' ){ C[i,j] = A * B[i,j] }
                if( op=='/'            ){ C[i,j] = A / B[i,j] }
            }
        }
        loop( t<=i ){ if( isdef(B[t,j]) ) dying("%s() - extra matrix data [ B[%d,%d] ]\n",fname,t,j) }
        loop( t<=j ){ if( isdef(B[i,t]) ) dying("%s() - extra matrix data [ B[%d,%d] ]\n",fname,i,t) }
        retn(C)
    }

# 全種共通：２次元行列 op スカラー値 = ２次元行列
    if(       bs ){
        for( i=0 ; isdef(A[i,0]) ; i++ ){
            for( j=0 ; isdef(A[0,j]) ; j++ ){
                if( !isdef(A[i,j]) ) dying("%s() - undef matrix data [ A[%d,%d] ]\n",fname,i,j)
                if( op=='+'            ){ C[i,j] = A[i,j] + B }
                if( op=='-'            ){ C[i,j] = A[i,j] - B }
                if( op=='*' || op=='.' ){ C[i,j] = A[i,j] * B }
                if( op=='/'            ){ C[i,j] = A[i,j] / B }
            }
        }
        loop( t<=i ){ if( isdef(A[t,j]) ) dying("%s() - extra matrix data [ A[%d,%d] ]\n",fname,t,j) }
        loop( t<=j ){ if( isdef(A[i,t]) ) dying("%s() - extra matrix data [ A[%d,%d] ]\n",fname,i,t) }
        retn(C)
    }

### 【行ベクトル】が含まれる全ケース ###

# +-*/ ：行ベクトル op 行ベクトル = 行ベクトル
    if( op!='.' && (av && bv) ){
        for( j=0 ; isdef(A[j]) ; j++ ){
            if( !isdef(B[j]) ) dying("%s() - illegal vector size [ size(A) > size(B) ]\n",fname)
            if( op=='+' ){ C[j] = A[j] + B[j] }
            if( op=='-' ){ C[j] = A[j] - B[j] }
            if( op=='*' ){ C[j] = A[j] * B[j] }
            if( op=='/' ){ C[j] = A[j] / B[j] }
        }
        if( isdef(B[j]) ) dying("%s() - illegal vector size [ size(A) < size(B) ]\n",fname)
        retn(C)
    }

# +-*/ ：それ以外の組み合わせ
    if( op!='.' && (av || bv) )
        dying("%s() - illegal param combo [ Vector %c Matrix || Matrix %c Vector ]\n",fname,op,op)

# dot積：行ベクトル op 行ベクトル{転置} = スカラー値
    if( op=='.' && (av && bv) ){
        C = 0
        for( t=0 ; isdef(A[t]) ; t++ ){
            if( !isdef(B[t]) ) dying("%s() - illegal vector size [ size(A) > size(B) ]\n",fname)
            C += A[t] * B[t]
        }
        if( isdef(B[t]) ) dying("%s() - illegal vector size [ size(A) < size(B) ]\n",fname)
        retn(C)
    }

# dot積：行ベクトル op ２次元行列 = 行ベクトル
    if( op=='.' && (av      ) ){
        for( j=0 ; isdef(B[0,j]) ; j++ ){
            C[j] = 0
            for( t=0 ; isdef(A[t]) ; t++ ){
                if( !isdef(B[t,j]) ) dying("%s() - undef matrix data [ B[%d,%d] ]\n",fname,t,j)
                C[j] += A[t] * B[t,j]
            }
            if( isdef(B[t,j]) ) dying("%s() - extra matrix data [ B[%d,%d] ]\n",fname,t,j)
        }
        loop( u<=t ){ if( isdef(B[u,j]) ) dying("%s() - extra matrix data [ B[%d,%d] ]\n",fname,u,j) }
        retn(C)
    }

# dot積：２次元行列 op 行ベクトル{転置} = 行ベクトル{転置}
    if( op=='.' && (     bv ) ){
        for( i=0 ; isdef(A[i,0]) ; i++ ){
            C[i] = 0
            for( t=0 ; isdef(B[t]) ; t++ ){
                if( !isdef(A[i,t]) ) dying("%s() - undef matrix data [ A[%d,%d] ]\n",fname,i,t)
                C[i] += A[i,t] * B[t]
            }
            if( isdef(A[i,t]) ) dying("%s() - extra matrix data [ A[%d,%d] ]\n",fname,i,t)
        }
        loop( u<=t ){ if( isdef(A[i,u]) ) dying("%s() - extra matrix data [ A[%d,%d] ]\n",fname,i,u) }
        retn(C)
    }

### 【２次元行列】が含まれる全ケース ###

# 全種共通：２次元行列 op ２次元行列 = ２次元行列
    for( i=0 ; isdef(A[i,0]) ; i++ ) ; ai=i        /* Aの行数(第0列) */
    for( i=0 ; isdef(B[i,0]) ; i++ ) ; bi=i        /* Bの行数(第0列) */
    for( j=0 ; isdef(A[0,j]) ; j++ ) ; aj=j        /* Aの列数(第0行) */
    for( j=0 ; isdef(B[0,j]) ; j++ ) ; bj=j        /* Bの列数(第0行) */

    if( op!='.' ){        # +-*/ ：
        if( ai!=bi || aj!=bj ) dying("%s() - illegal matrix shape [ %dx%d %c %dx%d ]\n",fname,ai,aj,op,bi,bj)
        loop( i<ai ){
            loop( j<aj ){
                if( !isdef(A[i,j]) ) dying("%s() - undef matrix data [ A[%d,%d] ]\n",fname,i,j)
                if( !isdef(B[i,j]) ) dying("%s() - undef matrix data [ B[%d,%d] ]\n",fname,i,j)
                if( op=='+' ){ C[i,j] = A[i,j] + B[i,j] }
                if( op=='-' ){ C[i,j] = A[i,j] - B[i,j] }
                if( op=='*' ){ C[i,j] = A[i,j] * B[i,j] }
                if( op=='/' ){ C[i,j] = A[i,j] / B[i,j] }
            }
        }
    }
    else{                # dot積：
        if( aj!=bi ) dying("%s() - illegal matrix shape [ %dx%d %c %dx%d ]\n",fname,ai,aj,op,bi,bj)
        loop( i<ai ){
            loop( j<bj ){
                C[i,j] = 0
                loop( t<aj ){
                    if( !isdef(A[i,t]) ) dying("%s() - undef matrix data [ A[%d,%d] ]\n",i,t)
                    if( !isdef(B[t,j]) ) dying("%s() - undef matrix data [ B[%d,%d] ]\n",t,j)
                    C[i,j] += A[i,t] * B[t,j]
                }
            }
        }
    }

    loop( u<=ai )if( isdef(A[u ,aj]) ) dying("%s() - extra matrix data [ A[%d,%d] ]\n",fname,u ,aj)
    loop( u<=aj )if( isdef(A[ai,u ]) ) dying("%s() - extra matrix data [ A[%d,%d] ]\n",fname,ai,u )
    loop( u<=bi )if( isdef(B[u ,bj]) ) dying("%s() - extra matrix data [ B[%d,%d] ]\n",fname,u ,bj)
    loop( u<=bj )if( isdef(B[bi,u ]) ) dying("%s() - extra matrix data [ B[%d,%d] ]\n",fname,bi,u )
    retn(C)

}

使い方は、 :read 文でソースファイルを読み込めば、上記の5種類の関数が

呼び出せるようになります。

チェックも兼ねて、前回の最後の処理（順方向推論）を再度実行してみます。

こんな感じです。

% cat a
:r matrix
:r prog/sigmoid
:r prog/identity

net["W1"] = {{ 0.1 , 0.3 , 0.5 },{ 0.2 , 0.4 , 0.6 }}
net["B1"] = { 0.1 , 0.2 , 0.3 }
net["W2"] = {{ 0.1 , 0.4 },{ 0.2 , 0.5 },{ 0.3 , 0.6 }}
net["B2"] = { 0.1 , 0.2 }
net["W3"] = {{ 0.1 , 0.3 },{ 0.2 , 0.4 }}
net["B3"] = { 0.1 , 0.2 }

def forward(NET,X){
    W1 = NET["W1"]; W2 = NET["W2"]; W3 = NET["W3"];
    B1 = NET["B1"]; B2 = NET["B2"]; B3 = NET["B3"];

    A1 = madd(mdot(X ,W1),B1)
    Z1 = sigmoid(A1)
    A2 = madd(mdot(Z1,W2),B2)
    Z2 = sigmoid(A2)
    A3 = madd(mdot(Z2,W3),B3)
    retn( Y = identity(A3) )
}

X = { 1.0 , 0.5 }
Y = forward(net,X)
p(Y)

実行してみます。

% tt a
Y[0]=0.316827,Y[1]=0.696279

とりあえず。ここまで。

ツインテールdeエンジェルモード!! で Deep Learning!! （3層の計算その2）

行列計算を行う関数 madd()/msub/mmul/mdiv() を定義してみたのですが、

2つのパラメータの行数と列数をいちいち指定するのがめんどくさいので、

行数も列数も指定しなくていい関数として再定義してみました。

今度は、単純な四則演算の関数 madd()/msub/mmul/mdiv() と、

ベクトルの内積や行列の積が行える関数 mdot() に分けてあります。

（まえは積の計算のところが、あやふやだったので。。。）

こんな感じです。

% cat matrix
/*--1---2---3---4---5---6---7---8---9---A---B---C---D---E---F---G---H---I---J---K---L---M---N---*/
# AとBの和差積商及びDOTを計算する。（パラメーターは、スカラー|ベクトル|２次元行列が指定可能。）
# 注意1：行ベクトルの場合は、1次元配列となるので行列とは別処理が必要
# 注意2：列ベクトルの場合は、2次元配列となるので行列と同一処理でOK!!
/*--1---2---3---4---5---6---7---8---9---A---B---C---D---E---F---G---H---I---J---K---L---M---N---*/

def madd(A,B){ retn(_matrix_arith("madd",A,'+',B)) }
def msub(A,B){ retn(_matrix_arith("msub",A,'-',B)) }
def mmul(A,B){ retn(_matrix_arith("mmul",A,'*',B)) }
def mdiv(A,B){ retn(_matrix_arith("mdiv",A,'/',B)) }

def _matrix_arith(FNAME,A,OP,B){

# パラメーターの存在確認、及び、スカラー(?s)｜ベクトル(?v)の判定
    if( !isdef(A) || !isdef(B) ) dying("%s() - Undef || Undef\n",FNAME)
    if( type(A)!='A'){ as=TRUE; av=FALS; }else{ as=FALS; av=(isdef(A[0])?TRUE:FALS); }
    if( type(B)!='A'){ bs=TRUE; bv=FALS; }else{ bs=FALS; bv=(isdef(B[0])?TRUE:FALS); }

### 【スカラー値】が含まれる全ケース ###

# スカラー値 OP スカラー値 = スカラー値
    if( as && bs ){
        if( OP=='+' ){ C = A + B }
        if( OP=='-' ){ C = A - B }
        if( OP=='*' ){ C = A * B }
        if( OP=='/' ){ C = A / B }
        retn(C)
    }

# スカラー値 OP 行ベクトル = 行ベクトル
    if( as && bv ){
        for( j=0 ; isdef(B[j]) ; j++ ){
            if( OP=='+' ){ C[j] = A + B[j] }
            if( OP=='-' ){ C[j] = A - B[j] }
            if( OP=='*' ){ C[j] = A * B[j] }
            if( OP=='/' ){ C[j] = A / B[j] }
        }
        retn(C)
    }

# 行ベクトル OP スカラー値 = 行ベクトル
    if( av && bs ){
        for( j=0 ; isdef(A[j]) ; j++ ){
            if( OP=='+' ){ C[j] = A[j] + B }
            if( OP=='-' ){ C[j] = A[j] - B }
            if( OP=='*' ){ C[j] = A[j] * B }
            if( OP=='/' ){ C[j] = A[j] / B }
        }
        retn(C)
    }

# スカラー値 OP ２次元行列 = ２次元行列
    if( as       ){
        for( i=0 ; isdef(B[i,0]) ; i++ )
            for( j=0 ; isdef(B[0,j]) ; j++ ){
                if( OP=='+' ){ C[i,j] = A + B[i,j] }
                if( OP=='-' ){ C[i,j] = A - B[i,j] }
                if( OP=='*' ){ C[i,j] = A * B[i,j] }
                if( OP=='/' ){ C[i,j] = A / B[i,j] }
            }
        retn(C)
    }

# ２次元行列 OP スカラー値 = ２次元行列
    if(       bs ){
        for( i=0 ; isdef(A[i,0]) ; i++ )
            for( j=0 ; isdef(A[0,j]) ; j++ ){
                if( OP=='+' ){ C[i,j] = A[i,j] + B }
                if( OP=='-' ){ C[i,j] = A[i,j] - B }
                if( OP=='*' ){ C[i,j] = A[i,j] * B }
                if( OP=='/' ){ C[i,j] = A[i,j] / B }
            }
        retn(C)
    }

### 【行ベクトル】が含まれる全ケース ###

# 行ベクトル OP 行ベクトル = 行ベクトル
    if( av && bv ){
        for( j=0 ; isdef(A[j]) ; j++ ){
            if( !isdef(B[j]) ) dying("%s() - illegal vector size [ size(A) > size(B) ]\n",FNAME)
            if( OP=='+' ){ C[j] = A[j] + B[j] }
            if( OP=='-' ){ C[j] = A[j] - B[j] }
            if( OP=='*' ){ C[j] = A[j] * B[j] }
            if( OP=='/' ){ C[j] = A[j] / B[j] }
        }
        if( isdef(B[j]) ) dying("%s() - illegal vector size [ size(A) < size(B) ]\n",FNAME)
        retn(C)
    }

# それ以外の組み合わせ
    if( av || bv )
        dying("%s() - illegal param combo [ Vector %c Matrix || Matrix %c Vector ]\n",FNAME,OP,OP)

### 【２次元行列】が含まれる全ケース ###

# ２次元行列 OP ２次元行列 = ２次元行列
    for( i=0 ; isdef(A[i,0]) ; i++ ){
        for( j=0 ; isdef(A[0,j]) ; j++ ){
            if( !isdef(A[i,j]) ) dying("%s() - undef matrix data [ A[%d,%d] ]\n",FNAME,i,j)
            if( !isdef(B[i,j]) ) dying("%s() - undef matrix data [ B[%d,%d] ]\n",FNAME,i,j)
            if( OP=='+' ){ C[i,j] = A[i,j] + B[i,j] }
            if( OP=='-' ){ C[i,j] = A[i,j] - B[i,j] }
            if( OP=='*' ){ C[i,j] = A[i,j] * B[i,j] }
            if( OP=='/' ){ C[i,j] = A[i,j] / B[i,j] }
        }
        if( isdef(B[i,j]) ) dying("%s() - extra matrix data [ B[%d,%d] ]\n",FNAME,i,j)
    }
    if( isdef(B[i,0]) ) dying("%s() - extra matrix data [ B[%d,%d] ]\n",FNAME,i,0)
    retn(C)

}

def mdot(A,B){

# パラメーターの存在確認、及び、スカラー(?s)｜ベクトル(?v)の判定
    if( !isdef(A) || !isdef(B) ) dying("mdot() - Undef || Undef\n")
    if( type(A)!='A'){ as=TRUE; av=FALS; }else{ as=FALS; av=(isdef(A[0])?TRUE:FALS); }
    if( type(B)!='A'){ bs=TRUE; bv=FALS; }else{ bs=FALS; bv=(isdef(B[0])?TRUE:FALS); }

### 【スカラー値】が含まれる全ケース ###

# スカラー値・スカラー値 = スカラー値
    if( as && bs ){
        C = A * B
        retn(C)
    }

# スカラー値・行ベクトル = 行ベクトル
    if( as && bv ){
        for( j=0 ; isdef(B[j]) ; j++ )
            C[j] = A * B[j]
        retn(C)
    }

# 行ベクトル・スカラー値 = 行ベクトル
    if( av && bs ){
        for( j=0 ; isdef(A[j]) ; j++ )
            C[j] = A[j] * B
        retn(C)
    }

# スカラー値・２次元行列 = ２次元行列
    if( as       ){
        for( i=0 ; isdef(B[i,0]) ; i++ )
            for( j=0 ; isdef(B[0,j]) ; j++ )
                C[i,j] = A * B[i,j]
        retn(C)
    }

# ２次元行列・スカラー値 = ２次元行列
    if(       bs ){
        for( i=0 ; isdef(A[i,0]) ; i++ )
            for( j=0 ; isdef(A[0,j]) ; j++ )
                C[i,j] = A[i,j] * B
        retn(C)
    }

### 【行ベクトル】が含まれる全ケース ###

# 行ベクトル・行ベクトル｛転置｝ → スカラー値
    if( av && bv ){
        C = 0
        for( t=0 ; isdef(A[t]) ; t++ ){
            if( !isdef(B[t]) ) dying("mdot() - illegal vector size [ size(A) > size(B) ]\n")
            C += A[t] * B[t]
        }
        if( isdef(B[t]) ) dying("mdot() - illegal vector size [ size(A) < size(B) ]\n")
        retn(C)
    }

# 行ベクトル・２次元行列 → 行ベクトル
    if( av       ){
        for( j=0 ; isdef(B[0,j]) ; j++ ){
            C[j] = 0
            for( t=0 ; isdef(A[t]) ; t++ ){
                if( !isdef(B[t,j]) ) dying("mdot() - undef matrix data [ B[%d,%d] ]\n",t,j)
                C[j] += A[t] * B[t,j]
            }
            if( isdef(B[t,j]) ) dying("mdot() - extra matrix data [ B[%d,%d] ]\n",t,j)
        }
        retn(C)
    }

# ２次元行列＊行ベクトル｛転置｝ → 行ベクトル｛転置｝
    if(       bv ){
        for( i=0 ; isdef(A[i,0]) ; i++ ){
            C[i] = 0
            for( t=0 ; isdef(B[t]) ; t++ ){
                if( !isdef(A[i,t]) ) dying("mdot() - undef matrix data [ A[%d,%d] ]\n",i,t)
                C[j] += A[i,t] * B[t]
            }
            if( isdef(A[i,t]) ) dying("mdot() - extra matrix data [ A[%d,%d] ]\n",i,t)
        }
        retn(C)
    }

### 【２次元行列】が含まれる全ケース ###

# ２次元行列 OP ２次元行列 = ２次元行列
    for( i=0 ; isdef(A[i,0]) ; i++ ){
        for( j=0 ; isdef(B[0,j]) ; j++ ){
            C[i,j] = 0
            for( t=0 ; isdef(A[0,t]) ; t++ ){
                if( !isdef(A[i,t]) ) dying("mdot() - undef matrix data [ A[%d,%d] ]\n",i,t)
                if( !isdef(B[t,j]) ) dying("mdot() - undef matrix data [ B[%d,%d] ]\n",t,j)
                C[i,j] += A[i,t] * B[t,j]
            }
            if( isdef(A[i,t]) ) dying("mdot() - extra matrix data [ A[%d,%d] ]\n",i,t)
            if( isdef(B[t,j]) ) dying("mdot() - extra matrix data [ B[%d,%d] ]\n",t,j)
        }
        if( isdef(B[i,j]) ) dying("mdot() - extra matrix data [ B[%d,%d] ]\n",i,j)
    }
    retn(C)

}

少し長いですが、大したことはしていません。

これをファイル matrix として保存しておきます。

これを使えば、先の演算例はもうすこしシンプルになります。

まず、第1の例です。

% cat x
:r prog/sigmoid
:r matrix

X  = { 1.0 , 0.5 }
W1 = {{ 0.1 , 0.3 , 0.5 },{ 0.2 , 0.4 , 0.6 }}
B1 = { 0.1 , 0.2 , 0.3 }

A1 = madd(mdot(X,W1),B1)
p(A1)

Z1 = sigmoid(A1)
p(Z1)

% tt x
A1[0]=0.300000,A1[1]=0.700000,A1[2]=1.100000
Z1[0]=0.574443,Z1[1]=0.668188,Z1[2]=0.750260

次に2番目の例です。

% cat y
:r matrix
:r prog/sigmoid

X  = { 1.0 , 0.5 }
W1 = {{ 0.1 , 0.3 , 0.5 },{ 0.2 , 0.4 , 0.6 }}
B1 = { 0.1 , 0.2 , 0.3 }

W2 = {{ 0.1 , 0.4 },{ 0.2 , 0.5 },{ 0.3 , 0.6 }}
B2 = { 0.1 , 0.2 }

A1 = madd(mdot(X,W1),B1)
p(A1)

Z1 = sigmoid(A1)
p(Z1)

A2 = madd(mdot(Z1,W2),B2)
p(A2)

Z2 = sigmoid(A2)
p(Z2)

% tt y
A1[0]=0.300000,A1[1]=0.700000,A1[2]=1.100000
Z1[0]=0.574443,Z1[1]=0.668188,Z1[2]=0.750260
A2[0]=0.516160,A2[1]=1.214027
Z2[0]=0.626249,Z2[1]=0.771011

そして、3番目の例です。

% cat z
:r matrix
:r prog/sigmoid
:r prog/identity

X  = { 1.0 , 0.5 }
W1 = {{ 0.1 , 0.3 , 0.5 },{ 0.2 , 0.4 , 0.6 }}
B1 = { 0.1 , 0.2 , 0.3 }

W2 = {{ 0.1 , 0.4 },{ 0.2 , 0.5 },{ 0.3 , 0.6 }}
B2 = { 0.1 , 0.2 }

W3 = {{ 0.1 , 0.3 },{ 0.2 , 0.4 }}
B3 = { 0.1 , 0.2 }

A1 = madd(mdot(X,W1),B1)
p(A1)

Z1 = sigmoid(A1)
p(Z1)

A2 = madd(mdot(Z1,W2),B2)
p(A2)

Z2 = sigmoid(A2)
p(Z2)

A3 = madd(mdot(Z2,W3),B3)
p(A3)

Y  = identity(A3)
p(Y)

% tt z
A1[0]=0.300000,A1[1]=0.700000,A1[2]=1.100000
Z1[0]=0.574443,Z1[1]=0.668188,Z1[2]=0.750260
A2[0]=0.516160,A2[1]=1.214027
Z2[0]=0.626249,Z2[1]=0.771011
A3[0]=0.316827,A3[1]=0.696279
Y[0]=0.316827,Y[1]=0.696279

スクリプトの記述が大部シンプルかつ直感的になりました。

最後は、再び処理を書籍の様に、まとめて書いてみます。

% cat a
:r matrix
:r prog/sigmoid
:r prog/identity

net["W1"] = {{ 0.1 , 0.3 , 0.5 },{ 0.2 , 0.4 , 0.6 }}
net["B1"] = { 0.1 , 0.2 , 0.3 }
net["W2"] = {{ 0.1 , 0.4 },{ 0.2 , 0.5 },{ 0.3 , 0.6 }}
net["B2"] = { 0.1 , 0.2 }
net["W3"] = {{ 0.1 , 0.3 },{ 0.2 , 0.4 }}
net["B3"] = { 0.1 , 0.2 }

def forward(NET,X){
    W1 = NET["W1"]; W2 = NET["W2"]; W3 = NET["W3"];
    B1 = NET["B1"]; B2 = NET["B2"]; B3 = NET["B3"];

    A1 = madd(mdot(X ,W1),B1)
    Z1 = sigmoid(A1)
    A2 = madd(mdot(Z1,W2),B2)
    Z2 = sigmoid(A2)
    A3 = madd(mdot(Z2,W3),B3)
    retn( Y = identity(A3) )
}

X = { 1.0 , 0.5 }
Y = forward(net,X)
p(Y)

% tt a
Y[0]=0.316827,Y[1]=0.696279

大部スクリプトがすっきりしました。

とりあえず、ここまで。

ツインテールdeエンジェルモード!! で Deep Learning!! （3層の計算）

ツインテールdeエンジェルモード!! で3層ニューラルネットワークの計算をしてみます。対象とする3層ニューラルネットワークは、書籍の図3.15をご参照下さい。

まずは、入力層から第1層までを計算してみますが、その前に、

関数 madd() のこの行が間違っていました。すいません。

# 行ベクトル＋行ベクトル → 行ベクトル
    if( (ai==1       ) && (bi==1       ) ){
        if(aj!=bj) retn(NULL)
        loop(j<aj){ C[j] = A[j] + B[j] }
        retn(C)
    }

関数 msub() のこの行が間違っていました。すいません。

# 行ベクトルー行ベクトル → 行ベクトル
    if( (ai==1       ) && (bi==1       ) ){
        if(aj!=bj) retn(NULL)
        loop(j<aj){ C[j] = A[j] - B[j] }
        retn(C)
    }

戻って、入力層から第1層までの計算です。

% cat x
:r prog/madd
:r prog/mmul
:r prog/sigmoid

X  = { 1.0 , 0.5 }
W1 = {{ 0.1 , 0.3 , 0.5 },{ 0.2 , 0.4 , 0.6 }}
B1 = { 0.1 , 0.2 , 0.3 }

A1 = madd( mmul(X,1,2,W1,2,3),1,3 , B1,1,3 )
p(A1)

Z1 = sigmoid( A1 )
p(Z1)

% tt x
A1[0]=0.300000,A1[1]=0.700000,A1[2]=1.100000
Z1[0]=0.574443,Z1[1]=0.668188,Z1[2]=0.750260

続いて、第1層から第2層までの計算です。

% cat y
:r prog/madd
:r prog/mmul
:r prog/sigmoid

X  = { 1.0 , 0.5 }
W1 = {{ 0.1 , 0.3 , 0.5 },{ 0.2 , 0.4 , 0.6 }}
B1 = { 0.1 , 0.2 , 0.3 }

W2 = {{ 0.1 , 0.4 },{ 0.2 , 0.5 },{ 0.3 , 0.6 }}
B2 = { 0.1 , 0.2 }

A1 = madd( mmul(X ,1,2,W1,2,3),1,3 , B1,1,3 )
p(A1)

Z1 = sigmoid( A1 )
p(Z1)

A2 = madd( mmul(Z1,1,3,W2,3,2),1,2 , B2,1,2 )
p(A2)

Z2 = sigmoid( A2 )
p(Z2)

% tt y
A1[0]=0.300000,A1[1]=0.700000,A1[2]=1.100000
Z1[0]=0.574443,Z1[1]=0.668188,Z1[2]=0.750260
A2[0]=0.516160,A2[1]=1.214027
Z2[0]=0.626249,Z2[1]=0.771011

続いて、第2層から出力層までの計算です。

% cat z
:r prog/madd
:r prog/mmul
:r prog/sigmoid
:r prog/identity

X  = { 1.0 , 0.5 }
W1 = {{ 0.1 , 0.3 , 0.5 },{ 0.2 , 0.4 , 0.6 }}
B1 = { 0.1 , 0.2 , 0.3 }

W2 = {{ 0.1 , 0.4 },{ 0.2 , 0.5 },{ 0.3 , 0.6 }}
B2 = { 0.1 , 0.2 }

W3 = {{ 0.1 , 0.3 },{ 0.2 , 0.4 }}
B3 = { 0.1 , 0.2 }

A1 = madd( mmul(X ,1,2,W1,2,3),1,3 , B1,1,3 )
p(A1)

Z1 = sigmoid( A1 )
p(Z1)

A2 = madd( mmul(Z1,1,3,W2,3,2),1,2 , B2,1,2 )
p(A2)

Z2 = sigmoid( A2 )
p(Z2)

A3 = madd( mmul(Z2,1,2,W3,2,2),1,2 , B3,1,2 )
p(A3)

Y  = identity( A3 )
p(Y)

% tt z
A1[0]=0.300000,A1[1]=0.700000,A1[2]=1.100000
Z1[0]=0.574443,Z1[1]=0.668188,Z1[2]=0.750260
A2[0]=0.516160,A2[1]=1.214027
Z2[0]=0.626249,Z2[1]=0.771011
A3[0]=0.316827,A3[1]=0.696279
Y[0]=0.316827,Y[1]=0.696279

ようやく答え（配列Y）が出ました。

次は、処理を書籍の様に、まとめて書いてみます。

% cat a
:r prog/madd
:r prog/mmul
:r prog/sigmoid
:r prog/identity

net["W1"] = {{ 0.1 , 0.3 , 0.5 },{ 0.2 , 0.4 , 0.6 }}
net["B1"] = { 0.1 , 0.2 , 0.3 }
net["W2"] = {{ 0.1 , 0.4 },{ 0.2 , 0.5 },{ 0.3 , 0.6 }}
net["B2"] = { 0.1 , 0.2 }
net["W3"] = {{ 0.1 , 0.3 },{ 0.2 , 0.4 }}
net["B3"] = { 0.1 , 0.2 }

def forward(NET,X){
    W1 = NET["W1"]; W2 = NET["W2"]; W3 = NET["W3"];
    B1 = NET["B1"]; B2 = NET["B2"]; B3 = NET["B3"];

    A1 = madd( mmul(X ,1,2,W1,2,3),1,3 , B1,1,3 )
    Z1 = sigmoid( A1 )
    A2 = madd( mmul(Z1,1,3,W2,3,2),1,2 , B2,1,2 )
    Z2 = sigmoid( A2 )
    A3 = madd( mmul(Z2,1,2,W3,2,2),1,2 , B3,1,2 )
    retn( Y = identity( A3 ) )
}

X = { 1.0 , 0.5 }
Y = forward(net,X)
p(Y)

% tt a
Y[0]=0.316827,Y[1]=0.696279

同じ結果が出ました。（よかった。）

なお、関数 identity() は、次の様な（入力をそのまま返す）関数です。

def identity(x){
    if( type(x)=='A' ){
        each( k=keys(x) )
            ans[k] = x[k]
    }
    else
            ans    = x
    retn(ans)
}

とりあえず、ここまで。

ツインテールdeエンジェルモード!! で多次元配列の計算（その2）

ツインテールdeエンジェルモード!! で、多次元配列計算用の関数を作ってみます。

2つのデータの和差積商を計算する関数を作成します。

なお、2つのデータは、スカラーでも、ベクトルでも、行列でもいいものとします。

まずは積。

/*--1---2---3---4---5---6---7---8---9---A---B---C---D---E---F---G---H---I---J---K---L---M---N---*/
# 行列A(ai行*aj列)と行列B(bi行*bj列)の積を計算する。
# 注意1：A又はBは、スカラー値や行ベクトルであっても良い。なお、スカラー値の場合は0行*0列、行ベクトルの場合は1行*N列と指定する。
# 注意2：行ベクトルの場合は、1次元配列となるので別処理が必要となる。
# 注意3：列ベクトルの場合は、2次元配列となるので別処理は不要である。（行列と同一処理でOK!!）
/*--1---2---3---4---5---6---7---8---9---A---B---C---D---E---F---G---H---I---J---K---L---M---N---*/

def mmul(A,ai,aj,B,bi,bj){

# パラメータの補完と確認
    if( type(A)!='A' ){ ai=aj=0 } elif( ai<1||aj<1 ){ retn(NULL) }
    if( type(B)!='A' ){ bi=bj=0 } elif( bi<1||bj<1 ){ retn(NULL) }

######################################
### 【スカラー値】が含まれるケース ###
######################################

# スカラー値＊スカラー値 → スカラー値
    if( (ai==0&&aj==0) && (bi==0&&bj==0) ){
        retn(C=A*B)
    }

# スカラー値＊行ベクトル → 行ベクトル
    if( (ai==0&&aj==0) && (bi==1       ) ){
        loop(j<bj)
            C[j] = A * B[j]
        retn(C)
    }

# 行ベクトル＊スカラー値 → 行ベクトル
    if( (ai==1       ) && (bi==0&&bj==0) ){
        loop(j<aj)
            C[j] = A[j] * B
        retn(C)
    }

# スカラー値＊２次元行列 → ２次元行列
    if( (ai==0&&aj==0)                   ){
        loop(i<bi)
            loop(j<bj)
                C[i,j] = A * B[i,j]
        retn(C)
    }

# ２次元行列＊スカラー値 → ２次元行列
    if(                   (bi==0&&bj==0) ){
        loop(i<ai)
            loop(j<aj)
                C[i,j] = A[i,j] * B
        retn(C)
    }

######################################
### 【行ベクトル】が含まれるケース ###
######################################

# 行ベクトル＊行ベクトル｛転置｝ → スカラー値
    if( (ai==1       ) && (bi==1       ) ){
        if(aj!=bj) retn(NULL)
        C=0; loop(t<bj){ C += A[t] * B[t] }
        retn(C)
    }

# 行ベクトル＊２次元行列 → 行ベクトル
    if( (ai==1       )                   ){
        if(aj!=bi) retn(NULL)
        loop(j<bj){
            C[j]=0; loop(t<bi){ C[j] += A[t] * B[t,j] }
        }
        retn(C)
    }

# ２次元行列＊行ベクトル｛転置｝ → 行ベクトル｛転置｝
    if(                   (bi==1       ) ){
        if(aj!=bj) retn(NULL)
        loop(i<ai){
            C[i]=0; loop(t<bj){ C[i] += A[i,t] * B[t] }
        }
        retn(C)
    }

######################################
### それ以外の全ケース
######################################

# ２次元行列＊２次元行列 → ２次元行列
    if(aj!=bi) retn(NULL)
    loop(i<ai){
        loop(j<bj){
            C[i,j]=0; loop(t<bi){ C[i,j] += A[i,t] * B[t,j] }
        }
    }
    retn(C)

}

次に、商。

/*--1---2---3---4---5---6---7---8---9---A---B---C---D---E---F---G---H---I---J---K---L---M---N---*/
# 行列A(ai行*aj列)と行列B(bi行*bj列)の商を計算する。
# 注意1：A又はBは、スカラー値や行ベクトルであっても良い。なお、スカラー値の場合は0行*0列、行ベクトルの場合は1行*N列と指定する。
# 注意2：行ベクトルの場合は、1次元配列となるので別処理が必要となる。
# 注意3：列ベクトルの場合は、2次元配列となるので別処理は不要である。（行列と同一処理でOK!!）
/*--1---2---3---4---5---6---7---8---9---A---B---C---D---E---F---G---H---I---J---K---L---M---N---*/

def mdiv(A,ai,aj,B,bi,bj){

# パラメータの補完と確認
    if( type(A)!='A' ){ ai=aj=0 } elif( ai<1||aj<1 ){ retn(NULL) }
    if( type(B)!='A' ){ bi=bj=0 } elif( bi<1||bj<1 ){ retn(NULL) }

######################################
### 【スカラー値】が含まれるケース ###
######################################

# スカラー値／スカラー値 → スカラー値
    if( (ai==0&&aj==0) && (bi==0&&bj==0) ){
        retn(C=A/B)
    }

# スカラー値／行ベクトル → 行ベクトル
    if( (ai==0&&aj==0) && (bi==1       ) ){
        loop(j<bj)
            C[j] = A / B[j]
        retn(C)
    }

# 行ベクトル／スカラー値 → 行ベクトル
    if( (ai==1       ) && (bi==0&&bj==0) ){
        loop(j<aj)
            C[j] = A[j] / B
        retn(C)
    }

# スカラー値／２次元行列 → ２次元行列
    if( (ai==0&&aj==0)                   ){
        loop(i<bi)
            loop(j<bj)
                C[i,j] = A / B[i,j]
        retn(C)
    }

# ２次元行列／スカラー値 → ２次元行列
    if(                   (bi==0&&bj==0) ){
        loop(i<ai)
            loop(j<aj)
                C[i,j] = A[i,j] / B
        retn(C)
    }

######################################
### それ以外の全ケース
######################################

    retn(NULL)

}

次に、和。

/*--1---2---3---4---5---6---7---8---9---A---B---C---D---E---F---G---H---I---J---K---L---M---N---*/
# 行列A(ai行*aj列)と行列B(bi行*bj列)の和を計算する。
# 注意1：A又はBは、スカラー値や行ベクトルであっても良い。なお、スカラー値の場合は0行*0列、行ベクトルの場合は1行*N列と指定する。
# 注意2：行ベクトルの場合は、1次元配列となるので別処理が必要となる。
# 注意3：列ベクトルの場合は、2次元配列となるので別処理は不要である。（行列と同一処理でOK!!）
/*--1---2---3---4---5---6---7---8---9---A---B---C---D---E---F---G---H---I---J---K---L---M---N---*/

def madd(A,ai,aj,B,bi,bj){

# パラメータの補完と確認
    if( type(A)!='A' ){ ai=aj=0 } elif( ai<1||aj<1 ){ retn(NULL) }
    if( type(B)!='A' ){ bi=bj=0 } elif( bi<1||bj<1 ){ retn(NULL) }

######################################
### 【スカラー値】が含まれるケース ###
######################################

# スカラー値＋スカラー値 → スカラー値
    if( (ai==0&&aj==0) && (bi==0&&bj==0) ){
        retn(C=A+B)
    }

# スカラー値＋行ベクトル → 行ベクトル
    if( (ai==0&&aj==0) && (bi==1       ) ){
        loop(j<bj)
            C[j] = A + B[j]
        retn(C)
    }

# 行ベクトル＋スカラー値 → 行ベクトル
    if( (ai==1       ) && (bi==0&&bj==0) ){
        loop(j<aj)
            C[j] = A[j] + B
        retn(C)
    }

# スカラー値＋２次元行列 → ２次元行列
    if( (ai==0&&aj==0)                   ){
        loop(i<bi)
            loop(j<bj)
                C[i,j] = A + B[i,j]
        retn(C)
    }

# ２次元行列＋スカラー値 → ２次元行列
    if(                   (bi==0&&bj==0) ){
        loop(i<ai)
            loop(j<aj)
                C[i,j] = A[i,j] + B
        retn(C)
    }

######################################
### 【行ベクトル】が含まれるケース ###
######################################

# 行ベクトル＋行ベクトル → 行ベクトル
    if( (ai==1       ) && (bi==1       ) ){
        if(aj!=bj) retn(NULL)
        loop(i<ai){ C[i] = A[i] + B[i] }
        retn(C)
    }

# それ以外の組み合わせ
    if( (ai==1       ) || (bi==1       ) )
        return(NULL);

######################################
### それ以外の全ケース
######################################

# ２次元行列＋２次元行列 → ２次元行列
    if(ai!=bi||aj!=bj) retn(NULL)
    loop(i<ai){
        loop(j<aj){
            C[i,j] = A[i,j] + B[i,j]
        }
    }
    retn(C)

}

次に、差。

/*--1---2---3---4---5---6---7---8---9---A---B---C---D---E---F---G---H---I---J---K---L---M---N---*/
# 行列A(ai行*aj列)と行列B(bi行*bj列)の差を計算する。
# 注意1：A又はBは、スカラー値や行ベクトルであっても良い。なお、スカラー値の場合は0行*0列、行ベクトルの場合は1行*N列と指定する。
# 注意2：行ベクトルの場合は、1次元配列となるので別処理が必要となる。
# 注意3：列ベクトルの場合は、2次元配列となるので別処理は不要である。（行列と同一処理でOK!!）
/*--1---2---3---4---5---6---7---8---9---A---B---C---D---E---F---G---H---I---J---K---L---M---N---*/

def msub(A,ai,aj,B,bi,bj){

# パラメータの補完と確認
    if( type(A)!='A' ){ ai=aj=0 } elif( ai<1||aj<1 ){ retn(NULL) }
    if( type(B)!='A' ){ bi=bj=0 } elif( bi<1||bj<1 ){ retn(NULL) }

######################################
### 【スカラー値】が含まれるケース ###
######################################

# スカラー値ースカラー値 → スカラー値
    if( (ai==0&&aj==0) && (bi==0&&bj==0) ){
        retn(C=A-B)
    }

# スカラー値ー行ベクトル → 行ベクトル
    if( (ai==0&&aj==0) && (bi==1       ) ){
        loop(j<bj)
            C[j] = A - B[j]
        retn(C)
    }

# 行ベクトルースカラー値 → 行ベクトル
    if( (ai==1       ) && (bi==0&&bj==0) ){
        loop(j<aj)
            C[j] = A[j] - B
        retn(C)
    }

# スカラー値ー２次元行列 → ２次元行列
    if( (ai==0&&aj==0)                   ){
        loop(i<bi)
            loop(j<bj)
                C[i,j] = A - B[i,j]
        retn(C)
    }

# ２次元行列ースカラー値 → ２次元行列
    if(                   (bi==0&&bj==0) ){
        loop(i<ai)
            loop(j<aj)
                C[i,j] = A[i,j] - B
        retn(C)
    }

######################################
### 【行ベクトル】が含まれるケース ###
######################################

# 行ベクトルー行ベクトル → 行ベクトル
    if( (ai==1       ) && (bi==1       ) ){
        if(aj!=bj) retn(NULL)
        loop(i<ai){ C[i] = A[i] - B[i] }
        retn(C)
    }

# それ以外の組み合わせ
    if( (ai==1       ) || (bi==1       ) )
        return(NULL);

######################################
### それ以外の全ケース
######################################

# ２次元行列ー２次元行列 → ２次元行列
    if(ai!=bi||aj!=bj) retn(NULL)
    loop(i<ai){
        loop(j<aj){
            C[i,j] = A[i,j] - B[i,j]
        }
    }
    retn(C)

}

そして、これらの関数をファイル mmul mdiv madd msub に、

それぞれ個別保存しておきます。ついでに、ファイルはディレクトリ prog に

入れておきます。

すると、先日やった演算は、こんな感じに単純化することができました。

まずは、1番目の演算です。

% cat x
:r prog/mmul

A={{1,2},{3,4}}
B={{5,6},{7,8}}

Y=mmul(A,2,2,B,2,2)
Y

% tt x
Y[0,0]=19,Y[0,1]=22,Y[1,0]=43,Y[1,1]=50

2番目の演算です。

% cat y
:r prog/mmul

A={{1,2,3},{4,5,6}}
B={{1,2},{3,4},{5,6}}

Y=mmul(A,2,3,B,3,2)
Y

% tt y
Y[0,0]=22,Y[0,1]=28,Y[1,0]=49,Y[1,1]=64

3番目の演算です。

% cat z
:r prog/mmul

A={{1,2},{3,4},{5,6}}
B={7,8}

Y=mmul(A,3,2,B,1,2)
Y

% tt z
Y[0]=23,Y[1]=53,Y[2]=83

最後の演算です。

% cat n
:r prog/mmul

X={1,2}
W={{1,3,5},{2,4,6}}

Y=mmul(X,1,2,W,2,3)
Y

% tt n
Y[0]=5,Y[1]=11,Y[2]=17

とりあえず、ここまで。

ツインテールdeエンジェルモード!! で Deep Learning!! （多次元配列の計算）

ツインテールdeエンジェルモード!! で多次元配列の計算をしてみます。

まずは、1次元配列リテラルの記述。

% tt
tt> a={1,2,3,4}
a[0]=1,a[1]=2,a[2]=3,a[3]=4

次に、2次元配列リテラルの記述。

% tt

tt> b={{1,2},{3,4},{5,6}}
b[2,0]=5,b[2,1]=6,b[0,0]=1,b[0,1]=2,b[1,0]=3,b[1,1]=4

次に、行列の積。

まずは、2X2行列 と 2X2行列の積。意外と難しい。

a = {{1,2},{3,4}}
b = {{5,6},{7,8}}

I=2 ; N=2 ; J=2

# (IxN)*(NxJ)
    loop(i<I){
        loop(j<J){
            c[i,j] = 0; loop(n<N){ c[i,j] += a[i,n]*b[n,j] }
        }
    }
p(c)

結果は、

c[0,0]=19,c[0,1]=22,c[1,0]=43,c[1,1]=50

もう1つ。

次は、2X3行列 と 3X2行列の積。

a = {{1,2,3},{4,5,6}}
b = {{1,2},{3,4},{5,6}}

I=2 ; N=3 ; J=2

# (IxN)*(NxJ)
    loop(i<I){
        loop(j<J){
            c[i,j] = 0; loop(n<N){ c[i,j] += a[i,n]*b[n,j] }
        }
    }
p(c)

結果は、

c[0,0]=22,c[0,1]=28,c[1,0]=49,c[1,1]=64

最後に、3X2行列 と 2X1行列の積。

a = {{1,2},{3,4},{5,6}}
b = {{7},{8}}

I=3 ; N=2 ; J=1

# (IxN)*(NxJ)
    loop(i<I){
        loop(j<J){
            c[i,j] = 0; loop(n<N){ c[i,j] += a[i,n]*b[n,j] }
        }
    }
p(c)

結果は、

c[2,0]=83,c[0,0]=23,c[1,0]=53

行列の積は、関数を組んだ方がいいかもしれない。

とりあえず、こんな感じの関数を組んでみました。

# 行列A(aiXaj)と行列B(biXbj)の積を計算する。
def dot(A,ai,aj,B,bi,bj){
    if(aj!=bi) retn(NULL)
    loop(i<ai){
        loop(j<bj){
            C[i,j] = 0; loop(n<aj){ C[i,j] += A[i,n]*B[n,j] }
        }
    }
    retn(C)
}

この関数があると、、、（ファイル dot にセーブ）

最初の例は、こんなに簡単になりました。

:r dot
a = {{1,2},{3,4}}
b = {{5,6},{7,8}}
c = dot(a,2,2,b,2,2)
p(c)

結果は、

c[0,0]=19,c[0,1]=22,c[1,0]=43,c[1,1]=50

2番目の例も、こんなに簡単になりました。

:r dot
a = {{1,2,3},{4,5,6}}
b = {{1,2},{3,4},{5,6}}
c = dot(a,2,3,b,3,2)
p(c)

結果は、

c[0,0]=22,c[0,1]=28,c[1,0]=49,c[1,1]=64

3番目の例も、こんなに簡単になりました。

:r dot
a = {{1,2},{3,4},{5,6}}
b = {{7},{8}}
c = dot(a,3,2,b,2,1)
p(c)

結果は、

c[2,0]=83,c[0,0]=23,c[1,0]=53

最後に、ニューラルネットワークの内積の例です。

これは、最初の行列Xがベクトルなので、

関数dot（）もそれに対応するように変更します。

変更後の関数dot()です。

# 行列A(aiXaj)と行列B(biXbj)の積を計算する。
def dot(A,ai,aj,B,bi,bj){
    if(aj!=bi) retn(NULL)
    if(ai==01){        # 行列Aはベクトル -> 結果もベクトル
        loop(j<bj){
            C[  j] = 0; loop(n<aj){ C[  j] += A[  n]*B[n,j] }
        }
        retn(C)
    }
    loop(i<ai){
        loop(j<bj){
            C[i,j] = 0; loop(n<aj){ C[i,j] += A[i,n]*B[n,j] }
        }
    }
    retn(C)
}

スクリプト本体です。

:r dot
X = {1,2}
W = {{1,3,5},{2,4,6}}
Y = dot(X,1,2,W,2,3)
p(Y)

結果です。

Y[0]=5,Y[1]=11,Y[2]=17

とりあえず、ここまで。