Matrix multiplication-different results depending on how you perform the multiplication

[Christophe Kolodziejczyk]

Dear Statalist users,
I am experiencing that multiplying two matrices can give different results depending on how you multiply those matrices.
For exemple, I have a matrix X and a matrix b, where X is N x k and b is a k x k, and I want to compute X*b.
I can just write xb = X*b in mata or I can loop over the columns of b, multiply X with each column of b and store the results in a new matrix xb (which is k x k).
Why I want to use the second method can be justified by the context of the problem. I would expect to find exactly the same results though. But it is apparently not the case in the current version of Stata.
I have written below an example which shows that there are differences. Although those differences are very small they may matter for your problem.
I am using Stata MP 18 on Windows 10 and Windows Server 2016. I ran the same program with an earlier version of Stata (12) and the issue did not show up. Running the program in version 18 with version control doesn't fix the problem.
Has someone any insight on why I obtain different results?

mata:
------------------------------------------------- mata (type end to exit) --------------------------------------------------------------
: N = 20

: k = 3

: rseed(101)

: X  = runiform(N,k)

: b  = runiform(k,k)

: 
: xb1 = X*b 

: xb0 = J(rows(xb1),cols(xb1),.)

: for (i = 1; i <= k ; i++)
> {
>         xb0[,i] = X*b[,i]
> 
> }

: 
: // xb0 = cross(X',b)
: xb0:!=xb1
        1   2   3
     +-------------+
   1 |  1   1   0  |
   2 |  1   0   0  |
   3 |  0   0   0  |
   4 |  0   0   0  |
   5 |  0   0   0  |
   6 |  1   0   0  |
   7 |  1   0   0  |
   8 |  0   0   0  |
   9 |  1   1   0  |
  10 |  0   1   0  |
  11 |  0   0   0  |
  12 |  0   0   1  |
  13 |  1   1   1  |
  14 |  0   0   1  |
  15 |  1   0   1  |
  16 |  1   0   0  |
  17 |  0   1   0  |
  18 |  0   0   0  |
  19 |  0   1   0  |
  20 |  0   0   0  |
     +-------------+

: mreldif(xb0,xb1)
  9.06980e-17

: 
: end 
----------------------------------------------------------------------------------------------------------------------------------------

[John Mullahy]

Curious.

Using this adaptation of your code:

mata:

N = 20

k = 3

rseed(101)

X  = runiform(N,k)

b  = runiform(k,k)

 
xb1 = X*b

xb0 = J(rows(xb1),cols(xb1),.)

for (i = 1; i <= k ; i++)
 {
         xb0[,i] = X*b[,i]
 
 }

 
 xb2 = cross(X',b)
 xb3 = quadcross(X',b)

  xb0:!=xb1
  xb2:!=xb1
  xb3:!=xb1
 
  mreldif(xb0,xb1)
  mreldif(xb2,xb1)
  mreldif(xb3,xb1)
end

I obtained the following results on two different Macs and two different Stata versions:

1) MacOS Sierra 10.12.6 and Stata v17

:   mreldif(xb0,xb1)
  9.06980e-17

:   mreldif(xb2,xb1)
  0

:   mreldif(xb3,xb1)
  1.02670e-16

2) MacOS Ventura 13.5.2 and Stata v18

:   mreldif(xb0,xb1)
  0

:   mreldif(xb2,xb1)
  0

:   mreldif(xb3,xb1)
  1.02670e-16

[Christophe Kolodziejczyk]

Thanks John. Curious, indeed. It also seems to depend on the OS. It is reassuring that X*b and cross(X',b) give the same result. I would expect small differences when using the quad precision version of cross() is used instead of the double precision.

In Stata 18 (and on Windows) I discovered that if I switch the Intel MKL LAPACK off, then the difference disappears. I think I will contanct Stata's technical support.

: set lapack_mkl off
...
: mreldif(xb0,xb1)
  0

[Mauricio Caceres]

Christophe Kolodziejczyk Very interesting. The mata manual p. 45 (https://www.stata.com/manuals/m.pdf) says the Intel MKL uses processor-specific code for different routines. Stata specifically says that the same code can cause differences across processors and operating systems. However, I also take it to mean that different routines that do the same thing may give slightly different results because they're basically implemented in slightly different ways (in this case, matrix times matrix vs matrix times vector would be different routines).

If you run

set lapack_mkl_cnr compatible

for instance (and restart Stata) then the differences go away (at least for me) even with

set lapack_mkl on

.

[Christophe Kolodziejczyk]

Mauricio Caceres Thanks for your insight and the reference to the mata manual, which I should have looked into after I mentioned the Intel MKL. It may well be the case that Mata uses different routines when it multiplies two matrices rather than a matrix with a vector. I will also look into the lapack_mkl_cnr settings. Thanks again.