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

Question

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

这是为什么？

## 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 

- Nairolf

---

我不确定，但我可以在一个非常相似的系统上重现那种行为。 - IRTFM

Answer 1

您最好使用.Machine$double.xmin而不是double.eps。这将使更少的数字归零并产生相同的效果。为了避免次标准数，您可能需要使用编译器标志重新编译BLAS，将这些数字设置为零而不是引发FP陷阱。

- Hilmar Berger

---

你是否碰巧知道如何使用这些标志重新编译BLAS？ - Nairolf

这个线程可能会很有趣：https://github.com/xianyi/OpenBLAS/issues/1237 - Avraham

Answer 2

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的值设为零的方法？ 手动检查每个矩阵中的小数似乎不太方便。

- Nairolf