Вот мы и добрались до написания программ.
В этой статье напишем скрипты для расчётов резистивного расстояния и для моделирования случайных блужданий. В качестве ЯП был выбран 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.

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

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