安全通论作业一

作业内容

用matlab/python等作图工具，数值分析攻击信道容量。

分析报告

基本概念与公式

在单挑盲攻防的背景下，将攻守双方过去 N 次对抗的“盲自评结果”记录下来，得到一个二维随机变量（X 和 Y 分别表示攻击方和守方盲自评的结果，0表示失败，1表示成功）：
$$(X,Y)=(X_1,Y_1),(X_2,Y_2),...,(X_N,Y_N)$$

当 N 趋于无穷大时，频率趋于概率Pr，所以，只要攻守双方足够长
时间对抗之后，得到随机变量X、Y的概率分布如下：

Pr(攻方盲自评为成功) = Pr(X=1) = p
Pr(攻方盲自评为失败) = Pr(X=0) = 1 - p，0 < p < 1
Pr(守方盲自评为成功) = Pr(Y=1) = q
Pr(守方盲自评为失败) = Pr(Y=0) = 1 - q，0 < q < 1
Pr(攻方盲自评为成功，守方盲自评为成功) = Pr(X=1，Y=1) = a，0 < a < 1
Pr(攻方盲自评为成功，守方盲自评为失败) = Pr(X=1，Y=0) = b，0 < b < 1
Pr(攻方盲自评为失败，守方盲自评为成功) = Pr(X=0，Y=1) = c，0 < c < 1
Pr(攻方盲自评为失败，守方盲自评为失败) = Pr(X=0，Y=0) = d，0 < d < 1

即，p 和 q 是攻方和守方盲自评成功的概率，a、b、c、d 表示攻方和守方盲自评结果的联合概率。X 和 Y 的联合概率与边缘概率表示如下：

这里，a、b、c、d、p、q 之间还满足下列关系：

• a + b + c + d = 1
• p = Pr(X=1) = Pr(X=1,Y=0) + Pr(X=1,Y=1) = a + b
• q = Pr(Y=1) = Pr(X=1,Y=1) + Pr(X=0,Y=1) = a + c
  由此可知，6 个变量 a、b、c、d、p、q 中，其实只有 3 个是独立的。我们设这三个标量分别是 a、p 和 q。

黑客攻击能力的极限

以 X 作为输入，Z 作为输出，构造出一个通信信道 F，称之为攻击信道。

如果黑客的某次攻击真正成功，那么攻击信道 F 就成功地传输 1 bit 到收端；如果有 1 bit 被成功地从攻击信道 F 的发端，传送到了收端，那么黑客 X 就获得了一次真正成功攻击。

引理 1：黑客获得一次真正成功的攻击，其实就对等于攻击信道 F 成功地传输了 1 bit。

定理 1（黑客攻击能力极限定理）：设由随机变量(X；Z)组成的攻击信道 F 的信道容量为 C。那么：如果黑客想真正成功地把红客打败 k 次，一定有某种技巧（对应于仙农编码），使得他能够在 k/C 次攻击中，以任意接近 1 的概率达到目的；如果黑客经过 n 次攻击，获得了 S 次真正成功的攻击，一定有 S ≤ nC。

根据定理 1，只要求出攻击信道 F 的信道容量 C，那么黑客的攻击能力极限就确定了。

攻击信道容量

根据补充章节 1，在单用户信道中，信道容量 = max 互信息。

随机变量 (X,Z) 的联合概率分布为：

Pr (X=0, Z=0) = Pr (X=0,Y=0) = d
Pr (X=0, Z=1) = Pr (X=0,Y=1) = c
Pr (X=1, Z=0) = Pr (X=1,Y=1) = a
Pr (X=1, Z=1) = Pr (X=1,Y=0) = b

结合 X 和 Z 构成的通信系统 F 的转移矩阵
A = [ A ( 0 , 0 ) A ( 0 , 1 ) A ( 1 , 0 ) A ( 1 , 1 ) ] = [ d 1 − p 1 − d 1 − p a p 1 − a p ] A =

$$\begin{bmatrix} A(0,0) & A(0, 1) \\ A(1,0) & A(1, 1) \\ \end{bmatrix}$$ = $$\begin{bmatrix} \frac{d}{1-p} & 1 - \frac{d}{1-p} \\ \frac{a}{p} & 1 - \frac{a}{p} \\ \end{bmatrix}$$ A=[A(0,0)A(1,0)​A(0,1)A(1,1)​]=[1−pd​pa​​1−1−pd​1−pa​​]

