Вычислить значение двойного суммирования функции по парам вершин графа

Ответ 1

Ожидаемый ответ R_inf[1,2] = 0.8 и R_inf[6,8] = 0.8 примера OP a является результатом алгоритма, применяемого к ориентированному графу a.T, а не неориентированному.

Здесь R_2 вычислено с кодом OP (с некоторыми незначительными исправлениями) и с моим кодом. И R_inf также вычисляется с помощью моего кода:

after 2 iterations with pp_loop:
0.8000 0      0      0      0      0      0      0.8000 0     
0      0.8000 0.8000 0      0      0.4000 0.3200 0      0.3200
0      0.8000 0.8000 0      0      0.4000 0.3200 0      0.3200
0      0      0      0.8000 0.4800 0      0      0      0     
0      0      0      0.4800 0.8000 0      0      0      0     
0      0.4000 0.4000 0      0      0.8000 0.1600 0      0.1600
0      0.3200 0.3200 0      0      0.1600 0.8000 0      0.8000
0.8000 0      0      0      0      0      0      0.8000 0     
0      0.3200 0.3200 0      0      0.1600 0.8000 0      0.8000
after 2 iterations with OP:
0.8000 0      0      0      0      0      0      0.8000 0     
0      0.8000 0.8000 0      0      0.4000 0.3200 0      0.3200
0      0.8000 0.8000 0      0      0.4000 0.3200 0      0.3200
0      0      0      0.8000 0.4800 0      0      0      0     
0      0      0      0.4800 0.8000 0      0      0      0     
0      0.4000 0.4000 0      0      0.8000 0.1600 0      0.1600
0      0.3200 0.3200 0      0      0.1600 0.8000 0      0.8000
0.8000 0      0      0      0      0      0      0.8000 0     
0      0.3200 0.3200 0      0      0.1600 0.8000 0      0.8000
in fixed point:
0.8000 0      0      0      0      0      0      0.8000 0     
0      0.8000 0.8000 0      0      0.4000 0.3200 0      0.3200
0      0.8000 0.8000 0      0      0.4000 0.3200 0      0.3200
0      0      0      0.8000 0.4800 0.0960 0.0256 0      0.0256
0      0      0      0.4800 0.8000 0.1280 0      0      0     
0      0.4000 0.4000 0.0960 0.1280 0.8000 0.1600 0      0.1600
0      0.3200 0.3200 0.0256 0      0.1600 0.8000 0      0.8000
0.8000 0      0      0      0      0      0      0.8000 0     
0      0.3200 0.3200 0.0256 0      0.1600 0.8000 0      0.8000

Мое решение быстрее и может вычислять фиксированную точку для сетей до ~1000 узлов.

Немного объяснения: стратегия состоит в том, чтобы выразить отношение рекуррентности по матричной алгебре, потому что тогда мы можем использовать высоко оптимизированную линейную алгебру blas, которая поставляется с numpy/scipy. Например, внутренние суммы - все они за один раз - могут быть выполнены путем умножения матрицы R_k текущих значений на транспонированную матрицу смежности. Это может показаться расточительным, учитывая все нули в матрице смежности, но blas является очень-очень быстрым. Аналогично, мы можем вычислить размеры пересечений, умножив матрицу смежности на ее транспонирование и т.д.

Предполагая, что конечной целью является не столько последовательность R, но и фиксированная точка, мы можем сэкономить еще некоторое время, используя правильный линейный решатель, потому что рекуррентность является линейной (+ смещение). Поскольку запись n^2xn^2 матрицы этой операции неэффективна, мы используем итеративный решатель, который может использовать вместо LinearOperator матрицу.

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

Код, включающий OP с исправлениями:

import numpy as np
from scipy import sparse
from scipy.sparse import linalg

def mock_data(n, p):
    data = (np.random.random((n, n)) < p).astype(int)
    np.einsum('ii->i', data)[:] = 0
    return data

import copy

import scipy
import scipy.sparse as sp
import time
import warnings

