Случайные блуждания: связь с резистивным расстоянием (часть 3)

Вот мы и добрались до написания программ.
В этой статье напишем скрипты для расчётов резистивного расстояния и для моделирования случайных блужданий. В качестве ЯП был выбран Octave (всё-таки математикой занимаемся).

И так, в прошлых двух частях мы познакомились с теорией данной темы. Настало время проверить теорию на практике.

Резистивное расстояние

Определим задачи для нашей программы:

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

Ввод и предварительная проверка данных

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

Здесь проверяется количество переданных параметров. В случае если их недостаточно, программа завершает работу с сообщением об ошибке. Далее происходит парсинг параметров: первый — имя файла, второй — начальный узел, третий — конечный узел. После происходит проверка введенных параметров на корректность. В случае неправильных данных вызывается ошибка с пояснением. Запускается таймер выполнения программы. В случае, если переданный файл уже обрабатывался, данные загрузятся из сохранённого файла.

args = argv(); if numel(args) < 3 error('Необходимо указать имя файла с матрицей смежности и',  'номера двух узлов.'); end filename = args{1};  a = str2double(args{2}); b = str2double(args{3}); if isnan(a) || isnan(b) || a <= 0 || b <= 0  error('Номера узлов должны быть положительными целыми ', 'числами.'); endif  tic; [~, name, ~] = fileparts(filename); pseudo_inverse_filename = [name, '_L_plus.mat'];  if exist(pseudo_inverse_filename, 'file') == 2     load(pseudo_inverse_filename, 'L_plus');     disp(['Псевдообратная матрица Лапласа загружена из файла ',  pseudo_inverse_filename, '.']); n = size(L_plus, 1);

Основные расчеты

Этот блок отвечает за тяжеловесные расчеты программы.

Считывается матрица смежности из указанного файла. Подсчитывается количество вершин. Вычисляются степени вершин, и создается диагональная матрица степеней. Затем вычисляется матрица Лапласиана. Происходит проверка корректности вычисления псевдообратной матрицы. Вычисляется псевдообратная матрица Лапласиана. Результат вычислений записывается в специальный файл. Считывание псевдообратной матрицы из этого файла позволит сократить время расчёта при очередном запуске программы с той же матрицей смежности.

else A = dlmread(filename); n = size(A, 1);   D = diag(sum(A, 2)); L = D - A; if rank(L) < n-1 error('Матрица Лапласиана не имеет полного ранга.'); end L_plus = pinv(L); save(pseudo_inverse_filename, 'L_plus'); disp(['Псевдообратная матрица Лапласа вычислена и ', 'сохранена в файл ', pseudo_inverse_filename, '.']); end

Расчет разницы потенциалов

Проверяется, что указанные узлы существуют в матрице. Создается результирующий вектор и производится расчет резистивного расстояния (эффективного сопротивления). Для оптимизации программы эти действия записаны в одну строчку.

if a > n || b > n error('Номера узлов выходят за пределы размерности ', 'матрицы.'); endif  answer = (L_plus(a, a) - L_plus(b, a)) - (L_plus(a, b) -  L_plus(b, b));

Вывод результата

Остановка таймера выполнения и вывод результата:

elapsed_time = toc; disp("==========================================================="); disp(['Резистивное расстояние между узлами ', num2str(a), ' и ',  num2str(b), ': ', num2str(answer)]); disp(['Время выполнения программы: ', num2str(elapsed_time),  ' секунд']);

Тестирование вычисления резистивного расстояния

Программа была протестирована на различных графах в три этапа:

Проверка правильности работы алгоритма на малых графах (до 5 вершин, заполненность 100%).
Тестирование на графах средних размеров (100-600 вершин, заполненность 90%).
Тестирование на больших графах (1000 и более вершин, заполненность 25-90%).

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

Результаты тестирования

