Optimise/replace this for loop used in calculation of mutual information

1 次查看（过去 30 天）
Ansh
Ansh 2018-1-16
评论： Greg 2018-1-20
So I previously posted on here with regards to some code given HERE. Thanks to the help of Jan Simons I have managed to eliminate one of the nested for loops to speed up the code pretty significantly. However, due to the dimensionality of the project and future projects I am working on I still need to optimise the code. I ran the profiler and the following block of code seems to be bottleneck:
Eps = zeros(nObs, 1);
Nn = zeros(nObs, 1);
nx1 = zeros(nObs, 1);
ny1 = zeros(nObs, 1);
nx2 = zeros(nObs, 1);
ny2 = zeros(nObs, 1);
for i = 1:nObs
dxSample = dx(i, :);
dxSample(i) = [];
dySample = dy(i, :);
dySample(i) = [];
dzSample = dz(i, :);
dzSample(i) = [];
[EpsSample, NnSample] = sort(dzSample, 'ascend');
Eps(i) = EpsSample(k);
Nn(i) = NnSample(k);
nx1(i) = sum(dxSample < Eps(i));
ny1(i) = sum(dySample < Eps(i));
nx2(i) = sum(dxSample <= Eps(i));
ny2(i) = sum(dySample <= Eps(i));
end
Note that the k is typically quite low: between 3-10, nObs is around 4364. Is there any way of speeding this up? Maybe even eliminating the for loop altogether if possible? Any help is very much appreciated!
  5 个评论
显示 3更早的评论
Walter Roberson
Walter Roberson 2018-1-16
I have seen several reports that implicit expansion is often slower than bsxfun. The details are likely to vary with release.
Ansh
Ansh 2018-1-17
Hi Greg, not sure if Jan has probably answered your question already but k is an input into the original function in the link in the original question. dx,dy and dz are already obtained from a different part of the code that Jan thankfully helped me to optimise with previously.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Greg
Greg 2018-1-16
编辑：Greg 2018-1-18
Taking some guesses:
eyeinds = eye(nObs,'logical');
dx2(eyeinds) = Inf;
dy2(eyeinds) = Inf;
dz2(eyeinds) = Inf;
% Edit to incorporate Jan's suggestions:
[EpsSample2, NnSample2] = mink(dz2,k,2);
Eps2 = EpsSample2(:,end);
Nn2 = NnSample2(:,end);
% [EpsSample2, NnSample2] = sort(dz2,2,'ascend');
% Eps2 = EpsSample2(:,k);
% Nn2 = NnSample2(:,k);
nx12 = sum(dx2 < Eps2,2);
ny12 = sum(dy2 < Eps2,2);
nx22 = sum(dx2 <= Eps2,2);
ny22 = sum(dy2 <= Eps2,2);
  8 个评论
显示 6更早的评论
Ansh
Ansh 2018-1-20
Yes the updated answer is now working and is giving about 30% speed up from the previous modification. I'm very happy with the speed up of the code overall and am happy to accept this answer. Thanks Greg and Jan (again)!

请先登录，再进行评论。

更多回答（1 个）

Jan
Jan 2018-1-16
Sorting the complete array is a waste of time if all you need is the k.th smallest element. Under Matlab 2017b see mink to solve this much faster.
A further idea: Omit setting the diagonal to Inf, but leave it at zero. Then the result of nx12 etc. is 1 to large.
[Eps2, Nn2] = mink(dz2, 2, k);
nx12 = sum(dx2 < Eps2, 2) - 1;
ny12 = sum(dy2 < Eps2, 2) - 1;
nx22 = sum(dx2 <= Eps2, 2) - 1;
ny22 = sum(dy2 <= Eps2, 2) - 1;
  2 个评论
Greg
Greg 2018-1-16
编辑：Greg 2018-1-17
Leaving at zero only works if all data is guaranteed positive. And that throws off the min counting in dz. I didn't want to make that assumption.
Just noticed your version leaves Eps2 as a non-vector for k > 1. I've edited my answer to incorporate the awesome suggestion to use mink.
Ansh
Ansh 2018-1-17
Thanks for the input again Jan and to Greg too! Unfortunately I don't have MATLAB 2017b - I'm using MATLAB 2016b currently so I'll download the newest version (might as well I guess) and re test to see how much faster it is (hopefully the speed is up is significant).
And yes Greg luckily all the data is positive since dx, dy and dz are essentially distances.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Reference Applications 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by