Why do these two operations behave so differently?

I used two different laptops to run the same code but got very different results (fac and fac2), could somebody tell me the reason? Thank you in advance. Please first download the data I uploaded. Then use the following code:

    load data.mat

    A=(stdclrdatas-stdclrdata)<0.0000001; % check whether two matrices are different
    sum(sum(A))
    abs(Rs-R)<0.00000001 % check whether two matrices are different
    abs(Ls-L)<0.0000001 % check whether two matrices are different

    n=30
    fac = stdclrdata * (R \ diag(ones(n,1))) * L;
    fac2 = stdclrdatas * (Rs \ diag(ones(30,1))) * Ls;

Answer 1

[Walter Roberson]

You did not define n for calculating fac. It appears it should be 30.

Both fac and fac2 report "Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = " . That is because both matrices are singular: rank() of them is 29.

diag(ones(30,1)) is the same as eye(30)

When you have a square matrix \ eye() of appropriate size, that is equivalent to inv() of the matrix except perhaps calculated slightly better.

stdclrdatas * inv(Rs) * Ls

is equivalent to stdclrdatas / Rs * Ls except that this version is calculated more accurately.

You should be replacing your calculations with

fac = stdclrdata / R  * L;
fac2 = stdclrdatas / Rs * Ls;

which should give you less trouble between versions.

The method you were using is quite unstable. Calculating the inverse tends to magnify numeric problems, especially singularities.

Using the forms I show will still have some variation between systems.

For sufficiently large systems, the high performance linear algebra toolboxes such as MKL get used. Those toolboxes are sometimes updated between releases. The toolboxes use multiple CPUs; different number of CPUs available on different systems can lead to subtle differences in the calculations due to chunks of operations being done in slightly different orders

[Jian Wang]

Hi Walter,

Thank you for your reply. But how can I go ahead facing such case where I need to compute inversion of nearly singular matrix?

Best, Jian