Soft Confidence-Weighted Learningを実装してみた

すごく久しぶりに書きます。いしたーです。

きょうは機械学習のおはなしです。


機械学習といえば、最近はDeep Learningが非常にアツいですね。

Googleがねこの画像を学習させたとかなんとか。

でも、Deep Learningはいろんなデメリットを抱えています。

  • 学習が遅い
  • パラメータの調整が難しい

などです。


Soft Confidence-Weighted Learning (SCW)という手法をご存知でしょうか？

SCWとは、逐次学習が可能な線形分類器の一種です。

かずー氏のブログが詳しいので詳細はそちらに譲るとします。(SCWについてはかずー氏にいろいろとお世話になりました。この場でお礼を述べておきます。)

ブログからそのまま引用すると、SCWは

  • 学習が高速
  • オンライン学習(逐次学習)が可能
  • ノイズに強い
  • 精度が良い

というメリットを持っています。

さらに、パラメータの調整も大雑把でいいそうです。Deep Learningは学習は遅いが複雑な分類が得意であるという特徴を持っているので、SCWとは正反対ですね。


では実装に移ります。実装と言っても、論文に書いてある数式をそのままコードに落としこむだけなので案外簡単です。

import numpy as np
from scipy.stats import norm

class SCW(object):
    def __init__(self, ndim, C=1.0, ETA=1.0):
        self.weights = np.zeros(ndim)
        self.covariance = np.ones(ndim)
        self.C = np.float64(C)
        self.cdf_values = self.calc_cdf_values(ETA)
    
    def sgn(self, x):
        t = np.dot(self.weights, x)
        if(t > 0):
            return 1
        if(t == 0):
            return 0
        if(t < 0):
            return -1

    def loss(self, x, teacher):
        t = teacher*np.dot(self.weights, x)
        if(t >= 1):
            return 0
        return 1-t

    def calc_cdf_values(self, ETA):
        phi = norm.cdf(ETA)
        psi = 1 + np.power(phi, 2)/2
        zeta = 1 + np.power(phi, 2)
        return (phi, psi, zeta)

    def calc_confidence(self, x, teacher):
        return np.dot(x, self.covariance*x)
    
    def calc_margin(self, x, teacher):
        return teacher*np.dot(self.weights, x)
    
    def calc_alpha(self, x, teacher):
        #calc in a child class
        pass

    def calc_beta(self, x, teacher):
        alpha = self.calc_alpha(x, teacher)
        v = self.calc_confidence(x, teacher)
        m = self.calc_margin(x, teacher) 
        phi, psi, zeta = self.cdf_values

        j = -alpha * v * phi
        k = np.sqrt(np.power(alpha*v*phi, 2) + 4*v)
        u = np.power(j+k, 2) / 4
        return (alpha * phi) / (np.sqrt(u) + v*alpha*phi)

    def update_covariance(self, x, teacher):
        beta = self.calc_beta(x, teacher)
        c = self.covariance
        self.covariance -= beta*c*c*x*x 

    def update_weights(self, x, teacher):
        alpha = self.calc_alpha(x, teacher)
        self.weights += alpha*teacher*self.covariance*x

    def update(self, x, teacher):
        y = self.sgn(x)
        if(self.loss(x, teacher) > 0):
            self.update_weights(x, teacher)
            self.update_covariance(x, teacher)

    def train(self, X, teachers):
        for x, teacher in zip(X, teachers):
            self.update(x, teacher)
        return self.weights, self.covariance


class SCW1(SCW):
    def calc_alpha(self, x, teacher):
        v = self.calc_confidence(x, teacher)
        m = self.calc_margin(x, teacher) 
        phi, psi, zeta = self.cdf_values
        
        j = np.power(m, 2) * np.power(phi, 4) / 4
        k = v * zeta * np.power(phi, 2)
        t = (-m*psi + np.sqrt(j+k)) / (v*zeta)
        return min(self.C, max(0, t))


class SCW2(SCW):
    def calc_alpha(self, x, teacher):
        v = self.calc_confidence(x, teacher)
        m = self.calc_margin(x, teacher) 
        phi, psi, zeta = self.cdf_values
        
        n = v+1/self.C
        a = np.power(phi*m*v, 2)
        b = 4*n*v * (n+v*np.power(phi, 2))
        gamma = phi * np.sqrt(a+b)
        
        c = -(2*m*n + m*v*np.power(phi, 2))
        d = np.power(n, 2) + n*v*np.power(phi, 2)
        t = (c+gamma)/(2*d)
        return max(0, t)

SCWにはSCW-ⅠとSCW-Ⅱがあるので両方実装してみました。両者の違いはαを計算する部分だけなので、αを計算している部分以外を親クラスとしてくくり出し、子クラスにαの計算を実装しています。

少し戸惑ったのは論文のこの部分。
f:id:sonicair:20140422150517p:plain
cumulative functionってナニ？と思ったのですが、調べたらすぐにわかりました。
Cumulative distribution function - Wikipedia, the free encyclopedia
ここではscipyのcdfという関数を使って計算しています。

