Behaviour of backslash operator for non-square matrices least-squares fitting

I am currently trying to reproduce a set of results from a pre-existing project and can't get to the bottom of the following difference.

Phi = design matrix of input data size(29507x97)

Ref = reference values to fit to, six data sets, size(29507x6)

Phi = rand(29507,97);
Ref = rand(29507,6);

c1 = Phi\Ref;
for ii = 1:6
    c2(:,ii) = Phi\Ref(:,ii);
end

all(all(c1==c2))

I would have expected c1 to give identical results to c2. I can't find detailed information on what the "\" algorithm is doing differently between these case.

Any help appreciated!

Answer 1

[the cyclist]

They are equal, to within floating-point precision. Notice that

max(abs(c1(:)-c2(:)))

is around 1e-16.

[Matt O'Donnell]

Thanks. Yes. I figured it is floating point related, but I'm wondering what the difference is between the two methods algorithmically. From what I can gather from the slash documentation it should be performing exactly the same calculations.

In my data set I'm getting larger deviation between the two sets of coefficients and want to be 100% sure I can ignore this.

[the cyclist]

At the level of floating-point error, even the order of operations on arrays matters. Take a look at this, for example:

rng default
A = rand(10000);
mean(mean(A,1)) - mean(mean(A,2))

ans =

     4.440892098500626e-16

So, even though the two calculations are "exactly the same calculations", there is error. So, if the matrix version is doing something row-wise, that could lead to a difference.

I don't know how big the difference that you are seeing is, but I expect this is the root cause.

[Matt J]

Also, internal parallelization of mldivide and other linear algebra operations is different depending on the size of the inputs. So, when you perform the mldivide column-by-column, the inputs are divided up into parallel chunks in a different way, leading to floating point differences.

[Matt O'Donnell]

@the cyclist Yes. I agree with that as the likely root cause.

I thought it would be doing the same least squares calculation on a column by column basis as each column should be treated independently from each other in terms of doing the expected calculation.i.e. column 1 should never know anything about the other 5 columns's calculations mathematically speaking. But as @Matt J has suggested it is doing something across the columns, or even splitting each column differently, which leads to the difference/calculation order effect as suggested.

Many thanks to both of you for helping to get my head around this.