根据梯度方程,矩阵乘法由
给出
其中@
和*
都需要。如果读者感兴趣,请阅读以下代码:
# parameters
beta = 0.98
alpha = 0.03
delta = 0.1
T = 1000
loop = 1
dif = 1
tol = 1e-8
kss = ((1 / beta - (1 - delta)) / alpha)**(1 / (alpha - 1))
k = np.linspace(0.5 * kss, 1.8 * kss, T)
k_reshaped = k.reshape(-1, 1)
c = k_reshaped ** alpha + (1 - delta) * k_reshaped - k
c[c<0] = 1e-11
c = np.log(c)
beta_square = beta**2
# multiplication
I = np.identity(T)
E = np.ones(T)[:,None]
Q2 = I
while np.any(dif > tol) and loop < 200:
J = beta * Q2
B = inv(I - J)
Q3 = np.zeros([T,T])
ini = np.argmax(c + (B @ (J * c) @ E).flatten(),axis=1)
Q3[np.arange(T),ini] = 1
gB = 2 * B @ (J * c @ E) @ (beta * Q2 * c @ E + B @ (np.linalg.matrix_power(I - J, 2) * c @ E)).T / beta_square
B += 0.1 * gB
dif = np.max(np.absolute(Q3 - Q2))
kcQ = k[ini]
Q2 = Q3
loop += 1
基本上,它遵循随机梯度下降算法,矩阵B
由B = inv(I - J)
初始化并由B += 0.1 * gB
演化,J
沿着Q2
变化,每次迭代都需要确定Q2
。然而Q2
是稀疏矩阵,每列只有一个1,其余为零,在代码中如下:
ini = np.argmax(c + (B @ (J * c) @ E).flatten(),axis=1)
Q3[np.arange(T),ini] = 1
...
Q2 = Q3
目前的代码演示了一个1000 × 1000的矩阵运算,这是否可以优化并运行得更好?
1条答案
按热度按时间svmlkihl1#
以下是一些改进:
beta * Q2
计算两次,第二次可以使用J
。J * c
也可以多次计算,而I - J
可以一次计算。B @ (J * c) @ E
和B @ (J * c @ E)
在数学上是等价的,但在您的情况下,后者更快,也可以计算一次。0.1 * (2 * Matrix / beta_square)
实际上计算了一个新矩阵M2 = 2 * Matrix
,然后是一个新矩阵M3 = M2 / beta_square
,然后是另一个M4 = 0.1 * M4
。创建许多这样的临时矩阵是昂贵的,因为它是一个内存限制的操作,并且内存带宽在现代机器上非常有限(与计算能力相比),更不用说填充新的临时数组通常比已经填充的要慢(因为虚拟内存,更具体地说是页面错误)。因此,执行(0.1 * 2 / beta_square) * Matrix
更快(因为乘以float比乘以大矩阵快得多)。np.argmax(c + tmp3.flatten(), axis=1)
或np.max(np.absolute(Q3 - Q2))
。out
参数(例如np.multiply(A, B, out=C)
),您可以使用它们。也就是说,这里的好处很小,因为inv
花费了大量时间。B
,您可以使用np.linalg.solve
代替,如Homer 512所述。对于大型矩阵,求解系统的速度明显更快(O(n**3)
与O(n**2)
),并且通常更准确。参见Don’t invert that matrix。例如,inv(I-J) @ b
可以被solve(I-J, b)
替换。使用solve
的好处并不是很大,尽管在您的特定用例中,因为I-J
矩阵是稀疏的。B
,那么这就有点复杂了。Numba可以帮助编写一个相对快速的矩阵求逆,专门针对稀疏矩阵,就像你的用例一样(因为Scipy的速度很慢)。np.linalg.matrix_power(tmp0, 2) * c
也可以在Numba中针对稀疏矩阵进行优化。下面是一个使用Numba的实现(几乎没有经过测试):
请注意,
invert_sparse_matrix
函数仍然需要很长时间(在我的机器上〉50%),但它比inv
快3倍,与solve
一样快。它是naive inversion algorithm的改进,对非常稀疏的矩阵进行了一些优化。它当然可以进一步优化(例如,使用平铺)但这肯定不是微不足道的事情(特别是对于新手程序员)。注意编译时间需要几秒钟。
总的来说,这个实现比我的i5- 9600 KF处理器(6核)上的初始实现快了大约4~5倍。