Улучшение кода построения матрицы MATLAB: или векторизация кода для начинающих

Как выполнять векторизацию, не зная, как это делать:

Во-первых, я расскажу только о векторизации первого двойного цикла, вы можете следовать той же логике для второго цикла.

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

Во-первых, я рекомендую «массировать» код в простых случаях, чтобы почувствовать его. Например, я использовал n=3 и v=1:6 и запустил первый цикл один раз, вот как выглядит B:

[N M]=size(B)
N =
     9
M =
    27

imagesc(B); 

Итак, вы можете видеть, что мы получаем матрицу, похожую на лестницу, которая довольно регулярна! Все, что нам нужно, это присвоить правильный индекс матрицы правильному значению v и все.

Есть много способов добиться этого, некоторые из них более элегантны, чем другие. Один из самых простых способов — использовать функцию find:

pos=[find(B==v(1)),find(B==v(2)),find(B==v(3)),...
     find(B==v(4)),find(B==v(5)),find(B==v(6))]

pos =
     1    10    19    28    37    46
    29    38    47    56    65    74
    57    66    75    84    93   102
    85    94   103   112   121   130
   113   122   131   140   149   158
   141   150   159   168   177   186
   169   178   187   196   205   214
   197   206   215   224   233   242

Приведенные выше значения являются линейными индексами матрицы B, где были найдены значения v. Каждый столбец представляет собой линейный индекс конкретное значение v в B. Например, все индексы [1 29 57 ...] содержат значение v(1) и т. д. Каждая строка полностью содержит v, поэтому индексы [29 38 47 56 65 74] содержат v=[v(1) v(2) ... v(6)]. Вы можете заметить, что для каждой строки разница между индексами равна 9, или каждый индекс разделен размером шага N, а их 6, что равно длине вектора v (тоже полученного numel(v)). Для столбцов разница между соседними элементами равна 28, или размер шага равен M+1.

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

B([1:N:numel(v)*N]+(M+1)*0)=v;
B([1:N:numel(v)*N]+(M+1)*1)=v;
...
B([1:N:numel(v)*N]+(M+1)*(N-2))=v;

Но это нецелесообразно для больших N-2, поэтому вы можете сделать это в цикле for, если хотите:

for kk=0:N-2;
     B([1:N:numel(v)*N]+(M+1)*kk)=v;
end

Matlab предлагает более эффективный способ получить все индексы сразу, используя bsxfun (это заменяет цикл for), например:

ind=bsxfun(@plus,1:N:N*numel(v),[0:(M+1):M*(N-1)+1]')

так что теперь мы можем использовать ind для присвоения v матрице N-1 раз. Для этого нам нужно «сгладить» ind в вектор-строку:

ind=reshape(ind.',1,[]);

и соедините v с самим собой N-1 раз (или сделайте еще N-1 копий v):

vec=repmat(v,[1 N-1]);

наконец получаем ответ:

B(ind)=vec;

Короче говоря, и записывая все это компактно, мы получаем решение из 2 строк (данный размер B уже известен: [N M]=size(B)):

ind=bsxfun(@plus,1:N:N*numel(v),[0:(M+1):M*(N-1)+1]');
B(reshape(ind.',1,[]))=repmat(v,[1 N-1]);

Для n=9 векторизованный код на моей машине на ~850 быстрее. (маленький n будет менее значимым)

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

N=3^(n-1);
M=3^n;
S=sparse([],[],[],N,M);
ind=bsxfun(@plus,1:N:N*numel(v),[0:(M+1):M*(N-1)+1]');
S(reshape(ind.',1,[]))=repmat(v,[1 N-1]);

Для n=10 я мог запустить только код разреженной матрицы (в противном случае нехватка памяти), и на моей машине это заняло ~ 6 секунд.

Теперь попробуйте применить это ко второму циклу...

Хотя ваша матрица имеет огромные размеры, она также очень «разреженная», что означает, что большинство ее элементов являются нулями. Чтобы повысить производительность, вы можете использовать поддержку разреженных матриц MATLAB, гарантируя, что вы работаете только с ненулевыми частями вашей матрицы.

Разреженные матрицы в MATLAB можно эффективно построить, создав координатную форму разреженной матрицы. Это означает, что должны быть определены три массива, определяющие строку, столбец и значение для каждой ненулевой записи в матрице. Это означает, что вместо того, чтобы присваивать значения с помощью обычного синтаксиса A(i,j) = x, мы вместо этого добавили бы эту ненулевую запись в нашу разреженную структуру индексации:

row(pos+1) = i;
col(pos+1) = j;
val(pos+1) = x;
% pos is the current position within the sparse indexing arrays!

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

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

rows = 3^(n-1);
cols = 3^(n+0);

% setup the sparse indexing arrays for non-
% zero elements of matrix B
row = zeros(rows*6,1);
col = zeros(rows*6,1);
val = zeros(rows*6,1);
pos = +0;

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

Я оставил логику для последней строки для вас!

for j = 1 : 3^(n-1)
    if (j < 3^(n-1))

% add entries for a general row
    for k = max(1,3*(j-1)+1) : min(3^n,3*(j-1)+6)             
        pos = pos+1;
        row(pos) = j;
        col(pos) = k;
        val(pos) = v(k-3*(j-1));                
    end

    else

% add entries for final row - todo!!

    end
end

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

% only keep the sparse indexing that we've used
row = row(1:pos);
col = col(1:pos);
val = val(1:pos);

Окончательную матрицу теперь можно построить с помощью команды sparse.

% build the actual sparse matrix
B = sparse(row,col,val,rows,cols);

Код можно запустить, сопоставив приведенные выше фрагменты. Для n = 9 мы получаем следующие результаты (для сравнения я также включил результаты для подхода на основе bsxfun, предложенного natan):

Elapsed time is 2.770617 seconds. (original)
Elapsed time is 0.006305 seconds. (sparse indexing)
Elapsed time is 0.261078 seconds. (bsxfun)

Исходному коду не хватает памяти для n = 10, хотя оба разреженных подхода все еще можно использовать:

Elapsed time is 0.019846 seconds. (sparse indexing)
Elapsed time is 2.133946 seconds. (bsxfun)

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

>> v=[-0,117377016134830 0,54433105395181 -0,0187057473531300 ...
-0,699119564792890 -0,136082763487960 0,426954037816980];
>> lendiff=length(v)-3;
>> B=repmat([v zeros(1,3^n-lendiff)],3^(n-1),1);
>> B=изменить форму(B',3^(n),3^(n-1)+1);
>> B(:,end-1)=B(:,end-1)+ B(:,end);
>> B=B(:,1:end-1)';

Здесь lendiff используется для создания 3^{n-1} копий строки с v, за которыми следуют нули, которые имеют длину 3^n+3, поэтому матрица размера [3^{n-1} 3^n+ 3].

Эта матрица преобразуется в размер [3^n 3^{n-1}+1] для создания сдвигов. Дополнительный столбец необходимо добавить к последнему, а B необходимо транспонировать.

Хотя должно быть намного быстрее.

ИЗМЕНИТЬ

Увидев решение Даррена и поняв, что reshape работает и с разреженными матрицами, я пришел к этому — без циклов for (не закодировал исходное решение).

Сначала значения для начала:

>> v=[-0.117377016134830  ...
       0.54433105395181   ...
      -0.0187057473531300 ...
      -0.699119564792890  ...
      -0.136082763487960  ...
       0.426954037816980];    
>> rows = 3^(n-1);                  % same number of rows
>> cols = 3^(n)+3;                  % add 3 cols to implement the shifts    

Затем сделайте матрицу с 3 дополнительными столбцами в строке

>> row=(1:rows)'*ones(1,length(v)); % row number where each copy of v is stored'
>> col=ones(rows,1)*(1:length(v));  % place v at the start columns of each row
>> val=ones(rows,1)*v;              % fill in the values of v at those positions
>> B=sparse(row,col,val,rows,cols); % make the matrix B[rows cols+3], but now sparse

Затем измените форму, чтобы реализовать сдвиги (дополнительная строка, правильное количество столбцов).

>> B=reshape(B',3^(n),rows+1);      % reshape into B[3^n rows+1], shifted v per row'
>> B(1:3,end-1)=B(1:3,end);         % the extra column contains last 3 values of v
>> B=B(:,1:end-1)';                 % delete extra column after copying, transpose

Для n=4,5,6,7 это приводит к процессорному времени в с:

n    original    new version
4    0.033       0.000
5    0.206       0.000
6    1.906       0.000
7    16.311      0.000

измеряется профилировщиком. Для исходной версии я не могу запустить n>7, но новая версия дает

n    new version
8    0.002
9    0.009
10   0.022
11   0.062
12   0.187
13   0.540
14   1.529
15   4.210

и вот как далеко идет моя оперативная память :)