ベクトル空間の元と成分表示について

ここ最近,ベクトル空間というものによく出くわす.
そして,ここ最近ずっと悩まされている.つらい.が,楽しくもある.
出来の悪い僕は,誰もが簡単に通過できる部分にひどく苦しめられる.(決してベクトル空間が簡単な概念であると言いたいわけではない.)

この記事はベクトル空間の元について,掘り下げて考えてみたことをまとめたものである.
(正確には,まとめたevernoteのノートをちょっと変えてポストしているわけだが.)
ベクトル空間の定義はベクトル空間 – Wikipediaに任せるとしよう.
ここで考えたいのは,その元についてである.
僕はベクトル空間の元は,その名前から,ベクトルであるとなんとなく思っていたが,少し考えればそうとは限らないということがわかる.
院試のための勉強のときに,N次正方行列の全体はベクトル空間である,数列全体はベクトル空間である,などといった問題に遭遇したことを思い出したためである.
そうすると,ベクトル空間の元は,ベクトルであるとは限らないという結論に達する.
しかしながら,多くの本では,ベクトル空間の元は必ずベクトルであるかのように書かれている.
ここに,僕の誤解釈による大きな勘違いが生じていたことに気がつく.

さて,ベクトル空間の元とはなんであろうか?
結論から言えば,大きな勘違いは,元と成分表示を同一視していたために生じたと言える.
例として,2次正方行列の全体を考えると,これはベクトル空間である.
そうすると,この は2次正方行列である.当たり前の話である.
ここで,基底が登場する.基底については基底 (線型代数学) – Wikipediaに任せるとしよう.
基底に,{[1 0; 0 0], [0 1; 0 0], [0 0; 1 0], [0 0; 0 1]}をとれば,元の成分表示は[x1, x2, x3, x4]で表現できる.
そうすると,以下の結論が得られると思うのだが,どうだろうか.

ベクトル空間の元はベクトルでなくてもよいが,ある基底に対する元の成分表示はベクトルで書ける.

これは僕にとっては,頭の中をすっきりさせてくれる文である.
ここで,ある本にこう書いてあったのを思い出す.

ベクトル空間Vnの基底にBを固定する

僕はこれをすっと通り過ぎたために,勘違いが発生したのだ.
今やこの文の意味するところは,ベクトル空間Vnの基底にBを固定すると,その元の成分表示はベクトルで書けるので,Vnの元の成分表示全体の元はベクトルである,という意味だと勝手に考えている.
ここからもう一つ言えるのは,

基底を構成する要素はベクトルでなくても良い

ということだ.上の例では基底を構成する要素は2次正方行列である.
しかしながら,基底の条件を見ると,線形独立性と,全域性が挙げられている.
そうすると,これら2つの性質はベクトルに限ったことではないと考えられる.

整理するとベクトル空間Vnを考える.
ここに基底Bを固定する(もちろんBの構成要素はベクトルでなくてよい).
Vnの元の基底Bによる成分表示の全体をB(Vn)とすると,B(Vn)の元はベクトルである.
この段階で,”Vnの元”と, B(Vn)の元 を同一視して考えることにすると,暗に仮定しているのではないか.
なので,Vnの元は必ずベクトルで表されるのだ(なぜならこの段階でVnの元,というとB(Vn)の元を指しているのだから).

もう一つ見えてくる.
僕は,”何かをベクトルで表現する”という言葉を好んで使うが,この時僕は暗にベクトル空間Vnを仮定し,ある基底Bによる成分表示の全体B(Vn)を考えていたということに気がついた.
何も考えずに1段階だと思っていたステップが,実は2段階だったのだ.
これには驚いた.ささやかな喜びだ.

さて,B(Vn)の元は必ずベクトルであり,n次元ベクトル全体はベクトル空間であるので,B(Vn)はベクトル空間である.
B(Vn)の元の要素数(次元数)はBの要素数,すなわちVnの次元,すなわちnである.
今,実数のみを考えることにすると,n次元実数ベクトル全体はRnと書かれることが多い.
そうすると,B(Vn)はRnの部分集合であるとも言えそうだ.
直感的には部分空間である,と言えそうだが,B(Vn)が和とスカラー倍に対して閉じているか検証していないのでここでは控えておく.

以上で終わる.
つっこみどころは満載かもしれない.
最初から元がベクトルであるベクトル空間を考えているのだ,と見なせば楽に読めるのかもしれない.
しかしながら,ベクトル空間としていろいろなことを表現できるのに,上のように見なすと,少しむなしい気がする.
以後もう少し掘り下げてみる予定だ.

広告

簡単なclassification

Python-nltkを使って簡単なclassificationをやっていきたいと思います.
今回使うのはロイターコーパスです(nltk.corpus.reuters).
このコーパスにはクラス(カテゴリ)がいくつかありますが,とりあえず2クラスのclassificationをやってみたいので,適当に2クラス選びます.

