Как я вектурирую этот код?

Я написал рекурсивную функцию, однако это занимает много времени. Поэтому он векторизовал его, но он не дает того же результата, что и рекурсивная функция. Это мой не-векторизованный код:

function visited = procedure_explore( u, adj_mat, visited )
visited(u) = 1;
neighbours = find(adj_mat(u,:));
for ii = 1:length(neighbours)
    if (visited(neighbours(ii)) == 0)
        visited = procedure_explore( neighbours(ii), adj_mat, visited );
    end
end
end

Это мой векторизованный код:

function visited = procedure_explore_vec( u, adj_mat, visited )
visited(u) = 1;
neighbours = find(adj_mat(u,:));
len_neighbours=length(neighbours);
visited_neighbours_zero=visited(neighbours(1:len_neighbours)) == 0;
if(~isempty(visited_neighbours_zero))
    visited = procedure_explore_vec( neighbours(visited_neighbours_zero), adj_mat, visited );
end
end

Это тестовый код

function main
    adj_mat=[0 0 0 0;
             1 0 1 1;
             1 0 0 0;
             1 0 0 1];
    u=2;
    visited=zeros(size(adj_mat,1));
    tic
    visited = procedure_explore( u, adj_mat, visited )
    toc
    visited=zeros(size(adj_mat,1));
    tic
    visited = procedure_explore_vec( u, adj_mat, visited )
    toc
end

Это алгоритм, который я пытаюсь реализовать:

Если векторизация невозможна, также будет полезно решение mex.

Тест производительности: Этот тест основан на MATLAB 2017a. Это показывает, что исходный код быстрее других методов.

Speed up between original and logical methods is 0.39672
Speed up between original and nearest methods is 0.0042583

Полный код

function main_recersive
    adj_mat=[0 0 0 0;
             1 0 1 1;
             1 0 0 0;
             1 0 0 1];
    u=2;
    visited=zeros(size(adj_mat,1));
    [email protected]()(procedure_explore( u, adj_mat, visited ));
    t_original=timeit(f_original);

    [email protected]()(procedure_explore_logical( u, adj_mat ));
    t_logical=timeit(f_logical);

    [email protected]()(procedure_explore_nearest( u, adj_mat,visited ));
    t_nearest=timeit(f_nearest);

    disp(['Speed up between original and logical methods is ',num2str(t_original/t_logical)])
    disp(['Speed up between original and nearest methods is ',num2str(t_original/t_nearest)])    

end

function visited = procedure_explore( u, adj_mat, visited )
    visited(u) = 1;
    neighbours = find(adj_mat(u,:));
    for ii = 1:length(neighbours)
        if (visited(neighbours(ii)) == 0)
            visited = procedure_explore( neighbours(ii), adj_mat, visited );
        end
    end
end

function visited = procedure_explore_nearest( u, adj_mat, visited )
    % add u since your function also includes it.
    nodeIDs = [nearest(digraph(adj_mat),u,inf) ; u];
    % transform to output format of your function
    visited = zeros(size(adj_mat,1));
    visited(nodeIDs) = 1;

end 

function visited = procedure_explore_logical( u, adj_mat )
   visited = false(1, size(adj_mat, 1));
   visited(u) = true;
   new_visited = visited;
   while any(new_visited)
      visited = any([visited; new_visited], 1);
      new_visited = any(adj_mat(new_visited, :), 1);
      new_visited = and(new_visited, ~visited);
   end
end

Ответ 1

Здесь интересная маленькая функция, которая выполняет нерекурсивный поиск по ширине на графе.

function visited = procedure_explore_logical( u, adj_mat )
   visited = false(1, size(adj_mat, 1));
   visited(u) = true;
   new_visited = visited;

   while any(new_visited)
      visited = any([visited; new_visited], 1);
      new_visited = any(adj_mat(new_visited, :), 1);
      new_visited = and(new_visited, ~visited);
   end
end

В Octave это работает примерно в 50 раз быстрее, чем ваша рекурсивная версия на матрице смежности 100x100. Вам нужно будет проверить его на MATLAB, чтобы узнать, что вы получаете.

Ответ 2

Вы можете представить свою матрицу смежности как список путей длины ровно один. Вы можете генерировать пути других длин n, перенося их на n-ю степень до ранга вашей матрицы. (adj_mat ^ 0 - единичная матрица)

В графе с n узлами самый длинный путь не может быть больше n-1, поэтому вы можете суммировать все мощности для анализа достижимости:

