Strange numerical error when transposing a matrix and solving a linear system in Matlab

Question

I stumbled upon rather strange behaviour in MATLAB. The operator for solving a system of linear equations, \, sometimes produces different results, though the only thing that is changed is the place of transpose operator.

Take a look at this example:

A0 = rand(4);
b = rand(4,1);
A1 = A0';
x0 = A0\b;
x1 = A1'\b;
x2 = linsolve(A0,b);
x3 = linsolve(A1',b);
x4 = mldivide(A0,b);
x5 = mldivide(A1',b);

x0 = x2 = x3 = x4 = x5 but x0 != x1 (they differ in magnitude of 10^-15)

A0 in the above example is:

   0.781931966588002   0.530872257027928   0.112283962156027   0.964422667215901
   0.100606322362422   0.091498731339412   0.784427890743913   0.432484993970361
   0.294066333758628   0.405315419880591   0.291570317906931   0.694752194617940
   0.237373019705579   0.104846247115757   0.603533438750887   0.758099275289454

and b is:

   0.432642326147101
   0.655498039803537
   0.109755050723052
   0.933759848385332

I know this probably does not cause any practical problems but still I am curious for an explanation why this might happen.

To add to the strangeness: `A0'\b` differs not only from `A1\b`, but also from `ctranspose(A0)\b` or `transpose(A0)\b` (the latter three are equal) — Luis Mendo, Jan 13 '22 at 21:14
My guess: Matlab's interpreter is being smart and replacing `A0'\b` by `(b'/A0)'` (it is the same for non-complex numbers; I didn't test for complex). Why would it do that? Because transposing two vectors is faster than transposing a matrix. Why would the results differ? Because ```\``` and ```/``` perhaps use different algorithms, leading to different, small numerical errors. Does the Matlab interpreter really apply that sort of optimization? Maybe, see [this](https://stackoverflow.com/questions/57954599/does-matlabs-implicit-broadcasting-optimise-based-on-surrounding-code) — Luis Mendo, Jan 13 '22 at 21:23

score 4 · Accepted Answer · edited Jan 14 '22 at 07:37

I suspect it is the parser and how it feeds the matrices to the LAPACK library routines. E.g., in the matrix multiplication case of A'*B where A and B are matrices, the transpose operation isn't explicitly done. Rather, MATLAB calls the appropriate BLAS routine (e.g., DGEMM) with appropriate flags so that the equivalent operation is done, but may result in a different order of operations than if you had explicitly done the transpose first. I suspect this might be the case with your example, and that the transpose isn't explicitly done but flags are passed to the LAPACK library routines in the background to have a mathematically equivalent operation done but the actual order of operations is different resulting in a slightly different answer.

Strange numerical error when transposing a matrix and solving a linear system in Matlab

1 Answers1