Incorrect result and non-deterministic hardware specific result with threaded OpenBLAS on Intel Xeon 8488C for matrix with sparse column

[jwasteneys-ccl]

Summary
We are observing incorrect and non-deterministic results when performing matrix multiplication of the form t(X) %*% X in R when using OpenBLAS on Intel Xeon 8488C CPUs. The issue appears only when multi-threading. The bug does not occur on other CPU architectures (e.g., AMD EPYC, Intel 8259CL) and is not present when limiting OpenBLAS to a single thread or using R’s internal matrix multiplication routines.

System Information:

Architecture: x86_64
CPU op-mode(s): 32-bit, 64-bit
Address sizes: 46 bits physical, 48 bits virtual
Byte Order: Little Endian
CPU(s): 8
On-line CPU(s) list: 0-7
Vendor ID: GenuineIntel
Model name: Intel(R) Xeon(R) Platinum 8488C
CPU family: 6
Model: 143
Thread(s) per core: 2
Core(s) per socket: 4
Socket(s): 1
Stepping: 8
BogoMIPS: 4800.00

OS: Debian Linux

BLAS: OpenBLAS (/usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 )

Threading: pthreads

LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0

R version 4.4.1 (2024-06-14)

Reproducible Example:

Xsmall.csv

In this sample data, column 2449 is all exactly zero. When calculating X'X the element at [2448,2449] should also be exactly zero:

R
X <- as.matrix(read.csv("Xsmall.csv"))
options(matprod = "blas")
(t(X) %*% X)[2448, 2449]  # returns 0.002292634

# issue fixed with crossprod:
crossprod(X)[2448, 2449]  # returns exactly 0

# Manual sum:
sum(X[, 2448] * X[, 2449])  # returns exactly 0

# issue fixed with internal R implementation
options(matprod = "internal") 
(t(X) %*% X)[2448, 2449]  # returns exactly 0

Only on affected hardware (Intel 8488C), when using more than 1 thread, we see an incorrect value (~0.00229). Setting OPENBLAS_NUM_THREADS=1 makes the issue go away.

We have only reproduced the issue so far with the attached matrix. Note this matrix is ill-conditioned. However, this should not affect the result to this magnitude.

Here is a chart of the result we get for element [2448,2449] with respect to the OPENBLAS_NUM_THREADS setting:

Observations

• The issue does not appear on AMD EPYC 7R13 or Intel 8259CL.
• Only observed with multithreaded OpenBLAS on AVX-512 capable hardware.
• Dot products involving entirely zero columns behave unexpectedly under multithreading.
• Manual dot product confirms result should be exactly 0.
• Disabling threading (OPENBLAS_NUM_THREADS=1) makes the issue disappear.

Hypothesis
This appears to be a bug in the multithreaded dgemm implementation when using AVX-512. Potential contributing factors include:

• Accumulation errors across SIMD lanes when computing many dot products in parallel.

• Internal tiling and block accumulation strategies may be introducing rounding errors even when one operand is all zero.

• The result (~0.0023) is too large to be explained by standard floating-point noise, suggesting a flaw in reduction logic.

[martin-frbg]

Both Intel systems would be using the same (AVX-512) DDOT kernel - are you sure you are using the same software versions on them ?

[jwasteneys-ccl]

All the jobs are running on the same docker image, and I've compared all the session/environment info between them. AVX was just a hypothesis. Could it be related to Cascade lake vs Saphire rapids?

[martin-frbg]

Same DGEMM kernel code and parameters between the two Intel systems (the main difference is in BFLOAT16). But given that you seem to be using version 0.3.20 from three years ago (did not notice this initially, sorry), it is possible that it may not have identified either cpu correctly, potentially falling back to much older code. I wonder if you could rerun with environment variable OPENBLAS_VERBOSE=2 to have it print out what set of kernels it is using ?

[jwasteneys-ccl]

When I turn on that variable I get:
Core: Cooperlake

I also tried setting OPENBLAS_CORETYPE=haswell, and I got the same (incorrect) results as on cooperlake

[jwasteneys-ccl]

Actually I take that back, with multiple threads it looks like the results are different on Haswell vs Cooperlake:

Core: Haswell
[1] -2.660589e-10

Core: Cooperlake
[1] 0.002292634

[martin-frbg]

Unfortunately I cannot reproduce this (with R-4.5.0) on any system I have access to (which sadly does not include 8488C specifically), and I do not see what would single out that particular cpu.

[jwasteneys-ccl]

Is there any other information that I could provide that would help to narrow it down?

[martin-frbg]

Just to be sure - are

BLAS: OpenBLAS (/usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 )
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;

both pointing to the same library ?

[jwasteneys-ccl]

This is what I currently have set up (reproduces issue):

Core: Cooperlake

> print((t(Xsmall) %*% Xsmall)[2448, 2449] )
[1] 0.002292634
> sessionInfo()
R version 4.4.1 (2024-06-14)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.5 LTS

Matrix products: default
BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: Etc/UTC
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
 [1] httr_1.4.7          compiler_4.4.1      R6_2.5.1           
 [4] cli_3.6.3           tools_4.4.1         base64enc_0.1-3    
 [7] curl_5.2.3          aws.signature_0.6.0 xml2_1.3.6         
[10] aws.s3_0.3.21       digest_0.6.37       rlang_1.1.4 

[jwasteneys-ccl]

Sorry I made a mistake for coretype=Haswell, the issue goes away.

I had a piece of test code that set Xsmall[,2449] <- 1e-10 to see if a small number behaved differently than zero which gave me this result:

Core: Haswell
[1] -2.660589e-10

When really the result should be:

Core: Haswell
[1] 0

And for Cooperlake:

Core: Cooperlake
[1] 0.002292634

I also rebuilt my docker container to be sure:

> sessionInfo()
R version 4.4.1 (2024-06-14)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.5 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0

I'm not sure what the difference between libblas.so.3 and libopenblasp-r0.3.20.so and why it switched from before.

So it seems like it is related to Cooperlake in particular.

I also tried with Skylake:

Core: SkylakeX
[1] 0

So it seems like the issue is might be what you originally thought, and it's identifying the wrong kernel?

[martin-frbg]

I wonder if you could install the current OpenBLAS in your docker container to check ? The libblas.so.3 and libopenblas... in your sessioninfo output are (or should be) links to the same shared library, it is just that stock R has them as libblas.so.3 and liblapack.so.3 (and corresponding to the respective unoptimized library from the Reference-LAPACK/BLAS by default)

[jwasteneys-ccl]

I think the version of ubuntu I'm using (ubuntu:jammy-20240125) pins the blas version, so I haven't been able to easily update it