Statsmodels GLM串行运行出现多重共线性，批量提交结果正常

同一代码、conda环境、数据及库版本条件下，使用sm.GLM结合自然三次样条拟合数据时，串行运行得到异常的超大系数结果，但通过Slurm批量提交运行结果正常。

所用代码 #

import numpy as np
import pandas as pd
import statsmodels.api as sm
from patsy import dmatrix

size = 5000
rounded_data = np.loadtxt('test.txt')
c1 = np.repeat(np.arange(rounded_data.shape[0]), rounded_data.shape[0])
c2 = np.tile(np.arange(rounded_data.shape[0]), rounded_data.shape[0])
c3 = rounded_data.flatten()
dist = np.abs(c1 - c2) * size
# 注：代码中bin1、bin2未定义，需确认其来源
dfs = pd.DataFrame({'c1': bin1, 'c2': bin2, 'c3': c3, 'distance': dist})
train_x = dfs["distance"]
train_y = dfs["c3"]
# 使用3阶自然三次样条
n_knots = 12

knots = np.percentile(train_x, np.linspace(0, 100, n_knots + 2)[1:-1])

transformed_x = dmatrix("cr(train_x, df=n_knots + 2)", {"train_x": train_x})

y_mean = np.mean(train_y)
y_var = np.var(train_y, ddof = 1)
alpha = y_mean ** 2 / (y_var - y_mean)

fit1 = sm.GLM(dfs["c3"], transformed_x, family=sm.families.NegativeBinomial(alpha=alpha))
result = fit1.fit()
print(result.summary())

串行运行异常结果 #

Generalized Linear Model Regression Results                  
==============================================================================
Dep. Variable:                  c3   No. Observations:               131627
Model:                            GLM   Df Residuals:                   131613
Model Family:        NegativeBinomial   Df Model:                           13
Link Function:                    Log   Scale:                          1.0000
Method:                          IRLS   Log-Likelihood:            -4.9962e+05
Date:                Thu, 11 May 2023   Deviance:                   1.4181e+05
Time:                        17:31:42   Pearson chi2:                 1.50e+05
No. Iterations:                   100   Pseudo R-squ. (CS):             0.9999
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const        8.07e+10   3.09e+10      2.616      0.009    2.02e+10    1.41e+11
x1          -8.07e+10   3.09e+10     -2.616      0.009   -1.41e+11   -2.02e+10
x2          -8.07e+10   3.09e+10     -2.616      0.009   -1.41e+11   -2.02e+10
x3          -8.07e+10   3.09e+10     -2.616      0.009   -1.41e+11   -2.02e+10
x4          -8.07e+10   3.09e+10     -2.616      0.009   -1.41e+11   -2.02e+10
...

Slurm批量运行正常结果 #

Generalized Linear Model Regression Results                  
==============================================================================
Dep. Variable:                  c3   No. Observations:               133083
Model:                            GLM   Df Residuals:                   133069
Model Family:        NegativeBinomial   Df Model:                           13
Link Function:                    Log   Scale:                          1.0000
Method:                          IRLS   Log-Likelihood:            -4.2261e+05
Date:                Thu, 11 May 2023   Deviance:                   1.4775e+05
Time:                        17:12:12   Pearson chi2:                 1.55e+05
No. Iterations:                   100   Pseudo R-squ. (CS):             0.9998
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const          1.9531      0.004    435.520      0.000       1.944       1.962
x1             3.7154      0.007    547.396      0.000       3.702       3.729
x2             1.4755      0.005    268.910      0.000       1.465       1.486
x3             1.1598      0.006    209.027      0.000       1.149       1.171
x4             0.6174      0.006    105.013      0.000       0.606       0.629
...

排查方向 #

• 确认数据一致性：两次运行的观测数不同（131627 vs 133083），说明数据加载或处理环节存在差异。需检查：

  • test.txt的路径是否为绝对路径，串行和Slurm环境下是否加载了同一文件
  • 代码中未定义的bin1、bin2变量，在两种环境下的值是否一致，这直接影响dfs的行数
  • 输出rounded_data.shape、dfs.shape确认数据维度

• 检查IRLS迭代收敛性：异常结果迭代达到100次（默认最大迭代次数），可能未收敛。可尝试：

  • 增加迭代次数：result = fit1.fit(max_iter=200)
  • 使用Slurm的正常系数作为初始值：result = fit1.fit(start_params=slurm_coef_array)

• 核实环境细节：

  • 确认串行和Slurm节点的numpy、scipy、statsmodels版本完全一致，运行print(np.__version__, sm.__version__)对比
  • 检查线性代数后端（如MKL/OpenBLAS）的差异，数值计算的细微差别可能影响迭代路径

• 验证样条基构造：对比两种环境下transformed_x的矩阵特征（如求和、特征值），确认样条基函数构造无差异

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