飯島泰蔵さんの本を読んだ

思ったことなどをメモがてら書く．

今回読んだのはこの本．

視覚情報の基礎理論―パターン認識問題の源流

飯島さんは部分空間法とか複合類似度法とかOCRとかその関連の人．

この本ではまず，画像とは何か？から入ってそこから，視覚情報の基礎方程式（電磁気学でのマクスウェル方程式のようなもの）を導出し，それを礎に様々な方向に発展させている．基礎方程式はアフィン変換に不変ってことを証明したり，最適な視野を求める方法（顕著性マップのようなもの？）だったり，空間エントロピーを求めたり，などなど．それらの中で文字認識のために具体的に実装したものが複合類似度法．

・基礎方程式は何を言っているかさっぱりで，現在これを理解している人はいるのだろうか？
・複合類似度の他にも様々な類似度を提案していて，その中に相互部分空間法のような，集合間類似度なども書いてあった．説明が全く書いてなかったが，結構ためになるのでは？などと思った．主成分分析の解析解についても述べていて，固有値の予測式も導出しており，これを使ったら面白いことできそうな気がした．
・どのように引用していたか分からなかったが，参考文献には大津展之先生の博士論文が書いてあった，やはり読んでみたい．
・パターン認識は認識の一番簡単な問題といっており，そろそろパターン認識から脱出して認識論に立ち戻ってもいい時期じゃないかなーなんて思ったり．

このような哲学から入って方程式を組み立てて実装までやる研究とか憧れる．

ちなみに，時系列の研究もしていたらしく，下のような本も書いているので，また機会があったら読んでみたい．

自然観測法の理論―瞬時性に着目した新しい波形解析法

ディジタル自然観測法―時系列解析のための新しい理論

行列の行/列の定数倍の処理と速度

ノルムの正規化などをする際に，行列Xの行や列単位で定数倍したい時がありますが，MATLABでどう書くのが一番早いかを測定したので，メモ．

実験コード

N = 3000;
D = 1000;
A = rand(D,N);
B = (1:N);

%% 実装1
for i = 1:100
    x = bsxfun(@times, A, B);
end

%% 実装2
for i = 1:100
    x = A.*repmat(B,D,1);
end

%% 実装3
for i = 1:100
    x = A*sparse(1:N,1:N, B);
end

%% 実装4
for i = 1:100
    x = A*diag(B);
end

時間の測定は，MATLABのプロファイラーを使用．

結果

• 実装1 : 1.19s
• 実装2 : 1.48s
• 実装3 : 1.36s
• 実装4 : 90.47s以上（終了しなかったので，途中で停止）

これまでは実装2をよく使っていたのですが，実装1が僅差ながらも早いという結果に．

渡辺慧先生の本を読んだ。

醜いアヒルの子定理や部分空間法で有名な渡辺慧先生。
渡辺慧 - Wikipedia

部分空間法は知っていて、よく使っているが、
この人がどのような人か、何をしたかなど全く知らなかったので、いくつか本を読んでみた。

知るということ 認識学序説 (ちくま学芸文庫)

認識とパタン (1978年) (岩波新書)

これらは、認識論について非常に分かりやすく書かれており、面白く読めた。

M系列符号の生成

トランジスタ技術(2012/8)のKinect特集でM系列符号が使われていて、面白そうなので、MATLABで実装してみた。

コード

%% M系列符号
clear;
N = 3;
M = 5;
R = randi([0 1], M,1);
nP = 2^M-1;
B = zeros(nP*30, 1);

% 生成
for iP = 1:nP*30
    R = [xor(R(N),R(end)); R(1:end-1)];    
    B(iP) = R(1);
end

% 結果
figure(1);
subplot(2,1,1);
plot(-40:40,xcorr(B(1:1:end),40,'coeff'));
subplot(2,1,2);
plot(-40:40,xcorr(B(1:10:end),40,'coeff'));

結果

自己相関のプロット。

理論では周期は2^M-1=31で、実験結果もちゃんとその周期となった。
（キレイに相関が出ていないのは、MATLABのxcorrの処理が原因？）
生成した符号の部分集合でも、周期が同じであった。
面白い。

参考

トランジスタ技術 2012年 08月号 [雑誌]

http://www.cqpub.co.jp/dwm/contents/0007/dwm000701150.pdf