随机变量 X 和 Z之间的互信息为
I ( X , Z ) = ∑ x ∑ z p ( x , z ) log ⁡ p ( x , z ) p ( x ) p ( z ) = d log ⁡ d ( 1 − p ) ( a + d ) + c log ⁡ c ( 1 − p ) ( b + c ) + a log ⁡ a p ( a + d ) + b log ⁡ b p ( b + c )

$$\begin{aligned} I(X, Z) &= \sum_x\sum_z p(x, z)\log\frac{p(x, z)}{p(x)p(z)} \\ &= d\log\frac{d}{(1-p)(a+d)} + c\log\frac{c}{(1-p)(b+c)} \\ &+ a\log\frac{a}{p(a+d)} + b\log\frac{b}{p(b+c)} \end{aligned}$$ I(X,Z)​=x∑​z∑​p(x,z)logp(x)p(z)p(x,z)​=dlog(1−p)(a+d)d​+clog(1−p)(b+c)c​+alogp(a+d)a​+blogp(b+c)b​​

又由于 a + b + c + d = 1，p = a + b，q = a + c，0 < a, b, c, d, p, q < 1，可以进一步转化为只与变量a和p有关的如下公式（注意：此时q已不再是变量，而是确定值了）
I ( X , Z ) = [ 1 + a − p − q ] log ⁡ 1 + a − ( p + q ) ( 1 − p ) ( 1 + 2 a − p − q ) + ( q − a ) log ⁡ q − a ( 1 − p ) ( p + q − 2 a ) + a log ⁡ a p ( 1 + 2 a − p − q ) + ( p − a ) log ⁡ p − a p ( p + q − 2 a )

$$\begin{aligned} I(X, Z) &= [1+a-p-q]\log\frac{1+a-(p+q)}{(1-p)(1+2a-p-q)} \\ &+ (q-a)\log\frac{q-a}{(1-p)(p+q-2a)} \\ &+ a\log\frac{a}{p(1+2a-p-q)} + (p-a)\log\frac{p-a}{p(p+q-2a)} \end{aligned}$$ I(X,Z)​=[1+a−p−q]log(1−p)(1+2a−p−q)1+a−(p+q)​+(q−a)log(1−p)(p+q−2a)q−a​+alogp(1+2a−p−q)a​+(p−a)logp(p+q−2a)p−a​​

利用此 $I(X,Z)$ 就可知：以 X 为输入，Z 为输出的信道 F 的信道容量 $C=max[I(X,Z)]$（这里最大值是针对 X 为所有可能的二元离散随机变量来计算的）。
在此处一对一的信道中，其容量 $C=ma{x}_{0<a,p<1}[I(X,Z)]$，这里的最大值是对仅仅两个变量 a 和 p 在条件 0 < a，p < 1 下之取的，也就是说 攻击信道容量 C，其实是 q 的函数。

攻击信道容量的数值分析

根据攻击信道容量公式

C = m a x 0 < a , p < 1 [ I ( X , Z ) ] = m a x 0 < a , p < 1 [ ( 1 + a − p − q ) log ⁡ 1 + a − ( p + q ) ( 1 − p ) ( 1 + 2 a − p − q ) + ( q − a ) log ⁡ q − a ( 1 − p ) ( p + q − 2 a ) + a log ⁡ a p ( 1 + 2 a − p − q ) + ( p − a ) log ⁡ p − a p ( p + q − 2 a ) ]

$$\begin{aligned} C &= max_{0<a,p<1}[I(X,Z)] \\ &= max_{0<a,p<1} [ (1+a-p-q)\log\frac{1+a-(p+q)}{(1-p)(1+2a-p-q)} \\ &+ (q-a)\log\frac{q-a}{(1-p)(p+q-2a)} \\ &+a\log\frac{a}{p(1+2a-p-q)} + (p-a)\log\frac{p-a}{p(p+q-2a)} ] \end{aligned}$$ C​=max0<a,p<1​[I(X,Z)]=max0<a,p<1​[(1+a−p−q)log(1−p)(1+2a−p−q)1+a−(p+q)​+(q−a)log(1−p)(p+q−2a)q−a​+alogp(1+2a−p−q)a​+(p−a)logp(p+q−2a)p−a​]​

下面给出以 q 为自变量的攻击信道容量的数值分析函数：

