Search code examples
matlabnumericlinear-equation

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


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.


Solution

  • 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.