ソースコードはGitHubに置いておきました。付属のtest.pyを実行すると動作を確認することができます。

論文に擬似コードが書いてあるぐらいなので実装はけっこう簡単です。他の言語で実装してみたり、実際にSCWを応用したアプリケーションを作ってみたりしてはいかがでしょうか。

HaskellでRSA暗号の鍵生成器作った

学校の課題でなんか作ってこいって言われたんで 

RSA暗号の鍵生成器を作ってみました。

ソースコードはGithubに置いてあります。

 

大まかな構造はこんなかんじです

乱数生成 -> 乱数をもとに素数を生成 -> 素数をもとに鍵ペアを生成 -> ファイルに出力

 

ではモジュールごとにソースコードを解説していきます 

まず乱数生成から。

module Random
(genSeed
,genInfiniteRandomList
)where

import System.Entropy
import System.Random
import Data.Time.Clock.POSIX

genSeed :: IO Integer
genSeed = do
    mbstr <- getEntropy 128
	let str' = show mbstr	
	let str = map toEnum $ filter (\x -> 48<=x && x<58) $ map fromEnum $ str':: [Char]
	let randomNum = read str :: Integer
	return randomNum

genInfiniteRandomList :: (Integer, Integer) -> IO [Integer]
genInfiniteRandomList (n1, n2) = do
 g <- getPOSIXTime
 let randNum = randomRs (n1, n2) (mkStdGen (round g)) :: [Integer]
 return $ map toInteger randNum

genSeedはpとqのseedを生成します。getEntropyで取得した文字列から数字だけを抜き出してIntegerに変換するというちょっと強引な方法を使ってます。

getEntropyは/dev/urandomから乱数を取得しているので暗号用の乱数としては脆弱です。誰か/dev/randomから取ってこられるように書きなおしてくださいおねがいします。

genInfiniteRandomListはその名の通り乱数の無限リストを生成します。こちらは素数判定に使うのでセキュリティは考慮する必要がありません。

 

次は素数判定ですがその前に素数判定と鍵生成に使う計算用モジュールを載せます。ほぼ全部コピペです。extEuclidに至っては動作原理さえよくわかってません。

module CalcAlgorithm
(power
,powMod
,inverse
) where

power :: (Integral a1, Num a) => a -> a1 -> a
power x n
    | n == 0    = 1
    | even n    = power (x*x) (n `div` 2)
    | otherwise = x * power x (n - 1)
 
powMod :: Integer -> Integer -> Integer -> Integer
powMod base ex m
  | ex == 0   = 1
  | even ex   = square (powMod base (ex `div` 2) m) `mod` m
  | otherwise = (base * (powMod base (ex - 1) m)) `mod` m
 where
   square x = x * x

extEuclid :: Integral t => t -> t -> (t, t)
extEuclid a b = recurse a b 1 0 0 1
  where recurse a 0 x0 y0 x1 y1 = (x0, y0)
        recurse a b x0 y0 x1 y1 = let q = a `div` b in 
                                  let r = a `mod` b in
                                  recurse b r x1 y1 (x0 - q * x1) (y0 - q * y1)

inverse :: Integer -> Integer -> Integer
inverse x n = (fst (extEuclid x n)) `mod` n

 

素数判定です。鍵の素となるpとqを生成します。

module PrimalityTest
(isMillerRabinPrime
,isProbablePrime
,getPrime
) where

import Control.Applicative
import Control.Monad
import System.Random
import Data.Time.Clock.POSIX
import CalcAlgorithm
import Random

factor2 :: (Integral t, Num t1) => t -> (t, t1)
factor2 x = factor2' (x,0)
 where
  factor2' (n,s)
   | even n = factor2' (n `div` 2, s+1)
   | otherwise = (n,s)

isProbablePrime :: Integer -> Integer -> Bool
isProbablePrime n a =
 let (d,s) = factor2 (n-1) in
 let p = powMod a d n /= 1 in
 let q = not $ elem (n-1) $ map (\r -> powMod a ((power 2 r)*d) n) [0..s-1] in
 not (p && q)


isMillerRabinPrime :: Integer -> IO Bool
isMillerRabinPrime n = do
 ns <- genInfiniteRandomList (1,n-1) 
 return . and $ map (isProbablePrime n) (take 12 ns)

getPrime :: Integer -> IO Integer
getPrime s = do
 if odd s then getPrime' s else getPrime (s+1)
  where getPrime' s = do b <- isMillerRabinPrime s
                         if b then return s else getPrime' (s+2)

MillerRabin素数判定法を使っています。こちらとWikipediaの記事を参考にしました。

MillerRabin素数判定法をざっと解説します。

ある数をMillerさんに渡します。

ある数が素数なら、Millerさんは「そいつは素数だね！！」と言います。

