Suy luận kiểu Eigen và C++ 11 không cho Cholesky của sản phẩm ma trận

Question

Tôi đang cố gắng phân tích cholesky của sản phẩm của ma trận với chuyển vị của nó, sử dụng loại Eigen và C++ 11 "tự động". Sự cố xảy ra khi tôi cố gắng làm

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

Tôi đang sử dụng XCode 6.1, Eigen 3.2.2. Lỗi loại tôi nhận được là here.

Ví dụ tối thiểu này cho biết sự cố trên máy của tôi. Thay đổi loại c từ auto thành MatrixXd để xem nó hoạt động.

#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; 
} 

Có cách nào để thực hiện công việc này trong khi vẫn sử dụng tự động không?

Là một câu hỏi phụ, có lý do hiệu suất nào để không khẳng định ma trận là MatrixXd ở mỗi giai đoạn không? Sử dụng tự động sẽ cho phép Eigen giữ câu trả lời như bất kỳ biểu thức mẫu lạ nào mà nó thích. Tôi không chắc chắn nếu gõ nó như MatrixXd sẽ gây ra vấn đề hay không.

Answer 1

Hãy để tôi tóm tắt những gì đang diễn ra và tại sao nó sai. Trước hết, chúng ta hãy nhanh chóng auto từ khóa với các loại họ đang dùng:

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

Lưu ý rằng Eigen là một thư viện biểu mẫu.Ở đây, GeneralProduct<> là một loại trừu tượng đại diện cho sản phẩm. Không có tính toán được thực hiện. Do đó, nếu bạn sao chép cTc đến một MatrixXd như:

MatrixXd d = cTc; 

tương đương với:

MatrixXd d = c.transpose() * c; 

sau đó sản phẩm a*b sẽ được thực hiện hai lần! Vì vậy, trong mọi trường hợp nó là nhiều lợi thế để đánh giá a*b trong một tạm thời rõ ràng, và tương tự cho c^T*c:

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

Dòng cuối cùng:

auto chol = cTc.llt(); 

cũng là khá sai. Nếu cTc là một loại sản phẩm trừu tượng, thì nó sẽ cố gắng tạo ra một hệ số Cholesky làm việc trên một loại sản phẩm trừu tượng không thể thực hiện được. Bây giờ, nếu cTc là một MatrixXd, sau đó mã của bạn nên làm việc nhưng điều này vẫn không phải là cách ưa thích như là phương pháp llt() khá để thực hiện biểu hiện một liner như:

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

Nếu bạn muốn có một đối tượng Cholesky đặt tên, sau đó chứ không phải sử dụng constructor của nó:

LLT<MatrixXd> chol(cTc); 

hoặc thậm chí:

LLT chol (c.transpose() * c);

tương đương trừ khi bạn phải sử dụng c.transpose() * c trong các tính toán khác.

Cuối cùng, phụ thuộc vào kích thước của a và b, nó có thể được thích hợp hơn để tính toán cTc như:

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

Trong tương lai (ví dụ, Eigen 3.3), Eigen sẽ có thể thấy:

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

làm sản phẩm của bốn ma trận m0.transpose() * m1.transpose() * m2 * m3 và đặt dấu ngoặc đơn ở đúng vị trí. Tuy nhiên, nó không thể biết rằng m0==m3 và m1==m2, và do đó nếu cách ưa thích là để đánh giá a*b tạm thời, thì bạn sẽ vẫn phải tự làm điều đó.

Answer 2

Vấn đề là phép nhân thứ nhất trả về Eigen::GeneralProduct thay vì MatrixXd và auto là chọn loại trả về. Bạn có thể tạo một cách rõ ràng MatrixXd từ một số Eigen::GeneralProduct để khi bạn khai báo rõ ràng loại nó hoạt động chính xác. Xem http://eigen.tuxfamily.org/dox/classEigen_1_1GeneralProduct.html

EDIT: Tôi không phải là chuyên gia về sản phẩm Eigen hoặc đặc tính hiệu suất của việc truyền. Tôi chỉ phỏng đoán câu trả lời từ thông báo lỗi và được xác nhận từ tài liệu trực tuyến. Hồ sơ luôn là lựa chọn tốt nhất của bạn để kiểm tra hiệu suất của các phần khác nhau trong mã của bạn.

Answer 3

Tôi không phải là chuyên gia tại Eigen, nhưng các thư viện như thế này thường trả về các đối tượng proxy từ các hoạt động và sau đó sử dụng chuyển đổi hoặc nhà thầu ngầm để ép buộc công việc thực tế. (Biểu thức mẫu là một ví dụ cực đoan về điều này.) Điều này tránh sao chép không cần thiết của ma trận đầy đủ của dữ liệu trong nhiều tình huống. Thật không may, auto khá vui khi chỉ tạo một đối tượng của loại proxy, mà người dùng thông thường của thư viện sẽ không bao giờ khai báo rõ ràng. Vì bạn cần phải tính toán số liệu cuối cùng, nên không có một lần truy cập hiệu suất nào từ việc truyền đến một số MatrixXd. (Scott Meyers, trong "hiệu quả Modern C++", đưa ra ví dụ có liên quan của việc sử dụng auto với vector<bool>, nơi operator[](size_t i) trả về một proxy.)

Answer 4

KHÔNG sử dụng auto với các biểu thức Eigen. Tôi tình cờ gặp thậm chí nhiều vấn đề "ấn tượng" với điều này trước đây, thấy

eigen auto type deduction in general product

và được tư vấn bởi một trong những người sáng tạo Eigen (Gael) không sử dụng auto với các biểu thức Eigen.

Việc truyền từ biểu thức sang loại cụ thể như MatrixXd phải cực kỳ nhanh, trừ khi bạn muốn đánh giá lười (vì khi thực hiện kết quả được đánh giá).