Vectorizing gives slightly different solution

I'm writing some code to implement a numerical procedure in evaluating blocking probabilities for a circuit-switched network using the Erlang-B formula, and after vectorising some of the code I end up receiving a slightly different answer. More specifically, in a 1 by 4 array that is the output, the two pieces of code produce the same result as the first element, but the second third and fourth all differ by about 0.01-0.03. The code involves a while loop that exits upon the difference between iterations being less than some user specified tolerance, but the tolerance is currently set at 1e-15, so I can't imagine this makes any difference.

Basically, I was hoping someone would be able to give me some insight as to why/when this may occur. If need be I can provide the code - it involves a calling script and three functions (one is iterative).

Thanks

[madhan ravi]

Upload your code and datas

[James McCusker]

Here is my code - the issue is in the calcy.m function, the commented code is my vectorised code and the for loop is the code is replaces. Notice that the solution output from the script.m code changes from

B = 0.1079 0.2270 0.2321 0.1180

with the vectorised code, to

B = 0.1079 0.2448 0.2555 0.1256

with the for loop. Any help is greatly appreciated!

Answer 1

[Guillaume]

There are three problems with your vectorised version:

1. you repmat alpha at each step of the j loop. Not critical, just unnecessary. The repmat only needs to occur once (in fact it's not needed at all once 2) is fixed)
2. your alpha(tempi~=0) results in a vector (by necessity it can't be a matrix) and thus your prod goes awry.
3. the dimension you take the product along is wrong. You want the product across the rows, so it's dimension 1.

The fix would be to replace your two vectorised lines by:

tempProduct = relevantLoads.*prod(1-alpha.*temp, 1);

However, you can vectorise the whole lot (not for the faint hearted)

mask = permute(~eye(size(A, 1)), [1 3 2]);  %mask to ignore the current links for each row
routealpha = permute(A .* alpha, [3 2 1]);
y = sum(prod(1 - routealpha .* mask, 3) .* loads .* A, 2);

As a bonus the above creates a column vector y as opposed to the row vector of your original code. Since the elements of y corresponds to the rows of A, it makes sense for y to be a column vector.

Note that find was completely unnecessary in your original code. It would have worked just fine if you'd use

    validroutes = A(j,:)~=0; 
    temp = A(:,validroutes); temp(j,:) = 0; % store these routes, but make the link we are on 0
    relevantLoads = loads(validroutes); % loads for these routes
    tempProduct = zeros(1,nnz(validroutes);

Finally, it's good that you documented the function and the code. I'd recommend that you add to any function documentation a description of the inputs and outputs, including size and type requirements, e.g.:

function y = calcy(loads,A,alpha)
% calcy is a function that calculates the vector of yj values for a given
% vector of blocking probabilities alpha, the links, loads and link/route
% matrix A.
%   A: a MxN matrix of double. row ares links, columns are routes.
%   loads: a 1xN vector of double. loads of each route'
%   alpha: a Mx1 vector of double. blocking probability of link.
%   y: a Mx1 vector of double. ???