R中构建大维度对角矩阵W内存不足，求Fischer Scoring替代方案

我有一个包含67856个元素的大型数据集，想要构造X^TWX，但尝试构建维度为67856×67856的对角矩阵W时，出现错误：vector memory exhausted (limit reached?)。请问有没有替代方法构造该对角矩阵W，以支持Fisher Scoring算法运行？

代码如下：

initial_model1=(lm(ClaimNb~VehValue+VehAge+DrivAge,data=final_data))
r=coef(summary(mleNB))[1,1]
beta0=initial_model1$coefficients[1]
beta1=initial_model1$coefficients[2]
beta2=initial_model1$coefficients[3]
beta3=initial_model1$coefficients[4]
y=c(final_data$ClaimNb)
xi1=c(final_data$VehValue)
xi2=c(final_data$VehAge)
xi3=c(final_data$DrivAge)

##Column vector of regressors
regressors=matrix(c(beta0,beta1,beta2,beta3),ncol=1)
n=nrow(final_data)
beta0=regressors[1]
beta1=regressors[2]
beta2=regressors[3]
beta3=regressors[4]
intercept=rep(1,nrow(final_data))
X=matrix(c(intercept,xi1,xi2,xi3),ncol=4)
XT=t(X)
W=matrix(rep(0),n,n)
# Error: vector memory exhausted (limit reached?)

问题根源在于显式构造67856×67856的矩阵需要约36GB内存（按双精度浮点数计算），远超出常规内存限制。完全不需要构造完整的W矩阵，利用对角矩阵的运算性质，用向量替代矩阵即可完成计算：

• 用向量存储W的对角元素
  对角矩阵W的核心信息就是其对角线上的n个元素，直接用长度为n的向量w_vec存储即可，内存需求仅约527KB，完全可控。
  针对你的负二项模型Fisher Scoring场景，W的对角元素为μ_i*(1+μ_i/r)，其中μ_i = exp(Xβ)，计算代码如下：

  mu = exp(X %*% regressors)
  w_vec = mu * (1 + mu / r)
  

• 利用向量运算替代矩阵乘法
  计算X^TWX时，无需构造W矩阵，直接通过行加权后的X矩阵与X的转置相乘即可。在R中用sweep函数实现行加权：

  # 对X的每一行乘以对应的w_vec元素
  X_weighted = sweep(X, 1, w_vec, `*`)
  # 计算X^TWX
  XTWX = t(X) %*% X_weighted
  

• 完整修改后的代码片段
  替换原代码中构造W的部分，整合后的关键代码如下：

  initial_model1=(lm(ClaimNb~VehValue+VehAge+DrivAge,data=final_data))
  r=coef(summary(mleNB))[1,1]
  regressors=initial_model1$coefficients
  # 转成列向量
  regressors=matrix(regressors, ncol=1)
  y=final_data$ClaimNb
  # 简化X矩阵构造
  X=cbind(1, final_data$VehValue, final_data$VehAge, final_data$DrivAge)
  n=nrow(final_data)
  
  # 计算W的对角元素向量
  mu = exp(X %*% regressors)
  w_vec = mu * (1 + mu / r)
  
  # 计算X^TWX
  X_weighted = sweep(X, 1, w_vec, `*`)
  XTWX = t(X) %*% X_weighted
  

这样既避免了内存耗尽的问题，又完全符合Fisher Scoring算法的计算要求。

内容的提问来源于stack exchange，提问作者mikaelaa