本文使用 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) > 0error("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 == 1D(1:end) = log(1 ./ p0);elseindices0 = (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為等長度的arrayD = 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))); elseerror("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 : Mresults((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 : Mexp_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多倍，爽爆了。

總結

以上是生活随笔為你收集整理的matlab repmat函数_Matlab向量化编程在二级劝退学科中的一个应用例子的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。