长话短说，我将一个函数应用于多个不同的时间间隔，然后将结果数组存储在ndarray中的不同索引处。目前，我正在通过使用一个for循环来实现这一点，该循环的numpy等价于enumerate函数。据我所知，这消除了numpy的主要优势：矢量化。这是不是我的rountine的一种特殊实现方式，可以保留这一优势？
这是我的代码：
大多数是功能psi_t的工作部件
import numpy as np
# Number of Walks and Number of Positions
N = 100
P = 2*N +1
hopping_rate = 0.5
psi_t0 = np.zeros(P)
psi_t0[N] = 1
#creates the line upon which the particle moves
#index N is the central position
def hamiltonian(line_length, hopping_rate):
'''
creates the simple non time dependent hamiltonian for H = γA
where A is the adjancency matrix
'''
return hopping_rate * line_adjacency_matrix(line_length)
def measurement_operator(positions,finished_quantum_state):
'''
Converts the finished quantum state into an array of probabilities for
being in each position.
Uses the measurement operator from Susan Blog
https://susan-stepney.blogspot.com/2014/02/mathjax.html
Improved on by this guy
https://github.com/Driminary/python-randomwalk-project/blob/master/quantum-2D.py
Apart from the fact that the measurement operator drops the extra dimensions of the spin space,
which isn't present in the continuous walk.
'''
probabilities = np.empty(P)
#M_hat = np.zeros((2*P,2*P,2*P))
for k in range(P):
posn = np.zeros(P) # values of positions to nought ..
posn[k] = 1 #except for the value we're interested in
#M_hat = np.kron(np.outer(posn,posn)) #perform measurement at the current pos
M_hat = np.outer(posn,posn)
proj = M_hat.dot(finished_quantum_state) #find the state the system is in
probabilities[k] = proj.dot(proj.conjugate()).real #Calculate Prob of Particle being there
return probabilities
def psi_t(initial_wave_function,positions,hopping_rate,time):
'''
Acts upon the initial state to give the 'position' of the quantum particle at time t. Applies the measurement operator
to return the probability of being at any position at time t.
'''
psi_t = np.matmul((LA.expm(-1j*hamiltonian(positions,hopping_rate)*time)),initial_wave_function) #state after the continuous walk after time evolution
probablities = measurement_operator(P, psi_t)
return probablities
time_evolution = 150 #how many 'seconds' the wavefunction is evolved for
time_interval = 0.5
number_of_intervals =int(time_evolution / time_interval )
number_of_positions = P
probabilities_at_t =np.ndarray((number_of_intervals,number_of_positions)) #creates the empty ndarray ready for the probabilites at time t
array_of_times = np.linspace(0,time_evolution,number_of_intervals) #produces the individual times at which psi_t is calculated,
for idx,time in np.ndenumerate(array_of_times):
probabilities_at_t[idx] = psi_t(psi_t0,P,hopping_rate,time) #the array probabillites_at_t is filled at index idx with the array of probabilities produced by psi_t.
#This is the step I am trying to vectorise
函数psi_t在for循环上被调用以单独作用于array_of_times中的每个time。有没有办法psi_t可以像对数组x执行x**2一样对数组array_of_times执行操作？这能一蹴而就吗？
p.S Eagle Eyed Overflowers会注意到，无论如何，在measurement_operator中都有一个for循环。然而，我认为没有办法摆脱这种情况！
这个问题实际上是不可重复的，因为一些被调用的函数丢失了，但这里是我对measurement_operator的矢量化实现。这是假设finished_quantum_state的形状为(P, )(不确定是否是这样，因为直到那个部分才能复制(。
def measurement_operator_vectorized(positions, finished_quantum_state):
M_hat = np.zeros((P, P, P))
M_hat[np.arange(0, P), np.arange(0, P), np.arange(0, P)] = 1
proj = np.tensordot(M_hat, finished_quantum_state, axes=((2), (0)))
probabilities = (proj * proj.conjugate()).sum(axis=1).real
return probabilities
以下是一些基准标记-
P = 1000
a = np.random.rand(P)
b = np.random.rand(P)
%timeit c1 = measurement_operator(a, b)
%timeit c2 = measurement_operator_vectorized(a, b)
%memit c1 = measurement_operator(a, b)
%memit c2 = measurement_operator_vectorized(a, b)
print(np.allclose(c1, c2))
给予-
1.18 s ± 46.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
308 ms ± 6.76 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
peak memory: 86.43 MiB, increment: 0.00 MiB
peak memory: 90.34 MiB, increment: 3.91 MiB
True
矢量化版本比p~1000的内存使用速度更快。请注意，对于非常高的P值，矢量化版本的内存使用量将大大增加。
这并不是OP所要求的，但要将另一个循环矢量化，一个更完整的代码会很有帮助。
但是，只有当finished_quantum_state是真实的时，这个基准才有效。对于复杂的值，tensordot操作非常缓慢且效率低下(在内存中(，因此使用非矢量化版本可能会更好。
P = 1000
a = np.random.rand(P) + np.random.rand(P)*1j
b = np.random.rand(P) + np.random.rand(P)*1j
%timeit -n1 -r1 c1 = measurement_operator(a, b)
%timeit -n1 -r1 c2 = measurement_operator_vectorized(a, b)
%memit c1 = measurement_operator(a, b)
%memit c2 = measurement_operator_vectorized(a, b)
np.allclose(c1, c2)
2.97 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)
3.49 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)
peak memory: 102.69 MiB, increment: 0.03 MiB
peak memory: 15365.38 MiB, increment: 15262.69 MiB
然而，如果你真的想要最好的性能，你最好暂时忘记测量等物理细节，只做
def measurement_operator_fastest(positions, finished_quantum_state):
return (finished_quantum_state * finished_quantum_state.conjugate()).real
P = 1000
a = np.random.rand(P) + np.random.rand(P)*1j
b = np.random.rand(P) + np.random.rand(P)*1j
%timeit -n1 -r1 c1 = measurement_operator(a, b)
%timeit -n1 -r1 c2 = measurement_operator_vectorized(a, b)
%timeit -n1 -r1 c3 = measurement_operator_fastest(a, b)
%memit c1 = measurement_operator(a, b)
%memit c2 = measurement_operator_vectorized(a, b)
%memit c3 = measurement_operator_fastest(a, b)
print(np.allclose(c1, c2))
print(np.allclose(c1, c3))
2.87 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)
3.48 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)
16.6 µs ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)
peak memory: 102.70 MiB, increment: 0.00 MiB
peak memory: 15365.39 MiB, increment: 15262.69 MiB
peak memory: 102.69 MiB, increment: -0.01 MiB
True
True
直接取内积，可以使函数的速度提高10^6倍左右。当然，这是假设测量操作员是定义的。