matlab repmat函数_Matlab向量化编程在二级劝退学科中的一个应用例子-CSDN博客

本文使用 Zhihu On VSCode 创作并发布

向量化编程对于用过Matlab或者numpy做过稍微复杂一点的数值运算的人来说并不陌生，甚至说刚入门Matlab或者numpy的小白都知道要用z = x .* y来代替一个for k = 1 : N {z(k) = x(k) * y(k);}，进阶一点的人也明白迭代算法最好把更新函数写成向量化的形式。以前觉得向量化编程就是个小技巧，对着入门教程中那些梯度下降例子无脑抄几遍就能掌握向量化编程的主要精髓，直到前两天因为要对一个二级劝退学科的实验数据后处理程序向量化，才体会到要掌握向量化编程思想还是得动脑子。。。

实验数据后处理需求

我有在

个settings

下分别测得的数据

，其中每组数据

是一个随机离散时间序列

，其中

。然而理论合作者那边要求获得的时间序列要满足每个数据点所对应的测量setting必须是从那

个settings中按照给定概率分布随机抽取的，那我必须生成一条随机setting序列

，其中

，然后用我那测得的

组数据按照这个setting序列生成一条新的随机序列：也就是说假设这个setting序列里面有

个是

，那么我就用数据

中的前

个点往对应的位置填充，以此类推，那么生成的随机序列

就满足了理论合作者的要求。

在获得了这个随机序列

以后，还需要用它来获得两条新序列

和

用来实验结果的曲线图。第一条序列是这样的：对于第

个点，先求得前

个数据点的平均值

，然后求

，其中

是KL散度函数，

是一个固定不变的值。第二条序列是这样的：对于第

个点，同样先求得前

个数据点的平均值

，然后求

其中

是KL散度

在固定

后的逆函数

（这种条件下的逆函数是存在的，因为单调，只不过需要用数值方法求解）。

这还不行，单次实验结果方差太大，需要做

次

实验，也就是说上面过程重复个

次。

向量化过程

这个实验数据后处理程序初始版本还真的是严格按照上面过程编写的：每次产生一个随机setting，然后从原始数据抽取数据点填充，计算

，再计算

和

。当时这个程序写出来是为了仿真一个低维高保真度的情况，

不需要太大（5000左右吧，

取50），跑的时间还是可以接受的（不到一分钟吧），毕竟就画个图给老板看。结果实验数据维度高了保真度低了，

得取到50000左右，结果运行时间直接飙升到十几分钟，但我还要对着跑出来图debug，这来回几下也太费时间了。

我对源程序向量化的过程主要用到的技巧是数组的逻辑寻址，然后对满足不同逻辑的数据点分别进行不同的运算。首先生成指定分布的随机序列是不需要一个一个来的，假设我在

个settings上的概率分布为

，那这个随机序列生成函数可以这么写：

function randseq = getRandSeq(len, probs)
cutpoints = cumsum(probs);	% 给不同结果划分区间
randseq = rand(1, len);
idx = randseq < cutpoints(1);	% 找到落在第一个区间的所有点的下标
randseq(idx) = 2;	% 给第一个区间上的点赋予对应的值
for k = 2 : max(size(probs))	% 对后面K-1个区间重复上面过程
    idx = randseq < cutpoints(k) & randseq >= cutpoints(k-1);
    randseq(idx) = k+1; 
end
randseq = randseq - 1;
end

随机setting序列生成之后，用原始数据填充过程也是可以向量化的：

results = zeros(1, N);       % 给生成的随机序列预分配空间
% 生成随机setting序列
setting_seq = getRandSeq(N, probs);
% 用原始数据填充setting序列
for k = 1 : setting_num    
	idx = (setting_seq == k);	% 在setting序列中找出所有S_k对应的下标
	results(idx) = raw_data{k}(1:sum(idx));	% 用对应原始数据往对应位置填充     
end

计算

序列同样可以向量化：

results = cumsum(results) ./ (1:N);	% 直接在原地更新，不新建数组了

计算KL散度的函数

也给它向量化了：

function D = klDiv(p1, p0)
% 计算两个伯努利分布之间的KL散度
% p1和p2要么是两个长度相等的array，要么其中一个是数，另一个是array

% 检查输入的范围
if sum(p0 < 0) > 0 || sum(p0 > 1) > 0 || sum(p1 > 1) > 0 || sum(p1 < 0) > 0
    error("illegal input.")
