Решение нескольких линейных систем с использованием векторизации

Question

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

Я пытаюсь решить линейные системы вида Ax = B с A матрицей 4x4 и B вектором 4x1.

Я знаю, что для одной системы я могу использовать mldivide, чтобы получить x : x=A\B.

Однако я пытаюсь решить множество систем (возможно,> 10000), и я неохотно использую цикл for, потому что мне сказали, что он заметно медленнее, чем матричная формулировка, во многих задачах MATLAB.

Тогда у меня вопрос: есть ли способ решить Ax = B , используя векторизацию с A 4x4x N и B a матрица 4x N ?

PS: я не знаю, важно ли это, но вектор B одинаков для всех систем.

Answer 1

Я не думаю, что вы можете оптимизировать это дальше. Как объясняет @Tom, так как A меняется, нет никакой пользы в том, чтобы заранее учитывать различные A ...

Кроме того, зацикленное решение довольно быстрое, учитывая указанные вами размеры:

A = rand(4,4,10000);
B = rand(4,1);          %# same for all linear systems

tic
X = zeros(4,size(A,3));
for i=1:size(A,3)
    X(:,i) = A(:,:,i)\B;
end
toc

Прошедшее время составляет 0,168101 секунды.

Answer 2

Вы должны использовать цикл for. Может быть полезным предварительное вычисление факторизации и ее повторное использование, если A остается неизменным и B изменяется. Но для вашей проблемы, когда A изменяется и B остается неизменным, альтернативы решению N линейных систем не существует.

Вам также не стоит слишком беспокоиться о производительности циклов: компилятор MATLAB JIT означает, что циклы часто могут быть такими же быстрыми в последних версиях MATLAB.

Answer 3

Вот довольно эзотерическое решение, использующее преимущества специфической оптимизации MATLAB. Построить огромную матрицу 4k x 4k разреженную с вашими блоками 4x4 по диагонали. Тогда решай все одновременно.

На моей машине получается то же решение с точностью до одинарной точности, что и у @ Amro / Tom's for-loop, но быстрее.

n = size(A,1);
k = size(A,3);
AS = reshape(permute(A,[1 3 2]),n*k,n);
S = sparse( ...
  repmat(1:n*k,n,1)', ...
  bsxfun(@plus,reshape(repmat(1:n:n*k,n,1),[],1),0:n-1), ...
  AS, ...
  n*k,n*k);
X = reshape(S\repmat(B,k,1),n,k);

для случайного примера:

For k = 10000
For loop: 0.122570 seconds.
Giant sparse system: 0.032287 seconds.

Если вы знаете, что ваши матрицы 4x4 являются положительно определенными, то вы можете использовать chol на S для повышения точности.

Это глупо . Но так же, как медленно в Matlab для циклов все еще в 2015 году, даже с JIT. Это решение, кажется, находит хорошее место, когда k не слишком велико, поэтому все по-прежнему помещается в память.

Answer 4

Вот проблема: вы пытаетесь выполнить 2D-операцию ( mldivide ) на 3D-матрице.Независимо от того, как вы смотрите на это, вам нужно ссылаться на матрицу по индексу, в котором начинается штраф времени ... Проблема не в цикле for, а в том, как люди его используют.Если вы можете по-разному структурировать вашу проблему, то, возможно, вы сможете найти лучший вариант, но сейчас у вас есть несколько вариантов:

1 - mex

2 - параллельная обработка (написать цикл parfor)

3 - CUDA

Answer 5

Я знаю, что этому посту уже лет, но я все равно внесу свои два цента. Вы МОЖЕТЕ поместить все свои матрицы A в большую диагональную матрицу блоков, где на диагонали большой матрицы будет 4x4 блока. В правой части будут все ваши b-векторы, сложенные друг на друга снова и снова. После того, как вы это настроите, она будет представлена ​​как разреженная система и может быть эффективно решена с помощью алгоритмов, выбранных mldivide. Блоки численно разделены, поэтому даже если там есть единичные блоки, ответы на неособые блоки должны быть правильными, когда вы используете mldivide. Существует код, который использовал этот подход в MATLAB Central:

http://www.mathworks.com/matlabcentral/fileexchange/24260-multiple-same-size-linear-solver

Я предлагаю поэкспериментировать, чтобы увидеть, является ли подход более быстрым, чем зацикливание. Я подозреваю, что это может быть более эффективным, особенно для большого количества небольших систем. В частности, если есть хорошие формулы для коэффициентов A в N матрицах, вы можете построить полную левую часть, используя векторные операции MATLAB (без циклов), что может дать вам дополнительную экономию затрат. Как уже отмечали другие, векторизованные операции не всегда выполняются быстрее, но, как мне кажется, они часто происходят.