Как рассчитать это в Matlab без цикла For

Question

Так как matlab работает медленно при выполнении цикла for, я обычно избегаю цикла for для всех моих кодов и превращаю их в матричный расчет, который был бы быстрым. Но вот проблема, которую я не могу найти умным способом:

У меня тревожная матрица

A=[a1,a2,a3,...,an], 

здесь a1,a2,a3....an - столбцы матрицы.

Еще одна матрица nxn

B=[b1,b2,b3,...,bn], 

аналогично b1,b2,b3... также являются столбцами B.

А также тревожная матрица М.

Я хочу вычислить матрицу NXN

C=[c1,c2,c3,...,cn], 

thus (M+diag(ai))*ci = bi.

а именно

ci = (M+diag(ai))\bi.

Я знаю один способ без цикла for:

C(:)=( blkdiag(M)+diag(A(:)) )\B(:).

Но это будет делать слишком много вычислений, чем необходимо.

Какие-нибудь умные решения? Вы можете предположить, что в вычислении нет проблемы сингулярности.

Answer 1

Утверждение "циклы медленны в Matlab" больше не соответствует действительности, так как Matlab...euhm, R2008a? (кто-то, пожалуйста, сообщите мне об этом:)

В любом случае, попробуйте это:

clc
clear all

M = rand(50);
a = rand(50);
b = rand(50);

% simple loop approach
tic
c = zeros(size(b));
for ii = 1:size(b,2)
    c(:,ii) = ( M+diag(a(:,ii)) ) \ b(:,ii);
end
toc

% not-so-simple vectorized approach
tic
MM = repmat({M}, 50,1);
c2 = (blkdiag(MM{:})+diag(a(:)))\b(:);
toc

norm(c(:)-c2(:))

Результаты:

Elapsed time is 0.011226 seconds.  % loop
Elapsed time is 1.049130 seconds.  % no-loop
ans =
     5.091221148787843e-10         % results are indeed "equal"

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

Некоторые проблемы просто не подходят для векторизации. Я думаю, что это один.