Портирование инверсии 1000x1000 матрицы комплексных чисел из Matlab на C++

Question

код Matlab:

invD = inv(D); 
Dew=2*invD-E; 

D является 1000x1000 матрицей комплексных чисел. Он инвертируется через 0,5 сек или меньше.

Я пробовал несколько способов портировать этот код на C++, используя разные подходы, но он всегда медленнее, чем 10 секунд. Что я делаю не так? Как возможно, Matlab оптимизирует этот код?

Я не написал этот код, я просто должен его портировать. Извините, мне гораздо лучше в кодировании, что в математике.

Answer 1

Я уверен, что вы используете Eigen неправильно. Я попытался инвертировать случайную матрицу 1000x1000 с октавой (бесплатная альтернатива MATLAB, возможно, немного менее оптимизированная) и со следующим кодом Eigen. Октав занимает около 1 секунды, Эйген о 1.5s (Я использую boost::timer вместо вашего решения синхронизации, но только потому, что это меньше хлопот):

#include <iostream> 
#include <Eigen/Dense> 
#include <boost/timer/timer.hpp> 

using Eigen::MatrixXcd; 

int main() { 
    MatrixXcd m,mi; 
    { 
    boost::timer::auto_cpu_timer t; 
    m = MatrixXcd::Random(1000,1000); 
    } 
    { 
    boost::timer::auto_cpu_timer t; 
    mi = m.inverse(); 
    } 
    std::cout << mi(4,4) << std::endl; 
} 

Первый возможный виновник я могу думать являются опциями компилятора. Вы скомпилировали свой код на C++ с -O2?

По-прежнему существует вероятность того, что ваши матрицы имеют особую структуру, так что численные соображения замедляют алгоритм Eigen в вашем конкретном случае, но прежде чем рассматривать это, я бы посмотрел в другом месте проблемы.

Answer 2

Вместо того, чтобы изобретать колесо, вы можете просто использовать boost::ublas и следовать за these instructions.

Answer 3

Мне непонятно, что вы делаете с этим обратным. Моя первая реакция заключается в том, что большинство людей не инвертируют матрицы. Например, декомпозиция LU намного эффективнее и стабильнее, чем полная инверсия.

Я предлагаю найти linear algebra library written in C++. Другие люди, намного превосходящие вас по математике и программированию, решили эту проблему. Пользователь их код; не пишите свои собственные.

Я также рекомендую посмотреть на CUDA, если у вас есть GPU.

1000x1000 не является потрясающе большой матрицей. Решение проблем с миллионами степеней свободы не является неслыханным.