import nltk
len(nltk.corpus.reuters.categories())

でクラス数を調べると, 90 クラスあるみたいです.
次に,

from nltk.corpus import reuters
for category in sorted(reuters.categories(), key=lambda x:len(reuters.fileids(categories=x)), reverse=True):
    print category,
    print len(reuters.fileids(categories=category))

で各クラスに属する文書数を調べます.どれを選びましょうか.笑
なるべく2つのクラスが関係なさそうで,かつ文書数が同じくらいのを選びましょうか.
“ship”(船)は286件,”wheat”(小麦)は283件で,互いにあまり関係なさそうなのでこれにします.
さて,クラスが決まりました.次に教師ありか教師なしかをどうするかですね.今回はとりあえず教師ありにしましょう.
というわけでデータを訓練データとテストデータに分けます.
ロイターコーパスではfileidを使って分けることができます(もっと良い方法を知っている方は教えてください).
以下で訓練データとテストデータを分けます(正確には訓練データのfileidたちとテストデータのfileidたち).

from nltk.corpus import reuters
class1, class2 = "ship", "wheat"

train_fileids1 = [fileid for fileid in reuters.fileids(categories=class1) if "train" in fileid]
train_fileids2 = [fileid for fileid in reuters.fileids(categories=class2) if "train" in fileid]
test_fileids1 = [fileid for fileid in reuters.fileids(categories=class1) if "test" in fileid]
test_fileids2 = [fileid for fileid in reuters.fileids(categories=class2) if "test" in fileid]

print len(train_fileids1), len(test_fileids1)
print len(train_fileids2), len(test_fileids2)

こんな感じですかね.
shipでは訓練データが197件,テストデータが89件,
wheatでは訓練データが212件,テストデータが71件となりました.

さて,次はこれらをベクトルで表現します.
どう表現するかは重要ですが,言い出すときりがないし,今回の趣旨から外れるので,単純に単語の頻度とします.
ここで言う頻度は,ある文書での単語iの出現回数を,総単語数で除したものです.
ここから先はコードを観て頂いた方が理解しやすいと思います.説明が下手なので….
ベクトルの次元をそろえるために,共通の単語集合Vを用いて各文書をベクトルで表現することとします.
というわけでVを求めます.なお,Vはトレーニングデータに出てくる単語とします.

V = set()
for fileid in train_fileids1+train_fileids2:
    V = V.union(set(reuters.words(fileid)))
print len(V)

これでVを求めると,Vは8127単語からなる集合となりました.
ここで,ステミングや見出し語化をするかは皆さんにおまかせするとして,ストップワードくらいは外しておこうと思います.
それと,大文字小文字を区別しないようにします.

sws = stopwords.words("english")
V = set()
for fileid in train_fileids1+train_fileids2:
    V = V.union(set([word.lower() for word in reuters.words(fileid) if word.lower() not in sws]))
print len(V)

最終的にはこんな感じですかね.以下でもいいですね.

sws = stopwords.words("english")
V = set(word.lower() for fileid in train_fileids1+train_fileids2 for word in reuters.words(fileid) if word.lower() not in sws)
print len(V)

そうすると,Vは6459単語からなる集合となりました.
つまり,これから扱う各ベクトルは6459次元になります.
featuresという関数を用いて,reuters.words(fileid)で与えられる単語のリスト(以下words)をベクトルに変換することとします.
というわけでfeaturesを書いていきます.
featuresはwordsとVを引数に取り,Vの基数次元のベクトル(list)を返す関数とします.

def features(words, V):
    v = dict.fromkeys(V, 0.0)
    for word in words:
        try:
            v[word.lower()] += 1
        except:
            pass
    return [val/len(words) for val in v.values()]

こんな感じになりました.この関数featuresを用いて,各文書をベクトルで表現します.
もちろん,これでできたベクトルは疎なベクトルです.

さあ,文書をベクトルで表現できたので,いよいよ今回の本題であるclassificationです.
今回の記事では,最も単純な方法を実装したいと思います.
優れたアルゴリズムを実感するために,あえて今回は簡単な方法を実装するのです.
今回実装する方法は非常にシンプルで,各クラスの代表ベクトルを求めて,それに近い方に分類するという方法です.
どんなベクトルを代表ベクトルにするかは,いろいろありますが,今回は単純に平均ベクトルにしたいと思います.
なお,以下からpython-numpyを使います.

import numpy
train_vectors1 = numpy.array([numpy.array(features(reuters.words(fileid), V)) for fileid in train_fileids1])
train_vectors2 = numpy.array([numpy.array(features(reuters.words(fileid), V)) for fileid in train_fileids2])
rep1 = train_vectors1.mean(axis=0)
rep2 = train_vectors2.mean(axis=0)

