向量化编程思路(矩阵计算)
Author: Labyrinthine Leo
Init_time: 2020.11.24
Index Words: Matrix Compuation
、python
、for循环
原文链接: 向量化编程小结
公众号:Leo的博客城堡
1、DTLZ测试函数
1.1、DTLZ1
公式:
- LaTeX语法
**Definition**
$$
\begin{align}
Minimize \quad & f_1(X) = \frac{1}{2}x_1x_2 \dots x_{M-1}(1+g(X_M)), \\
Minimize \quad & f_2(X) = \frac{1}{2}x_1x_2 \dots (1-x_{M-1})(1+g(X_M)),\\
\vdots \quad & \vdots \\
Minimize \quad & f_{M-1}(X) = \frac{1}{2}x_1(1-x_2)(1+g(X_M)),\\
Minimize \quad & f_M(X) = \frac{1}{2}(1-x_1)(1+g(X_M)),\\
subject \; to \quad & 0 \leq x_i \leq 1, \quad for \; i = 1,2,\cdots,n.
\end{align}
$$
The last $k = (n - M + 1)$ variables are represented as $X_M$. The functional $g(X_M)$ requires $\lvert X_M \rvert = k$ variables and must take any function with $g \geq 0$.
We suggest the following:
$$
g(X_M) = 100 \Big[ \lvert X_M \rvert + \sum_{x_i \in X_M}(x_i - 0.5)^2 - cos(20\pi(x_i - 0.5)) \Big]
$$
**Optimum**
The Pareto-optimal solution corresponds to $x_i = 0.5$(for all $x_i \in X_M$)and the objective function values lie on the linear hyper-plane:$\sum_{m=1}^M f_m^*=0.5$.
- LaTeX公式显示
Definition
M
i
n
i
m
i
z
e
f
1
(
X
)
=
1
2
x
1
x
2
…
x
M
−
1
(
1
+
g
(
X
M
)
)
,
M
i
n
i
m
i
z
e
f
2
(
X
)
=
1
2
x
1
x
2
…
(
1
−
x
M
−
1
)
(
1
+
g
(
X
M
)
)
,
⋮
⋮
M
i
n
i
m
i
z
e
f
M
−
1
(
X
)
=
1
2
x
1
(
1
−
x
2
)
(
1
+
g
(
X
M
)
)
,
M
i
n
i
m
i
z
e
f
M
(
X
)
=
1
2
(
1
−
x
1
)
(
1
+
g
(
X
M
)
)
,
s
u
b
j
e
c
t
t
o
0
≤
x
i
≤
1
,
f
o
r
i
=
1
,
2
,
⋯
,
n
.
\begin{aligned} Minimize \quad & f_1(X) = \frac{1}{2}x_1x_2 \dots x_{M-1}(1+g(X_M)), \\ Minimize \quad & f_2(X) = \frac{1}{2}x_1x_2 \dots (1-x_{M-1})(1+g(X_M)),\\ \vdots \quad & \vdots \\ Minimize \quad & f_{M-1}(X) = \frac{1}{2}x_1(1-x_2)(1+g(X_M)),\\ Minimize \quad & f_M(X) = \frac{1}{2}(1-x_1)(1+g(X_M)),\\ subject \; to \quad & 0 \leq x_i \leq 1, \quad for \; i = 1,2,\cdots,n. \end{aligned}
MinimizeMinimize⋮MinimizeMinimizesubjecttof1(X)=21x1x2…xM−1(1+g(XM)),f2(X)=21x1x2…(1−xM−1)(1+g(XM)),⋮fM−1(X)=21x1(1−x2)(1+g(XM)),fM(X)=21(1−x1)(1+g(XM)),0≤xi≤1,fori=1,2,⋯,n.
The last
k
=
(
n
−
M
+
1
)
k = (n - M + 1)
k=(n−M+1) variables are represented as
X
M
X_M
XM. The functional
g
(
X
M
)
g(X_M)
g(XM) requires
∣
X
M
∣
=
k
\lvert X_M \rvert = k
∣XM∣=k variables and must take any function with
g
≥
0
g \geq 0
g≥0.
We suggest the following:
g
(
X
M
)
=
100
[
∣
X
M
∣
+
∑
x
i
∈
X
M
(
x
i
−
0.5
)
2
−
c
o
s
(
20
π
(
x
i
−
0.5
)
)
]
g(X_M) = 100 \Big[ \lvert X_M \rvert + \sum_{x_i \in X_M}(x_i - 0.5)^2 - cos(20\pi(x_i - 0.5)) \Big]
g(XM)=100[∣XM∣+xi∈XM∑(xi−0.5)2−cos(20π(xi−0.5))]
Optimum
The Pareto-optimal solution corresponds to
x
i
=
0.5
x_i = 0.5
xi=0.5(for all
x
i
∈
X
M
x_i \in X_M
xi∈XM)and the objective function values lie on the linear hyper-plane:
∑
m
=
1
M
f
m
∗
=
0.5
\sum_{m=1}^M f_m^*=0.5
∑m=1Mfm∗=0.5.
- 解析
首先简单看一下这个测试函数,该测试函数的决策向量为n
维,其中决策变量取值范围为[0,1]
;包含M
个目标,其中每个目标函数中均含有1+g(XM)
项,这里的XM
表示X
向量的后n-M+1
个元素组成的尾部向量。
笔者在编写该测试函数时,惯有的python
思维一定是直接使用显式的for
循环用每个目标函数对决策变量进行求解,然后对每个目标函数值进行整合为最终的目标向量,这里笔者给出最低效的显示for
循环写法:
import numpy as np
def DTLZ1(X: np.ndarray, M: int = 3) -> np.ndarray:
"""
测试函数DTLZ1
:param X: 决策向量
:param M: 目标数
:return: 目标向量
"""
# 目标数大于等于3,假设决策向量X是2维ndarray,可接受多个个体
assert X.shape[1] >= M+4, 'Dimension is invalid!'
n = X.shape[1] # 决策变量维度
objvs_list = [] # 目标值向量
for x in X: # 遍历X向量
objv_list = [] # 单个个体目标向量
XM = x[M-1:] # XM向量
for i in range(M):
fx = 0.5
for j in range(M-i-1):
fx *= x[j]
if i != 0:
fx *= (1 - x[M-i-1]) # 除了第一个目标函数均有该项
fx *= (1 + g(XM))
objv_list.append(fx)
objvs_list.append(objv_list)
objvs_list = np.array(objvs_list)
return objvs_list
def g(X: np.ndarray) -> float:
"""
计算gx函数值
:param X: n维向量
:return: 标量值
"""
lens = X.shape[0] # 获取维度
gx = 0 # 对返回值进行初始化
for i in range(lens):
gx += ((X[i] - 0.5)**2 - np.cos(20*np.pi*(X[i] - 0.5)))
gx = 100 * (lens + gx)
return gx
由上述代码可以看出,特别是在求和求平均时,使用for
循环虽然简单易懂且易用,但是都是以时间为代价的,因此我们需要转换思维,将matlab
中矩阵计算或并行计算的技巧运用起来,这在后期如果使用GPU
加速上是很方便的(笔者也是刚开始刻意培养该编码思维,因此书此博客以交流)。
于是,这里使用numpy
的矩阵计算进行改写代码:
import numpy as np
def DTLZ1(X: np.ndarray, M: int = 3) -> np.ndarray:
"""
测试函数DTLZ1
:param X: 决策向量,可以是多维,接收多个个体
:param M: 目标数
:return: 目标向量
"""
assert X.shape[1] >= M+4, 'Dimension is invalid!'
N = X.shape[0] # 获取个体个数
XM = X[:, (M-1):] # 截取矩阵中每个个体的XM部分
g = 100 * (XM[0].shape[0] + np.sum((XM - 0.5)**2 - np.cos(20*np.pi*(XM - 0.5)), axis=1, keepdims=True)) # 计算gx函数值,其中axis=1表示按列取元素,keepdims表示保持维度不变
ones_matrix = np.ones((N, 1)) # 创建(N, 1)的1矩阵
# 矩阵计算目标函数值
f = 0.5 * np.fliplr(np.cumprod(np.hstack((ones_matrix, X[:, :M-1])), axis=1)) * np.hstack((ones_matrix, 1-X[:, range(M-2, -1, -1)])) * np.tile(1+g, (1, M))
return f
可以发现一个问题,就是在for
循环中求解gx
函数的10行代码和求解目标函数的20行代码,在矩阵计算中均可用1行代码替代,这会带来一个问题就是可读性很差(矩阵计算的源码是组里平台Plat-EMO:https://github.com/BIMK/PlatEMO
中的,说实话笔者当时读那求解多维目标函数矩阵的1行代码读懂花了好久,读懂后发现运用的太巧妙了,简直为作者(大牛老师)折服,但是对初学者不太友好(还是我太菜了),因此这里讲解一下我的不成熟的理解)
我们先观察M个目标函数的公式,目的就是最后能将目标函数公式的计算分解为多个矩阵的乘积,因此首先要做的是将每个公式的共同项找到,方便后面直接同维矩阵点乘;看到有
1
2
\frac{1}{2}
21,
(
1
+
g
(
X
M
)
)
(1+g(X_M))
(1+g(XM)),这个可以直接等维化进行计算,其中还有1项
(
1
−
X
i
)
i
=
1
,
2
,
⋯
,
M
−
1
(1 - X_i) \quad i=1,2,\cdots,M-1
(1−Xi)i=1,2,⋯,M−1,除了 f1
其他目标函数均有,那我们思考是否在创建这一项的矩阵时在此矩阵的第一个元素上添加1;不同的发现f1
有
x
1
x
2
…
x
M
−
1
x_1x_2 \dots x_{M-1}
x1x2…xM−1,f2
有
x
1
x
2
…
x
M
−
2
x_1x_2 \dots x_{M-2}
x1x2…xM−2,fM
没有单个
x
i
x_i
xi的乘积项,与前面同理我们可以在构建该项矩阵时使用1进行替代。既然已经分解出各项的对应矩阵了,那就可以编写矩阵计算代码,在此之前,笔者简单介绍一下numpy
中的几个矩阵计算函数。
- np.rondom.random()
该函数根据传入的size
参数随机生成[0.0,1.0)
之间的值。>>>np.random.random((2,3)) # 产生两行三列的矩阵 [[0.90340192 0.1374747 0.13927635] [0.80739129 0.39767684 0.1653542 ]] >>>1-np.random.random((3,)) # 取值范围变为(0.0,1.0] [0.07249142 0.65223414 0.2491879 ]
- np.random.uniform()
该函数有三个参数取值的下界、上界和矩阵维度,同上,左闭右开(貌似numpy
中均是如此),然后在此范围中均匀的取size
大小的矩阵。>>>np.random.uniform(1, 10, (2, 3)) [[7.72409478 6.00616211 2.22809703] [1.53925921 2.0920911 1.40096691]]
- np.random.randint()
该函数只能取整数,注意当前两个定界的参数有一个缺失时,则默认取到[0, low)
,但此时对维度必须使用size
关键字指定,且界值必须大于0。>>>np.random.randint(-5, 10, (2, 6)) [[ 3 -4 5 -1 7 9] [ 6 9 -2 -2 6 5]] >>>np.random.randint(6, size=(2, 5)) # 必须使用size关键字,且界值大于0 [[5 4 3 0 0] [0 4 5 1 5]]
- np.sum()
注意里面的axis
参数,表示对某个维度的值进行求和,如果是0,则表示从行选取元素那就是对列进行求和,同理如果是1,则表示从列中选取元素那就是对行进行求和。我们要保持输出结果的维度和输入矩阵相同,则设置参数keepdims
为True
即可。>>>x = np.random.random((3, 10)) [[0.99732285 0.17234051 0.13713575 0.93259546 0.69681816 0.06600017 0.75546305 0.75387619 0.92302454 0.71152476] [0.12427096 0.01988013 0.02621099 0.02830649 0.24621107 0.86002795 0.53883106 0.55282198 0.84203089 0.12417332] [0.27918368 0.58575927 0.96959575 0.56103022 0.01864729 0.80063267 0.23297427 0.8071052 0.38786064 0.86354185]] >>>np.sum(x, axis=1, keepdims=True) [[6.14610144] [3.36276484] [5.50633085]] >>>np.sum(x, axis=0, keepdims=True) [[1.40077749 0.77797991 1.13294248 1.52193217 0.96167652 1.72666079 1.52726839 2.11380336 2.15291607 1.69923993]]
- np.tile()
该函数传入两个参数,一个A
矩阵,一个rpes
表示对矩阵A
复制的形状,比如reps=(3,2)
,从外到内进行复制,即按行复制三次,按列复制两次,如:>>>np.tile([1,2,3,4,5],(3, 2)) [[1 2 3 4 5 1 2 3 4 5] [1 2 3 4 5 1 2 3 4 5] [1 2 3 4 5 1 2 3 4 5]]
- np.hstack()
该函数对两个矩阵进行合并,在最里层的元素进行堆叠:
- np.fliplr()
该函数对矩阵进行按行翻转,其实等同于a[:, ::-1]
- np.cumprod()
该函数就是按照给定的axis
对元素进行求累积乘,比如[1,2,3]
生成[1 2 6]
,表示[1*1 1*2 1*2*3]
,之所以介绍该函数,正是符合上述目标函数中不断变化的 x 1 , x 2 , … x M − 1 x_1,x_2, \dots x_{M-1} x1,x2,…xM−1的乘积。
了解完这几个函数,那就可以对目标函数公式集进行拆解:
1、首先对
x
1
x
2
…
x
M
−
1
x_1x_2 \dots x_{M-1}
x1x2…xM−1部分进行构建,因为这部分是由
f
1
f_1
f1到
f
M
f_M
fM不断变化的,即项数递减,因此我们使用刚才讲的np.cumprod()
函数,又由于
f
M
f_M
fM没有该项,所以补1,np.hstack(np.ones((N,1), X[:, :M-1])
生成矩阵:
1
x
11
x
12
…
x
1
M
−
1
1
x
21
x
22
…
x
2
M
−
1
⋮
⋮
⋮
⋱
⋮
1
x
N
1
x
N
2
…
x
N
M
−
1
\begin{array}{ccc} 1 & x_{11} & x_{12} & \dots & x_{1M-1} \\ 1 & x_{21} & x_{22} & \dots & x_{2M-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{N1} & x_{N2} & \dots & x_{NM-1} \\ \end{array}
11⋮1x11x21⋮xN1x12x22⋮xN2……⋱…x1M−1x2M−1⋮xNM−1
2、每一行即一个个体的参数,继续,使用np.cumprod()
对每行进行累乘,np.cumprod(np.hstack(np.ones((N,1), X[:, :M-1]))
,得到如下矩阵(注意:每一行都是M维):
1
1
∗
x
11
1
∗
x
11
∗
x
12
…
1
∗
x
11
∗
…
x
1
M
−
1
1
1
∗
x
21
1
∗
x
11
∗
x
22
…
1
∗
x
21
∗
…
x
2
M
−
1
⋮
⋮
⋮
⋱
⋮
1
1
∗
x
N
1
1
∗
x
N
1
∗
x
N
2
…
1
∗
x
N
1
∗
…
x
N
M
−
1
\begin{array}{ccc} 1 & 1*x_{11} & 1*x_{11}*x_{12} & \dots & 1*x_{11}* \dots x_{1M-1} \\ 1 & 1*x_{21} & 1*x_{11}*x_{22} & \dots & 1*x_{21}* \dots x_{2M-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & 1*x_{N1} & 1*x_{N1}*x_{N2} & \dots & 1*x_{N1}* \dots x_{NM-1} \\ \end{array}
11⋮11∗x111∗x21⋮1∗xN11∗x11∗x121∗x11∗x22⋮1∗xN1∗xN2……⋱…1∗x11∗…x1M−11∗x21∗…x2M−1⋮1∗xN1∗…xNM−1
3、是不是每一行的每一项就是目标函数中的元素累乘的形式,只不过和
f
1
f_1
f1到
f
M
f_M
fM顺序相反,所以我们要使用np.fliplr()
将每行翻转,np.fliplr(np.cumprod(np.hstack(np.ones((N,1), X[:, :M-1])))
,结果:
1
∗
x
11
∗
…
x
1
M
−
1
…
1
∗
x
11
∗
x
12
1
∗
x
11
1
1
∗
x
21
∗
…
x
2
M
−
1
…
1
∗
x
21
∗
x
22
1
∗
x
21
1
⋮
⋱
⋮
⋮
⋮
1
∗
x
N
1
∗
…
x
N
M
−
1
…
1
∗
x
N
1
∗
x
N
2
1
∗
x
N
1
1
\begin{array}{ccc} 1*x_{11}* \dots x_{1M-1} & \dots & 1*x_{11}*x_{12} & 1*x_{11} & 1 \\ 1*x_{21}* \dots x_{2M-1} & \dots & 1*x_{21}*x_{22} & 1*x_{21} & 1 \\ \vdots & \ddots & \vdots & \vdots & \vdots \\ 1*x_{N1}* \dots x_{NM-1} & \dots & 1*x_{N1}*x_{N2} & 1*x_{N1} & 1 \\ \end{array}
1∗x11∗…x1M−11∗x21∗…x2M−1⋮1∗xN1∗…xNM−1……⋱…1∗x11∗x121∗x21∗x22⋮1∗xN1∗xN21∗x111∗x21⋮1∗xN111⋮1
4、是不是已经很明显了,已经成型了吧,后面就很简单了,先要构建
(
1
−
x
i
)
i
∈
[
1
,
2
,
…
,
M
−
1
]
(1 - x_i) \; i \in [1,2,\dots ,M-1]
(1−xi)i∈[1,2,…,M−1]部分矩阵,与上同理,使用ones
矩阵凑一项,然后需要逆序,np.hstack((ones_matrix, 1-X[:, range(M-2, -1, -1)]))
,得到矩阵:
1
1
−
x
1
M
−
1
…
1
−
x
11
1
1
−
x
2
M
−
1
…
1
−
x
21
⋮
⋮
⋱
⋮
1
1
−
x
N
M
−
1
…
1
−
x
N
1
\begin{array}{ccc} 1 & 1-x_{1M-1} & \dots & 1-x_{11} \\ 1 & 1-x_{2M-1} & \dots & 1-x_{21} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & 1-x_{NM-1} & \dots & 1-x_{N1} \\ \end{array}
11⋮11−x1M−11−x2M−1⋮1−xNM−1……⋱…1−x111−x21⋮1−xN1
5、可以发现每一行即每个个体的目标向量仍然是M维(这是必须的),接着就是对共同项和共同系数乘即可,使用np.tile()
对
(
1
+
g
(
X
M
)
)
(1+g(X_M))
(1+g(XM))进行扩展(1, M),为何是1行,因为会自动广播。
6、使用f = 0.5 * np.fliplr(np.cumprod(np.hstack((ones_matrix, X[:, :M-1])), axis=1)) * np.hstack((ones_matrix, 1-X[:, range(M-2, -1, -1)])) * np.tile(1+g, (1, M))
即可同时获得多个个体的目标函数值。
至此,通过一层层抽丝剥茧,发现矩阵计算很巧妙,其实也不算过于难,当然这只是一个例子,小试牛刀当然无法看到庐山真面目,笔者还会继续分享更多的矩阵计算思路,因为笔者也是刚开始学习。
- 对比
既然对于for
循环和矩阵计算的代码我们都写好了,不妨我们使用大一点的数据进行比较一下,看看是否和我们预想的那样孰优孰劣。
import numpy as np
import time
def DTLZ1_0(X: np.ndarray, M: int = 3) -> np.ndarray:
"""
测试函数DTLZ1
:param X: 决策向量
:param M: 目标数
:return: 目标向量
"""
# 目标数大于等于3,假设决策向量X是2维ndarray,可接受多个个体
assert X.shape[1] >= M+4, 'Dimension is invalid!'
n = X.shape[1] # 决策变量维度
objvs_list = [] # 目标值向量
for x in X: # 遍历X向量
objv_list = [] # 单个个体目标向量
XM = x[M-1:] # XM向量
for i in range(M):
fx = 0.5
for j in range(M-i-1):
fx *= x[j]
if i != 0:
fx *= (1 - x[M-i-1]) # 除了第一个目标函数均有该项
fx *= (1 + g(XM))
objv_list.append(fx)
objvs_list.append(objv_list)
objvs_list = np.array(objvs_list)
return objvs_list
def g(X: np.ndarray) -> float:
"""
计算gx函数值
:param X: n维向量
:return: 标量值
"""
lens = X.shape[0] # 获取维度
gx = 0 # 对返回值进行初始化
for i in range(lens):
gx += ((X[i] - 0.5)**2 - np.cos(20*np.pi*(X[i] - 0.5)))
gx = 100 * (lens + gx)
return gx
def DTLZ1_1(X: np.ndarray, M: int = 3) -> np.ndarray:
"""
测试函数DTLZ1
:param X: 决策向量,可以是多维,接收多个个体
:param M: 目标数
:return: 目标向量
"""
assert X.shape[1] >= M+4, 'Dimension is invalid!'
N = X.shape[0] # 获取个体个数
XM = X[:, (M-1):] # 截取矩阵中每个个体的XM部分
g = 100 * (XM[0].shape[0] + np.sum((XM - 0.5)**2 - np.cos(20*np.pi*(XM - 0.5)), axis=1, keepdims=True)) # 计算gx函数值,其中axis=1表示按列取元素,keepdims表示保持维度不变
ones_matrix = np.ones((N, 1)) # 创建(N, 1)的1矩阵
# 矩阵计算目标函数值
f = 0.5 * np.fliplr(np.cumprod(np.hstack((ones_matrix, X[:, :M-1])), axis=1)) * np.hstack((ones_matrix, 1-X[:, range(M-2, -1, -1)])) * np.tile(1+g, (1, M))
return f
np.random.seed(1)
x = np.random.uniform(0, 1, (4,10))
print("-----------x--------------\n",x)
np.random.seed(1)
y = np.random.uniform(0, 1, (4,10))
print("-----------y--------------\n",y)
print("-----------for--------------\n",DTLZ1_0(x, 3))
print("-----------matrix--------------\n",DTLZ1_1(y, 3))
我们先新建简单的决策变量,(4,10),测试两个方法输出的值是否相同,如下所示,可以看到结果相同,证明编码方式没有问题,且由于变量维度很小,时间差感觉上无差。
-----------x--------------
[[4.17022005e-01 7.20324493e-01 1.14374817e-04 3.02332573e-01
1.46755891e-01 9.23385948e-02 1.86260211e-01 3.45560727e-01
3.96767474e-01 5.38816734e-01]
[4.19194514e-01 6.85219500e-01 2.04452250e-01 8.78117436e-01
2.73875932e-02 6.70467510e-01 4.17304802e-01 5.58689828e-01
1.40386939e-01 1.98101489e-01]
[8.00744569e-01 9.68261576e-01 3.13424178e-01 6.92322616e-01
8.76389152e-01 8.94606664e-01 8.50442114e-02 3.90547832e-02
1.69830420e-01 8.78142503e-01]
[9.83468338e-02 4.21107625e-01 9.57889530e-01 5.33165285e-01
6.91877114e-01 3.15515631e-01 6.86500928e-01 8.34625672e-01
1.82882773e-02 7.50144315e-01]]
-----------y--------------
[[4.17022005e-01 7.20324493e-01 1.14374817e-04 3.02332573e-01
1.46755891e-01 9.23385948e-02 1.86260211e-01 3.45560727e-01
3.96767474e-01 5.38816734e-01]
[4.19194514e-01 6.85219500e-01 2.04452250e-01 8.78117436e-01
2.73875932e-02 6.70467510e-01 4.17304802e-01 5.58689828e-01
1.40386939e-01 1.98101489e-01]
[8.00744569e-01 9.68261576e-01 3.13424178e-01 6.92322616e-01
8.76389152e-01 8.94606664e-01 8.50442114e-02 3.90547832e-02
1.69830420e-01 8.78142503e-01]
[9.83468338e-02 4.21107625e-01 9.57889530e-01 5.33165285e-01
6.91877114e-01 3.15515631e-01 6.86500928e-01 8.34625672e-01
1.82882773e-02 7.50144315e-01]]
-----------for--------------
[[103.98274008 40.37267339 201.80237159]
[118.02765095 54.22029425 238.65424765]
[261.30214527 8.56516313 67.1531585 ]
[ 18.98340529 26.09629444 413.29499287]]
-----------matrix--------------
[[103.98274008 40.37267339 201.80237159]
[118.02765095 54.22029425 238.65424765]
[261.30214527 8.56516313 67.1531585 ]
[ 18.98340529 26.09629444 413.29499287]]
那现在就用大一点的决策变量进行计算比较,这里设置为(100,10000),M=10,当然可能不符合要求,没有关系,此处只是测试使用:
test_data = np.random.uniform(0, 1, (100, 10000))
print("-----------test_date--------------\n",test_data)
b_time = time.time()
print("-----------for--------------\n")
m = DTLZ1_0(test_data, 10)
e_time = time.time()
print("-----------for_time--------------\n","{:.5f}s".format(e_time-b_time))
b_time = time.time()
print("-----------matrix--------------\n")
n = DTLZ1_1(test_data, 10)
e_time = time.time()
print("-----------matrix_time--------------\n","{:.5f}s".format(e_time-b_time))
结果:
-----------test_date--------------
[[0.98886109 0.74816565 0.28044399 ... 0.64718918 0.0402554 0.38441391]
[0.68292767 0.69363631 0.87008649 ... 0.85004188 0.51956084 0.41380247]
[0.58685434 0.10883402 0.29372002 ... 0.84505411 0.85041685 0.23719826]
...
[0.94824474 0.23706907 0.75472375 ... 0.0654176 0.68374453 0.3275583 ]
[0.28864111 0.12874645 0.21830357 ... 0.1653693 0.69696769 0.431051 ]
[0.00305066 0.5161635 0.3510851 ... 0.83417877 0.22503453 0.75468994]]
-----------for--------------
-----------for_time--------------
22.97153s
-----------matrix--------------
-----------matrix_time--------------
0.03890s
可以看到显式for
循环的计算时间为22.97秒,而矩阵计算的时间为0.038秒,将近600倍的加速,如果在种群不断迭代的计算中,那耗时的倍数可能保持不变,但是时长的差距会拉的吓人,比如说最终矩阵计算耗时10分钟,那for
循环按600倍就需要将近100小时(当然这里有夸大的地方,不可能一直保持600倍),笔者想强调的是矩阵计算的真的具有极其迷人的高效性,而此也是以其神秘性(可读性)为代价。
人生到处知何似,应似飞鸿踏雪泥,至此笔者对于矩阵计算思维的简单理解到这里讲完了,当然这只是其中一个小小的例子,理论的有效性还需要在更多的时间和实践证明,后面笔者也会分享更多的矩阵计算例子方便大家共同学习,共同进步。
临渊羡鱼不如退而结网
创作不易,如果您觉得这篇文章对你有用,可以点个赞,算是对笔者的支持和激励!这里是Leo的博客城堡,以Python为核,ML&DL为主,泛之形形色色,输寥寥拙见,摄浮光掠影,讲三两故事。临渊羡鱼,不如退而结网,持续干货输出,有趣的灵魂值得你的关注!
原文可以去笔者的github
主页:https://github.com/LabyrinthineLeo/Yxs_Git_Learning_repos
查看(如果可以,点个star
也无妨呀,嘿嘿)。