【发布时间】：2021-08-04 13:46:24
【问题描述】：
问题：如何在一个实值的一维数组中找到（全局）最大值的索引，该数组在循环迭代期间以（相对）少量的已知值进行部分更新，但在 Python 中变化的索引最有效？特别是，如果一维数组的大小增加而更新值的数量保持有限，我如何保持每次迭代的计算时间不变？
上下文：我正在以一种有效的变体 (Clark, 1980) 实现 CLEAN 反卷积算法 (Högbom, 1974; Roberts+, 1987)，以去除功率谱中的频谱旁瓣（别名）作为频率和方位角阶 m 的函数。该算法使用两个嵌套循环，主要（外部）和次要（内部）循环。在次要迭代期间，我需要找到一维功率谱密度（PSD）数组psd_check_temp（复值傅里叶变换数组fu_check_temp的绝对平方）的最大值的索引ind_sel_max，即在索引ind_within_beam 的小迭代之间部分更新。索引ind_within_beam 的值和元素数量在次要迭代之间可能会有所不同，但与len(psd_check_temp) 相比，在给定足够多的迭代的情况下，最大元素数量会变小。示例性地，len(ind_check_temp) 最大可以为 ~2500，len(psd_check_temp) 最大可以增长到 >~300000。 len(psd_check_temp) 仅在主要迭代之间更改，但在次要迭代之间不更改。在给定的迭代过程中，我不知道psd_check_temp 的值究竟会如何变化（通常但并不总是它们会减小），但它可能会在下一次迭代中计算出来。但是，ind_sel_max 始终包含在 ind_with_beam 中，psd_check_temp[ind_sel_max] 将在迭代期间更新。
循环结构：
i_major = 0
while continue_major_loop:
i_minor = 0
# code of outer loop
...
while continue_minor_loop:
# code of inner loop
...
i_minor += 1
# code of outer loop
...
i_major += 1
代码示例 1（内部循环，慢）：
# compute the PSD
if i_minor == 0:
psd_check_temp = np.real(fu_check_temp)**2 + np.imag(fu_check_temp)**2
else:
psd_check_temp[ind_within_beam] = np.real(fu_check_temp[ind_within_beam])**2 + np.imag(fu_check_temp[ind_within_beam])**2
# find the index of the maximum PSD
ind_sel_max = np.argmax(psd_check_temp)
...
ind_within_beam = ...
fu_check_temp[ind_within_beam] -= ...
...
问题 1： 在每次迭代中通过 np.argmax() 查找最大值的索引（请参见上面的代码）将使用越来越多的计算时间，因为 len(psd_check_temp) 在每次主要迭代中变得更大，因此经过足够多的迭代后变得相当慢。
问题：由于在迭代过程中只有一小部分 PSD 数组被更新，直觉上我猜应该有一种方法可以将每次迭代的计算时间限制在被更新（它本身如上面所写的示例性地限制为 ~2500），而不是检查整个数组 psd_check_temp 的最大值。我怎样才能最好地做到这一点？
我的想法是在第一次小迭代期间对psd_check_temp 进行一次排序，并且在每次迭代中仅对更新后的 PSD 值进行排序，保留未更新的 PSD 值，找到将排序后的更新 PSD 值插入排序后的索引通过np.searchsorted() 插入未更新的 PSD 值，并通过np.insert() 插入排序后的更新 PSD 值。由于psd_check_temp 将保持排序，因此最大 PSD 的索引是最后一个。但是，我还需要相应地处理fu_check_temp 和几个一维（inu_check_temp、im_check_temp）和二维（ind_2d）索引数组。
代码示例 2（甚至更慢）：
...
# create a 2D index array, filled with index -1 (marking invalid points)
ind_2d = np.full((nnu_use,nm), fill_value=-1, dtype=int)
...
while continue_major_loop:
...
# get a 1D array of indices corresponding to the 1D frequency and m index arrays
ind_1d = np.arange(len(inu_check_temp))
# fill the index array with the indices corresponding to the 1D frequency and m index arrays
ind_2d[inu_check_temp, im_check_temp] = ind_1d
...
while continue_minor_loop:
if i_minor == 0:
psd_check_temp = np.real(fu_check_temp)**2 + np.imag(fu_check_temp)**2
# find the indices that sort the PSD array and sort the PSD and Fourier transform arrays as well
# as the frequency and m index arrays accordingly
ind_sort_temp = np.argsort(psd_check_temp)
psd_check_temp = psd_check_temp[ind_sort_temp]
fu_check_temp = fu_check_temp[ind_sort_temp]
inu_check_temp = inu_check_temp[ind_sort_temp]
im_check_temp = im_check_temp[ind_sort_temp]
# also update the 2D index array accordingly
ind_2d[inu_check_temp, im_check_temp] = ind_1d
else:
psd_updated_temp = np.real(fu_check_temp[ind_within_beam])**2 + np.imag(fu_check_temp[ind_within_beam])**2
# (1) get the indices to sort the PSD of the updated Fourier transform
ind_sort_temp = np.argsort(psd_updated_temp)
# (2) sort the PSD, Fourier transform, and frequency and m index arrays accordingly
psd_updated_temp = psd_updated_temp[ind_sort_temp]
fu_updated_temp = fu_check_temp[ind_within_beam][ind_sort_temp]
inu_updated_temp = inu_check_temp[ind_within_beam][ind_sort_temp]
im_updated_temp = im_check_temp[ind_within_beam][ind_sort_temp]
# (3) get the 1D indices of the array elements that have not been updated; np.setdiff1d() only returns
# unique elements, but that is ok for index arrays where there are no duplicates; also it returns
# the indices in order
ind_complement_temp = np.setdiff1d(ind_1d, ind_1d[ind_within_beam], assume_unique=True)
# (4) keep only the values that were not updated
psd_cut_temp = psd_check_temp[ind_complement_temp]
fu_cut_temp = fu_check_temp[ind_complement_temp]
inu_cut_temp = inu_check_temp[ind_complement_temp]
im_cut_temp = im_check_temp[ind_complement_temp]
# (5) find the indices where to insert the updated values in the sorted arrays
ind_insert_temp = np.searchsorted(psd_cut_temp, psd_updated_temp)
# (6) insert the updated values in the sorted arrays
psd_check_temp = np.insert(psd_cut_temp,ind_insert_temp,psd_updated_temp)
fu_check_temp = np.insert(fu_cut_temp,ind_insert_temp,fu_updated_temp)
inu_check_temp = np.insert(inu_cut_temp,ind_insert_temp,inu_updated_temp)
im_check_temp = np.insert(im_cut_temp,ind_insert_temp,im_updated_temp)
# (7) also update the 2D index array accordingly
ind_2d[inu_check_temp, im_check_temp] = ind_1d
# the index of the maximum corresponds to the last array element, since the array is sorted
ind_sel_max = len(psd_check_temp) - 1
...
ind_within_beam = ...
fu_check_temp[ind_within_beam] -= ...
...
问题 2： 这种实现甚至比代码示例 1 中的实现还要慢。对于 10000 次迭代，该算法对于代码示例 1 需要大约 0.2 秒，对于代码示例 2 需要 12 秒，而我需要 10^6 次或更多次迭代。步骤np.insert() (6)、np.setdiff1d() (3)、ind_2d update (7) 的计算时间，保持未更新的值 (4)、np.searchsorted() (5)、np.argsort(psd_check_temp)(1 ) 并对更新后的值 (2) 进行排序，分别为 4.5、2.7、2.0、1.2、0.6、0.5 和 0.2 秒。此外，此实现会多次请求内存。
非常感谢任何建议/帮助！
【问题讨论】：
-
你能展示修改
psd_check_temp数组大小的代码吗？有psd_check_temp[ind_within_beam] = ...但这不会改变数组的大小，所以我想知道为什么你需要处理这个大小改变问题？ -
@a_guest:
psd_check_temp增加了大小，因为fu_check_temp确实如此。fu_check_temp在主循环中通过fu_check_temp = fu_resid[:nnu_use][mask_check_temp]获得。mask_check_temp是 PSD 大于某个 PSD 阈值（在每次主要迭代中发生变化）的元素的布尔掩码，我用它来索引二维数组fu_resid。随着 PSD 阈值的变化（通常会变小），mask_check_temp中True元素的数量会增加，布尔掩码索引产生的一维数组fu_check_temp的大小也会增加。 -
@bproxauf 你不能保存你的最大值和它的索引，每次你更新数组时你也更新保存的最大值吗？这样做的好处是您只需查看所有内容一次，因为您可以将传入的条目与保存的最大值进行比较。
-
@user2640045：在给定的迭代期间，找到的最大值，即
psd_check_temp[ind_sel_max]将被更改，因为ind_sel_max始终包含在ind_with_beam中（我编辑了问题以澄清这一点）。当我们将更新后的值与之前保存/存储的最大值进行比较并且更新后的值更小（通常会如此）时，我们仍然不知道psd_check_temp的新最大值是在ind_with_beam内还是在未更新的值。
标签： python arrays numpy indexing max