Малые графы: алгоритм выдает решение менее чем за сотую долю секунды.
Средние графы: время выполнения менее полусекунды для графов с количеством вершин менее 500. При 500-600 вершинах расчет занимает 3.5 секунды.
Большие графы: для графа с 1000 вершинами и заполненностью 90% расчет занимает 6.5 секунд. Граф с 3000 вершинами и заполненностью 25% рассчитывается за 500 секунд.

Зависимость времени выполнения программы от количества вершин в графе

Основное время выполнения программы приходится на вычисление псевдообратной матрицы (сложность O(n^3) ). Поэтому, время выполнения увеличивается на три порядка при увеличении размера графа в 10 раз.

Реализация программы случайного блуждания

Аналогично использовался Octave. Но т.к. в этом скрипте достаточно большое количество переменных, думаю не лишним будет дать описание каждой.

Входные переменные

$\text{ adj}$ — матрица смежности графа.
$\text{start}$ — индекс начального узла.
$\text{end_}$ — индекс конечного узла.
$\text{num_sim}$ — количество тестов.

Выходные переменные

$\text{fht}$ — вектор, содержащий время первого попадания в конечный узел для каждого теста.
$\text{ct}$ — вектор, содержащий время полного обхода графа для каждого теста.
$\text{cmt}$ — вектор, содержащий время перемещения от начального узла к конечному и обратно для каждого теста.
$\text{sfh}$ — отсортированный вектор времени первого попадания и их количества.
$\text{sct}$ — отсортированный вектор времени полного обхода графа и их количества.
$\text{scmt}$ — отсортированный вектор времени перемещения и их количества.
$\text{mfht}$ — среднее время первого попадания в конечный узел.
$\text{mct}$ — среднее время полного обхода графа.
$\text{mcmt}$ — среднее время перемещения от начального узла к конечному и обратно.
$\text{eff_res}$ — эффективное сопротивление графа, рассчитанное как среднее время перемещения, деленное на удвоенное количество ребер.

Промежуточные переменные

$\text{n}$ — количество узлов в графе.
$\text{m}$ — количество ребер в графе, рассчитанное как половина суммы всех элементов матрицы смежности.
$\text{curr}$ — текущий узел в процессе блуждания.
$\text{visited}$ — булевый вектор, указывающий, были ли посещены узлы.
$\text{steps}$ — количество шагов, сделанных в текущем тесте.
$\text{first_hit}$ — булева переменная, указывающая, было ли достигнуто первое попадание в конечный узел.
$\text{visit_count}$ — количество посещенных узлов.
$\text{neighbors}$ — вектор соседних узлов для текущего узла.
$\text{next}$ — следующий узел для перехода.
$\text{cms}$ — количество шагов для перемещения от начального узла к конечному и обратно в текущем тесте.
$\text{file_path}$ — путь к файлу с матрицей смежности.
$\text{elapsed_time}$ — время, затраченное на выполнение тестов.

Скрипт случайных блужданий

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

function [fht, ct, mfht, mct, eff_res, mcmt] = random_walk(adj, start, end_, num_sim) n = size(adj, 1); m = sum(adj(:)) / 2; fht = zeros(num_sim, 1); ct = zeros(num_sim, 1); cmt = zeros(num_sim, 1);

Следующией код выполняет цикл по числу тестов для выполнения случайных блужданий.

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

for sim = 1:num_sim   curr = start;   visited = false(n, 1);   visited(start) = true;   steps = 0;   first_hit = false;   visit_count = 1;

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

while visit_count < n   neighbors = find(adj(curr, :));   next = neighbors(randi(length(neighbors)));   steps = steps + 1;   curr = next;

Если это первое попадание в конечный узел, сохраняется количество шагов до первого попадания.

if ~first_hit && curr == end_   fht(sim) = steps;   first_hit = true; end

Если узел еще не был посещен, он помечается как посещенный, и увеличивается счетчик посещений.

  if ~visited(curr)     visited(curr) = true;     visit_count = visit_count + 1;   end end

После завершения обхода всех узлов записывается количество шагов.