adj_mat + adj_mat^2 + adj_mat^3
ans =
   0   0   0   0
   4   0   1   3
   1   0   0   0
   3   0   0   3

Это количество (разных) способов, которые вы можете использовать для перехода от одного node к другому. Для простой доступности, проверьте, больше ли это значение:

visited(v) = ans(v, :) > 0;

В зависимости от вашего определения вам может потребоваться изменить столбцы и строки в результате (т.е. взять ans (:, v)).

Для производительности вы можете использовать более низкие мощности для создания более высоких. Например, было бы эффективно рассчитано что-то вроде A + A ^ 2 + A ^ 3 + A ^ 4 + A ^ 5:

A2 = A^2;
A3 = A2*A
A4 = A2^2;
A5 = A4*A;
allWalks= A + A2 + A3 + A4 + A5;

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

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

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

См. также этот ответ: fooobar.com/questions/459043/...

Ответ 3

Я не думаю, что вы можете правильно обозначить свою функцию: ваша оригинальная функция никогда не достигает одного и того же node несколько раз. С помощью векторизации вы передадите все непосредственно связанные узлы одновременно с следующей функцией. Поэтому тогда было бы возможно, что в следующих случаях один и тот же node получает несколько раз. Например. в вашем примере node 1 будет достигнуто 3 раза. Таким образом, в то время как у вас больше не будет цикла, функция может в зависимости от вашей сети быть вызвана более рекурсивно, что увеличит время вычислений.

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

Немного offtopic, но с Matlab 2016a вы можете использовать nearest(), чтобы найти все доступные узлы (без запуска node)., Он вызывает алгоритм ширины в отличие от вашего алгоритма глубины:

% add u since your function also includes it.
nodeIDs = [nearest(digraph(adj_mat),u,inf) ; u]; 

% transform to output format of your function
visited = zeros(size(adj_mat,1));
visited(nodeIDs) = 1;

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

Ответ 4

Проблема с рекурсивной функцией связана с visited(u) = 1;. Это связано с тем, что MATLAB использует метод copy-on-write для передачи/назначения переменных. Если вы не изменяете visited в теле функции, ее копия не производится, но когда она изменяется, копия ее создается и модификация применяется на ее копии. Чтобы предотвратить использование объекта handle object, необходимо передать ссылку на функцию.

Определите класс дескриптора (сохраните его на visited_class.m):

classdef visited_class < handle
    properties
        visited
    end
    methods
        function obj = visited_class(adj_mat)
            obj.visited = zeros(1, size(adj_mat,1));
        end
    end
end

Рекурсивная функция:

function procedure_explore_handle( u, adj_mat,visited_handle )
    visited_handle.visited(u) = 1;
    neighbours = find(adj_mat(u,:));
    for n = neighbours
        if (visited_handle.visited(n) == 0)
            procedure_explore_handle( n, adj_mat , visited_handle );
        end
    end
end

Инициализация переменных:

adj_mat=[0 0 0 0;
         1 0 1 1;
         1 0 0 0;
         1 0 0 1];
visited_handle = visited_class(adj_mat);
u = 2;

Назовите его как:

procedure_explore_handle( u, adj_mat,visited_handle );

Результат сохраняется в visited_handle:

disp(visited_handle.visited)

Ответ 5

Если вы хотите перейти от одной точки графика к другому, наиболее эффективным способом найти ее с точки зрения ресурсов является алгоритм Дейкстра. Алгоритм Флойда-Варшалла вычисляет все расстояния между всеми точками и может быть параллелизирован (начиная с нескольких точек).

Зачем нужна векторизация (или использование mex-программирования)? Если вы просто хотите максимально использовать быстрые алгоритмы умножения Matlab, то использование продуктов A должно быстро найти вас:

adj_mat2=adj_mat^2;               % allowed to use 2 steps
while (adj_mat2 ~= adj_mat)       % check if new points were reached
      adj_mat=adj_mat2;           % current set of reachable points
      adj_mat2=(adj_mat^2)>0;     % allowed all steps again: power method
end

Ответ 6

Этот ответ просто дает явную, векторизованную реализацию предложения DasKrümelmonster, который, я думаю, быстрее, чем код в вопросе (по крайней мере, если размеры матрицы не слишком велики). Он использует функцию polyvalm для оценки суммы степеней матрицы смежности.

function visited = procedure_explore_vec(u, adj_mat)
    connectivity_matrix = polyvalm(ones(size(adj_mat,1),1),adj_mat)>0;

    visited = connectivity_matrix(u,:);
end