徵和C++ 11的類型推斷失敗的矩陣產品

的喬萊斯基我試圖把其轉置矩陣的產品的Cholesky分解，利用特徵向量和C++ 11「自動」類型。我試圖做的時候出現問題徵和C++ 11的類型推斷失敗的矩陣產品

auto c = a * b 
auto cTc = c.tranpose() * c; 
auto chol = cTc.llt();

我正在使用XCode 6.1，Eigen 3.2.2。我得到的類型錯誤是here。

這個小例子，顯示了我的機器上的問題。將c的類型從auto更改爲MatrixXd以查看其工作。

#include <iostream> 
#include <Eigen/Eigen> 
using namespace std; 
using namespace Eigen; 


int main(int argc, const char * argv[]) { 
    MatrixXd a = MatrixXd::Random(100, 3); 
    MatrixXd b = MatrixXd::Random(3, 100); 
    auto c = a * b; 
    auto cTc = c.transpose() * c; 
    auto chol = cTc.llt(); 
    return 0; 
}

有沒有辦法讓這項工作，同時仍然使用汽車？

作爲一個方面的問題，是否有性能的理由不能斷言基質是在每個階段MatrixXd？使用自動將允許Eigen保留答案，無論它幻想哪種奇怪的模板表達。我不確定是否將它鍵入MatrixXd會導致問題。

來源

2014-11-24 c0g

讓我總結什麼是怎麼回事，爲什麼這是錯的。首先，讓我們來實例化auto關鍵字他們正在採取的類型：

typedef GeneralProduct<MatrixXd,MatrixXd> Prod; 
Prod c = a * b; 
GeneralProduct<Transpose<Prod>,Prod> cTc = c.transpose() * c;

注意徵是一個表達式模板庫。這裏，GeneralProduct<>是代表產品的抽象類型。不執行計算。因此，如果複製cTc到MatrixXd爲：

MatrixXd d = cTc;

這相當於：

MatrixXd d = c.transpose() * c;

那麼該產品a*b將進行兩次！因此，在任何情況下，更優選地將內明確的臨時評估a*b，和同爲c^T*c：

MatrixXd c = a * b; 
MatrixXd cTc = c.transpose() * c;

最後一行：

auto chol = cTc.llt();

也相當錯誤的。如果cTc是一個抽象產品類型，那麼它會嘗試實例化一個Cholesky因式分解，這是一個不可能的抽象產品類型。現在，如果cTc是一個MatrixXd，那麼你的代碼應該工作，但是這仍然沒有方法llt()的首選方法是相當實行一班輪表達，如：

VectorXd b = ...; 
VectorXd x = cTc.llt().solve(b);

如果你想要一個名爲喬萊斯基對象，然後而使用它的構造：

LLT<MatrixXd> chol(cTc);

甚至：

LLT CHOL（c.transpose（）* C）;

這是等效的，除非你必須在其他計算中使用c.transpose() * c。

最後，根據的a和b大小的，它可能是最好計算cTc爲：

MatrixXd cTc = b.transpose() * (a.transpose() * a) * b;

未來（即本徵3.3），本徵將能夠看到：

auto c = a * b; 
MatrixXd cTc = c.transpose() * c;

作爲四個矩陣m0.transpose() * m1.transpose() * m2 * m3的乘積，並將括號放在正確的位置。但是，它不知道m0==m3和m1==m2，因此如果首選的方法是臨時評估a*b，那麼您仍然必須自己做。

來源

2014-11-25 08:57:23 ggael

謝謝 - 從圖書館的開發人員那裏聽到真的很棒！我的理由是看Eigen是否可以在他們有用的形狀時正確地優化'm0.transpose（）* m1.transpose（）* m2 * m3' - 因此我想把所有東西放在表達式空間中，直到最後一刻。是否由於模板的工作原因，我無法對GeneralProduct進行cholesky分解，難道只有沒有人有足夠的關注將其添加到Eigen中嗎？還是有理由這麼做是愚蠢的？ – c0g 2014-11-25 13:00:00

問題是第一次乘法返回Eigen::GeneralProduct而不是MatrixXd而auto正在拾取返回類型。你可以從Eigen::GeneralProduct隱含地創建MatrixXd，所以當你顯式聲明它的類型時它能正常工作。見http://eigen.tuxfamily.org/dox/classEigen_1_1GeneralProduct.html

編輯：我不是做鑄件的本徵產品或性能方面的專家。我剛剛從錯誤信息中猜出了答案，並從在線文檔中得到了確認。性能分析始終是您檢查代碼不同部分性能的最佳選擇。

來源

2014-11-24 20:45:55

對於投射到「MatrixXd」有沒有性能上的提升？顯然，在這個受限制的例子中，它會很小，但在現實生活中呢？你會說什麼解決方案在這裏 - 使用ProductReturnType來做這件事似乎在執行'c.tranpose（）* c'步驟時會有些生氣。 – c0g 2014-11-24 21:02:42

我不是在Eigen專家，但像這樣的圖書館經常回運營代理對象，然後用隱式轉換或構造迫使實際工作。（表達式模板是這方面的一個極端例子。）這樣可以避免在許多情況下不必要地複製整個數據矩陣。不幸的是，auto很高興創建一個代理類型的對象，通常庫的用戶永遠不會顯式聲明。由於您需要最終計算出數字，因此從鑄造到MatrixXd本身不會造成任何性能下降。（斯科特邁爾斯，在「有效的現代C++」，給出了使用auto與vector<bool>，其中operator[](size_t i)返回代理的相關示例。）

來源

2014-11-24 21:21:59 sfjac

不要使用auto與本徵表達式。我碰上這個更「戲劇化」的問題之前，請參閱

eigen auto type deduction in general product

，並通過本徵的創作者（蓋爾）不使用auto與本徵的表現之一建議。

從表達式轉換爲MatrixXd這樣的特定類型應該非常快，除非您想懶惰評估（因爲在執行投射時評估結果）。

來源

2014-11-24 23:09:46 vsoftco

我認爲這與你的稍有不同。你從'id'返回一個對象，然後由'autoresult'引用 - 然後'Id'的對象消失，'autoresult'引用一些不會退出的東西。 'c.transpose（）'引用'c'，它仍然存在，所以我的'auto cTc'不應該引用任何不存在的東西 - 類似於ggael作爲一種解決方案提供的使用'auto id = Id（Foo， 4）'。我想測試Eigen是否可以正確地注意到'c.transpose（）* c = b.transpose（）* a.transpose（）* a * b'簡化了數學運算，所以理想情況下它將保留爲表達式。 – c0g 2014-11-24 23:35:29

是的，你是對的，但是，如果我沒有得到解釋，我完全不知道爲什麼代碼表現異常。由於Eigen有很多類型（基本上每個表達式都是不同的類型），加上間接引用，'auto'會讓事情變得複雜，因爲你並不真正意識到底層發生了什麼，以及如何評估這些表達式。 – vsoftco 2014-11-24 23:43:39

徵和C++ 11的類型推斷失敗的矩陣產品

回答

相關問題