Улучшение производительности многих операций подматрицы левого деления (mldivide, \) - PullRequest
Новая
       47

Улучшение производительности многих операций подматрицы левого деления (mldivide, \)

1 голос
/ 30 июня 2019

У меня есть две матрицы, a1 и a2. a1 - это 3x12000, а a2 - это 3x4000. Я хочу создать еще один массив размером 3x4000, который является левым матричным делением (mldivide, \) подматриц 3x3 a1 и подматриц 3x1 a2. Вы можете сделать это легко с помощью цикла for: 

for ii = 1:3:12000
  a = a1(:,ii:ii+2)\a2(:, ceil(ii/3));
end

Однако мне было интересно, есть ли более быстрый способ сделать это.

Редактировать : я знаю, что предварительное распределение увеличивает скорость, я просто показывал это для визуальных целей.

Edit2 : Убрано итеративное увеличение массива. Кажется, мои вопросы были немного неверно истолкованы. В основном мне было интересно, есть ли какие-нибудь матричные операции, которые я мог бы выполнить для достижения своей цели, поскольку это, скорее всего, было бы быстрее, чем цикл for, то есть изменение формы a1 на матрицу 3x3x4000 и a2 на матрицу 3x1x4000, а левая матрица делит каждый уровень тем не менее, за один раз вы не можете разделить матрицу с помощью трехмерных матриц.

Ответы [ 3 ]

4 голосов
/ 30 июня 2019

Вы можете создать одну систему уравнений, содержащую множество независимых «подсистем» уравнений, поместив подматрицы a1 в диагональ матрицы 12000x12000, например: 

a1(1,1) a1(1,2) a1(1,3)    0       0      0        0       0      0       
a1(2,1) a1(2,2) a1(2,3)    0       0      0        0       0      0       
a1(3,1) a1(3,2) a1(3,3)    0       0      0        0       0      0       
   0       0       0    a1(1,4) a1(1,5) a1(1,6)    0       0      0       
   0       0       0    a1(2,4) a1(2,5) a1(2,6)    0       0      0       
   0       0       0    a1(3,4) a1(3,5) a1(3,6)    0       0      0       
   0       0       0       0       0       0    a1(1,7) a1(1,8) a1(1,9)   
   0       0       0       0       0       0    a1(2,7) a1(2,8) a1(2,9)   
   0       0       0       0       0       0    a1(3,7) a1(3,8) a1(3,9)

и затем влево разделите его на a2(:).

. Это можно сделать, используя kron и разреженную матрицу, подобную этой ( source ): 

a1_kron = kron(speye(12000/3),ones(3));
a1_kron(logical(a1_kron)) = a1(:);
a = a1_kron\a2(:);
a = reshape(a, [3 12000/3]);

Advantage - Speed ​​: это примерно в 3-4 раза быстрее, чем цикл for с предварительным распределением на моем ПК.

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

Примечание : Как показано в Stefano M 's ответ , используя одну большую системуem уравнений (с использованием kron и разреженной матрицы) быстрее, чем использование цикла for (с предварительным распределением) только для очень небольшого размера подсистем уравнений (на моем ПК для числа уравнений <= 7) для больших размеровподсистем уравнений, используя цикл for быстрее. </p>

Сравнение различных методов

Я написал и запустил код для сравнения 4 различных методов решения этой проблемы:

  1. для цикла, без предварительного выделения
  2. для цикла, с предварительным распределением
  3. kron
  4. cellfun

Тест: 

n = 1200000;
a1 = rand(3,n);
a2 = rand(3,n/3);

disp('Method 1: for loop, no preallocation')
tic
a_method1 = [];
for ii = 1:3:n
  a_method1 = [a_method1 a1(:,ii:ii+2)\a2(:, ceil(ii/3))];
end
toc
disp(' ')

disp('Method 2: for loop, with preallocation')
tic
a1_reshape = reshape(a1, 3, 3, []);
a_method2 = zeros(size(a2));
for i = 1:size(a1_reshape,3)
    a_method2(:,i) = a1_reshape(:,:,i) \ a2(:,i);
end
toc
disp(' ')

disp('Method 3: kron')
tic
a1_kron = kron(speye(n/3),ones(3));
a1_kron(logical(a1_kron)) = a1(:);
a_method3 = a1_kron\a2(:);
a_method3 = reshape(a_method3, [3 n/3]);
toc
disp(' ')

disp('Method 4: cellfun')
tic
a1_cells = mat2cell(a1, size(a1, 1), repmat(3 ,1,size(a1, 2)/3));
a2_cells = mat2cell(a2, size(a2, 1), ones(1,size(a2, 2)));
a_cells = cellfun(@(x, y) x\y, a1_cells, a2_cells, 'UniformOutput', 0);
a_method4 = cell2mat(a_cells);
toc
disp(' ')

Результаты: 

Method 1: for loop, no preallocation
Elapsed time is 747.635280 seconds.

Method 2: for loop, with preallocation
Elapsed time is 1.426560 seconds.

Method 3: kron
Elapsed time is 0.357458 seconds.