ct(sim) = steps;

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

  cms = 0;   curr = start;   while curr ~= end_     neighbors = find(adj(curr, :));     next = neighbors(randi(length(neighbors)));     cms = cms + 1;     curr = next;   end   while curr ~= start     neighbors = find(adj(curr, :));     next = neighbors(randi(length(neighbors)));     cms = cms + 1;     curr = next;   end   cmt(sim) = cms; end

Здесь рассчитывается среднее время первого попадания, среднее время обхода, среднее время перемещения и эффективное сопротивление.

  mfht = mean(fht);   mct = mean(ct);   mcmt = mean(cmt);    eff_res = mcmt / (2 * m); end

Для запуска функции random walk и получения результатов был разработан скрипт run random walk. Этот скрипт задаёт параметры графа, запускает функцию random walk и выводит результаты на экран.

Указание пути к файлу с матрицей смежности, начального узла, конечного узла и числа тестов. Затем чтение матрицы смежности из указанного файла с помощью функции $\texttt{dlmread}$ и сохранение в переменную матрицы смежности.

args = argv(); if numel(args) < 3         error('Необходимо указать имя файла с матрицей смежности и',                 'номера двух узлов.'); end  file_path = args{1}; start = str2double(args{2}); end_ = str2double(args{3}); num_sim = 1000;  adj = dlmread(file_path);

Выполнение функции случайного блуждания с заданными параметрами и измерение затраченного времени с помощью $\texttt{toc}$ и $\texttt{toc}$ . Результаты сохраняются в соответствующих переменных, а время выполнения — в переменной $\texttt{elapsed_time}.$

tic; [fht, ct, mfht, mct, eff_res, mcmt] = random_walk(adj, start, end_, num_sim); elapsed_time = toc;

Этот блок кода выводит основные результаты вычислений, такие как: среднее время первого попадания, среднее время прохода, среднее время обхода всего графа, эффективное сопротивление, время выполнения программы

fprintf('Среднее время первого попадания в вершину %d из' , 'вершины %d: %f шага.\n', end_, start, mfht); fprintf('Среднее время прохода из вершины %d в вершину %d и' ,  'обратно: %f шага.\n', start, end_, mcmt); fprintf('Среднее время обхода всего графа: %f шага.\n', mct); fprintf('Эффективное сопротивление: %f.\n', eff_res); fprintf('Время выполнения программы: %f секунд.\n', elapsed_time);

Проведение тестов

В качестве примеров мы рассмотрим три графа. Первый граф будет тестовым для понимания методологии решения, вторым графом является Московское метро, третий — линейный граф.

Пример 1

$\text{adj_matrix} = \begin{pmatrix} 0 & 1 & 1 & 0 \\ 1 & 0 & 1 & 1 \\ 1 & 1 & 0 & 1 \\ 0 & 1 & 1 & 0 \end{pmatrix}$

Начальная вершина была выбрана как вершина 1, а конечная вершина как вершина 4. Количество тестов было установлено на 1000. Результаты тестов показали следующее:

Среднее время первого попадания в вершину 4 из вершины 1: 5.047200 шага.
Среднее время прохода из вершины 1 в вершину 4 и обратно: 10.004400 шага.
Среднее время обхода всего графа: 5.934400 шага.
Сопротивление: 1.000440.
Время выполнения программы: 7.940940 секунд.

Стоит заметить, что при стократном увеличении тестов, результаты изменятся незначительно.

Пример 2

В качестве второго примера была взята матрица Московского метро (442 вершины, 630 ребер), так как она позволяет проверить случайное блуждание в реальной транспортной сети. Начальная вершина была выбрана как вершина 1 (станция Новокосино), а конечная вершина как вершина 442 (станция Апрелевка). Количество тестов было установлено на 100000. Результаты тестов показали следующее:

Среднее время первого попадания в вершину 442 из вершины 1: 17159.75 шага.
Среднее время прохода из вершины 1 в вершину 442 и обратно: 20142.68 шага.
Среднее время обхода всего графа: 40777.67 шага.
Сопротивление: 15.99.
Время выполнения программы: 3194.03 секунд.