Woodbury行列反転公式を使って、逆行列演算を高速化。

回帰の問題などでは(X'X)^-1をよく使う。
Xが横長であるときX'Xが非常に大きな行列になるので、
逆行列を求めるのが大変で高速化が重要になる。

そこで、Woodbury行列反転公式(Woodbury matrix inversion formula)
(A+BCD)^-1 = A^-1 - A^-1 B(C^-1+DA^-1 B)^-1DA^-1
という超有名公式を使って、どれくらい高速化できるか実験した。

コード

%% setting
nLoop = 10;
h = 100;    %1000;
w = 2000; %1000;
X = rand(h,w);

%% normal
tic;
for i = 1:nLoop
    Y = X'*X + 0.05*eye(w,w);
    Yi = inv(Y);
end
t1 = (toc/nLoop);

%% woodbury matrix inversion formula
%(A+BCD)^-1 = A^-1 - A^-1 B(C^-1+DA^-1 B)^-1DA^-1
tic;
for i = 1:nLoop
    Ai = 1/0.05*eye(w,w);
    XAi = X*Ai;
    Yi = Ai - Ai*X'/(XAi*X')*XAi;
end
t2 = (toc/nLoop);

%%
disp([t1 t2]);

結果

h=100、w= 100のとき、通常 0.0009s 高速化 0.0012s
h=100、w= 500のとき、通常 0.0320s 高速化 0.0205s
h=100、w=1000のとき、通常 0.1977s 高速化 0.0625s
h=100、w=2000のとき、通常 1.3256s 高速化 0.1901s

通常の場合に逆行列を求める際にinv()を使用しているが、実際の回帰の問題などではバックスラッシュ(/)などの演算子を使うので、もう少し早くなると思われる。

結論

行列の縦横比が5以上であれば、高速化の価値があると思われる。

追記(2013-02-25)

コメントで指摘をもらったので，訂正．
実験に使ったコードの21行目は

    Yi = Ai - Ai*X'/(XAi*X')*XAi;

ではなく，正しくは

    Yi = Ai - Ai*X'/(eye(h) + XAi*X')*XAi;

結果自体にはそれほど影響はないかと．

MACでのmexの設定の注意点

MACのMATLABで

mex -O hogehoge.cpp

とやると，

error: stdio.h: No such file or directory

が出力され，対処に思いの外，時間がかかったので，メモ．

対処法

/Applications/MATLAB_R2011a_Student.app/bin/gccopts.shの

CC='gcc-4.2'
SDKROOT='/Developer/SDKs/MacOSX10.6.sdk'

を以下のとおりに書きかえる．

CC='gcc'          
SDKROOT='/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.7.sdk'

使うコンパイラによって変更するファイルが違うみたいなので注意．
GCCを使う場合はgccopts.sh，
system ANSI compilerを使う場合はmexopts.sh．

これは試行錯誤で動くのを確認した程度なので，なぜこれでいいのかは分からない．
どこかにちゃんと仕様が書いてあるのだろう．

参考

Matlabの設定（mex） - 豊橋の研究日記

JETの画像とカラー画像を重ねて描画

MATLABでJETの画像とカラー画像を重ねた画像をよく見るのだが，
どう作ればいいのか分からなかったので，調べたついでにメモ．

% 準備
clear;
close all;

img = im2double(imread('gundam.jpg'));
img = imresize(img, 0.1);
gimg = rgb2gray(img);
h = fspecial('log', 20, 7);
simg = filter2(h, gimg);

% jetの画像を生成 
mask = simg>0.001;
ssimg = simg;
ssimg(~mask) = 0;
jimg=ind2rgb(gray2ind(ssimg/max(ssimg(:)), 256),jet(256));
jimg(repmat(~mask,[1,1,3])) = img(repmat(~mask,[1,1,3]));

% 合成して描画
figure(1);
imshow(img);
hold on;
h=imshow(jimg);
set(h,'AlphaData', 0.5);
hold off;

figure(2);imshow(img);
figure(3);imagesc(simg);
figure(4);imshow(jimg);

結果はこのようになる．

JETをどのように描画したいかは場合によるが，基本は同じ．
アルファブレンディングのところは，

imshow(0.5*img+0.5*jimg);

と書いてもいいかもしれない．