2014-02-28

Understanding How cbind/rbind works in R

Rを使うようになってきた最大の動機はR自身で柔軟にデータを作成できることにあります。例えば、パッケージユーザーのための機械学習(1)：決定木 - 銀座で働くデータサイエンティストのブログでは以下のようなXORパターンのデータを用いて決定木の解説を行っています。

「少しランダムが入った市松模様」とでもいうのでしょうか。上述のものは「シンプルなXORパターン」と称されていました。少し点群の数を増やすと市松模様を基調としている様子がはっきりしてきます。

自分も是非、こういうデータを短時間でゴリゴリ作れるようになりたいと思いました。先のエントリでは、ありがたいことにデータを作成した際のRの実装を掲載してくださっていました。

# シンプルなXORパターン
> p11<-cbind(rnorm(n=25,mean=1,sd=0.5),rnorm(n=25,mean=1,sd=0.5))
> p12<-cbind(rnorm(n=25,mean=-1,sd=0.5),rnorm(n=25,mean=1,sd=0.5))
> p13<-cbind(rnorm(n=25,mean=-1,sd=0.5),rnorm(n=25,mean=-1,sd=0.5))
> p14<-cbind(rnorm(n=25,mean=1,sd=0.5),rnorm(n=25,mean=-1,sd=0.5))
> t<-as.factor(c(rep(0,50),rep(1,50)))
> d1<-as.data.frame(cbind(rbind(p11,p13,p12,p14),t))
> names(d1)<-c("x","y","label")

...これで出来るのか、という印象。初見では何をやっているか全く検討もつきませんでした。

部分的に実行してコードリーディングを行った結果、cbindとrbindという関数が何回も登場することに気付きました。

> help(cbind)
Combine R Objects by Rows or Columns

Description:

     Take a sequence of vector, matrix or data frames arguments and
     combine by _c_olumns or _r_ows, respectively.  These are generic
     functions with methods for other R classes.

Usage:

     cbind(..., deparse.level = 1)
     rbind(..., deparse.level = 1)

                                   :

なるほど、cbindもrbindもベクトルや行列を結合する関数のようです。cとrはそれぞれ、ColumnとRowを表すとのこと。要素数の少ないベクトルや行列で実験を行った結果、なかなか一貫性のある設計であることに気付きました。

まずベクトルのcbindとrbindから考えてみることにします。ベクトルaとbを以下のように置きます。

Rでもベクトルを定義しましょう。

> a = c(1, 2, 3)
> b = c(4, 5, 6)
> a
[1] 1 2 3
> b
[1] 4 5 6

これらのベクトルをcbindしたものは行列となります。これをM_cとします。M_cは以下のような、2つのベクトルを結合した行列となります。

> cbind(a, b)
     a b
[1,] 1 4
[2,] 2 5
[3,] 3 6

同様に、ベクトルをrbindしたものをM_rとします。M_rの結果は以下のようになりました。

> rbind(a, b)
  [,1] [,2] [,3]
a    1    2    3
b    4    5    6

次に行列のcbind/rbindを考えてみましょう。以下のような行列SとTを考えます。

Rでも行列を定義しましょう。

> S = matrix(c(1,  4, 2,  5, 3,  6), 2, 3)
> T = matrix(c(7, 10, 8, 11, 9, 12), 2, 3)
> S
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
> T
     [,1] [,2] [,3]
[1,]    7    8    9
[2,]   10   11   12

ベクトルの場合と同様に、行列をcbindしたものをM_cとします。M_cは以下のような行列となります。

> cbind(S, T)
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    2    3    7    8    9
[2,]    4    5    6   10   11   12

行列をrbindしたものをM_rとします。M_rは以下のような行列となります。

> rbind(S, T)
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
[3,]    7    8    9
[4,]   10   11   12

さて、どうでしょう? とても一貫していますね。引数がベクトルである場合でも行列である場合でもcbind(a, b)は以下のような行列を作成するのです。

同様にrbind(a, b)は以下のような行列を作成します。

ベクトルのcbind/rbindには少し注意が必要です - ベクトルが列ベクトルとして扱われるか行ベクトルとして扱われるかは関数によって異なります。また、cbindで以下のように合成することは出来ません(このような結合を行うにはc(a, b)とします)。

さて、ここまでのcbind/rbindの知識を元に、今一度XORパターンの実装を読んでみましょう(再掲)。

# シンプルなXORパターン
> p11<-cbind(rnorm(n=25,mean=1,sd=0.5),rnorm(n=25,mean=1,sd=0.5))
> p12<-cbind(rnorm(n=25,mean=-1,sd=0.5),rnorm(n=25,mean=1,sd=0.5))
> p13<-cbind(rnorm(n=25,mean=-1,sd=0.5),rnorm(n=25,mean=-1,sd=0.5))
> p14<-cbind(rnorm(n=25,mean=1,sd=0.5),rnorm(n=25,mean=-1,sd=0.5))
> t<-as.factor(c(rep(0,50),rep(1,50)))
> d1<-as.data.frame(cbind(rbind(p11,p13,p12,p14),t))
> names(d1)<-c("x","y","label")

