复杂网络与网络博弈
本文省略大量预备知识，并以极简化的方式介绍网络博弈和演化动力学。
常见三种复杂网络
复杂网络本质是一种高连接属性的拓扑结构。在网络上进行博弈过程，网络结构本身对博弈过程会产生影响。
常用的三种网络有：ER随机网络（Erdős Rényi Random）、BA无标度网络（Barabási Albert Scale-free）、小世界网络（Small-world）。
ER网络和SW网络属于均质网络，ER网络中每个节点的度均相同，收益信息在网络上的传导同样快，用于研究合作行为在“人人平等”的情况下演化的过程。而BA无标度网络属于异质网络，其特点是度分布满足幂律特征，hub节点所传导的信息比一般节点更多。SW网络充分体现“六度分隔”现象，常用于研究合作行为在小世界现象下如何演化。
通过特定参数下的模拟，这三种拓扑结构基本如下图：
所有图的规模N = 50，ER网络的连边概率p = 0.01；SW网络中每个节点的邻居数量nei = 3，重连接概率为p = 0.2；SF无标度网络的平均度<k> = 6。
网络结构上，视每个节点为一个博弈主体，而将每条边视为玩家间的直接联系，这样的种群被称为结构种群。
具有收益矩阵的常见网络博弈
网络博弈多种多样，我们能在不同的博弈类型中得到合作行为演化的不同轨迹。常见的博弈框架根据收益计算的不同基本分为具有收益矩阵的、非收益矩阵表达的。
这里讨论三种具有单边收益矩阵的两人博弈，和一种非矩阵表达的博弈，分别为：
- 囚徒困境（PDG）
- 雪堆博弈（SDG）
- 捐赠博弈（DG）
- 最后通牒（UG）
PDG、SDG、DG 的单边收益矩阵，UG 的收益计算形式分别表示为：
上图中的参数 b 和 c 通常用来控制困境强度，以便更好的研究不同程度的困境将如何影响合作行为的演化，通常将 c 设置为 1，控制 b 值从大到小逐步变化。
其中最后通牒博弈中的 (\(p_i\), \(q_i\)) 代表第 i 个 player 的策略，一般习惯在 0-1 上取值，当 i 与邻居j进行博弈时，\(p_i\) 表示个体 i 最多能提供给邻居的金额，\(q_i\) 表示个体 i 最少能接受邻居提供的金额。
上述四种博弈均只讨论两玩家的情形。
当重复博弈进行一轮时，个体 i 通过网络连边中的直接连接关系（邻接矩阵中有连边的部分），可以分别计算出它与所有邻居 j 博弈所得的收益，进而通过加总全部收益得到累计收益（Total Payoff）。
定义 M 为博弈矩阵，f 为博弈所带来的收益，S 表示个体选择的策略，这里个体间的策略我们抽象地概括为亲社会性的 “C” 策略、自私性的 “D” 策略。另外记 G 为结构种群。
则S可以表示为：
初始状态下，每个个体的策略随机给定 C or D，进一步的，定义在网络 G(V,E) 上全部与个体 i 的邻居集合为 \(\Omega_i\) 。则收益计算可以表示如下：
具体在计算时实际上是个体 i 的策略向量与个体 j 的策略向量通过收益矩阵 M 进行向量乘法运算，可以表达为：
然而分别对每一个体i反复迭代运算无形中增加计算负担，于是根据原理，可以利用邻接矩阵对上式重新定义，首先定义邻接矩阵 A：
再记 \(\Phi_{ij}\) 为全部 \(S_i^T M S_j\) 构成的矩阵，则直接将 \(\Phi_{ij}\) 与邻接矩阵做 Hadamard 积（对应位置相乘），i 与 j 不相连的项全部乘零消失，留下的项均为有连边的收益计算结果，收益矩阵表示为：
自然的，个体 i 的累计收益表示为 F 的第 i 行或列求和即可：
根据上述步骤进行一轮博弈，计算出全部个体的累计收益向量后需要更新每个个体的初始策略。常见的动力学更新规则有：
- 模仿更新
- 从众更新
- 期望更新
- 自反更新
这里的只采用 Fermi 更新，一种简单的模仿更新。该模仿规则能够根据个体收益差计算出个体 i 的策略互换概率 w，即：
其中 \(\kappa\) 可以理解为决策时产生的扰动，一般取很小的值，例如0.01。
Fermi 更新在通过概率 w 决定是否要与对手互换策略时，如果个体 i 有多个邻居之间相连，则在满足 w 条件时，需要先选择其中一个邻居进行策略互换，这里无论策略好坏（所带来的收益高低），每个邻居被抽到的概率均等。
Fermi 在更新过程中也可以遵循两种不同的动力学规则，即同步更新或者异步更新:
同步更新是指一轮博弈结束后，所有 N 个节点需要顺序遍历，也就是每个节点都需要进行且只进行一次更新。
异步更新指的是系统总共需要迭代 N 次（节点个数次），每次更新的对象从 N 个个体中等概率抽样，也就是最极端的情况下可能有节点从未更新过，或者有节点重复更新 N 次。
事实上进行博弈的过程中你会发现，随着博弈重复进行下去，在特定的更新机制下复杂系统中个体策略将不断变化，直到一定步长时，C 策略和 D 策略各占据一定比例不再改变（非理想状态下应该围绕某个比例小幅度波动）。
动力学数据生成实现
生成的数据主要是累计收益向量 \(Y = \sum{F_{N \cdot}}\)，以及每轮博弈的矩阵\(\Phi_{ij}\)，即有多少轮模拟就会有多少列Y，有多少个节点就会有多少个 \(\Phi_{ij}\)
原创，仅供自己回忆使用，请勿转载
For R
gen_dynamic<-function(N, L, aver_k=6, net_type="BA", game_type="PDG"){
# Calculate the payoff between i and j
# which need Payoff Matrix（for others except SUG）
payoff<-function(si, sj, M){
return(t(as.matrix(si)) %*% M %*% as.matrix(sj))
}
# which no need Payoff Matrix（for SUG）
payoff_sug<-function(s_i, s_j, S){
p_i<-S[,s_i][1]; q_j<-S[,s_j][2]; p_j<-S[,s_j][1]; q_i<-S[,s_i][2]
if ((p_i >= q_j) & (p_j >= q_i)) {p_sug<-p_j + 1 - p_i}
if ((p_i >= q_j) & (p_j < q_i)) {p_sug<-1 - p_i}
if ((p_i < q_j) & (p_j >= q_i)) {p_sug<-p_j}
if ((p_i < q_j) & (p_j < q_i)) {p_sug<-0}
return(p_sug)
}
# Generate different Game structure except SUG
gen_PDG<-function(m.R, m.S, m.T, m.P){
return(matrix(c(m.R, m.S, m.T, m.P), 2, 2, byrow=T)) # PDG Payoff Matrix
}
gen_SDG<-function(m.r){
return(matrix(c(1, 1-m.r, 1+m.r, 0), 2, 2, byrow=T)) # SDG Payoff Matrix
}
gen_DG<-function(m.b, m.c){
return(matrix(c(m.b-m.c, -m.c, m.b, 0), 2, 2, byrow=T)) # DG Payoff Matrix
}
################################### Simulation of Dynamics ############################################
# ER and SW networks are transformed into degree control
p_or_m<-function(N, k){
return(round(k / (N-1)))
}
# Generate different Complex Networks
if(net_type == "BA"){G<-barabasi.game(n=N, power=1, m=aver_k, directed = F)} # BA Scale-Free
if(net_type == "ER"){G<-erdos.renyi.game(n=N, p.or.m=p_or_m(N, aver_k), type="gnp", directed = F)} # ER Random
if(net_type == "SW"){G<-sample_smallworld(dim=1, size=N, nei=aver_k/2, p=.0001)} # SW Regular
# Init Payoff Matrix
if(game_type == "PDG"){M<-gen_PDG(1, 0, 1.2, 0)} # R、S、T、P参数随意设置
if(game_type == "SDG"){M<-gen_SDG(1.2)} # r参数随意设置
if(game_type == "DG"){M<-gen_DG(2, 1)} # b、c参数随意设置
# Init Strategy
if(game_type == "SUG"){
# for SUG
S<-as.data.frame(matrix(NA, 2, N))
for (i in 1:N) {
S[,i]<-c(runif(1), runif(1))
}
}else{
# For others
S<-as.data.frame(matrix(NA, 2, N))
for (i in 1:N) {
if(runif(1) >= .5){
S[,i]<-c(0, 1)
}else{
S[,i]<-c(1, 0)
}
}
}
# Record Adj
real_adj<-as.matrix(get.adjacency(G))
storage.mode(real_adj)<-"integer" # 因为Adj全是0-1，所以强行由double型储存为integer型
# init X, Y
Y<-data.frame(matrix(NA, L, N)) # 每一列为一个时间步的累计收益
X<-list() # 为包含N个矩阵的列表（即分块对角矩阵分块存储）
for (i in 1:N) {
X[[i]]<-matrix(NA, L, N)
}
# simulation process
for (l in 1:L) {
F_temp<-matrix(NA, N, N)
if(game_type == "SUG"){
for (u in 1:N) {
for (v in 1:N) {
F_temp[u, v]<-payoff_sug(u, v, S)
}
}
}else{
for (u in 1:N) {
for (v in 1:N) {
F_temp[u, v]<-payoff(S[,u], S[,v], M)
}
}
}
# Record Y and X
F_cumulate<-apply((real_adj*F_temp), 1, sum)
# F_cumulate<-apply((real_adj*F_temp), 2, sum)
Y[l, ]<-F_cumulate # The row of Y
for (n in 1:N) {
X[[n]][l,]<-F_temp[n,]
}
# Fermi Updating (Synchronization Update)
if(game_type == "SUG"){
for (i in 1:N) {
update_index<-sample(as.vector(neighbors(G, i, "all")), 1)
w<-1/(1 + exp((F_cumulate[i] - F_cumulate[update_index])/.1))
if(runif(1) >= w){
S[, i]<-S[, update_index]
S[, i]<-S[, i] + .05
}
}
}else{
for (i in 1:N) {
update_index<-sample(as.vector(neighbors(G, i, "all")), 1)
w<-1/(1 + exp((F_cumulate[i] - F_cumulate[update_index])/.1))
if(runif(1) >= w){
S[, i]<-S[, update_index]
}
}
}
}
# X block diag treatment
# X_diag_block<-as.matrix(bdiag(X)) #使用此代码将X保存为分块对角矩阵需要用到开头备注的 “bigmemory” 包
return(list(Y = Y, X = X, real_adj = real_adj))
}
# dyn_vars<-gen_dynamic(100, 15, 6, "BA", "SUG")
# dyn_vars[["Y"]]; dyn_vars[["X"]]; dyn_vars[["real_adj"]]
For Python
import numpy as np
import math, random
import networkx as nx
# Gen Network
class Gen_Net:
def basf(n: int, k: int):
G_sf = nx.barabasi_albert_graph(n, k)
return G_sf
def ER(n: int, p: float):
G_er = nx.random_graphs.erdos_renyi_graph(n, p)
return G_er
def SW(n: int, k: int, p: float):
G_sw = nx.random_graphs.watts_strogatz_graph(n, k, p)
return G_sw
# Define Useful Tools
class Use_Tools:
def chunk(lst, size): # make selected array sort with rule size
return list(map(lambda x: lst[x * size: x * size + size], list(range(0, math.ceil(len(lst) / size)))))
def payoff(si, M, sj) -> float: # si*M*sj' and input is np.array
phi = np.matmul(np.matmul(si, M), sj.T)
return phi
# Define Game Structure
class Game:
def PDG(b: float):
M = np.array([[1, 0], [b, 0]]) # b = 1.2; M[0][0]=1; M[0][1]=0; M[1][0]=1.2; M[1][1]=0
return M
def SDG(r: float):
M = np.array([[1, 1 - r], [1 + r, 0]])
return M
# Initalize Strategies Set
G = Gen_Net.basf(100, 6) # average degree <k> = 6: k = 3
node_list = np.array(list(G.nodes)) + 1 # all nodes
edge_list = np.array(G.edges) + 1 # linkage's edges
nebmat = nx.to_numpy_matrix(G) # Adjacent Matrix of G
# np.savetxt('Adjacent_Matrix.csv', nebmat, delimiter = ',')
s_set = [] # every node which choose C or D with equal prob
for _ in range(len(node_list)):
if np.random.rand(1) > 1 / 2: s_set.append((0, 1))
else: s_set.append((1, 0)) # strategies set of #N nodes
# Monte-Carlo Simulation
L = 300; Phi = []
M = Game.PDG(1.2) # MC Rounds
Y = np.zeros([L, G.number_of_nodes()])
for _ in range(G.number_of_nodes()): # Store Phi
Phi.append(np.zeros([L, G.number_of_nodes()]))
for l in range(L):
for f in range(G.number_of_nodes()): # asynchronous updating
temp_payoff, neb_temp_payoff = [], []
focal_player = random.sample(list(node_list), 1)[0]
# print("\n", focal_player)
# find focal's neb and calculate payoff of focal player
rdn_neb = []
for i in range(len(edge_list)):
if edge_list[i][0] != focal_player:
if edge_list[i][1] == focal_player:
rdn_neb.append(edge_list[i][0])
temp_payoff.append(Use_Tools.payoff(np.array(s_set[edge_list[i][1] - 1]), M, np.array(s_set[edge_list[i][0] - 1])))
else:
rdn_neb.append(edge_list[i][1])
temp_payoff.append(Use_Tools.payoff(np.array(s_set[edge_list[i][0] - 1]), M, np.array(s_set[edge_list[i][1] - 1])))
# print("\n", rdn_neb)
# print("\n", temp_payoff)
f_i = np.sum(temp_payoff)
Y[l][focal_player - 1] = f_i
# find neb's neb and calculate payoff of player's neb
select_neb = random.sample(rdn_neb, 1)[0]
for i in range(len(edge_list)):
if edge_list[i][0] != select_neb:
if edge_list[i][1] == select_neb:
neb_temp_payoff.append(Use_Tools.payoff(np.array(s_set[edge_list[i][1] - 1]), M, np.array(s_set[edge_list[i][0] - 1])))
else:
neb_temp_payoff.append(Use_Tools.payoff(np.array(s_set[edge_list[i][0] - 1]), M, np.array(s_set[edge_list[i][1] - 1])))
# print("\n", select_neb)
# print("\n", neb_temp_payoff)
f_j = np.sum(neb_temp_payoff)
# Fermi updating
# fi_lambda = np.sum(nebmat[focal_player - 1]) # normlization
# fj_lambda = np.sum(nebmat[select_neb - 1])
# w = 1 / (1 + math.exp(((f_i / fi_lambda) - (f_j / fj_lambda)) / 0.1)) # k = .1
w = 1 / (1 + math.exp((f_i - f_j) / 0.1)) # k = .1
if np.random.rand(1) <= w: s_set[focal_player - 1] = s_set[select_neb - 1]
print("\n th {} Rounds of Fermi updating:\n {}".format(f, s_set))
# find Phi
for m in range(G.number_of_nodes()): # #(m) matrix
for c in range(G.number_of_nodes()): # #(col) of #(m)
Phi[m][l][c] = Use_Tools.payoff(np.array(s_set[m]), M, np.array(s_set[c]))
print("\n\n th {} Rounds of MC Simulation".format(l))
print("\n\n Y = \n {} \n\n Dim(Y) = {}".format(Y, Y.shape))
# check weather Y and Phi are correct
Y1_hat = np.matmul(np.matrix(Phi[0]),nebmat[:,0])
Y1 = Y[:,0].reshape(L, 1)
Y1 == Y1_hat