2017-08-22 16 views
4

私は線形代数C++ライブラリEigen3でランダム対称行列を得ようとしています。私は次のようにしています:C++ Eigen3行列の不思議な振る舞い

Eigen::MatrixXd m(3, 3); 
m.setRandom(); 
m = 0.5 * (m + m.transpose()); 

しかし、再考は全く間違っています。しかし、変数mを書き換えずに単純にコンソールに出力すると、次のようになります。

Eigen::MatrixXd m(3, 3); 
m.setRandom(); 
cout << 0.5 * (m + m.transpose()) << endl; 

すべてが正しく動作しているようです。私はどこに問題があるのか​​理解できません。それは、transposeや*や+のようなメソッドは、怠惰なやり方で行列mを参照する代わりに、新しい行列を即座に作成しないからですか?しかし、どうすれば正式な文書からそれを知るべきですか?そして、圧倒的にエラーを起こしやすいこのような行動ではないのですか?

更新: はい、私は怠惰な計算についての私の推測は正しいと思います。だから今

/** \returns an expression of the transpose of *this. 
    * 
    * Example: \include MatrixBase_transpose.cpp 
    * Output: \verbinclude MatrixBase_transpose.out 
    * 
    * \warning If you want to replace a matrix by its own transpose, do \b NOT do this: 
    * \code 
    * m = m.transpose(); // bug!!! caused by aliasing effect 
    * \endcode 
    * Instead, use the transposeInPlace() method: 
    * \code 
    * m.transposeInPlace(); 
    * \endcode 
    * which gives Eigen good opportunities for optimization, or alternatively you can also do: 
    * \code 
    * m = m.transpose().eval(); 
    * \endcode 
    * 
    * \sa transposeInPlace(), adjoint() */ 

を私は計算の長鎖を実行するときに私が使用すべきかのパターン疑問に思って:それはtranspose方法のdocummentationに記載されていますか?どこでも.eval()を書く?正直言って、それはかなり醜いですが、まだ誤りがちです。

+0

@Someprogrammerdude m.transpose()の場所に行列Mを転置が、結果を返すので、それは問題ないはずはありません –

+0

@Someprogrammerdudeだから私はそれが重複しているとは思わない、少なくともリンクされた質問ではない –

+0

'm =(0.5 +(m + m.transpose()))。eval();'式の最後に単一の 'eval()'をddします。あなたの場合など、エイリアシングの問題が発生した場合は、 'eval()'を追加する必要があります。 –

答えて

3

「私は怠惰な計算についての私の推測は正しいと思う」

はい、あなたは正しいです。遅延評価の規則はhereと記載されています。私は以下の点のいくつかを抽出しました:

Eigenは、それを一時変数に評価するかどうかを、各部分式ごとに自動的に決定します。 [...]

エクスプレッションテンプレートベースのライブラリでは、サブ表現を一時的に評価することを避けることができ、多くの場合、速度が大幅に向上します。これは、式が直ちにではなく、できるだけ遅く評価されるため、遅延評価と呼ばれます。しかし、ほとんどの他のエクスプレッションテンプレートベースのライブラリは、常にレイジー評価を選択します。それには2つの問題があります。まず、遅延評価が必ずしもパフォーマンスの良い選択とは限りません。第2に、例えばマトリックス製品の場合のように、遅延評価は非常に危険です。matrix = matrix*matrixを実行すると、マトリックス製品が動作する方法のために、マトリックス製品が遅延評価された場合に間違った結果が生じます。

"だから私は長い計算チェーンを実行するときにどのようなパターンを使用すべきですか?

一般的に、表現テンプレートは、直ちに/遅れて評価するかどうかの問題を解決するはずですが、時にはすべてのエッジケースを見つけられない場合があるので、matrix.transpose()は1つです。

  • フォース遅延評価: -

  • 行列1と行列2の間にはエイリアシング matrix1.noalias() = matrix2 * matrix2;遅延評価はOKですが、あなたは他がそうでなければ選択されていた遅延または即時評価を、強制的に .eval().noalias()を追加することができます

    フォース即時評価:matrix1 = matrix1.transpose().eval()遅延評価がちょうど.eval()を追加し、あなたの転置ケースについては

OKではない:

m = 0.5 * (m + m.transpose()).eval(); 
1

pingulの答えがあなたのソリューションが失敗する理由を説明しています。一時的なものを避けたい場合は、できることを追加したかっただけです。基本的に、トリックは自分の行列の1三角形の半分に式を割り当てることです:

#include <iostream> 
#include <Eigen/Core> 

int main() { 
    Eigen::Matrix4f M; 
    M.setRandom(); 
    std::cout << "M:\n" << M << '\n'; 
    // for comparison, evaluate into other matrix (similar to using .eval()): 
    Eigen::Matrix4f M2 = 0.5f*(M+M.transpose()); 
    // in place operation, only assing to lower half, then to upper half: 
    M.triangularView<Eigen::StrictlyLower>() = 0.5f*(M + M.transpose()); 
    M.triangularView<Eigen::StrictlyUpper>() = M.transpose(); 
    std::cout << "M2:\n" << M2 << "\nM:\n" << M << '\n'; 
} 

実際には、あなただけの「ランダムな対称行列」と与えられた行列Mの必ずしも対称一部をしたい場合は、あなたは、単に下の部分(またはその逆)に上部をコピーすることができます。

M.triangularView<Eigen::StrictlyLower>() = M.transpose(); 
関連する問題