EDIT: after the code provided by OP, the problem became clear. Trick - convert it to data frame:
> x = array(rnorm(437216*8*3), dim = c(437216,8,3))
> system.time(apply(x,1:2,mean))
user system elapsed
107.06 0.18 107.34
> Tmp <- data.frame(V1=as.vector(x[,,1]),
+ V2=as.vector(x[,,2]),
+ V3= as.vector(x[,,3]))
> system.time({
+ Means <- rowMeans(Tmp)
+ Sd <- sqrt(rowSums((Tmp-Means)^2)/(3-1))
+ })
user system elapsed
6.72 0.40 7.12
To get the results in the correct matrix:
Means <- matrix(Means,ncol=8)
Sd <- matrix(Sd,ncol=8)
Proof of concept:
x = array(rnorm(10*8*3), dim = c(10,8,3))
m1 <- apply(x,1:2,mean)
sd1 <- apply(x,1:2,sd)
Tmp <- data.frame(V1=as.vector(x[,,1]),
V2=as.vector(x[,,2]),
V3= as.vector(x[,,3]))
m2 <- rowMeans(Tmp)
sd2 <- sqrt(rowSums((Tmp-m2)^2)/2)
m2 <-matrix(m2,ncol=8)
sd2 <- matrix(sd2,ncol=8)
> all.equal(m1,m2)
[1] TRUE
> all.equal(sd1,sd2)
[1] TRUE