Random Projectionを試す

はじめに

言語処理を行う場合、単語数を考えると高次元スパースなベクトルを扱うことが多い。
次元削減を行える手法の一つである、Random Projectionを試してみる。

Random Projectionとは

  • 乱数を要素に持ち、各列ベクトルの大きさが1である行列Rを用意して、行列Xをかけることで次元を落とすことができる
    • X_rp = R * X
  • また、このRの各要素がN(0,1)の正規乱数の場合、各列ベクトル間のユークリッド距離をできるだけ保ったまま、次元削減できることが証明されている
  • この乱数行列Rの作り方として、以下が提案されている
    • Rの各要素r_ijについて、以下の近似を用いる
    • 1/6の確率で、r_ij = sqrt(3)
    • 2/3の確率で、r_ij = 0
    • 1/6の確率で、r_ij = -sqrt(3)

準備

ドキュメント群からcos類似度の近い文書を検索するということを、次元削減した行列に対してできるかどうか試してみる。
D : d*N行列 (dは単語種類数, Nはドキュメント数)
R : k*d行列 (kは次元削減後の次元数)

#形態素に分割されたドキュメント部分だけ取り出す
$ cut -f2- -d" " matome.data > data 

#単語の種類だけ取り出す
$ sed 's/ /\
/g' data|LC_ALL="C" sort |LC_ALL="C" uniq > wordlist

#下記で扱う形式(単語数 単語1 単語2...)に直す
$ echo 1987 > docs
$ cat data| perl -e 'while(<>){chomp; $t=$_; @sp=split(/\s/); print scalar(@sp)," ",$t,"\n"; }' >> docs

コード

#include <iostream>
#include <vector>
#include <cmath>
#include <map>
#include <fstream>

class Matrix {
  int m, n;
  std::vector<double> val;

  void splice(int i, int j){
    int N = val.size();
    if(i+(j-1)<N){
      std::vector<double> vl(N-j, 0.0);
      int vl_idx = 0, val_idx = 0;
      while(val_idx<N){
        if(val_idx<i||i+j-1<val_idx){
          vl[vl_idx] = val[val_idx];
          vl_idx++;
        }
        val_idx++;
      }
      val = vl;
    }
  }
  Matrix clone(){
    Matrix ret(m, n);
    for(int i=0; i<m; i++){
      for(int j=0; j<n; j++){
        ret.setVal(i, j, val[n*i+j]);
      }
    }
    return ret;
  }
  
public:
  Matrix():m(0),n(0){}
  Matrix(int m, int n):m(m),n(n){
    if(m>0 && n>0){
      for(int i=0; i<m*n; i++){
        val.push_back(0.0);
      }
    }
  }

  int getM() const { return m; }
  int getN() const { return n; }
  double getVal(int i, int j) const { return val[n*i+j]; }
  void setVal(int i, int j, double x){ val[n*i+j] = x; }
  bool isSquare(){ return n==m; }

  Matrix add(const Matrix& mat){
    if(m == mat.m && n == mat.n){
      Matrix ret(m, n);
      for(int i=0; i<m; i++){
        for(int j=0; j<n; j++){
          ret.setVal(i, j, val[n*i+j] + mat.getVal(i,j));
        }
      }
      return ret;
    }
    return Matrix(-1, -1);
  }
  Matrix operator+(const Matrix& mat){ return add(mat); }

  Matrix sub(const Matrix& mat){
    if(m == mat.m && n == mat.n){
      Matrix ret(m, n);
      for(int i=0; i<m; i++){
        for(int j=0; j<n; j++){
          ret.setVal(i, j, val[n*i+j] - mat.getVal(i,j));
        }
      }
      return ret;
    }
    return Matrix(-1, -1);
  }
  Matrix operator-(const Matrix& mat){ return sub(mat); }

  Matrix prod(const Matrix& mat){
    if(n == mat.m){
      Matrix ret(m, mat.n);
      for(int i=0; i<m; i++){
        for(int j=0; j<mat.n; j++){
          double d = 0.0;
          for(int k=0; k<n; k++){
            d += val[n*i+k] * mat.getVal(k,j);
          }
          ret.setVal(i, j, d);
        }
      }
      return ret;
    }
    return Matrix(-1, -1);
  }
  Matrix operator*(const Matrix& mat){ return prod(mat); }