class Algo():
    def __init__(self, c=0.8, max_iter=5, is_verbose=False, eps=1e-3):
        self.c = c
        self.max_iter = max_iter
        self.is_verbose = is_verbose
        self.eps = eps

    def compute_similarity(self, a, R0=None, use_undirected=False,
                           no_neighbor_policy='10'):
        '''

        Note: The adjacency matrix is a (0,1)-matrix with zeros on its 
        diagonal.
        '''                                 
        a_transpose=np.transpose(a)
        a_undirected=np.add(a,a_transpose)
        a_use = a_undirected if use_undirected else a

        self.num_nodes = np.shape(a)[0]
        current_sim = np.eye(self.num_nodes) if R0 is None else R0
        prev_sim = np.zeros(shape=(self.num_nodes, self.num_nodes))

        #Determine the set of Neighbors for each node
        nbs = []

        for i in range(self.num_nodes):

            nbs.append(np.nonzero(a_use[i])[0])
        R = []

        for itr in range(self.max_iter):

            if scipy.allclose(current_sim, prev_sim, atol=self.eps):
                break

            prev_sim = copy.deepcopy(current_sim)
            R.append(prev_sim)
            for i in range(self.num_nodes):
                for j in range(self.num_nodes):
                    if len(nbs[i]) == 0 and len(nbs[j]) == 0:
                        current_sim[i][j] = (
                            0 if (no_neighbor_policy=='00' or
                                  (no_neighbor_policy=='10' and i!=j))
                            else self.c)
                    elif i == j:
                        current_sim[i][j] = self.c
                    elif len(nbs[i]) == 0 or len(nbs[j]) == 0:
                        current_sim[i][j] = 0
                    else:
                        #commonSet=0.0
                        differenceSetofp_ij=0.0
                        differenceSetofq_ij=0.0
                        commonSet=len(np.intersect1d(nbs[i],nbs[j]))/np.maximum(len(np.union1d(nbs[i],nbs[j])), 1)


                        for x in np.setdiff1d(nbs[i],nbs[j]):
                            for y in nbs[j]:

                                differenceSetofp_ij+=prev_sim[x][y]
                        differenceSetofp_ij*=1/np.maximum(len(np.union1d(nbs[i],nbs[j]))*len(nbs[j]), 1)
                                #print differenceSetofp



                        for x in np.setdiff1d(nbs[j],nbs[i]):
                            for y in nbs[i]:

                                differenceSetofq_ij+=prev_sim[x][y]
                        differenceSetofq_ij*=1/np.maximum(len(np.union1d(nbs[j],nbs[i]))*len(nbs[i]), 1) 

                        current_sim[i][j] = self.c*(commonSet+differenceSetofp_ij+differenceSetofq_ij)


        if self.is_verbose:

            print('SimRank: converged after {0} iterations'.format(itr))
        R.append(current_sim)
        return np.array(R)

def f_OP(adj, c, max_n, R0=None, use_undirected=False, no_neighbor_policy='10'):
    return Algo(c, max_n).compute_similarity(adj, R0, use_undirected,
                                             no_neighbor_policy)

def f_pp(adj, c, max_n, R0=None, use_undirected=False,
         no_neighbor_policy='10', solver=linalg.gcrotmk):
    n, n = adj.shape
    if use_undirected:
         adj = adj | adj.T
    aind = np.where(adj.T)
    n_neigh = adj.sum(axis=1)
    n_inter = adj @ adj.T
    n_union = np.add.outer(n_neigh, n_neigh) - n_inter
    r_union = c / np.maximum(n_union, 1)
    if no_neighbor_policy == '11':
        n_inter[n_union==0] = 1
        r_union[n_union==0] = c
    elif no_neighbor_policy == '10':
        np.einsum('ii->i', n_inter)[np.einsum('ii->i', n_union)==0] = 1
        np.einsum('ii->i', r_union)[np.einsum('ii->i', n_union)==0] = c
    elif no_neighbor_policy != '00':
        raise ValueError
    avg_ngh = adj.T / np.maximum(n_neigh, 1)
    if solver is None:
        R = np.empty((max_n+1, n, n))
        R[0] = 0 if R0 is None else R0
        if R0 is None:
            np.einsum('ii->i', R[0])[:] = 1       
        for i in range(max_n):
            rowsums = R[i] @ avg_ngh
            rowsums[aind] = 0
            rec_sum = adj @ rowsums
            R[i+1] = r_union * (n_inter + rec_sum + rec_sum.T)
            if np.allclose(R[i+1], R[i], rtol=1e-7, atol=1e-10):
                return R[:i+2]
        return R
    else:
        def LO(x):
            R = x.reshape(n, n)
            rowsums = R @ avg_ngh
            rowsums[aind] = 0
            rec_sum = adj @ rowsums
            return (r_union * (rec_sum + rec_sum.T) - R).ravel()
        R0 = np.identity(n) if R0 == None else R0
        LO = linalg.LinearOperator((n*n, n*n), matvec=LO)#, rmatvec=LO)
        res = solver(LO, (-r_union * n_inter).ravel(), R0.ravel(), tol=1e-9)
        assert res[1]==0
        return res[0].reshape(n, n)