こんな感じになりました.rep1は”ship”の代表ベクトルで,rep2は”wheat”の代表ベクトルです.
さて,いよいよテストデータを分類していきますが,
rep1,rep2のどちらに近いか,あるいは似ているかをどのように決めたらよいでしょうか.
シンプルなのはユークリッド距離やコサイン類似度を使うことでしょう.
今回はコサイン類似度を用います.1に近いほど2文書は似ているということになります.
具体的には,predictという関数を用いて,ベクトルで表現された文書を0か1に写像しようと思います.
predictはvector, rep1, rep2を引数に取り,rep2よりrep1に近ければ1を,rep1よりrep2に近ければ0を返します.

def predict(v, rep1, rep2):
    def simcos(v1, v2):
        return v1.dot(v2) / numpy.sqrt(v1.dot(v1))*numpy.sqrt(v2.dot(v2))
    sim1 = simcos(v, rep1)
    sim2 = simcos(v, rep2)
    if sim1 >= sim2:
        return 1
    return 0

こんな感じになりました.これを使うと,

test_vectors1 = numpy.array([numpy.array(features(reuters.words(fileid), V)) for fileid in test_fileids1])
test_vectors2 = numpy.array([numpy.array(features(reuters.words(fileid), V)) for fileid in test_fileids2])
predicted = [predict(v, rep1, rep2) for v in test_vectors1] + [predict(v, rep1, rep2) for v in test_vectors2]
print predicted

sol = [1 for fileid in test_fileids1] + [0 for fileid in test_fileids2]
print sol

このように書けまして,結果は以下のようになりました.

