将无法向量化的“for”循环转换为稀疏矩阵答案

【问题标题】：Converting `for` loop that can't be vectorized to sparse matrix将无法向量化的“for”循环转换为稀疏矩阵
【发布时间】：2018-05-10 14:25:12
【问题描述】：

有 2 个盒子和一个小间隙，每秒允许 1 个粒子从一个盒子进入另一个盒子。一个粒子是从 A 到 B，还是从 B 到 A 取决于 Pa/Ptot 的比值（Pa：A 盒中的粒子数，Ptot：两个盒子中的总粒子数）。

为了让它更快，I need to get rid of the for loops，但是我找不到将它们向量化或将它们转换为代表我的for 循环的稀疏矩阵的方法：

对于无法矢量化的 for 循环怎么办？迭代 n 的结果取决于您在迭代 n-1、n-2 等中计算的结果。您可以定义一个代表您的 for 循环的稀疏矩阵，然后进行稀疏矩阵求解。

但我不知道如何从中定义稀疏矩阵。模拟归结为计算：

在哪里

是在尝试表达我的问题时给我带来麻烦的部分，如here 所述。（注意：括号内的内容是布尔运算）

问题：

我可以矢量化for 循环吗？
如果不是，如何定义稀疏矩阵？
（额外问题）为什么 Python 中的执行时间（0.027 秒）比 Octave（0.75 秒）快 x27？

_{注意：我在 Python 和 Octave 中都实现了模拟，很快就会在 Matlab 上进行，因此标签是正确的。}

八度码

1; % starting with `function` causes errors

function arr = Px_simulation (Pa_init, Ptot, t_arr)
  t_size = size(t_arr);
  arr = zeros(t_size);   % fixed size array is better than arr = [] 
  rand_arr = rand(t_size);  % create all rand values at once
  _Pa = Pa_init;
  for _j=t_arr()
    if (rand_arr(_j) * Ptot > _Pa)
      _Pa += 1;
    else
      _Pa -= 1;
    endif
    arr(_j) = _Pa;
  endfor
endfunction


t = 1:10^5;

for _i=1:3
  Ptot = 100*10^_i;
  tic()
  Pa_simulation = Px_simulation(Ptot, Ptot, t);
  toc()
  subplot(2,2,_i);
  plot(t, Pa_simulation, "-2;simulation;")
  title(strcat("{P}_{a0}=", num2str(Ptot), ',P=', num2str(Ptot)))
endfor

Python

import numpy
import matplotlib.pyplot as plt
import timeit
import cpuinfo

from random import random

print('\nCPU: {}'.format(cpuinfo.get_cpu_info()['brand']))


PARTICLES_COUNT_LST = [1000, 10000, 100000]
DURATION = 10**5

t_vals = numpy.linspace(0, DURATION, DURATION)


def simulation(na_initial, ntotal, tvals):
    shape = numpy.shape(tvals)
    arr = numpy.zeros(shape)
    na_current = na_initial

    for i in range(len(tvals)):
        if random() > (na_current/ntotal):
            na_current += 1
        else:
            na_current -= 1
        arr[i] = na_current
    return arr


plot_lst = []
for i in PARTICLES_COUNT_LST:
    start_t = timeit.default_timer()
    n_a_simulation = simulation(na_initial=i, ntotal=i, tvals=t_vals)
    execution_time = (timeit.default_timer() - start_t)
    print('Execution time: {:.6}'.format(execution_time))
    plot_lst.append(n_a_simulation)


for i in range(len(PARTICLES_COUNT_LST)):
    plt.subplot('22{}'.format(i))
    plt.plot(t_vals, plot_lst[i], 'r')

    plt.grid(linestyle='dotted')
    plt.xlabel("time [s]")
    plt.ylabel("Particles in box A")

plt.show()

【问题讨论】：

你能添加一个小的可重现数据集和你想要的数据集吗？
我认为 Octave 链接中的稀疏矩阵建议是制作一个线性方程，M*x=f，其中M 是一个稀疏矩阵。有限元和有限差分使用此求解 PDE。稀疏矩阵被广泛用于表示此类线性方程的刚度矩阵。
@MaxU（我是向量化和稀疏矩阵的新手，所以这听起来可能很傻......但是）你指的是什么数据集？我不应该先把它变成像here 这样的Ax=B 吗？我不知道如何创建矩阵A，所以我无法提供数据集。
@hpaulj 这是否意味着除非我可以从p_n 的公式中创建一个线性方程，否则我不能对给定问题使用稀疏矩阵？

标签： python matlab numpy octave sparse-matrix

【解决方案1】：

IIUC 你可以在Octave 和Numpy 中使用cumsum()：

八度：

>> p = rand(1, 5);
>> r = rand(1, 5);
>> p
p =

   0.43804   0.37906   0.18445   0.88555   0.58913

>> r
r =

   0.70735   0.41619   0.37457   0.72841   0.27605

>> cumsum (2*(p<(r+0.03)) - 1)
ans =

   1   2   3   2   1

>> (2*(p<(r+0.03)) - 1)
ans =

   1   1   1  -1  -1

还要注意以下函数将返回值（[-1, 1]）：

【讨论】：