def f_pp_loop(adj, c, max_n, R0=None, use_undirected=False,
              no_neighbor_policy='10'):
    return f_pp(adj, c, max_n, R0, use_undirected, no_neighbor_policy, None)


# test on OP example:
a = np.array([ [0, 1, 1, 0, 0, 1, 0, 0, 0],
               [0, 0, 0, 0, 1, 0, 0, 0, 0],
               [0, 0, 0, 1, 0, 0, 0, 0, 0],
               [0, 0, 0, 0, 0, 0, 0, 0, 0],
               [0, 0, 0, 0, 0, 0, 1, 0, 1],
               [0, 0, 0, 1, 0, 0, 0, 0, 0],
               [0, 0, 0, 0, 0, 1, 0, 0, 0],
               [0, 0, 0, 0, 0, 0, 1, 0, 1],
               [0, 0, 0, 0, 0, 0, 0, 0, 0]])

a = a.T
m, c, nnp = 2, 0.8, '11'

sol = {}
for f in f_pp_loop, f_OP:
    name = f.__name__[2:]
    sol[name] = f(a, c, m, None, False, nnp)
    print(f"after {len(sol[name])-1} iterations with {name}:")
    print("\n".join(" ".join("0     " if i == 0 else f"{i:6.4f}" for i in row) for row in sol[name][-1]))
print("in fixed point:")
print("\n".join(" ".join("0     " if i == 0 else f"{i:6.4f}" for i in row) for row in f_pp(a, c, 64, None, False, nnp)))

print()

# demo large network
n, m, p, c, nnp = 1000, 256, 0.02, 0.9, '11'
adj = mock_data(n, p)
print(f'size {n}, maxiter {m}, connection prob {p}, c {c}')
from time import time
t = time()
spp_loop = f_pp_loop(adj, c, m, None, False, nnp)
print('pp_loop: {:8.5f} sec'.format(time()-t))
print('squared rel err after {} iterations: {}'.format(len(spp_loop)-1, np.nanmax(((spp_loop[-1]-spp_loop[-2])/(spp_loop[-1]+spp_loop[-2]))**2)))
t = time()
spp = f_pp(adj, c, m, None, False, nnp)
print('pp: {:8.5f} sec'.format(time()-t))
def comp(a, b):
    print(np.allclose(a, b))
    print('abserr', np.nanmax(np.abs(a-b)), 'relerr', np.nanmax(2 * np.abs(a-b)/(a+b)))
print('fixed points equal:', end=' ')
comp(spp_loop[-1], spp)

Ответ 2

Я использую массив NumPy для хранения вывода (для упрощения создания и индексирования) и networkx для создания примера графика и предоставления наборов соседей. Ни один из модулей не является существенным, поскольку я не использую никаких дополнительных функций. Вот настройка:

import numpy as np
import networkx as nx 
G = nx.petersen_graph()    # example graph
vlist = list(G.nodes)      # list of vertices
V = len(vlist)             # number of vertices
kmax = 4                   # maximal k
C = 2                      # mysterious constant C
R = np.zeros((kmax, V, V)) # array of results, R[k, i, j]
np.fill_diagonal(R[0], 1)  # initialize R[0] by putting 1 on its diagonal

И вот цикл, вычисляющий R [k + 1, i, j]. Основными символами являются Lp и Lq, множества соседей p и q. Наборы Python пересекаются с Lp & Lq, соединены с Lp | Lq и вычитаются с помощью Lp - Lq. Строки с присвоением Rij представляют собой довольно много дословных копий данных формул.

for k in range(kmax-1):
  for i in range(V):
    for j in range(V):
      Lp, Lq = set(G[vlist[i]]), set(G[vlist[j]])
      Rij = len(Lp & Lq) / len(Lp | Lq)
      Rij += sum([R[k, p1, q1] for q1 in Lq for p1 in Lp - Lq]) / (len(Lp | Lq) * len(Lq))
      Rij += sum([R[k, p1, q1] for p1 in Lp for q1 in Lq - Lp]) / (len(Lp | Lq) * len(Lp))
      R[k+1, i, j] = C*Rij