# Compute mutual information
def get_i(a, p ,q):
	c1 = 1 + a - p - q
	c2 = q - a
	c3 = p - a
	c4 = 1 + 2 * a - p - q
	c5 = p + q - 2 * a

	i1 = c1 * np.log(c1 / ((1 - p) * c4))
	i2 = c2 * np.log(c2 / ((1 - p) * c5))
	i3 = a * np.log(a / (p * c4))
	i4 = c3 * np.log(c3 / (p * c5))

	i = i1 + i2 + i3 + i4
	return i

# def get_max_i(i):
def get_max_i(a, p ,q):
	i = get_i(a, p ,q)
	i = np.nan_to_num(i)
	max_i = np.max(i)
	# print(max_i)
	# print(i)
	return max_i
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

其中 get_i 函数是根据公式 I(X, Z) 计算互信息，get_max_i 函数是根据 $C=ma{x}_{0<a,p<1}[I(X,Z)]$ 计算信道容量，即最大的互信息。画出攻击信道容量的数值分析图像如下：

import numpy as np
import matplotlib.pyplot as plt

_p = np.arange(0.01, 1.00, 0.01)
_q = np.arange(0.01, 1.00, 0.01)

# Meshgrid for 2D array creation
p, q = np.meshgrid(_p, _q)

# 2D array outer product
a = np.outer(_p, _q)

# Compute mutual information
def get_i(a, p ,q):
	c1 = 1 + a - p - q
	c2 = q - a
	c4 = p - a
	c5 = 1 + 2 * a - p - q
	c6 = p + q - 2 * a

	i1 = c1 * np.log(c1 / ((1 - p) * c5))
	i2 = c2 * np.log(c2 / ((1 - p) * c6))
	i3 = a * np.log(a / (p * c5))
	i4 = c4 * np.log(c4 / (p * c6))

	i = i1 + i2 + i3 + i4
	return i

i = get_i(a, p ,q)

# Plotting
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(p, q, i, cmap='jet')

# Axis labels
ax.set_xlabel('p')
ax.set_ylabel('q')
ax.set_zlabel('I(X,Z)')

# Title
ax.set_title('Image of I(X,Z) when a=pq')

# plt.show()
plt.savefig('Ivspq.png', dpi=200)

# print(np.argmin(i), np.max(i))
# Reshape the 2D array i into a 1D array
i_flat = i.flatten()

# Find the index of the maximum value in the 1D array
max_idx = np.argmax(i_flat)

# Find the corresponding p and q values using the index
max_p_idx, max_q_idx = np.unravel_index(max_idx, i.shape)
max_p = p[max_p_idx, max_q_idx]
max_q = q[max_p_idx, max_q_idx]

print("Maximum I(X,Z) value:", i_flat[max_idx])
print("Corresponding p value:", max_p)
print("Corresponding q value:", max_q)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61

完整绘图代码见 hw1-Cvsq.py 文件。从图中可以看到：攻击信道容量 C 随 q 的增大呈现先减小后增大的趋势，即攻击信道容量 C 是一个凹函数。在 q = 0.5 附近时，攻击信道容量 C 达到最小值 0.069。

a = pq (p 和 q 互相独立) 时，互信息与 p 和 q 的关系

完整绘图代码见 hw1-Ivspq.py 文件。从图中可以看到：

• 当 q 趋近于 0 或 1 时，互信息最大；
• 当 q 趋近于 0.5 时，互信息最小；
• 当 q 不变，p 趋近于 0.5 时，互信息相对较大。

import numpy as np
import matplotlib.pyplot as plt

a_range = np.arange(0.01, 1.00, 0.001)
p_range = np.arange(0.01, 1.00, 0.001)
q_range = np.arange(0.02, 0.99, 0.001)

# Meshgrid for 2D array creation
a, p = np.meshgrid(a_range, p_range)

# Compute mutual information
def get_i(a, p ,q):
	c1 = 1 + a - p - q
	c2 = q - a
	c3 = p - a
	c4 = 1 + 2 * a - p - q
	c5 = p + q - 2 * a

	i1 = c1 * np.log(c1 / ((1 - p) * c4))
	i2 = c2 * np.log(c2 / ((1 - p) * c5))
	i3 = a * np.log(a / (p * c4))
	i4 = c3 * np.log(c3 / (p * c5))

	i = i1 + i2 + i3 + i4
	return i

# def get_max_i(i):
def get_max_i(a, p ,q):
	i = get_i(a, p ,q)
	i = np.nan_to_num(i)
	max_i = np.max(i)
	# print(max_i)
	# print(i)
	return max_i

