Самый быстрый способ добавить несколько разреженных матриц в цикле в MATLAB

Question

У меня есть код, который многократно вычисляет разреженную матрицу в цикле (он выполняет этот расчет 13472 раза, если быть точным). Каждая из этих разреженных матриц уникальна.

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

Когда все матрицы 13742 добавлены, код выходит из цикла, и программа завершается.

Недостаток кода возникает при добавлении разреженных матриц. Я сделал фиктивную версию кода, который демонстрирует то же поведение, что и мой настоящий код. Он состоит из функции MATLAB и сценария, приведенного ниже.

(1) Функция, которая генерирует разреженную матрицу:

function out = test_evaluate_stiffness(n) 
ind = randi([1 n*n],300,1); 
val = rand(300,1); 
[I,J] = ind2sub([n,n],ind); 
out = sparse(I,J,val,n,n); 
end 

(2) Основной скрипт (программа)

% Calculate the stiffness matrix 
n=1000; 
K=sparse([],[],[],n,n,n^2); 
tic 
for i=1:13472 
    temp=rand(1)*test_evaluate_stiffness(n); 
    K=K+temp; 
end 
fprintf('Stiffness Calculation Complete\nTime taken = %f s\n',toc) 

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

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

Отчет профайлер также приводится ниже:

Answer 1

Если вам нужна только сумма этих матриц, вместо того, чтобы строить все из них в отдельности, а затем суммируя их, просто сцепить векторы I, J и vals и звоните sparse только один раз. Если есть повторяющиеся строки [i,j] в [I,J], соответствующие значения S(i,j) будут суммированы автоматически, поэтому код абсолютно эквивалентен. Поскольку вызов sparse предполагает внутренний вызов алгоритма сортировки, вы сохраняете промежуточные сорта 13742-1 и можете уйти с одного.

---

Это предполагает изменение подписи test_evaluate_stiffness для вывода [I,J,val]:

function [I,J,val] = test_evaluate_stiffness(n) 

и удаление линии out = sparse(I,J,val,n,n);.

Вы измените свою другую функцию:

n = 1000; 
[I,J,V] = deal([]); 
tic; 
for i = 1:13472 
    [I_i, J_i, V_i] = test_evaluate_stiffness(n); 
    nE = numel(I_i); 
    I(end+(1:nE)) = I_i; 
    J(end+(1:nE)) = J_i; 
    V(end+(1:nE)) = rand(1)*V_i; 
end 
K = sparse(I,J,V,n,n); 
fprintf('Stiffness Calculation Complete\nTime taken = %f s\n',toc); 

Если вы знаете длину выхода test_evaluate_stiffness раньше времени, возможно, вы можете сэкономить время, предварительное выделение массивов I, J и V с соответствующего размера zeros матрицы и установить их, используя нечто вроде:

I((i-1)*nE + (1:nE)) = ... 
J((i-1)*nE + (1:nE)) = ... 
V((i-1)*nE + (1:nE)) = ... 

Answer 2

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

Почти ... но один заключительный трюк: если вы можете создать векторы так, что J сортируется в порядке возрастания, то вы будете значительно улучшить скорость sparse вызова, о факторе 4 в моем опыте.

(Если это легче иметь I сортируют, а затем создать транспонированную матрицу sparse(J,I,V) и снимите транспонировать его впоследствии.)