計算矩陣元素的最快方法

Question

我正在研究一些需要我重複計算大方形矩陣元素的東西。該過程涉及讀取存儲在另一個矩陣中的數據，然後計算矩陣元素。目前我正在使用循環的雙重來做到這一點。

library(matrixcalc) 

data <- matrix(nrow=3,ncol=1000) 

for(x in 1:ncol(data)){ 
    for(y in 1:ncol(data)){ 
     matrix[x,y]=exp(-entrywise.norm(data[,x]-data[,y],2)) 
    } 
} 

問題是，這是非常慢，因爲我的矩陣非常大。這個程序的最快選擇是什麼？

Answer 1

短而且速度非常快：

mat <- exp(-as.matrix(dist(t(data)))) 

我也想建議fields::rdist功能作爲一個更快的替代dist計算歐幾里得距離的矩陣，那麼如果裝載包不是一個問題，可以考慮：

library(fields) 
mat <- exp(-rdist(t(data))) 

---

爲了讓你的速度提高了一個想法：

data <- matrix(runif(3000), nrow=3, ncol=1000) 

OP <- function(data) { 
    require(matrixcalc) 
    mat <- matrix(0, ncol(data), ncol(data)) 
    for(x in 1:ncol(data)){ 
    for(y in 1:ncol(data)){ 
     mat[x,y]=exp(-entrywise.norm(data[,x]-data[,y],2)) 
    } 
    } 
    mat 
} 

flodel1 <- function(data) exp(-as.matrix(dist(t(data)))) 
flodel2 <- function(data) { 
    require(fields) 
    exp(-rdist(t(data))) 
} 

system.time(res1 <- OP(data)) 
# user system elapsed 
# 22.708 2.080 24.602 
system.time(res2 <- flodel1(data)) 
# user system elapsed 
# 0.112 0.025 0.136 
system.time(res3 <- flodel2(data)) 
# user system elapsed 
# 0.048 0.000 0.049 

（注意，在OP和flodel2的情況下，這些運行時間不包括包的負載，因爲他們已經被之前的測試加載。 ）

Answer 2

R語言使用column-major-order數組。更改for循環順序可以提高性能。因爲這樣，您可以更加連續地訪問內存，從而實現CPU緩存的優勢。

for(y in 1:dim) //outer is y now 
{ 
    for(x in 1:dim) //now x is count inside 
    { 
     matrix[x,y]=exp(-entrywise.norm(data[,x]-data[,y],2)) 
    } 
} 

你的 「矩陣」 是二維數組吧？

如果你需要更多的速度，你可以展開一些內循環，以減少cpu的分支負載和更好的緩存/預取。

for(y in 1:dim) 
{ 
    for(x in 1:(dim/8)) //lets imagine dimension is a multiple of 8 
    { 
     matrix[x,y]=exp(-entrywise.norm(data[,x]-data[,y],2)) 
     matrix[x+1,y]=exp(-entrywise.norm(data[,x+1]-data[,y],2)) 
     matrix[x+2,y]=exp(-entrywise.norm(data[,x+2]-data[,y],2)) 
     matrix[x+3,y]=exp(-entrywise.norm(data[,x+3]-data[,y],2)) 
     matrix[x+4,y]=exp(-entrywise.norm(data[,x+4]-data[,y],2)) 
     matrix[x+5,y]=exp(-entrywise.norm(data[,x+5]-data[,y],2)) 
     matrix[x+6,y]=exp(-entrywise.norm(data[,x+6]-data[,y],2)) 
     matrix[x+7,y]=exp(-entrywise.norm(data[,x+7]-data[,y],2)) 
    } 
} 

Answer 3

這應該是相當快：

nc <- ncol(data) 

mat <- diag(nc) 

for(x in 2:nc){ 
    for(y in 1:x){ 
     mat[x, y] <- exp(-(sum((data[ , x] - data[ , y])^2)^.5)) 
    } 
} 

mat[upper.tri(mat)] <- t(mat)[upper.tri(mat)] 

Answer 4

您可以使用colSums代替內環的。根據答案通過@Sven海恩斯坦：

nc <- ncol(data) 

mat <- diag(nc) 

for(x in 2:nc){ 
    mat[x, 1:x] <- exp(-(colSums((data[ , 1:x] - data[ ,x])^2)^.5)) 
} 

mat[upper.tri(mat)] <- t(mat)[upper.tri(mat)]