end
if max(size(p0)) == 1	% p1为array, p0为数值
    D = zeros(size(p1));	% 预分配内存
    indices0 = (p1 <= p0);	% 特殊情况0点下标
    indices1 = (p1 == 1);	% 特殊情况1点下标
    indices2 = (p1 > p0) & (p1 < 1);	% 正常情况点下标
    D(indices0) = 0;
    D(indices1) = log(1 / p0);
    D(indices2) = p1(indices2) .* log(p1(indices2) / p0) +...
                            (1 - p1(indices2)) .* log((1 - p1(indices2)) / (1 - p0));
elseif max(size(p1)) == 1	% p0为array, p1为数值
    D = zeros(size(p0));	% 预分配内存
    if p1 == 1
        D(1:end) = log(1 ./ p0);
    else
        indices0 = (p0 >= p1);	% 特殊情况0点下标
        indices1 = (p0 < p1);	% 正常情况点下标
        D(indices0) = 0;
        D(indices1) = p1 * log(p1 ./ p0) +...
                            (1 - p1) * log((1 - p1) ./ (1 - p0));
    end
elseif max(size(p1)) == max(size(p0))	% p0和p1为等长度的array
    D = zeros(size(p1));	% 预分配内存
    indices0 = (p1 <= p0);	% 特殊情况0点下标
    indices1 = (p1 == 1);	% 特殊情况1点下标
    indices2 = (p1 > p0) & (p1 < 1);	% 正常情况点下标
    D(indices0) = 0;
    D(indices1) = log(1 ./ p0(indices1));
    D(indices2) = p1(indices2) .* log(p1(indices2) ./ p0(indices2)) +...
                            (1 - p1(indices2)) .* log((1 - p1(indices2)) ./ (1 - p0(indices3)));
else
    error("illegal input.")
end
end

那么计算

序列就可以一口气完成了：

Delta = exp(-klDiv(results, p0) .* (1:N));

最头疼的是

序列的计算，因为每个点都需要用数值方法求解（我这用的是二分法，因为函数是单调的），当时我的思路是能不能想出什么近似方法把KL散度的逆函数显式表达出来，这样就可以向量化计算了，结果没想出来。。。后来就直接把二分法这个过程给向量化了：

function x = invKLDiv(ps, y)
threshold = 1e-4;	% 计算精度
left = zeros(size(ps));	
right = ps;
x = (left + right) / 2;
res = klDiv(ps, x) - y;
idx = (right - left) > threshold;	% 找出精度仍然不达标的所有下标
idx1 = (res > 0) & idx;	% 找出idx中需要更新区间左边界的所有下标
idx2 = (res <= 0) & idx;	% % 找出idx中需要更新区间右边界的所有下标
while sum(idx) > 0	% 循环直到所有点都达到精度要求
    left(idx1) = x(idx1);
    right(idx2) = x(idx2);
    x(idx) = (left(idx) + right(idx)) / 2;
    res(idx) = klDiv(ps(idx), x(idx)) - y(idx);
    idx = (right - left) > threshold;
    idx1 = (res > 0) & idx;
    idx2 = (res <= 0) & idx;
end
end

那现在

序列也可以一口气计算了：

eps_A = 16 * (1 - invKLDiv(results, -log(delta)./ (1:N))) / 3;

最后就是把那

次

实验数据后处理给向量化了，很简单，就是把上面生成的随机setting序列的长度从

拓展到

就是了，最后整个过程长这样：

% start post-processing the raw data
one_to_N = 1 : N;
one_to_N_M = repmat(one_to_N, 1, M);
results = zeros(1, N*M);       % test results sequence
% generate a random setting sequence
setting_seq = getRandSeq(N*M, probs);
% fill the results sequence by raw data according to setting_seq
for k = 1 : setting_num    
	idx = (setting_seq == k);
	results(idx) = raw_data{k}(1:sum(idx));       
end
% compute the significance and infidelity sequence
for k = 1 : M
	results((k-1)*N+1 : k*N) = cumsum(results((k-1)*N+1 : k*N)) ./ one_to_N;
end
Delta = exp(-klDiv(results, p0) .* one_to_N_M);
eps_A = d * (1 - invKLDiv(results, -log(delta) ./ one_to_N_M)) / ((d + 1) * nu);
% store the M i.i.d. exp data
for k = 1 : M
    exp_data{4*(k - 1) + 1} = results(k*N);
    exp_data{4*(k - 1) + 2} = Delta((k-1)*N+1 : k*N);
    exp_data{4*(k - 1) + 3} = eps_A((k-1)*N+1 : k*N);
    exp_data{4*(k - 1) + 4} = results((k-1)*N+1 : k*N);
end

整个运行时常缩短为11秒，比之前的十几分钟快了60多倍，爽爆了。