[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

非常におもしろいですね笑.なんと1つだけしか1がない,つまり”ship”に属する文書は1つしかないと言っているわけです.
一応精度を出しておきます.

sol = numpy.array(sol)
accuracy =  float(len([i for i in predicted-sol if i == 0])) / len(sol)
print accuracy

精度は, 0.45 でした.
「半分くらいはうまく分類できているのか」という印象を受けますが,もし,正例と負例の数が同じくらいなら全部片方に分類してもこのくらいの精度が出てしまうんですね.
そこで,もう1つの評価指標,F値を導入します.

def F_measure(predicted, solution):
    tp, tn, fp, fn = 0.0, 0.0, 0.0, 0.0
    for p, s in zip(predicted, solution):
        if p == 1:
            if p == s:
                tp += 1
            else:
                fp += 1
        else:
            if p == s:
                tn += 1
            else:
                fn += 1
    precision = tp / (tp + fp)
    recall = tp / (tp + fn)
    return 2 * precision * recall / (precision + recall)

これで評価すると,今回の分類結果は 0.02 ととても低くなります.
F値は横着者を許してはくれないんですね笑.

Summary

今回の分類では,ほとんど”wheat”に分類しちゃったんですね.
まったく横着者です笑.
次回からこの結果を改善していくことで,各アルゴリズムのすごさを実感したいと思います.

「正しい」証明とは何か

最近「正しい」証明とは何か,というのを考えてしまっています. これは数学基礎論とかで研究されている分野なのだろうか. 例えば,本を読んでいて,こうこうこうだからこうなる. なんて簡単に書かれていますが,本当にそうなの?って思っちゃうんです. どういう根拠でそうなの?って思っちゃうんですよね. どうしようもないですね笑.

まだいろいろ調べている段階ですが,普通の本なんかでは公理系っていうのを暗に仮定している気がします. 公理系の中には,公理があって導出規則なんかもあって,どっちかというと論理学に近いのでしょうか. 例えば,pならばqっていう形の定理を証明している本では,もちろんpを仮定してqが出て来たらpならばqを導出して良いという演繹という導出規則を仮定しているでしょう. で,そんなことを書いているスペースはないから察してくれよ,と. まぁそうですよね.だってその本は数学基礎論の本ではないし,論理学の本でもないし,そんな当たり前のところに割くページなんかないですよね. というわけでいつまでたってもぼんや〜り雲に覆われたような数学観が養われてきた僕です. しかしですね.これに立ち入るのはかなり勇気いると思うんですよね. 少し考えただけでノイローゼになりそうです.笑 なので,どういう公理系があってどういう導出規則があるのか,そういうところだけ拾って来て, 「ああ,この本のここではこういう規則を使っているんだな」,というのを意識しながら本を読むことにしようかなと思っています. 実際何の抵抗もなく受け入れられる導出規則なんかは,本の中で文章として使われていても,多分読んでて違和感を感じないんですよね. ん?って思うのはだいたい「任意の」とか「存在する」とか背理法とか. この辺は「全称例化」「全称汎化」「存在例化」「存在汎化」を抑えれば良いと思うんですよね. で,証明を読むときには,そこに書いてある事実が,何を「仮定」して成り立っているのかを見極めるようにする. 言い換えればどんな仮定の上にその事実は成り立っているのか,を見極めながら読む. 数学のできる人って,これらを全て頭の中でやっているんだろうけど,僕は多分できないので,どんな仮定の上に成り立っているのか,きちんと書こうと思います. そうして全てを追えれば,最初に書いた疑問なんかは浮かんでこないんだろうと思います.

\epsilon - \delta論法なんかそのいい例ですよね. そんなわけでここに飛び込んだんですよね. そうするとまたいろんなことが気になって来て,最終的に解析学を復習しよう,なんてことになってきました. なんですかね. 機械学習とか勉強してても,フワフワしてる感じがすごいんですよね. フワフワしてると不安じゃないですか. だからここらで時間かけて復習しようと思ったんですね. そんなわけでそういうポストも増えるかもしれない,いやそういうポストもするぞ!という意思表示のポストですこれは笑.

Convex Optimization / Stephen Boyd & Lieven Vandenberghe

タイトルの本を読んでいます.
意味不明なところやわからないところをポストすると思われます. この本は以下で入手できます.

発端は凸計画問題の論文を読むにあたって,知識が足りなさすぎてわけがわからなかったからです. 専門用語が多すぎるのもあって英語の本を探したらこれが無料だったので,読み始めました.
目標は凸計画問題はどんなものかを知ることです. あと,ラグランジュの未定乗数法で解を求めることができますが, 「何をやっているんだろう」という感じが半端ないので,その辺の謎も解きたかったです. 凸性とは何か,何がうれしいのか,絶えない疑問をどんどん解消していきたいです.

最近思ったのですが,「これ,なぜだ?」と思ったとき,賢い人は「ひとまず置いておき,次に進む」ということができるのだろうけど, 僕みたいに賢くもないのにその疑問が頭から離れない人は「さっさと解消してしまう」のが結局一番早い気がしています.

半空間(half spaces)

半空間(Half spaces)はアフィン集合でないが,凸集合であることを示す.PDFファイルで出しました.

不等式の扱いって大変ですね.悩みました. 自分なりに書いてみたのですが,おかしいところがないか不安です.

最近数学ばっかりやってるからプログラミングが恋しいなぁ.

最急降下法(gradient descent)

勾配法の下るバージョンの最急降下法(gradient descent)を実装してみた. ちなみに上るバージョンは最急上昇法(gradient ascent)と呼ばれる. 最近書いている微分を近似するという内容の記事は実はこれの準備のためなのでした.笑
導関数や勾配を明示的に与えることの無い実装を,したくてしたくてたまらなかったのです.
このアルゴリズムにより極小値を求めることができる. 通常勾配を明示的に与える必要があるが,最近書いている勾配の近似で代用した. つまり明示的に勾配を与える必要がない.
実装は以下.

 1: # f(x1, ..., xn)のxiによる偏導関数を返す
 2: def partial(f, n, i, dx=0.001):
 3:     I = np.identity(n)
 4:     def part(x):
 5:         return (f(x+I[i]*dx) - f(x)) / dx
 6:         #return (f(x+I[i]*2*dx) - f(x)) / 2*dx
 7:     return part
 8: 
 9: # f(x1, ..., xn)の勾配を返す
10: def Gradient(f, n, dx=0.01):
11:     def g(x):
12:         return np.array([G[i](x) for i in np.arange(n)], dtype=float)
13:     G = [partial(f, n, i, dx=dx) for i in np.arange(n)]
14:     return g
15: 
16: # n変数関数fを最小にするxを勾配法により求める. 
17: def gradient_descent(f, n, h0=0.03, dx=1./10**2):
18:     fx = Gradient(f, n, dx=dx)
19:     pp = np.random.random_sample(n)
20:     g = fx(pp)
21:     p = pp - h0*g
22:     h = h0
23:     while (p-pp).dot(p-pp) > 1./10**12:
24:         ph = h; h = 2 * ph
25:         while f(p-h*g) < f(p-ph*g):
26:             ph = h; h *= 2
27:         h /= 2; pp = p
28:         g = fx(p); p = p - h*g; h = h0
29:     return p

実は,この本を読んでいます.

これなら分かる最適化数学―基礎原理から計算手法まで

この本の例題3.2に出てくる関数

\displaystyle  f(x,y) = x^{3} + y^{3} - 9xy + 27

に上記のアルゴリズムを適用します.

n = 2
f = lambda x: x[0]**3 + x[1]**3 - 9*x[0]*x[1] + 27
solution = gradient_descent(f, n)
print solution, f(solution)

実行結果は以下のようになりました.

[ 2.98999076  2.98999076] 0.000899659139463

解は(x,y) = (3, 3)なので,良い感じだと思われます.
この例題はニュートン法の問題ですが,実はニュートン法も実装してみたのでまた記事を書きたいと思います. いやあ,楽しいなぁ.