производительность — Почему MATLAB / Octave стирает слово с C ++ в задачах на собственные значения?

Я надеюсь, что ответ на вопрос в названии заключается в том, что я делаю что-то глупое!

Здесь проблема. Я хочу вычислить все собственные значения и собственные векторы реальной, симметричной матрицы. Я реализовал код в MATLAB (на самом деле, я запускаю его с помощью Octave) и C ++, используя Научная библиотека GNU. Я предоставляю мой полный код ниже для обеих реализаций.

Насколько я понимаю, GSL поставляется с собственной реализацией API BLAS (в дальнейшем я называю это GSLCBLAS), и для использования этой библиотеки я компилирую, используя:

g++ -O3 -lgsl -lgslcblas

GSL предлагает Вот использовать альтернативную библиотеку BLAS, такую ​​как самооптимизирующаяся АТЛАС библиотека, для повышения производительности. Я использую Ubuntu 12.04 и установил пакеты ATLAS из Ubuntu репозиторий. В этом случае я компилирую, используя:

g++ -O3 -lgsl -lcblas -latlas -lm

Для всех трех случаев я провел эксперименты со случайно сгенерированными матрицами размеров от 100 до 1000 с шагом 100. Для каждого размера я выполнил 10 собственных разложений с различными матрицами и усреднил время, затраченное на это. Результаты таковы:

График результатов

Разница в производительности просто смешная. Для матрицы размером 1000 Октава выполняет разложение менее чем за секунду; GSLCBLAS и ATLAS занимают около 25 секунд.

Я подозреваю, что я могу использовать библиотеку ATLAS неправильно. Любые объяснения приветствуются; заранее спасибо.

Некоторые примечания по коду:

  • В реализации C ++ нет необходимости делать матрицу
    симметричный, потому что функция использует только нижнюю треугольную часть
    из этого
    .

  • В Октаве линия triu(A) + triu(A, 1)' заставляет матрицу быть симметричной.

  • Если вы хотите скомпилировать код C ++ на своем собственном компьютере с Linux, вам также необходимо добавить флаг -lrtиз-за clock_gettime функция.

  • К сожалению я не думаю clock_gettime выходы на другие платформы. Рассмотрите возможность изменения на gettimeofday,

Октавный код

K = 10;

fileID = fopen('octave_out.txt','w');

for N = 100:100:1000
AverageTime = 0.0;

for k = 1:K
A = randn(N, N);
A = triu(A) + triu(A, 1)';
tic;
eig(A);
AverageTime = AverageTime + toc/K;
end

disp([num2str(N), " ", num2str(AverageTime), "\n"]);
fprintf(fileID, '%d %f\n', N, AverageTime);
end

fclose(fileID);

Код C ++

#include <iostream>
#include <fstream>
#include <time.h>

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>

int main()
{
const int K = 10;

gsl_rng * RandomNumberGenerator = gsl_rng_alloc(gsl_rng_default);
gsl_rng_set(RandomNumberGenerator, 0);

std::ofstream OutputFile("atlas.txt", std::ios::trunc);

for (int N = 100; N <= 1000; N += 100)
{
gsl_matrix* A = gsl_matrix_alloc(N, N);
gsl_eigen_symmv_workspace* EigendecompositionWorkspace = gsl_eigen_symmv_alloc(N);
gsl_vector* Eigenvalues = gsl_vector_alloc(N);
gsl_matrix* Eigenvectors = gsl_matrix_alloc(N, N);

double AverageTime = 0.0;
for (int k = 0; k < K; k++)
{
for (int i = 0; i < N; i++)
{
for (int j = 0; j < N; j++)
{
gsl_matrix_set(A, i, j, gsl_ran_gaussian(RandomNumberGenerator, 1.0));
}
}

timespec start, end;
clock_gettime(CLOCK_MONOTONIC_RAW, &start);

gsl_eigen_symmv(A, Eigenvalues, Eigenvectors, EigendecompositionWorkspace);

clock_gettime(CLOCK_MONOTONIC_RAW, &end);
double TimeElapsed = (double) ((1e9*end.tv_sec + end.tv_nsec) - (1e9*start.tv_sec + start.tv_nsec))/1.0e9;
AverageTime += TimeElapsed/K;
std::cout << "N = " << N << ", k = " << k << ", Time = " << TimeElapsed << std::endl;
}
OutputFile << N << " " << AverageTime << std::endl;

gsl_matrix_free(A);
gsl_eigen_symmv_free(EigendecompositionWorkspace);
gsl_vector_free(Eigenvalues);
gsl_matrix_free(Eigenvectors);
}

return 0;
}

15

Решение

Я не согласен с предыдущим постом. Это не проблема многопоточности, это проблема алгоритма. Причина, по которой matlab, R и octave стирают слово с помощью библиотек C ++, заключается в том, что их библиотеки C ++ используют более сложные и лучшие алгоритмы. Если вы прочитаете страницу октавы, вы можете узнать, что Oни делать [1]:

