2017-03-10 8 views
2

私はいくつかの行列ライブラリにアクセスできますが、このプロジェクトではコンパイル時の定義とSVDの包含のために私はEigenを使用しています。効率的なマトリックス転置行列の固有値化

今、私は、次の操作をしています:

私が理解したよう
Eigen::Matrix<double,M,N> A;  // populated in the code 

Eigen::Matrix<double,N,N> B = A.transpose() * A; 

が、これはAのコピーを作成し、再び乗算され転置を形成します。この操作は、比較的小さな行列(M = 20〜30、N = 3)で実行されますが、1秒間に何百万回も実行されるため、できるだけ高速でなければなりません。

私は次を使用して高速であることを読んで:

B.noalias() = A.transpose() * A; 

私は、入力として受け入れ、Bを満たす自分のサブルーチンを書くことができますが、使用して効率的に、既存の実装がある場合、私は思っていましたサイクルの最小量。

+0

これを見ることを検討してください:http://scicomp.stackexchange.com/questions/25283/beating-typical-blas-libraries-matrix-multiplication-performance –

+0

これは役に立ちますか? http://stackoverflow.com/questions/39606224/does-eigen-have-self-transpose-multiply-optimization-like-h-transposeh – kennytm

答えて

1

まず、Eigenはテンプレート式に依存しているため、A.transpose()は一時的に評価されません。

第二に、中:

Matrix<double,N,N> B = A.transpose() * A; 

固有したがって、(ここでは、コンパイラはBのコンストラクタを呼び出しているため)、およびBは、式の右側に表示することはできません何も一時的には全く作成されないことを知っています。これは同等です:

Matrix<double,N,N> B;    // declare first 
B.noalias() = A.transpose() * A; // eval later 

最後に、このような小さな行列のために、私はB.selfadjointView().rankUpdate(A)の使用は(kennytmコメントで提案されているように)役立つことを期待しないでください。 otherhandで

、N = 3で、それは怠惰な実装を試みる価値があるかもしれません:

B = A.transpose().lazyProduct(A)

は、念のために。 Eigen'sには最高の製品実装を選択するヒューリスティックが組み込まれていますが、ヒューリスティックは簡単で評価が速いため、100%正しいとは限りません。

+0

ありがとうございました。遅延プロジェクトのヒントは非常にです。さて、エイゲンはGPU上のcudaで動作しないことがわかったので、私は別のものを完成させた。私は図書館が好きです。さらに、Aをまったく構築しないのが最も効率的です。これは私が行ったことです。 –

関連する問題