さて、もう読めるようになりましたね...と簡単に行けばよいのですが、やはりまだ難解ですね。

この実装では結合した行列から最終的にデータフレームを生成していますが、行列を作成することを目標としてみましょう。as.data.frame(...)の中身を作成するところまで、と考えてください。

少し分解しましょう。まず、p11〜p14までのデータを考えてみます。これらはX座標、Y座標共に標準偏差0.5の正規分布に従ったランダムな点群です。X座標、Y座標のそれぞれのベクトルとして用意し、cbindを用いて合成をしています。以下のようなイメージです。

meanを変えることで分布の中心を変えながら、p12、p13、p14を同様に作成します。二つ目の添字が象限を表しています。例えば、p13は第三象限に分布した点群です。そのため、ほとんどの値は負の値となります。

さて、ベクトルを合成した結果は行列となるのでした。これをrbindで合成します。市松模様とするために少し工夫を施し、p11、p13、p12、p14の順番で合成しましょう。

これとは別に0と1のラベルからなるベクトルを作成します。先ほど少し触れましたが、ベクトルとベクトルを直列に合成するにはc()を用います。

最後にcbindを用いて行列とベクトルを合成します。

これで目的の行列を作成することができました。データフレームに変換しなくても行列の段階でプロットすることも可能です。

plot(M[,1:2], col=M[,3], pch=16, xlab='x', ylab='y')

まとめ

ベクトルや行列等を合成して新たな行列を作成する場合はcbind/rbindを使うと便利です。Rにはまだ慣れていないためどんな場合でもデータフレームにまとめようとする傾向にありますが、今後はcbind/rbindをどんどん活用していきたいと考えています。

もちろん好みによるとは思いますが、個人的にはmatrixを使うよりもrbindを使ったほうが直感的に行列が記述出来るよう感じています。

> S = rbind(c(1, 2, 3), c(4, 5, 6))
> S
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6

以下にmatrixを使った行列の作成を再掲します。是非比較してみてください。

> S = matrix(c(1, 4, 2, 5, 3, 6), 2, 3)

2014-02-28

Understanding How cbind/rbind works in R

agw2014-02-28

2014-02-27

How matrix and vector are multiplied in R?

構文的な一貫性をなかなか感じられないのでRを敬遠していました。しかし時勢には敵わず、使用頻度が増えてきました。

どうせやるのであれば手が増えたほうがよいので、スニペットを多数作成し検証しています。気付いたことをブログに少しずつまとめようかと思っています。

さて、今日はvector型とmatrix型のオブジェクトの乗算を検証しました。行列とベクトルの演算は空気のように扱えないとストレスが貯まるので、どの処理系を学ぶときにも比較的早い段階で検証することにしています。

例えばベクトルの向きを調べることは重要です。vector型が列ベクトルとして解釈されるのであるか、行ベクトルとして解釈されるかで話は大きく違ってきます。

具体例を挙げてみましょう。3 × 3の行列をM、要素数3のベクトルをvとします。行列Mとベクトルvを乗じたとき、ベクトルが列ベクトルであれば、結果も列ベクトルになります。列ベクトルは行列の右側からかけます。

この演算方法を基調とする処理系をPre-multiplyな処理系などと呼びます。右辺がベクトルであることを簡潔に示すためにv'を導入しました。v'は以下のように計算されます。

これに対しベクトルが行ベクトルであれば行列の左側からかけないといけませんし、結果も行ベクトルになります。

この演算方法を基調としている処理系をPost-multiplyな処理系と呼びます。v'の要素は以下のように計算します。そのまま記述すると水平方向に長くなってしまうため右辺のベクトルを転置していることに注意してください。

Pre-multiplyな処理系とPost-multiplyな処理系で同様の結果を得るにはベクトルの方向だけでなく、行列を転置しなければいけないことに注意してください。

例えば同時座標系での移動行列を考えてみましょう。簡単に、移動行列とはあるベクトルの要素を平行移動する行列である、と考えてください。まずPre-multiplyな処理系での移動行列を考えてみることにします。

この行列は二次元ベクトル(x, y)を(T_x, T_y)だけ平行移動します。

さて、移動行列を適用することによって(x, y)が(x + T_x, y + T_y)に射影されました。これをPost-multiplyな処理系でもやってみましょう。Post-multiplyな処理系では行ベクトルを用いますが、要素の結果は等しくならなければなりません。

試しにPre-multiplyな処理系の移動行列をPost-multiplyな処理系で使ってみましょう。

何かめちゃくちゃなことになってしまいました。Pre-multiplyな処理系の行列をPost-multiplyな処理系で用いるには転置をしなければならなかったのです。

今回はうまくいっているようですね。つまり、Post-multiplyな処理系での移動行列は以下のようなものなのでした。Pre-multiplyな処理系のものと見比べて転置されていることを確認するのもよいでしょう。

