Eigen: как ускорить a = = coeffs * coeffs.transpose()

Question

Мне нужно вычислить множество (около 400 тыс.) Решений небольших линейных наименьших квадратов. Каждая проблема содержит 10-300 уравнений с 7 переменными. Для решения этих проблем я использую собственную библиотеку. Прямое решение занимает слишком много времени, и я превращаю каждую проблему в решение системы линейных уравнений 7x7, производя мои производные.

Я получаю приятное ускорение, но я хочу увеличить производительность снова.

Я использую vagrind для профилирования своей программы, и я обнаружил, что операция с самой высокой стоимостью самообслуживания является оператором + = собственной матрицы. Эта операция занимает более десяти вызовов a.ldlt(). Solve (b);

Я использую этот оператор, чтобы составить вектор матрицы и B каждой системы уравнений

//I cal these code to solve each problem 
const int nVars = 7; 
//i really need double precision 
Eigen::Matrix<double, nVars, nVars> a = Eigen::Matrix<double, nVars, nVars>::Zero(); 
Eigen::Matrix<double, nVars, 1> b = Eigen::Matrix<double, nVars, 1>::Zero(); 
Eigen::Matrix<double, nVars, 1> equationCoeffs; 
//............................ 
//Somewhere in big cycle. 
//equationCoeffs and z are updated on each iteration 
a += equationCoeffs * equationCoeffs.transpose(); 
b += equationCoeffs * z; 

Где г некоторых скалярной

Так что мой вопрос: Как я могу улучшить производительность этих операций?

PS Извините за мой плохой английский

Answer 1

У меня есть одна проблема Ах = Ь с 480K флоат переменных. Матрица А разрежена и ее решение с Eigen BiCGSTAB заняло 4,8 секунды.

Я также работал с ViennaCL раньше, поэтому я попытался решить ту же проблему, и потребовалось всего 1,2 секунды. Увеличение spead осуществляется обработкой на GPU.

Answer 2

Вместо того, чтобы составлять матрицу и векторные компоненты нормального уравнения вручную, по одному уравнению за раз, вы можете попытаться выделить достаточно большую матрицу один раз (например, 300 х 7) для хранения всех коэффициентов, а затем позволить оптимизировать Eigen ядра продукта матрицы-матрица сделать работу для вас:

Matrix<double,Dynamic,nbVars> D(300,nbVars); 
VectorXd f(300); 
for(...) 
{ 
    int nb_equations = ...; 
    for(i=0..nb_equations-1) 
    { 
    D.row(i) = equationCoeffs; 
    f(i) = z; 
    } 
    a = D.topRows(nb_equations).transpose() * D.topRows(nb_equations); 
    b = D.topRows(nb_equations).transpose() * f.head(nb_equations); 
    // solve ax=b 
} 

Вы могли бы скамейку как с столбцами и строками основных хранением матрицы D, чтобы увидеть, какой из них лучше.

Другим возможным подходом было бы объявить a, equationCoeffs и b, как 8х8 или 8х1 матрицы или векторы, убедившись, что equationCoeffs(7)==0. Таким образом, вы максимизируете использование SIMD. Затем при вызове LDLT используйте a.topLeftCorners<7,7>(), b.head<7>(). Вы даже можете объединить эту стратегию с предыдущей.

И, наконец, если ваш процессор поддерживает AVX или FMA, вы можете использовать ветвь devel и скомпилировать с -mavx или -mfma, чтобы получить значительное ускорение.

Answer 3

Если вы можете использовать g ++ 5.1, вы можете взглянуть на OpenMP (http://openmp.org/mp-documents/OpenMP4.0.0.Examples.pdf). G ++ 5.1 (или gcc5.1 для C) также имеет базовую поддержку OpenACC, вы также можете попробовать это. В будущем должна быть больше реализации OpenACC.

Также, если у вас есть доступ к компилятору intel (icc, icpc), он ускорил мой код, даже используя его.

Если вы можете использовать NVCC NVidia, вы можете использовать библиотеку упорную (http://docs.nvidia.com/cuda/thrust/#axzz3g8xJPGHe), есть много примеров кода на их GitHub, а (https://github.com/thrust/thrust).Однако использование тяги не так прямолинейно и требует некоторого реального мышления.

РЕДАКТИРОВАТЬ: Укороченный также требует Nvidia GPU. Для карт AMD я считаю, что есть библиотека под названием ArrayFire, которая очень похожа на Thrust (я еще не пробовал этот).