  void time(double x){
    for(int i=0; i<m; i++){
      for(int j=0; j<n; j++){
        val[n*i+j] *= x;
      }
    }
  }

  Matrix transpose(){
    Matrix ret(n, m);
    for(int i=0; i<m; i++){
      for(int j=0; j<n; j++){
        ret.setVal(j, i, val[n*i+j]);
      }
    }
    return ret;
  }
  
  double cofactor(int i, int j){
    Matrix mat = clone();
    mat.splice(i*mat.n, mat.m);
    mat.m -= 1;
    for(int k=mat.m-1; k>=0; k--){
      mat.splice(k*mat.n+j, 1);
    }
    mat.n -= 1;
    return mat.det() * ( ((i+j+2)%2==0) ? 1 : -1);
  }
  double det(){
    if(isSquare()){
      if(m == 2){
        return val[0]*val[3]-val[1]*val[2];
      }else{
        double d = 0;
        for(int k=0; k<n; k++){
          d += val[k] * cofactor(0, k);
        }
        return d;
      }
    }else{
      return 0.0;
    }
  }
  Matrix cofactorMatrix(){
    Matrix mat(m, n);
    for(int i=0; i<m; i++){
      for(int j=0; j<n; j++){
        mat.setVal(j, i, cofactor(i, j));
      }
    }
    return mat;
  }
  Matrix inverse(){
    if(isSquare()){
      double d = det();
      if(d != 0){
        Matrix mat;
        if(m>2){
          mat = cofactorMatrix();
        } else {
          mat = Matrix(2, 2);
          mat.setVal(0, 0, val[3]);
          mat.setVal(0, 1, -val[1]);
          mat.setVal(1, 0, -val[2]);
          mat.setVal(1, 1, val[0]);
        }
        mat.time(1 / d);
        return mat;
      }
    }else{
      return Matrix(-1, -1);
    }
    return Matrix(-1, -1);
  }
};
std::ostream& operator<<(std::ostream& os, const Matrix& mat){
  for(int i=0; i<mat.getM(); i++){
    for(int j=0; j<mat.getN(); j++){
      os << mat.getVal(i,j) << " ";
    }
    os << std::endl;
  }
  return os;
}

//乱数
// 注意: longではなくint(32bit)にすべき
unsigned long xor128(){
  static unsigned long x=123456789,y=362436069,z=521288629,w=88675123;
  unsigned long t;
  t=(x^(x<<11));x=y;y=z;z=w; return( w=(w^(w>>19))^(t^(t>>8)) );
}

//Random Projection
class RandomProjection {
  int d, N, k; //d:次元数  N:ドキュメント数  k:削減後の次元数
  Matrix R;
  Matrix X;
public:
  RandomProjection(int d, int N, int k):d(d),N(N),k(k),R(k,d){
    //R行列の作成
    for(int i=0; i<k; i++){
      for(int j=0; j<d; j++){
        int r = xor128()%60000;
        if(r<10000) R.setVal(i, j, sqrt(3.0));
        else if(r<50000) ;
        else R.setVal(i, j, -sqrt(3.0));
      }
    }
  }
  
  //行列matの次元削減
  void build(const Matrix& mat){
    X = R * mat;
  }

  //クエリ検索(1番cos類似度の高いものを1つ返す)
  int search(const Matrix& query){
    Matrix Q = R * query;
    int ret = -1;
    double val = 0;
    
    for(int i=0; i<N; i++){
      double tmp = 0, sz1 = 0, sz2 = 0;
      for(int j=0; j<k; j++){
        tmp += Q.getVal(j, 0) * X.getVal(j, i);
        sz1 += Q.getVal(j, 0) * Q.getVal(j, 0);
        sz2 += X.getVal(j, i) * X.getVal(j, i);
      }
      tmp /= sqrt(sz1) * sqrt(sz2);
      if(val < tmp){
        ret = i;
        val = tmp;
      }
    }

    return ret;
  }

};

