xiangze's sparse blog

機械学習、ベイズ統計、コンピュータビジョンと関連する数学について

redsvdのwrapper RRedSVDとRパッケージ作成、公開についての覚え書き

疎行列に対する特異値分解の高速な近似実装 redsvdのR wrapperである RRedSVDを作成、公開しました。

https://github.com/xiangze/RRedsvd

ライセンスはBSDです。

 

TokyoRでは特異値分解とは何か、どんな手法があるかなどの基本的説明、RとC++を橋渡しするRcpp, RcppEigen packageの紹介などを解説しました(全ては説明しきれませんでした)。

 

アルゴリズムと性能に関してはredsvdのプロジェクトページ(https://code.google.com/p/redsvd/ 日本語、英語), 製作者のhillbigさんのブログ(http://hillbig.cocolog-nifty.com/do/2010/08/redsvd-aa59.html 日本語)、NIPS tutorial(http://amath.colorado.edu/faculty/martinss/Talks/2009_NIPS_tutorial.pdf 英語)が詳しいです。

 

以下はRcppを用いたパッケージ開発に関する覚書です。

続きを読む

Python, Sympyを用いたVanishing Component Analysisの実装(と動作)

未だ挙動が怪しいですが、Pythonの数式処理ライブラリSympyを用いて Vanishing Component Analysis(論文pdf)を実装しました。
Sympyは可読性が高いのが利点ですが、速度を考慮した実装ではないため特に次数の高い多項式まで求めようとした場合にどんどん動作が遅くなってしまいます。
また動作に未だ不審な点があります。

既知の問題点

  • 速度が遅い
    • pure pythonで実装されており、抽象度の高いことも実現できるSympyの問題と思われますが、
  • Errorを吐いて落ちる場合がある。
  • x^2+y^2=0のような単純な多項式+noiseのデータを入力してもかなり複雑なVanishing componentとしてかなり複雑な多項式が返ってきてしまう。
  • Sympyではlimit関数で変数に代入ができますが、多変数多項式の変数に一括して代入する方法がすぐに見つけられなかったので独自実装してしまいました(一方でsumが使えます)。

今後

デバッグとともに得られたVanishing componentのプロットをsympy単体またはmatplotlibで行いたいです。

link

どんなデータでも(※)線形分離可能にしてしまう技術,Vanishing Component Analysis(ICML 2013)を紹介してきました
Rのほうが1行が長いせいか行数が少なく、またライブラリに依存せずに書けるようです。
R で Vanishing Component Analysis - Mi manca qualche giovedi`?
Vanishing Component Analysis を試作してみました 塚原裕史のブログ matlabです。componentのプロットがされています。

Pythonの数式処理ライブラリについてまとめました。C++で実装されたより速そうながあるようです。
Groebner bases in Python
Groebner基底の紹介です。
駆け足で読む『What is … a Grobner Basis?』

おまけ?

論文pdf擬似コードの間違い

  1. FindRangeNull(Figure 2.)の添字 t->i
  2. D_ij -> D_ii (そもそもDの非対角要素は0)

CVPRの論文タイトルの頻出単語と可視化

コンピュータビジョン分野の査読つき国際会議であるCVPR(Conference on Computer Vision and Pattern Recognition)の論文題名を
http://www.cvpapers.com/
から取得し、その変遷をRのtmモジュールを使い分析、可視化してみました。
CVPR2007からCVPR2013までの論文タイトルが対象となります。


2014/12/9追記
http://colinlea.com/guides/cvpr2013.html
では本文の単語を用いたトピックモデルで論文が分類された結果を見ることが出来ます。
NIPS,ICMLのものもあります。

続きを読む

QsysとverilogのTIPS

Altera社のツールQSysについての細かい事項メモ

  • generate文の2重ループが使える。
  • ただしgenerate文にはラベルが必要 例:
   genvar   i,j;

   generate
      for(i=0;i<16;i=i+1)begin :rep_tmpx
     assign tmpx[i]=(idx==i)?(i<<1):0;
      end
      for(j=0;j<16;j=j+1)begin :rep_tmp
         if(j==0)assign tmp[j]=tmpx[j];
         else assign tmp[j]=tmp[j-1]|tmpx[j];
      end
   endgenerate
  • generate文のある独自回路を含んだQsysのSimulationモデルは作成できないらしい。
  • 独自回路で定義したparameterはQSys上のGUIで設定できる。
  • 独自回路をQsysに追加した場合、独自回路内部の修正であればQSysのregenerate(再合成)は必要ない。
  • 複数clockを使う場合はそれぞれにclockモジュールを用意する必要がある?(要追加調査)

おそらくはQSys,Quartusに限らないこと(要調査)

  • register,wireのビット範囲指定に定数以外を使うことはできない。register,wireの配列の添字ならOKなのでassignを使って切り替える。

QsysでのAvalon busへの独自回路追加

AlteraのSoC開発ツールQsysではAvalon busに独自に作成した回路モジュールを接続することができます。
その際にAvalon busの仕様にあるinput,output信号をモジュールでは定義する必要があります。Qsysでの設定手順と一緒にまとめます。
この記事は基本的にはmarseeさんのAvalon-MMスレーブペリフェラルを参考にさせていただきました。
Qsys追加した独自回路モジュールをAvalon busを介してNIOSで操作する構成は同じですが、今回はoutput信号が直接7segの各セグメントに接続される普通(?)の構成なので、追加モジュールはより単純になります。
追加モジュール My7seg.v
avs_s1_*がavalon busの信号です。

module My7seg(/*AUTOARG*/
   // Outputs
   avs_s1_readdata, hex0, hex1, hex2, hex3, hex4, hex5, hex6, hex7,
   // Inputs
   csi_global_reset, csi_global_clk, avs_s1_address, avs_s1_read,
   avs_s1_write, avs_s1_writedata
   );
   input         csi_global_reset;   
   input         csi_global_clk;     

   input [2:0] 	 avs_s1_address;     
   input         avs_s1_read;        
   output [31:0] avs_s1_readdata;    
   input         avs_s1_write;       
   input [31:0]  avs_s1_writedata;   

   output [6:0]  hex0;
   output [6:0]  hex1;
   output [6:0]  hex2;
   output [6:0]  hex3;

   output [6:0]  hex4;
   output [6:0]  hex5;   

   output [6:0]  hex6;
   output [6:0]  hex7;   

 reg [31:0]		avs_s1_readdata;
   reg [3:0] 		iDIG[0:7];

   // avalon bus write
   always @(posedge csi_global_clk, posedge csi_global_reset) begin : SEVEN_VALUE_PROCESS
      integer i;
      if (csi_global_reset) begin
         for (i=0; i<8 ; i=i+1) begin
            iDIG[i] <= 8'd0;
         end
      end else begin
         if (avs_s1_write) begin
            case(avs_s1_address)
              3'b000 :
                iDIG[0] <= avs_s1_writedata[3:0];
              3'b001 :
                iDIG[1] <= avs_s1_writedata[3:0];
              3'b010 :
                iDIG[2] <= avs_s1_writedata[3:0];
              3'b011 :
                iDIG[3] <= avs_s1_writedata[3:0];
              3'b100 :
                iDIG[4] <= avs_s1_writedata[3:0];
              3'b101 :
                iDIG[5] <= avs_s1_writedata[3:0];
              3'b110 :
                iDIG[6] <= avs_s1_writedata[3:0];
              3'b111 :
                iDIG[7] <= avs_s1_writedata[3:0];
            endcase
         end
      end
   end
    
    // avalon bus read
    always @* begin
        avs_s1_readdata[31:8] = 24'd0;
        case (avs_s1_address)
            3'b000 :
                avs_s1_readdata[7:0] <= {4'b0,iDIG[0]};
            3'b001 :
                avs_s1_readdata[7:0] <= {4'b0,iDIG[1]};
            3'b010 :
                avs_s1_readdata[7:0] <= {4'b0,iDIG[2]};
            3'b011 :
                avs_s1_readdata[7:0] <= {4'b0,iDIG[3]};
            3'b100 :
                avs_s1_readdata[7:0] <= {4'b0,iDIG[4]};
            3'b101 :
                avs_s1_readdata[7:0] <= {4'b0,iDIG[5]};
            3'b110 :
                avs_s1_readdata[7:0] <= {4'b0,iDIG[6]};
            default : // 3'b111
                avs_s1_readdata[7:0] <= {4'b0,iDIG[7]};
        endcase
    end

   
/*   SEG7_LUT AUTO_TEMPLATE(
 .oSEG			(hex@[]),
 .iDIG			(iDIG[@]));
  );
 */
  
   SEG7_LUT seq0(/*AUTOINST*/
		 // Outputs
		 .oSEG			(hex0[6:0]),		 // Templated
		 // Inputs
		 .iDIG			(iDIG[0]));		 // Templated
   SEG7_LUT seq1(/*AUTOINST*/
		 // Outputs
		 .oSEG			(hex1[6:0]),		 // Templated
		 // Inputs
		 .iDIG			(iDIG[1]));		 // Templated
   SEG7_LUT seq2(/*AUTOINST*/
		 // Outputs
		 .oSEG			(hex2[6:0]),		 // Templated
		 // Inputs
		 .iDIG			(iDIG[2]));		 // Templated
   SEG7_LUT seq3(/*AUTOINST*/
		 // Outputs
		 .oSEG			(hex3[6:0]),		 // Templated
		 // Inputs
		 .iDIG			(iDIG[3]));		 // Templated
   SEG7_LUT seq4(/*AUTOINST*/
		 // Outputs
		 .oSEG			(hex4[6:0]),		 // Templated
		 // Inputs
		 .iDIG			(iDIG[4]));		 // Templated
   SEG7_LUT seq5(/*AUTOINST*/
		 // Outputs
		 .oSEG			(hex5[6:0]),		 // Templated
		 // Inputs
		 .iDIG			(iDIG[5]));		 // Templated
   SEG7_LUT seq6(/*AUTOINST*/
		 // Outputs
		 .oSEG			(hex6[6:0]),		 // Templated
		 // Inputs
		 .iDIG			(iDIG[6]));		 // Templated
   SEG7_LUT seq7(/*AUTOINST*/
		 // Outputs
		 .oSEG			(hex7[6:0]),		 // Templated
		 // Inputs
		 .iDIG			(iDIG[7]));		 // Templated
   
endmodule // My7seg

アドレスの各値には表示した数値を書き込む形になっていて、それを
数値を7segへの信号へ変換するモジュールSEG7_LUTとしては以下のようなものを使いました。
http://users.ece.gatech.edu/~hamblen/DE2/DE2_demonstrations/DE2_Default/SEG7_LUT.v
SEG7_LUTモジュールが多数になる場合にはgenerate文を使ったほうがよさそうです(参考:Verilogのgenerate文でのマルチプレクサの記述)。

Qsysを開き、左側のTreeのProject->New Componentを選択し、出てきたウィンドウで最初に表示されるComponent typeタブでName,Display nameを設定します。ここではMy7segにしました。
FilesタブのSynthesis filesの部分の+ボタンをクリックしてこの記述のファイルMy7seg.vを選択し、Analyze synthesis fileを押します。Verilog HDLの記述に問題がなければ成功のメッセージが表示されます。

しかしComponent Editorのウィンドウには大概の場合エラーが出てきてしまいます。これはQsysがAvalonバスの各信号を名前から推測するのに失敗している場合がほとんどなので、手でそれを直してあげます。具体的にはComponent EditorのSignalsタブで各信号の設定を行います。ここではavalon busの信号としてInterfaceでavalon_slave0を、Signal typeではaddress,readdata,readdatavalid,writedata,writedatavalidなどを選択します。
7segへ出力されることとなるhex*の信号ではnew Conduitを選択し、Signal typeはexportにします(1つしかありませんが)。reset信号はglobal_reset,clockはglobalを選択します。

その後Interfaceタブでも設定をしてあげなけなければいけません。まずRemove interface With No Signal ボタンを押して、余計なInterfaceを消します。
その結果としてSignalsタブで選択したInterface(バスのような信号の集合)が残ります。ここで各Interfaceのclock,resetの設定をします。
"global"の囲みの部分にclockがあることを確認し、"global_reset"のassociated clockをglobalに、avalon_slave0のassociated resetをglobal_reset,associated clockをglobalに選択します。

Component Editorの下部に表示されるErrorがなくなったらFinishとします。今回はhex*はすべてoutput信号なのでErrorが出ませんでしたが、avalon以外のinput信号がある場合には同様にclock,resetを設定する必要があります。
Finish後に左側のTreeのNew Componentの下にMy7segが表示されるのでこれをクリックで選択し,Addボタンを押してsystem componentsに追加します(下図)。

その後他のQsys内のモジュールと同様にbus信号、clock,resetを白丸をクリックして接続し、Base adressを他のモジュールとかぶらないように設定します(ここでは0x5000)。7seg信号への出力となるピンも忘れないようにクリックし名前を設定します。これで準備ができたのでQsysの合成を行います。
QuartusのほうではMy7seg.v, SEG7_LUT.vをプロジェクトに追加します。7segへの出力信号hex*を追加したので、Qsysでインスタンス記述をコピーし、TOPファイルにペーストして、該当する信号に忘れずに接続し、合成します。

NIOSのソフトウェアは以下のようになります。

#include <stdint.h>
#include "altera_avalon_pio_regs.h"
#include "system.h"

static void SevenSegCount( void ){
	int count;
	int i;
	int c;

	for (count = 0; count < 16; count++){
	for(i=0;i<8;i++){
	  IOWR_ALTERA_AVALON_PIO_DATA(MY7SEG_0_BASE+i*4, (count+i)&0xf);
	}
	usleep(500000);
  }
}

int main(){ 
  /* Event loop never exits. */
  while (1){
	  SevenSegCount();
  }
  return 0;
}

IOWR_ALTERA_AVALON_PIO_DATAというマクロはaltera_avalon_pio_regsに定義されており、直接該当アドレスに書き込むことはできません(タイマー割り込みの場合の参考:NIOS2奮闘記【22】タイマー割り込み②)
MY7SEG_0_BASEはQsysで生成されたsystem.h内に記述があります。Qsysで設定したアドレスを直接叩くよりはこちらの記述のほうが汎用性が若干高いです。
わざわざヘッダファイルからマクロを探してくるのが面倒な場合はOSを入れる必要があるかもしれません。

今回はSlaveとなるモジュールの場合でしたが、Master,割り込み信号を持ったモジュールも同様に設定できるはずです。。。

へびの補題の写経

ブログ書き初めです。巳年なのでへびの補題の特殊な場合であるコホモロジー群の長い完全列の導出の証明の写経です。

意義、応用範囲

完全列とは \psi \cdot \phi =0
となる準同形写像ψ、Φを用いて作られた列

のことでコホモロジー群とは
Mの開被覆Uと可換環R上に値をとる関数(より正確には局所的な関数の張り合わせである"切断")とを対応づける写像であるコチェインの集合
 C^q(U,R)

 Z^q(U,R)=\{f \in C^q(U,R) | \delta f =0 \}
に対して
 H^q(U,R)=Z^q(U,R) / \delta C^{q-1}(U,R)
と定義される群です。ここでコバウンダリー作用素δは f \in C^{q}(U,R)に対して
 (\delta f)(U_0,U_1, \cdots, U_{q+1}) =\sum_{i=0}^{q+1} (-1)^i f(U_0,U_1,  \cdots, \hat{U}_i \cdots,U_{q+1} )
と定義されます。( \hat{U_i}はUiをのぞいたものという意味)
 \delta  \cdot \delta=0であることから
 H^q(U,R)=I
であるのが自明なようですが、Mの形状(トポロジー)によっては自明な形にならない場合があります。
このコホモロジー群の考えと以下で示す証明は1次元複素多様体の形状とその上での正則関数の極、零点の数とを関係づけるリーマン・ロッホの定理の証明などで用いられます。リーマン・ロッホの定理を多次元の場合に一般化したものが指数定理であり、物理学においてはある種の対称性の破れを説明する場合に用いられる(らしい)です。
以下が長い証明です。

主張

M: パラコンパクトなハウスドルフ空間
U: M上の開被覆
R,S,T M上の可換環
短完全列

Φ :R->S 単射
ψ :S->T 全射
があるとき、長いコホモロジ−群の完全列

が存在する。

証明の流れ

  1. 用語の定義と前提
  2.  Im(\phi^*) \subseteq Ker(\psi^*)
  3.  Ker(\psi^*) \subseteq Im(\phi^*)
  4. δ*の定義と代表元の取り方によらないこと
  5. 帰納的極限とパラコンパクト性を用いた拡張
用語の定義と前提

上記のようにコチェイン C^q(U,R)
 Z^q(U,R)=\{ f \in C^q(U,R) | \delta f =0 \}
に対してq次のコホモロジー群は
 H^q(U,R)=Z^q(U,R) / \delta  C^{q-1}(U,R)
と定義される。
また以下での証明のために
 \bar{ C}^q(U,T) =\psi(C^q(U,S) \subset C^q(U,T) }
を定義する。これによってコチェインの短完全列
 0 -> C^q(U,R) -> C^q(U,S) -> \bar{ C}^q(U,T)  -> 0
ができる。
コサイクルの境界写像 \delta: C^q(U,*) -> C^{q+1}(U,*) とψ、Φは可換であり、

のような可換図式が得られる。
またコホモロジー群間の写像Φ*,ψ*が引き戻しで定義できて、完全列

が得られる。この完全性の証明( Im(\phi^*) \subseteq Ker(\psi^*) , Ker(\psi^*) \subseteq Im(\phi^*) )をまず以下で行う。

 Im(\phi*) \subseteq Ker(\psi*) となること

 \forall g \in Im(\phi^*), \exists h \in H^q
 [g]=[ \phi (h) ]
 \phi \dot \psi=0 , \phi ( \psi (h) )= \phi (g) =0
なので
 g \in Ker(\phi^*)
といえる。

 Ker(\psi*) \subseteq Im(\phi*) となること

 \forall g \in Ker(\psi^*)
 \exists f \in C^{q-1} (U,T) s.t. \psi  g=\delta f
となる。
ψが全射であるので
 \exists x \in C^{q-1}(U,S) s.t. \psi x=f
 \psi g= \delta f =\delta \psi x= \psi \delta x
 \psi(g - \delta x)=0

Im(Φ)=Ker(ψ) (仮定)から
 \exists h \in C^q(U,R) s.t. \phi h = g - \delta x
 \delta \phi h = \delta (g - \delta x)=\delta g - \delta \delta x=\delta g
 g \in Z^q(U,S)から \delta g=0よって
\delta \phi h =0
 \phi単射なのでKer(Φ)=0
 \delta \phi h = \phi \delta h =0
なので
 \delta h =0 \Rightarrow h \in Z^q(U,R)
 H^q(U,S)=Z^q(U,S)/ \delta C^{q-1}(U,S)
なので H^q(U,R),H^q(U,S)におけるh, gの代表元をそれぞれ[h],[g]とすると
 \phi^*[h]=[g]
すなわち
 Ker(\phi^*) \subset Im(\psi^*)
 Ker(\phi^*) = Im(\psi^*)

は完全列

δ*の定義と代表元の取り方によらないこと

 \delta^* : \bar{H}^q(U,T) \rightarrow H^q(U,R) の構成

 \bar{H}^q(U,T)=\bar{Z}^q(U,T)/\delta \bar{C}^{q-1}(U,T)なので
 \forall f \in \bar{C}^q(U,T) \delta f=0 \exists g \in C^q(U,S) s.t. \psi g =f
というgが存在し
 \psi \delta g = \delta \psi g = \delta f=0
とコチェインの列の完全性 (ψが全射、Φが単射)であることから
 \phi h=\delta g
となる h \in C^{q+1}(U,R)が存在する。

同値類
 [f] \in \bar{H}^q(U,T)
 [h] \in \bar{H}^{q+1}(U,R)
に対して \delta^*
 \delta^* [ f ] = [ g ]
と定義する。
ここで
 \phi \delta h = \delta \phi h = \delta \delta g=0
Φは単射なのでKer(Φ)=0, よって
 \phi (\delta h) \Rightarrow \delta h=0

  • δ*が同値類の代表元の取り方によらないこと

  f,f' \in \bar{Z}^q(U,T)に対して
 \psi g=f, \psi g'=f', \phi h=\delta g, \phi h'=\delta g'
 f \sim f'
すなわち
 \exists c \in \bar{C}^{q-1}(U,T) s.t f-f'=\delta c
また \bar{C}^{q-1}(U,T)の定義から
 \exists a s.t. \phi a=c
このとき
 \psi(-\delta a+(g-g'))= \delta \phi a+(\phi g-\psi g')=\delta c +f-f'=0
よって

の完全性から
 \psi(-\delta a+(g-g'))=0 \Rightarrow \exists b s.t. \phi b = -\delta a+(g-g')
 \phi \delta b=\delta \phi  b = - \delta \delta a +(\delta g - \delta g') = \phi h - \phi h' = \phi(h-h')
Φは単射なので
 \delta b=h-h'
よって
 h \sim h'
代表元fの取り方によらないことが示された。

帰納的極限とパラコンパクト性を用いた拡張

帰納的極限の概念を使うと開被覆UをMに広げることができる。
さらにMがパラコンパクトであることから
 \bar{H}^q(M,T)=H^q(M,T)
とできる(らしい)ので長い完全列

を得ることができる。

links

曲面上の関数論 の3章page 53 ~ 60を写経しました。

曲面上の関数論―リーマン‐ロッホの定理へのいざない

曲面上の関数論―リーマン‐ロッホの定理へのいざない

蛇の補題について
latexのパッケージamscdを使って可換図式を書きました。
http://www.jmilne.org/not/Mamscd.pdf