max_i = np.zeros_like(q_range)
for i, q in enumerate(q_range):
	# i = get_i(a, p ,q)
	max_i[i] = get_max_i(a, p ,q)

# Plotting
fig = plt.figure()
ax = fig.add_subplot()
ax.plot(q_range, max_i)

# Axis labels
ax.set_xlabel('q')
ax.set_ylabel('max I(X,Z)')

# Title
ax.set_title('Image of C = max I(X,Z) vs q')

plt.show()
# plt.savefig('Cvsq.png', dpi=100)

print(q_range[np.argmin(max_i)], np.min(max_i))
print(q_range[np.argmax(max_i)], np.max(max_i))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57

固定 q 时，互信息与 p 和 a 的关系

以下为 q = 0.1 到 q = 0.9 时， 互信息 I(X, Z) 与 p 和 a 的关系。（代码见 hw1-Ivsap.py）

import numpy as np
import matplotlib.pyplot as plt

q = 0.1

# Compute mutual information
def get_i(a, p ,q):
	c1 = 1 + a - p - q
	c2 = q - a
	c3 = p - a
	c4 = 1 + 2 * a - p - q
	c5 = p + q - 2 * a

	i1 = c1 * np.log(c1 / ((1 - p) * c4))
	i2 = c2 * np.log(c2 / ((1 - p) * c5))
	i3 = a * np.log(a / (p * c4))
	i4 = c3 * np.log(c3 / (p * c5))

	i = i1 + i2 + i3 + i4
	return i

# def get_max_i(i):
def get_max_i(a, p ,q):
	i = get_i(a, p ,q)
	i = np.nan_to_num(i)
	max_i = np.max(i)
	# print(max_i)
	# print(i)
	return max_i

# while (q <= 0.9):
# 	p_range = np.arange(0, 1, 0.01)
# 	a_range = np.arange(0, q, q/100)
# 	p, a = np.meshgrid(p_range, a_range)
	
# 	i = get_i(a, p, q)

# 	fig = plt.figure()
# 	ax = fig.add_subplot(projection='3d')
# 	ax.plot_surface(p, a, i, cmap='jet')
# 	string = 'Image of I(X,Z) when q = %.2f' % q
# 	ax.set_zlim3d(zmin=0, zmax=0.5)
# 	ax.set_title(string)
# 	ax.set_xlabel('p')
# 	ax.set_ylabel('a')
# 	ax.set_zlabel('I(X,Z)')
# 	# plt.savefig('q=%.2f.png' % q, dpi=200)
# 	plt.show()

# 	q = q + 0.1

fig, axs = plt.subplots(3, 3, figsize=(15, 15))

for i, ax in enumerate(axs.flatten()):
	q = i / 10.0 + 0.1
	p_range = np.arange(0, 1, 0.001)
	a_range = np.arange(0, q, q/1000)
	p, a = np.meshgrid(p_range, a_range)

	i_val = get_i(a, p, q)

	ax.axis("off")
	ax = fig.add_subplot(3, 3, i + 1, projection='3d')
	string = 'Image of I(X,Z) when q = %.2f' % q
	ax.plot_surface(p, a, i_val, cmap='jet')
	ax.set_title(string)
	ax.set_xlabel('p')
	ax.set_ylabel('a')
	ax.set_zlabel('I(X,Z)')
	ax.set_zlim3d(zmin=0, zmax=0.5)

plt.tight_layout()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73

在 q 从 0.1 变到 0.9 时，可以发现：

• 在局部图像上，图像中的曲面随着 q 的增大，呈现先凹后凸的趋势，中间部分相对于便于两边时凸起的，随后开始下陷，即两边升起；
• 在 q 增大的全过程中，互信息随着 p 先增大后减小的，即从 p 的角度看，互信息始终是一个凸函数；但凸起的程度则是先减小后增大的，即凸起的程度是一个凹函数，这个分界点在 p = 0.5 附近；
• 在 q 增大的全过程中，互信息随着 a 的变化则则是先凹函数，后凸函数，即随着 q 的增大，互信息随着 a 的变化趋势也在变化；
• 固定 a，随着 q 的增大，互信息与 p 的关系由凸函数变为凹函数；
• 固定 p，随着 q 的增大，互信息与 a 的关系始终是凹函数。
• 结合以上的第 2 和 4 条、第 3 和 5 条，进一步可以发现，互信息随着 p 和 a 的变化，在宏观和微观上有着互补的变化趋势。