ここまでで書いたように、処理系でベクトルが列ベクトルで扱われるか行ベクトルで扱われるかで行列の意味合いが全く異ります。そのため、ベクトルが列ベクトルなのか行ベクトルなのか予め調べておくことはとても重要です。言葉を変えれば、処理系でベクトルが列ベクトルで扱われるか行ベクトルで扱われるかを調べることによって、その処理系がPre-multiplyを基調として処理系かPost-multiplyを基調とした処理系かを見極めることができます。

さて、Rのvector型オブジェクトは以下のように記述します。

> v <- c(1, 2, 3)
> v
[1] 1 2 3

同様に、Rのmatrix型オブジェクトは以下のように記述します。

> M <- matrix(c(1, 0, 0, 0, 1, 0, 3, 5, 1), 3, 3)
> M
     [,1] [,2] [,3]
[1,]    1    0    3
[2,]    0    1    5
[3,]    0    0    1

Rはvector型とmatrix型の*演算子による乗算を許していますが、これは望む乗算ではありませんでした。調べてみたところ、%*%演算子を使えば望む行列ベクトル演算をするようです。早速やってみましょう。Rの処理系はPre-multiplyな処理系であると仮定し、(x, y) = (0, 0)を(T_x, T_y) = (3, 5)で移動してみましょう。

> M
     [,1] [,2] [,3]
[1,]    1    0    3
[2,]    0    1    5
[3,]    0    0    1
> v
[1] 0 0 1
> M %*% v
     [,1]
[1,]    3
[2,]    5
[3,]    1

数式で表すと以下のようになります。

さて、うまく移動出来たようですね。どうやらRの行列は列ベクトルを基調としていると言えそうです。また、戻り値がmatrix型であることにも注意しましょう。

念のため、Post-multiplyな処理系での演算も試してみましょう。数式で表すと以下です。

> M <- t(M)
> M
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    3    5    1
> v %*% M
     [,1] [,2] [,3]
[1,]    3    5    1

驚くことに、これもPost-multiply系であるかのように動作してしまいました。この結果から考えるに、Rで行列ベクトル演算を行う際、ベクトルは以下のように扱われることが分かります。

行列 %*% ベクトル
- ベクトルは列ベクトルであると解釈される
ベクトル %*% 行列
- ベクトルは行ベクトルであると解釈される

とても柔軟な処理であるようにも見えますが少し注意が必要だな、と改めて感じました。

まとめ

行列とベクトルの乗算を行うには%*%演算子を用いる
Rでのベクトルは文脈によって行ベクトルとも列ベクトルとも扱われる
RはPre-multiply前提の処理系でもPost-multiply前提の処理系でもなさそうだ

ベクトルの向きは、例えばLaTeXでベクトルを記述する際にも影響します。結果の通りRの処理系は予想以上に柔軟であり決定打には欠けますが、今後ベクトルをLaTeXで表記するときは列ベクトルを基調とし、行列演算はPre-multiplyな処理系を前提にしようと考えています。その理由はmatrixオブジェクトを生成するのに列ベクトルを並べるように要素を記述しなければならなかったことにあります。もっとも、Rの言語全体の一貫性のなさはよく話題になるようですので、最終的には列ベクトル + Pre-multiplyな処理系であるように考えるのが好みであることも少なからず影響しています。

2014-02-27

How matrix and vector are multiplied in R?

agw2014-02-27

2014-02-04

How to Find Palindromic Substrings?

Algorithm !

回文に関連する実装が苦手です。

回文について考えなければならないときはいつも気が重かったのですが、与えられた文字列から回文となる部分文字列を列挙する方法を@nico_shindanninさんがニコニコ生放送で解説していらっしゃいました。とても分かりやすい解説だったので、メモとしてまとめました。

次の文字列を考えてみます。

例えば先頭から5文字の部分文字列、CBABCは回文になっています。同様に、どのくらいの部分文字列が回文になっているでしょう? まずナイーブなO(N³)を実装して確かめてみました。

def palindromes(s):
    size = len(s)
    for i in range(size):
        for j in range(i, size):
            l, u = i, j
            while l <= u:
                if s[l] != s[u]:
                    break
                l += 1
                u -= 1
            else:
                yield s[i:j+1]

これはとても直感的な実装です。

部分文字列を列挙し
回文であるかどうか

を調べています。

この実装の動作を図示すると次のようになります。回文となる文字列の数は意外と多く、9個ありました。

さて、続いては@nico_shindanninさんが紹介していた方法です。回文の性質をうまく利用して、部分文字列の中央から広げて行くようにチェックを行うことによって計算量をO(N²)に落としています。

def palindromes(s):
    size = len(s)
    for i in range(size):
        l = u = i
        while 0 <= l and u < size:
            if s[l] != s[u]:
                break
            yield s[l:u+1]
            l -= 1
            u += 1

        l, u = i, i + 1
        while 0 <= l and u < size:
            if s[l] != s[u]:
                break
            yield s[l:u+1]
            l -= 1
            u += 1