Собственные значения вычисляются в несколько этапов процесса, который начинается с разложения Гессенберга, за которым следует разложение Шура, из которого очевидны собственные значения. Собственные векторы, при желании, вычисляются путем дальнейших манипуляций с разложением Шура.

Решение задач на собственные значения / собственные векторы нетривиально. На самом деле это одна из немногих вещей, которые «Численные рецепты в Си» рекомендует вам не реализовать себя. (P461). GSL часто медленный, что было моим первоначальным ответом. ALGLIB также медленен для своей стандартной реализации (я получаю около 12 секунд!):

#include <iostream>
#include <iomanip>
#include <ctime>

#include <linalg.h>

using std::cout;
using std::setw;
using std::endl;

const int VERBOSE = false;

int main(int argc, char** argv)
{

int size = 0;
if(argc != 2) {
cout << "Please provide a size of input" << endl;
return -1;
} else {
size = atoi(argv[1]);
cout << "Array Size: " << size << endl;
}

alglib::real_2d_array mat;
alglib::hqrndstate state;
alglib::hqrndrandomize(state);
mat.setlength(size, size);
for(int rr = 0 ; rr < mat.rows(); rr++) {
for(int cc = 0 ; cc < mat.cols(); cc++) {
mat[rr][cc] = mat[cc][rr] = alglib::hqrndnormal(state);
}
}

if(VERBOSE) {
cout << "Matrix: " << endl;
for(int rr = 0 ; rr < mat.rows(); rr++) {
for(int cc = 0 ; cc < mat.cols(); cc++) {
cout << setw(10) << mat[rr][cc];
}
cout << endl;
}
cout << endl;
}

alglib::real_1d_array d;
alglib::real_2d_array z;
auto t = clock();
alglib::smatrixevd(mat, mat.rows(), 1, 0, d, z);
t = clock() - t;

cout << (double)t/CLOCKS_PER_SEC << "s" << endl;

if(VERBOSE) {
for(int cc = 0 ; cc < mat.cols(); cc++) {
cout << "lambda: " << d[cc] << endl;
cout << "V: ";
for(int rr = 0 ; rr < mat.rows(); rr++) {
cout << setw(10) << z[rr][cc];
}
cout << endl;
}
}
}

Если вам действительно нужна быстрая библиотека, вероятно, нужно заняться настоящей охотой.

[1] http://www.gnu.org/software/octave/doc/interpreter/Basic-Matrix-Functions.html

6

Другие решения

Я также столкнулся с проблемой. Настоящая причина в том, что eig () в matlab не вычисляет собственные векторы, а код версии C выше. Разница во времени может быть больше, чем на один порядок, как показано на рисунке ниже. Так что сравнение не честно.

В Matlab, в зависимости от возвращаемого значения, фактическая вызванная функция будет отличаться. Для принудительного расчета собственных векторов [V, D] = EIG (A) следует использовать (см. коды ниже).

Фактическое время для вычисления задачи на собственные значения сильно зависит от свойств матрицы и желаемых результатов, таких как

  • Реальный или сложный
  • Эрмитово / симметрично или нет
  • Плотный или редкий
  • Только собственные значения, Собственные векторы, Только максимальное собственное значение и т. Д.
  • Последовательный или параллельный

Существуют алгоритмы, оптимизированные для каждого из приведенных выше случаев. В gsl эти алгоритмы выбрано вручную, поэтому неправильный выбор значительно снизит производительность. Некоторые классы-оболочки C ++ или некоторые языки, такие как matlab и mathematica, выберут оптимизированную версию с помощью некоторых методов.

Кроме того, Matlab и Mathematica использовали распараллеливание. Это еще больше увеличивает зазор, который вы видите в несколько раз, в зависимости от машины. Разумно сказать, что вычисление собственных значений и собственных векторов общего комплекса 1000×1000 составляет около секунды и десять секунд без распараллеливания.

Сравните Matlab и C
Рис. Сравните Matlab и C. «+ vec» означает, что коды включали вычисления собственных векторов. CPU% — это очень грубое наблюдение за использованием CPU при N = 1000, которое ограничено верхним пределом на 800%, хотя предполагается, что они полностью используют все 8 ядер. Разрыв между Matlab и C меньше, чем в 8 раз.

Сравните другой тип матрицы в Mathematica
Рис. Сравните различные типы матриц в Mathematica. Алгоритмы автоматически выбираются программой.

Matlab (с вычислением собственных векторов)

K = 10;

fileID = fopen('octave_out.txt','w');

for N = 100:100:1000
AverageTime = 0.0;

for k = 1:K
A = randn(N, N);
A = triu(A) + triu(A, 1)';
tic;
[V,D] = eig(A);
AverageTime = AverageTime + toc/K;
end

disp([num2str(N), ' ', num2str(AverageTime), '\n']);
fprintf(fileID, '%d %f\n', N, AverageTime);
end

fclose(fileID);

C ++ (БЕЗ расчета собственных векторов)

#include <iostream>
#include <fstream>
#include <time.h>

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>

