Trying to avoid for cycles

Hello world!

I have to obtain a matrix T starting from a matrix A and a vector w. Specifically:

• matrix A is a reciprocal (i.e. A(i,j) = 1/A(j,i)) square matrix of dimension N composed by all positive elements between 0 and 9;
• w is a column vector of dimension N whose entries are between 0 and 1.

T is defined element-wise as T(i,j) = A(i,j)*w(i)/w(j), so it is still a square matrix of dimension n. I used two ways: a double for cycle and element-wise moltiplications. Before going for the element-wise operations, I tested if the two methods provide the same result using logical operations (T==T2), and it happened to be incorrect. An example follows. It occurs that matrices T and T2 are not completely identical (some elements in the logical variables ans are zeros), eventhough if I open T and T2 from workspace or display the critical elements in command window the two seem identical.

I'd say that the problem depends on the well-known problems of matlab with decimal digits, and most probably the bias is negligible, but here's where my doubt emerges: which one of the two is more correct (assuming that both are correct)?

A = [1 2 3 1 3 8;
    1/2 1 7 3 1/5 1;
    1/3 1/7 1 1/5 1/5 1/2;
    1 1/3 5 1 1 1/3;
    1/3 5 5 1 1 3;
    1/8 1 2 3 1/3 1];
n = size(A,1);
[w, lambda] = eig(A);
[Lmax, indL] = max(diag(lambda));
w = w(:,indL)/(sum(w(:,indL)));
T = A.*(w'./w);
T2 = zeros(size(A));
for i = 1:n
   for j = 1:n
     T2(i,j) = A(i,j)*w(j)/w(i);
   end
end
T == T2

Thank you in advance,

Riccardo

[dpb]

The problem depends on the well-known problems of matlab with decimal digits, ..."

The "problem" has nothing whatever to do with MATLAB, it's totally the result of IEEE floating point precision. See <What Every Computer Scientist Should Know>

For your particular example,

>> max(abs(T(:)-T2(:)))
ans =
   4.4409e-16
>> 

As far as the Q?, "Neither". Both will be within rounding as close as can get; the difference simply is whether is rounding up or down or the sequence of operations changing an intermediary rounding result.

[Fabio Freschi]

I completely agree with dpb. You are not comparing the same sequence of operation. In fact, if you just put the w(j)/w(i) operation in brackets in the double loop, you have all ones in your T == T2 test.

     T2(i,j) = A(i,j)*(w(j)/w(i));

[Riccardo Sorvillo]

Firstly, thank you both for the answers.

"The "problem" has nothing whatever to do with MATLAB, it's totally the result of IEEE floating point precision."

I didn't know that, I thought it was a problem of MATLAB only, my mistake.

Then, thank you Fabio Freschi for let me noticing the very silly mistake in the double for cycle.

I appreciated your helps.

[Bjorn Gustavsson]

That is not really a "mistake" in the general sense - only here for exact comparison between two construction-methods is it a "mistake". The general lesson to learn is that you have to "accept" floating-point accuracy limitations - even your algorithms should be stable to such small numerical errors.