Резистивное расстояние по результатам выполнения программы для расчёта сопротивления с помощью законов Киргхофа равно 16.
Тогда из формулы (\ref{C_xy}) следует, что среднее время коммутирования равно:

$C^t_{(1,442)} = 2mR_{1,442} = 2 * 630 * 16 = 20160$

Среднее время коммутирования по результатам случайных блужданий:

$C^r_{(1,442)} \approx 20142$

Результаты вычислений показывают, что среднее время прохождения из вершины 1 в вершину 442 и обратно почти точно удовлетворяет соотношению формулы (\ref{C_xy}).

$C^t_{(1,442)} \approx C^r_{(1,442)}$

Ошибка составлятет 0,01%.

Пример 3

Третьим примером стал линейный граф с количеством вершин 442 и количеством ребер 441. Количество тестов было установлено на 10000.

Среднее время первого попадания в вершину 442 из вершины 1: 191980.58 шага.
Среднее время прохода из вершины 1 в вершину 442 и обратно: 389058.45 шага.
Среднее время обхода всего графа: 191980.58.
Сопротивление: 441.11.
Время выполнения программы: 32729.89 секунд.

Резистивное расстояние по результатам выполнения программы для расчёта сопротивления с помощью законов Киргхофа равно 441.
Тогда из формулы (\ref{C_xy}) следует, что среднее время коммутирования равно:

$C^t_{(1,442)} = 2mR_{1,442} = 2 * 441 * 441 = 388962$

Среднее время коммутирования по результатам случайных блужданий:

$C^r_{(1,442)} \approx 389058.45$

Результаты вычислений показывают, что среднее время прохождения из вершины 1 в вершину 442 и обратно достаточно точно удовлетворяет соотношению формулы (\ref{C_xy}).

$C^t_{(1,442)} \approx C^r_{(1,442)}$

Ошибка составлятет 0,03%.

Теоретические расчеты

Пример 1

Теперь с помощью теории рассчитаем количество шагов до попадания в вершину 4 из вершины 1. Составим матрицу степеней из матрицы смежности, используемой выше:

$\Delta = \begin{pmatrix} 2 & 0 & 0 & 0 \\ 0 & 3 & 0 & 0 \\ 0 & 0 & 3 & 0 \\ 0 & 0 & 0 & 2 \end{pmatrix}$

Из этих двух матриц составим матрицу Лапласа (Лапласиан), используя данную формулу: ${L = \Delta - A}$

$\begin{pmatrix} 2 & 0 & 0 & 0 \\ 0 & 3 & 0 & 0 \\ 0 & 0 & 3 & 0 \\ 0 & 0 & 0 & 2 \end{pmatrix} - \begin{pmatrix} 0 & 1 & 1 & 0 \\ 1 & 0 & 1 & 1 \\ 1 & 1 & 0 & 1 \\ 0 & 1 & 1 & 0 \end{pmatrix} = \begin{pmatrix} 2 & -1 & -1 & 0 \\ -1 & 3 & -1 & -1 \\ -1 & -1 & 3 & -1 \\ 0 & -1 & -1 & 2 \end{pmatrix} = \text{L}$

Далее мы получаем собственные значения и собственные вектора матрицы Лапласа, чтобы их получить я воспользуюсь Octave:

laplacian_matrix = [2 -1 -1 0;                                    -1 3 -1 -1;                    -1 -1 3 -1;                     0 -1 -1 2];   [eigenvectors, eigenvalues_matrix] = eig(laplacian_matrix);  eigenvalues = diag(eigenvalues_matrix);  disp('Собственные значения:'); disp(eigenvalues);  disp('Собственные вектора:'); disp(eigenvectors);

Где eigenvalues – это $\lambda_k$ , а eigenvectors – это $u_{i}$

