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

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

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

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

Я использую этот оператор для составления матрицы A и вектора 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;

Где z - некоторый скаляр

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

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


c++ eigen performance matrix eigen3
person Dark_Daiver    schedule 17.07.2015    source источник

Ответы (3)


arrow_upward
3
arrow_downward

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

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

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

person ggael    schedule 17.07.2015
comment
Спасибо за ответ! Выглядит очень интересно. Я попробую подход с огромными матрицами. Также я пытался использовать матрицы 8x8/8x1 раньше, но это не дает мне ускорения. - person Dark_Daiver; 17.07.2015

arrow_upward
1
arrow_downward

Если вы можете использовать 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). Однако использование тяги не так просто и требует серьезного обдумывания.

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

person kn0t3k    schedule 17.07.2015
comment
Спасибо за ответ. У меня нет доступа к продвинутым компиляторам, таким как icc. Но, может быть, я попробую использовать openmp. Также я думаю, что можно использовать openmp со старой версией gcc. - person Dark_Daiver; 17.07.2015

arrow_upward
0
arrow_downward

У меня есть единственная проблема Ax=b с 480k переменных с плавающей запятой. Матрица A разреженная, и ее решение с помощью Eigen BiCGSTAB заняло 4,8 секунды.

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

person Deepfreeze    schedule 17.07.2015
comment
Спасибо за ответ. Но у меня много очень маленьких плотных проблем, а не одна большая и редкая. Также я не могу использовать GPU сейчас =( - person Dark_Daiver; 17.07.2015