Method 4: cellfun
Elapsed time is 3.390576 seconds.

Сравнивая результаты четырех методов, вы можете увидеть, что использование метода 3 - kron дает несколько разные результаты: 

disp(['sumabs(a_method1(:) - a_method2(:)): ' num2str(sumabs(a_method1(:)-a_method2(:)))])
disp(['sumabs(a_method1(:) - a_method3(:)): ' num2str(sumabs(a_method1(:)-a_method3(:)))])
disp(['sumabs(a_method1(:) - a_method4(:)): ' num2str(sumabs(a_method1(:)-a_method4(:)))])

Результат: 

sumabs(a_method1(:) - a_method2(:)): 0
sumabs(a_method1(:) - a_method3(:)): 8.9793e-05
sumabs(a_method1(:) - a_method4(:)): 0
1 голос
/ 30 июня 2019

Вы решаете серию из N систем с m линейными уравнениями в каждой, N систем имеют вид 

Ax = b

Вы можете преобразовать их в единую систему из Nm линейных уравнений: 

|A1 0  0  ... 0 | |x1|   |b1|
|0  A2 0  ... 0 | |x2|   |b2|
|0  0  A3 ... 0 | |x3| = |b3|
|.  .  .  ... . | |. |   |. |
|0  0  0  ... AN| |xN|   |bN|

Однако решение этой единой системы уравнений намного дороже, чем решение всех маленьких.Обычно стоимость составляет O (n ^ 3), поэтому вы переходите от O (N m ^ 3) к O ((Nm) ^ 3).Огромная пессимизация.( Элиаху доказал, что я здесь не прав , очевидно, можно использовать разреженность матрицы.)

Можно уменьшить вычислительные затраты, но вам необходимо предоставить гарантии относительно данных.Например, если матрицы A положительно определены, системы могут быть решены дешевле.Тем не менее, учитывая, что вы имеете дело с матрицами 3х3, выигрыш там будет небольшим, так как это довольно простые системы для решения.

Если вы спрашиваете об этом, потому что считаете, что циклы в MATLAB по своей сути медленные, выследует знать, что это больше не так, и не было, так как MATLAB получил JIT 15 лет назад.Сегодня многие попытки векторизации приводят к одинаково быстрому коду и зачастую к более медленному коду, особенно для больших данных.(Я мог бы использовать некоторые моменты времени, которые я разместил здесь на SO, чтобы доказать это при необходимости.)

Я бы подумал, что решение всех систем за один раз может сократить количество проверок, которые MATLAB делает каждый раз, когда оператор\ называется.То есть жесткое кодирование размера и типа проблемы может улучшиться повсюду.Но единственный способ сделать это - написать MEX-файл.

1 голос
/ 30 июня 2019

Marginal Самым большим улучшением было бы предварительное распределение матрицы вывода вместо ее увеличения: 

A1 = reshape(A1, 3, 3, []);
a = zeros(size(A2));
for i = 1:size(A1,3) 
    a(:,i) = A1(:,:,i) \ A2(:,i);
end

В массиве preallocate, если доступен Parallel Toolbox, вы можете попробовать parfor

Редактировать

Этот ответ более не актуален, поскольку ОП перефразировал вопрос, чтобы избежать увеличения массива результатов, который был первоначальным основным узким местом.

Проблема здесь в том, что нужно решить 4000 независимых линейных систем 3х3. Матрица настолько мала, что может представлять интерес решение ad hoc , особенно если имеется некоторая информация о свойствах матрицы (симметричная или нет, номер условия и т. Д.). Тем не менее, придерживаясь оператора \ matlab, лучший способ ускорить вычисления - это явно использовать параллелизм, например, по команде parfor.

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

С помощью этой функции вы можете исследовать различные размеры проблем: 

function [t2, t3] = sotest(n, n2)

a1 = rand(n,n*n2);
a2 = rand(n,n2);

tic
a1_reshape = reshape(a1, n, n, []);
a_method2 = zeros(size(a2));
for i = 1:size(a1_reshape,3)
    a_method2(:,i) = a1_reshape(:,:,i) \ a2(:,i);
end
t2 = toc;

tic
a1_kron = kron(speye(n2),ones(n));
a1_kron(logical(a1_kron)) = a1(:);
a_method3 = a1_kron\a2(:);
a_method3 = reshape(a_method3, [n n2]);
t3 = toc;

assert ( norm(a_method2 - a_method3, 1) / norm(a_method2, 1) < 1e-8)

Действительно, для n=3 метод разреженной матрицы явно лучше, но для увеличения n он становится менее конкурентоспособным

solution times

Приведенный выше показатель был получен с 

>> for i=1:20; [t(i,1), t(i,2)] = sotest(i, 50000); end
>> loglog(1:20, t, '*-')

Мой последний комментарий заключается в том, что явный цикл с оператором \ действительно быстрый; формулировка разреженной матрицы немного менее точна и может стать проблематичной в крайних случаях; и наверняка разреженное матричное решение не очень читабельно. Если число n2 систем для решения очень и очень велико (> 1e6), то, возможно, следует изучить ad hoc решений.

...