$\lambda_1 = 0, \quad \lambda_2 = 2, \quad \lambda_3 = 4, \quad \lambda_4 = 4$ $\text{U} = \begin{pmatrix} -0.50 & -0.71 & 0.49 & 0.09 \\ -0.50 & 0.00 & -0.62 & 0.60 \\ -0.50 & 0.00 & -0.36 & -0.79 \\ -0.50 & 0.71 & 0.49 & 0.09 \\ \end{pmatrix}$

У собственных значений убираем нулевое значения и берем нулевую и третью строки из матрицы $\text{U}$ для расчета $R_{14}$ . Получаем следующее:

$\lambda_2 = 2, \quad \lambda_3 = 4, \quad \lambda_4 = 4$

Теперь нам нужны значения и . У нас уже известна, она равна 4, поэтому осталось найти . Для этого находим сумму всех элементов матрицы смежности и делим их на 2.

$\sum_{i,j} A_{ij} = 0 + 1 + 1 + 0 + 1 + 0 + 1 + 1 + 1 + 1 + 0 + 1 + 0 + 1 + 1 + 0 = 10$ $E = \frac{10}{2} = 5$

В итоге получаем, что = 5 и = 4. Получив эти значения, мы можем посчитать сопротивление

$R_{14} = \sum_{k=1}^{n-1} \frac{(u_{k1} - u_{k4})^2}{\lambda_k} = \frac{(0.71 - (-0.71))^2}{2} +$ $+ \frac{(0.49 - 0.49)^2}{4} + \frac{(0.09 - 0.09)^2}{4} \approx 1.0082$

Теперь нужно посчитать commute time, для этого воспользуемся формулой из предыдущей статьи (см часть 2):

$C_{14} = 2mR_{14} = 2 * 5 * 1.0083 \approx 10.0082$

Нам осталось посчитать только hitting time, но так как наш граф является маленьким и симметричным, то у нас есть возможность воспользоваться данным вариантом формулы:

$C_{14} = H_{14} + H_{41}$

Для нашего маленького и симметричного графа, это равносильно что:

$C_{14} \approx 2H_{14}$

Следовательно $H_{14}$ равен:

$H_{14} \approx \frac{C_{14} }{2} \approx \frac{10.0082}{2} \approx 5.0041$

Пример 2

Пример 1 был приведен для демонстрации методологии. Поскольку нахождение вручную данных графа Московского метро практически невозможно, воспользуемся программным методом. Для этого была написана программа, которая воспроизводит вычисления формул: $R_{xy}$ , $H_{xy}$ , $C_{xy}$ .

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

args = argv(); if numel(args) < 3         error('Необходимо указать имя файла с матрицей смежности и',                  'номера двух узлов.'); end file_path = args{1}; adj_matrix = dlmread(file_path);  n = size(adj_matrix, 1); m = sum(adj_matrix(:)) / 2;

Здесь мы вычисляем по формулам собственные значения и собственные векторы. Далее извлекаем собственные значения в вектор. И отбрасываем первое нулевое собственное значение и соответствующий ему собственный вектор (они не нужны для дальнейших вычислений).

degree_matrix = sum(adj_matrix, 2);  laplacian_matrix = diag(degree_matrix) - adj_matrix;  [eigenvectors, eigenvalues_matrix] = eig(laplacian_matrix); eigenvalues = diag(eigenvalues_matrix);  eigenvalues_nonzero = eigenvalues(2:end); eigenvectors_nonzero = eigenvectors(:, 2:end);

В этой части код инициализируем матрицу H нулями. Задаем максимальное количество итераций и допустимую ошибку.

H = zeros(n, n);  max_iter = 200000; tol = 1e-16;

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

tic;  neighbors = cell(n, 1); for i = 1:n     neighbors{i} = find(adj_matrix(i, :) == 1); end

В этом блоке в каждой итерации обновляем значения H для всех вершин, используя их соседей. Затем каждые 10000 итераций выводим текущую ошибку. Если ошибка становится меньше заданного порога (tol), процесс останавливается. И в конце записываем время выполнения итерационного процесса.