int main(){
  //次元削減後の次元
  int k = 500;

  //make wordlist
  std::map<std::string, int> wordmap;
  int d = 0; //単語の種類数
  std::string word;
  std::ifstream wordlist("wordlist");
  while(std::getline(wordlist, word)){
    wordmap[word] = d++;
  }
  std::cout << "dim: " << d << "->" << k << std::endl;


  //make Document Matrix
  std::ifstream docs("docs");
  int N; //ドキュメント数
  docs >> N;
  std::cout << "doc: " << N << std::endl;
  std::cout << "documents load..." << std::flush;;

  Matrix Doc(d, N);
  std::vector<std::string> doc_text;
  int cnt;
  for(int i=0; i<N; i++){
    std::string text;
    docs >> cnt;
    for(int j=0; j<cnt; j++){
      docs >> word;
      text += word;
      int id = wordmap[word];
      Doc.setVal(id, i, 1+Doc.getVal(id, i));
    }
    doc_text.push_back(text);
  }
  std::cout << "finish!" << std::endl;


  //make Random Projection Matrix
  RandomProjection rp(d,N,k);
  std::cout << "dimension reduction..." << std::flush;
  rp.build(Doc);  
  std::cout << "finish!" << std::endl;



  //query search
  // Input:  単語数 単語1 単語2 ...
  while(std::cin >> cnt){
    Matrix query(d, 1);
    for(int i=0; i<cnt; i++){
      std::cin >> word;
      if(wordmap.count(word)>0){
        int id = wordmap[word];
        query.setVal(id, 0, 1+query.getVal(id, 0));
      }
    }

    int ret = rp.search(query);
    if(ret != -1){
      std::cout << doc_text[ret] << std::endl;
    }
  }


  return 0;
}

結果

とりあえず、4608次元を500次元まで落としてみた。
いくつか検索してみる。

$ ./a.out
dim: 4608->500
doc: 1987
documents load...finish!
dimension reduction...finish!
2 Google メガネ
Googleメガネを規制騒動過熱
3 全国 託児 警察
全国の警察に託児所整備へ
1 サル
サル、まばたきで情報伝達か
2 練馬 区
練馬区で4棟焼く火災3人死亡

次元削減後の行列Xを使っても検索できてるっぽい。

2次元まで次元削減してプロット

ちょっと気になったので、k=2にして、各要素をプロットしてみた。

80個に分かれてる!?いくつか見てみる。

$ paste map.dat matome.data| grep " 0 0"| head
17 0 0	6 白 鵬 25 度目 優勝 朝 青龍 と 並ぶ
25 0 0	4 小平 住民 投票 5 万 票 開票 なし
45 0 0	1 いじめ 対策 法案 どう 評価 する
60 0 0	4 就寝 中 に 胸 刺さ れ 、 男性 死亡
71 0 0	1 九州 - 中 四国 で 早め の 梅雨入り
81 0 0	2 1 円 刻み 運賃 導入 に 温度 差
98 0 0	7 シリア 双方 が 和平 会議 出席 へ
101 0 0	4 2 児童 菓子 食べ 一 晩 過ごす
109 0 0	2 形骸 化 も ノー 残業 デー の 実態
150 0 0	1 「 待機 児 ゼロ 」 宣言 本気 度 は

$ paste map.dat matome.data| grep " 0 -1.73205"| head
1 0 -1.73205	3 エガ ちゃん 全裸 で 客席 ダイブ
8 0 -1.73205	6 イチロー 、 サヨナラ 危機 救う
20 0 -1.73205	0 Google メガネ を 規制 騒動 過熱
46 0 -1.73205	2 就活 後ろ 倒し 効果 と 問題 点 は
83 0 -1.73205	1 無意味 な 正座 は 体罰 事例 示す
89 0 -1.73205	5 ダニ 媒介 疾患 国内 死者 9 人 に
113 0 -1.73205	6 広島 ・ 堂林 の 特権 剥奪 を 示唆
120 0 -1.73205	2 女性 版 クールビズ 支持 得る か
128 0 -1.73205	2 「 標準時 2 時間 早める 」 課題 は
138 0 -1.73205	4 20 キロ 圏 警戒 区域 すべて 再編

ちょっと期待して見てみたけど、特に何かピンとくるようなまとまり方ではなかった。。。
詳しくみたら何かあるかも。(本当は、正当にk-meansとかを試すべきだろうけど...)