Reputation: 332
I've been trying to vectorize some operations from code a bachelor student wrote for us. He used ugly nested loops which I tried to change into bsxfun calls. This works when I use this function as is, but when I try to mex compile it using coder, I get an error.
Edit: I realize I should clarify some parts of the code. S, W, and N are nxn matrices that indicate preference of a pairwise comparison ij. mu is a row vector containing the (predicted) scores of the absolute rating reconstruction.
To check what the code does, you might wanna check the very bottom of the page.
The original code reads as follows:
for i = 1:n
for j = i:n
if i~= j
mu = x(i) - x(j);
md1 = normcdf(-delta1, mu); %minus d0
md0 = normcdf(-delta0, mu);
d0 = normcdf( delta0, mu);
d1 = normcdf( delta1, mu);
Y = Y + S(i,j) .* log(1 - d1) + ...
W(i,j) .* log(d1 - d0) + ...
N(i,j) .* log(d0 - md0) + ...
W(j,i) .* log(md0 - md1) + ...
S(j,i) .* log(md1);
end
end
end
Y = -Y; %this will be used in fminsearch, therefore the negative sign
I wrote a quick and dirty vectorized solution for this:
diffs = bsxfun(@minus,x,x');
% create the deltas array augmented by -Inf and Inf for easy diff calculation
deltass = sort([-Inf, Inf, deltas, -1*deltas]);
deltass = reshape(deltass, [1, 1, length(deltass)]);
normcdfs = diff(1-bsxfun(@normcdfupper, diffs, deltass), 1, 3);
normcdfs(repmat(logical(tril(ones(size(normcdfs(:,:,1))))),[1,1,length(deltass)-1])) = 1;
Y = -sum(sum(sum(normcdfs_input.*log(normcdfs))));
Edit: normcdfupper is a wrapper class I wrote that calls normcdf(diff, delta, 'upper') since I was running into the problem, that tail ends could cause an error.
Now like I said, this works fine in my day to day use. But since mexing the old nested for loop version brought an immense speedup, I tried doing the same with this sped up version, only to encounter the following error in coder:
Expansion is only supported along dimensions where one input argument or the other has a fixed length of 1.
Error in lik_ml (line 27)
normcdfs = diff(1-bsxfun(@normcdfupper,diffs,deltass),1,3);
diffs is a :infx:inf double, and deltass is a 1x1x:inf double. Coder tells me so:
If anyone could help me out with this error message, I'd greatly appreciate it.
Edit: Some code to test it on:
clear all;
% rng(322);
delta0 = 1;
delta1 = 2;
deltas = [delta0, delta1];
sigma = 0.5;
options = 5;
maxVotes = 10000;
voteStep = 3;
initialVoteStep = 3;
raw_mu = rand(1,options);
mu = sort(zscore(raw_mu));
S = zeros(options);
W = zeros(options);
N = zeros(options);
for i = 1:options
for j = i:options
if i ~= j
S(i,j) = 1 - normcdf(delta1, mu(i)-mu(j));
W(i,j) = normcdf( delta1, mu(i)-mu(j)) - normcdf( delta0, mu(i)-mu(j));
N(i,j) = normcdf( delta0, mu(i)-mu(j)) - normcdf(-delta0, mu(i)-mu(j));
W(j,i) = normcdf(-delta0, mu(i)-mu(j)) - normcdf(-delta1, mu(i)-mu(j));
S(j,i) = normcdf(-delta1, mu(i)-mu(j));
end
end
end
diffs = bsxfun(@minus,mu,mu');
deltass = sort([-Inf, Inf, deltas, -1*deltas]);
deltass = reshape(deltass,[1,1,length(deltass)]);
normcdfs = diff(bsxfun(@(x,y) normcdf(y,x),diffs,deltass),1,3);
Upvotes: 0
Views: 94
Reputation: 15837
Here is a vectorized version of for loop
[I,J] = meshgrid(1:n);
idx = I<J;
I = I(idx);
J = J(idx);
IJ = sub2ind([n n],I,J);
JI = sub2ind([n n],J,I);
mu = x(I)-x(J);
md1 = normcdf(-delta1, mu).';
md0 = normcdf(-delta0, mu).';
d0 = normcdf( delta0, mu).';
d1 = normcdf( delta1, mu).';
Y = sum(S(IJ) .* log(max(0,1 - d1)) + ...
W(IJ) .* log(max(0,d1 - d0)) + ...
N(IJ) .* log(max(0,d0 - md0)) + ...
W(JI) .* log(max(0,md0 - md1)) + ...
S(JI) .* log(max(0,md1)));
Upvotes: 1