int main()
{
const int K = 10;

gsl_rng * RandomNumberGenerator = gsl_rng_alloc(gsl_rng_default);
gsl_rng_set(RandomNumberGenerator, 0);

std::ofstream OutputFile("atlas.txt", std::ios::trunc);

for (int N = 100; N <= 1000; N += 100)
{
gsl_matrix* A = gsl_matrix_alloc(N, N);
gsl_eigen_symm_workspace* EigendecompositionWorkspace = gsl_eigen_symm_alloc(N);
gsl_vector* Eigenvalues = gsl_vector_alloc(N);

double AverageTime = 0.0;
for (int k = 0; k < K; k++)
{
for (int i = 0; i < N; i++)
{
for (int j = i; j < N; j++)
{
double rn = gsl_ran_gaussian(RandomNumberGenerator, 1.0);
gsl_matrix_set(A, i, j, rn);
gsl_matrix_set(A, j, i, rn);
}
}

timespec start, end;
clock_gettime(CLOCK_MONOTONIC_RAW, &start);

gsl_eigen_symm(A, Eigenvalues, EigendecompositionWorkspace);

clock_gettime(CLOCK_MONOTONIC_RAW, &end);
double TimeElapsed = (double) ((1e9*end.tv_sec + end.tv_nsec) - (1e9*start.tv_sec + start.tv_nsec))/1.0e9;
AverageTime += TimeElapsed/K;
std::cout << "N = " << N << ", k = " << k << ", Time = " << TimeElapsed << std::endl;
}
OutputFile << N << " " << AverageTime << std::endl;

gsl_matrix_free(A);
gsl_eigen_symm_free(EigendecompositionWorkspace);
gsl_vector_free(Eigenvalues);
}

return 0;
}

Mathematica

(* Symmetric real matrix + eigenvectors *)
Table[{NN, Mean[Table[(
M = Table[Random[], {i, NN}, {j, NN}];
M = M + Transpose[Conjugate[M]];
AbsoluteTiming[Eigensystem[M]][[1]]
), {K, 10}]]
}, {NN, Range[100, 1000, 100]}]

(* Symmetric real matrix *)
Table[{NN, Mean[Table[(
M = Table[Random[], {i, NN}, {j, NN}];
M = M + Transpose[Conjugate[M]];
AbsoluteTiming[Eigenvalues[M]][[1]]
), {K, 10}]]
}, {NN, Range[100, 1000, 100]}]

(* Asymmetric real matrix *)
Table[{NN, Mean[Table[(
M = Table[Random[], {i, NN}, {j, NN}];
AbsoluteTiming[Eigenvalues[M]][[1]]
), {K, 10}]]
}, {NN, Range[100, 1000, 100]}]

(* Hermitian matrix *)
Table[{NN, Mean[Table[(
M = Table[Random[] + I Random[], {i, NN}, {j, NN}];
M = M + Transpose[Conjugate[M]];
AbsoluteTiming[Eigenvalues[M]][[1]]
), {K, 10}]]
}, {NN, Range[100, 1000, 100]}]

(* Random complex matrix *)
Table[{NN, Mean[Table[(
M = Table[Random[] + I Random[], {i, NN}, {j, NN}];
AbsoluteTiming[Eigenvalues[M]][[1]]
), {K, 10}]]
}, {NN, Range[100, 1000, 100]}]
3

В реализации C ++ нет необходимости делать матрицу
симметричный, потому что функция использует только нижнюю треугольную часть
Это.

Это может быть не так. в ссылка, заявлено, что:

int gsl_eigen_symmv (gsl_matrix * A, gsl_vector * eval, gsl_matrix * evec, gsl_eigen_symmv_workspace * w)

Эта функция вычисляет собственные значения и собственные векторы реальная симметричная матрица
. Дополнительное рабочее пространство соответствующего размера должно быть предоставлено в w.
Диагональная и нижняя треугольная часть А разрушаются во время
вычисление, но строгая верхняя треугольная часть не упоминается.
Собственные значения хранятся в векторе eval и неупорядочены.
соответствующие собственные векторы хранятся в столбцах матрицы
evec. Например, собственный вектор в первом столбце соответствует
первое собственное значение. Собственные векторы гарантированно будут взаимно
ортогональный и нормированный на единицу величины.

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

На стороне MATLAB декомпозиция собственных значений может быть быстрее из-за многопоточного выполнения, как указано в эта ссылка:

Встроенная многопоточность

Линейная алгебра и числовые функции, такие как fft, \ (mldivide), eig,
svd и sort многопоточные в MATLAB. Многопоточные вычисления
был включен по умолчанию в MATLAB с момента выпуска 2008a. Эти
функции автоматически выполняются в нескольких вычислительных потоках в
один сеанс MATLAB, что позволяет им выполняться быстрее на
многоядерные машины. Кроме того, многие функции в изображении
Processing Toolbox ™ являются многопоточными.

Чтобы протестировать производительность MATLAB для одного ядра, вы можете отключить многопоточность:

Файл> Настройки> Общие> Многопоточность

в R2007a или новее, как указано Вот.

1