Known issue or bug in solving system of equations

Is it known that MATLAB gives different results when solving a single system of equations versus multiple ones or is this a bug? Here is my sample code:

num = 4;
A = reshape(mod(1:25, 20), 5, 5)+eye(5);
x = ones(5, num);
b = A*x;
x_batch = A\b;
x_seq = [];
for i=1:num,
x_seq(:, i) = A\b(:, i);
end
assert(norm(x_batch(:)-x_seq(:))==0)

If num is less than 4, then the solutions are exactly the same. I would have thought solving multiple related systems would yield the same results but faster, since MATLAB would only have to decompose the A matrix once for whatever solver it decides to use. However, x_batch is not always exactly equal to x_seq. When I run the above code, I get a difference of 1.7e-14. While this number is small, the difference accumulates in my algorithm because I'm using the solution from the previous iteration to solve a related problem.

Other than just running in batches of 3, does anyone know how to work around this issue?

[Patrik Ek]

Interesting observation. That should be an independent operation. I guess it may depend on whether matlab have done some trick to speed up the operation and so induced a dependence.

I tried your code and found that norm(x_batch(:)-x_seq(:)) gave a norm 1.3323e-15 for me for num = 4. However, looking at the maximum error max(abs(x_batch(:)-x_seq(:))) I found it to be 6.6613e-16, (which is of course still quite good). This is exactly 3*eps. I also tried to do the same operation for A=rand(5). I observed that the error per element always gave integer or half intger fractions of eps

So to conclude, I have no final answer, but some interesting observations. There seems to be an error that appears when b has enough number of columns and most likely due to some matlab optimization process. The relation to eps gives some indication that this could be related to some operation using eps. This error also gets larger when the sizes of the matrices increases. You may need to contact the tecnical support about this issue.