好累啊，绕那么大个弯……是问题没问好，还是都太不专业了？

a <- array(rnorm(20),dim=c(4,5,100))

这样就得到结果了，何必非要用matrix？
好吧，楼主这确实是个好问题，涉及到存储方式的问题，说一下自己的想法：
首先，矩阵实际上是有一个附加属性（维数dim）的向量，维数即为一个长度为2的向量，用来指定矩阵的行数和列数。
其次，当a被赋值NULL，如果a下一步被matrix赋值，实质上是用matrix生成的矩阵整体覆盖a的NULL值，a整体属性发生改变
最后，由于语句中，a是被循环赋值，matrix看起来只是覆盖a中的一个标量，而由于标量无法单独改变属性，因此matrix生成的矩阵赋值失败，a中的那个标量仅仅被matrix生成矩阵中的第一个值所覆盖，可以通过下列语句来判断其过程：
a <- NULL
for(i in 1:100) a[i] <- matrix(2:21,nrow = 4, ncol = 5)
此时可以发现a全部被matrix生成矩阵中的第一个值“2”所覆盖

存成3为数组
a初始化为4*5*100，然后给赋值a[,,i]
另外，不推荐用for循环，试试apply族函数

多维可以
list也可选择

这样都不给悬赏，感觉不会再爱了……