ある数が合成数なら、Millerさんは「そいつは合成数だよ」と言います。

でもMillerさんはドジっ子なので、合成数を持っていったときにたまに間違えて「そいつは素数だね！！」と言っちゃうことがあります。

そんなちょっとおっちょこちょいなアルゴリズムです。

 

Millerさんはたまに間違えるので、

「ねぇ、この数って素数？」

「うん、多分ね」

「ほんとに素数？こいつが合成数だったらどうなるかわかってんでしょうね」

「きっと素数だよ...」

 

といった感じに、Millerさんに何回も聞くことによって精度を上げています。別にいじめているわけじゃありません。

このコードではMillerさんに12回問い詰めています。Millerさんがかわいそうですね。

 

最後に、鍵を計算して出力する部分なんですが、はてなブログのエラーのせいでうまく表示できません。すみませんが代わりにGithubを見てください。

鍵のファイルは他のプログラムから読みやすいようにCSV形式にしています。

鍵を読み込んで暗号化・復号化するPython製のプログラムもGithubに上げてあります。Pythonはよくわからないのでググりながら適当に書きました。上記のコードを含めおかしいところがあったら言ってください。

 

あと、このプログラムを使って暗号化したら解読されて秘密が漏れた！なんてことがあっても僕は責任を負えません。「俺ら今暗号使って通信してんだぜ〜」って感じで遊び程度に使ってくれたら幸いです。

 

7/23

getEntropyで取ってきた乱数をハッシュ関数に通すようにしました。アップグレード済みのソースコードはGithubにあります。(1行と5文字しか追加してない)

P≠NP予想とは何かをまとめた

プログラミングを始めてある程度の時間が経つとどうしても「P≠NP予想」という言葉を目にする。しかし実際にどういうものなのかを調べてみても全くわからない。

Wikipediaを見ても

http://ja.wikipedia.org/wiki/P%E2%89%A0NP%E4%BA%88%E6%83%B3

クラスPって何？Javaのクラスと違うの？

とか

多項式時間って何だよ

とか思ってしまって...

とにかく全くわからない。

 

んで

学校の図書館に数学ガールがあったので手にとってみると

この問題についてある程度詳しく書かれていたのでまとめてみる。

 

 

P≠NP予想とは何か

Wikipediaをそのまま引用すると

http://ja.wikipedia.org/wiki/P≠NP予想

P≠NP予想（P≠NPよそう、英: P is not NP）は、計算複雑性理論（計算量理論）におけるクラスPとクラスNPが等しくないという予想である。

だそうだ

 

ではクラスPとは何かというと

「多項式時間で解が求められる問題の集合のこと」

 

では多項式時間とは何か

多項式時間で解ける問題とは

解くのにかかる時間がO(n^c)である問題のこと

それだけ

 

改めてクラスPとは何かというと

「O(n^c)で解が求められる問題の集合」のことである。

 

 

 

では今度はクラスNPとは何かというと

「yes となる証拠が与えられたとき、その証拠が本当に正しいかどうかを多項式時間で判定できる問題の集合」

のことだ。

http://ja.wikipedia.org/wiki/NP

要は解の候補が与えられた時にそれが正しい解なのかどうかを多項式時間で”判定できる”問題のこと。

これについて注意すべきことがある。

解が”求められる”のではなく、解の候補が正しいかどうかを多項式時間で”判定できる”問題であるということと、

NPはNon-deterministic Polynomial time（非決定性多項式時間）の略であり、決して「解が求められない問題」ではないということだ。

 

 

 

P問題がNP問題に含まれることはすでに証明されている。

言い換えると

「多項式時間で解を求められる問題は、解の候補が与えられたときその解が正しいかどうかを多項式時間で判定可能である」

ということがすでに証明済みであるということだ。

 

ではその逆はどうか。

NP問題はP問題に含まれるか。

これがP≠NP予想である。

 

P問題がNP問題に含まれることは証明済みなので、

仮にNP問題がP問題に含まれた場合、P=NPだといえる。

そうでないのならP≠NPだといえる。

 

 

P≠NP予想を解く鍵の一つに

NP完全問題

というものがある。

NP完全問題のうち一つでもP問題であるものが存在すれば、その時点でP=NPなんだそうだ。

つまり、NP完全問題のうち一つでもO(n^c)で解を求められるものが存在すればその時点でP=NPだということ。

NP完全問題は世界中にたくさんある。

充足可能性問題やハミルトン閉路問題などがそうだ。

このうち一つでもO(n^c)で解ければコンピュータ世界や数学界などいろいろな分野での大発見なのだが

世界中のコンピュータ科学者たちがいくら挑んでもNP問題をO(n^c)で解くことができないので

P≠NPだと予想されているということである。

 

この記事のほとんどは数学ガールの第4弾、「乱択アルゴリズム」で説明されているので読んでみるといい。

 

NP完全問題の定義を私は理解できていないのでこれから勉強する。理解でき次第記事にする。

今日はここまで。