for iteration = 1:max_iter     H_prev = H;     for u = 1:n         if degree_matrix(u) > 0             H(u, :) = 1 + (1 / degree_matrix(u)) * sum(H(neighbors{u}, :), 1);         end         H(u, u) = 0;  % Условие H_ii = 0     end     if mod(iteration, 10000) == 0         disp(['Итерация: ', num2str(iteration), ' H - ', num2str(max(max(abs(H - H_prev))))]);     end     if max(max(abs(H - H_prev))) < tol         disp(['Сошлось после ', num2str(iteration), ' итераций']);         break;     end end elapsed_time = toc;

Здесь мы рассчитаем среднее время прохождения, время возвращения и сопротивление между заданными вершинами на основе полученных данных.

i = str2double(args{2}); j = str2double(args{3});  H_ij = H(i, j);  R_ij = sum((eigenvectors_nonzero(i, :) - eigenvectors_nonzero(j, :)).^2 ./ eigenvalues_nonzero');  C_ij = 2 * m * R_ij;

Выводим вычисленные результаты.

disp(['Среднее время прохода из вершины ', num2str(i), ' в вершину ', num2str(j), ': ', num2str(H_ij)]); disp(['Среднее время прохода из вершины ', num2str(i), ' в вершину ', num2str(j), ' и обратно: ', num2str(C_ij)]); disp(['Сопротивление: ', num2str(R_ij)]); disp(['Время выполнения программы: ', num2str(elapsed_time), ' секунд']);

Результаты данных вычислений:

Сошлось после 293492 итераций
Среднее время прохождения из вершин 1 в вершину 442: 17085.752085552034
Среднее время прохождения из вершины 1 в вершину 442 и обратно: 20199.377213668093
Сопротивление: 16.03125175687944
Время выполнения программы: 2841.8106 секунд

Пример 3

Результаты вычислений для линейного графа:

Сошлось после 2016601 итераций
Среднее время прохождения из вершин 1 в вершину 442: 194480.99999830942
Среднее время прохождения из вершины 1 в вершину 442 и обратно: 388961.9999977657
Сопротивление: 440.9999999974668
Время выполнения программы: 8452.6199 секунд

Заключение

Результаты вычислительных экспериментов и теоретических расчетов для графа демонстрируют высокую степень согласованности, что свидетельствует о правильности применяемых моделей и алгоритмов. Тестирование проводилось на графе с четырьмя вершинами, на котором было установлено, что теоретические предсказания совпадают с практическими результатами. Кроме того, применение данных методов на лиенйном графе и на графе Московского метро показало, что алгоритмы могут работать не только в теоретических условиях, но и с реальной транспортной сетью. Среднее время первого попадания, прохождения и обхода всего графа, а также сопротивление, полученные в результате тестов, находятся в хорошем согласии с теоретическими расчетами. Эти три случая подтверждают, что предложенные модели и алгоритмы могут быть успешно применены для анализа и прогнозирования поведения более сложных и масштабных транспортных сетей.

Все скрипты доступны на GitHub.

Вот и написаны программы. Теперь можно применить статистический аппарат для анализа полученных данных. Но эта уже тема следующей главы.

ссылка на оригинал статьи https://habr.com/ru/articles/831490/

Случайные блуждания: связь с резистивным расстоянием (часть 3)

Резистивное расстояние

Ввод и предварительная проверка данных

Основные расчеты

Расчет разницы потенциалов

Вывод результата

Тестирование вычисления резистивного расстояния

Результаты тестирования

Реализация программы случайного блуждания

Входные переменные

Выходные переменные

Промежуточные переменные

Скрипт случайных блужданий

Проведение тестов

Пример 1

Пример 2

Пример 3

Теоретические расчеты

Пример 1

Пример 2

Пример 3

Заключение

Комментарии

Добавить комментарий Отменить ответ