当矩阵的值很小时，为什么矩阵乘积会变慢？

Question

我创建两个矩阵A和B相同的尺寸。A包含比更大的值B。矩阵乘法的A %*% A速度大约是的10倍B %*% B。

为什么是这样？

## disable openMP
library(RhpcBLASctl); blas_set_num_threads(1); omp_set_num_threads(1)

A <- exp(-as.matrix(dist(expand.grid(1:60, 1:60))))
summary(c(A))
#     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
# 0.000000 0.000000 0.000000 0.001738 0.000000 1.000000 

B <- exp(-as.matrix(dist(expand.grid(1:60, 1:60)))*10)
summary(c(B))
#      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
# 0.0000000 0.0000000 0.0000000 0.0002778 0.0000000 1.0000000 

identical(dim(A), dim(B))
## [1] TRUE

system.time(A %*% A)
#    user  system elapsed 
#   2.387   0.001   2.389 
system.time(B %*% B)
#    user  system elapsed 
#  21.285   0.020  21.310

sessionInfo()
# R version 3.6.1 (2019-07-05)
# Platform: x86_64-pc-linux-gnu (64-bit)
# Running under: Linux Mint 19.2

# Matrix products: default
# BLAS:   /usr/lib/x86_64-linux-gnu/openblas/libblas.so.3
# LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so

当矩阵包含许多小条目时，问题可能与base :: chol（）变慢了。

编辑：有一些小数字，这似乎减慢了计算速度。其他人没有。

slow <-  6.41135533887904e-164
fast1 <- 6.41135533887904e-150
fast2 <- 6.41135533887904e-170

Mslow <- array(slow, c(1000, 1000)); system.time(Mslow %*% Mslow)
#   user  system elapsed 
# 10.165   0.000  10.168 

Mfast1 <- array(fast1, c(1000, 1000)); system.time(Mfast1 %*% Mfast1)
#   user  system elapsed 
#  0.058   0.000   0.057 

Mfast2 <- array(fast2, c(1000, 1000)); system.time(Mfast2 %*% Mfast2)
#   user  system elapsed 
#  0.056   0.000   0.055 

Answer 1

R-devel 邮件列表的回复表明，这可能是非正规数的问题，或者 openBLAS 处理小数的速度可能会变慢。

来自https://en.wikipedia.org/wiki/Denormal_number：

在计算机科学中，非正规数或非正规化数（现在通常称为次正规数）填补了浮点运算中零附近的下溢间隙。任何大小小于最小正规数的非零数都是“次正规数”。[...] 在极端情况下，涉及非正规操作数的指令运行速度可能会慢 100 倍。

事实上，B包含非常小的数字：

sum(B<.Machine$double.eps)
[1] 12832980
sort(unique(B[B>0]))[10^(0:3)]
[1] 4.940656e-324 2.280607e-320 6.302966e-295 2.185410e-141

如果较小的数字设置为零，则计算具有预期的计算时间：

C <- B; C[abs(C)<.Machine$double.eps] <- 0
system.time(C %*% C)
   user  system elapsed 
  2.266   0.032   2.298 

有没有办法自动将以下值设置.Machine$double.eps为零？ 手动检查每个矩阵中的小数字似乎并不方便。