应避免直接调用 arma::pinv()或 arma::inv(),因 GLM 需 IRLS 迭代求解加权最小二乘,显式构造 XᵀWX 易放大误差;推荐用 arma::solve(Xw, yw, arma::solve_opts::no_approx)基于 QR 分解稳健求解。

为什么 不用现成库直接调用 arma::pinv() 或 arma::solve()
因为广义线性回归(GLM)不是简单最小二乘;它需要迭代重加权最小二乘(IRLS),每轮都要解一个带权重的加权最小二乘问题:w * X 和 w * y。直接对设计矩阵 X 求逆(比如用 arma::inv())既不稳定也不必要——尤其当 X 列满秩都不满足时,inv() 会崩溃或返回垃圾值。更可靠的做法是用 QR 分解或 SVD 解加权系统,而 Armadillo 的 arma::solve(A, b, arma::solve_opts::fast) 默认走 QR,已足够稳健。
如何用 Armadillo 正确实现加权最小二乘更新步
IRLS 每轮需解:argmin_β || W^(1/2) (X β − y) ||²,等价于求解:(Xᵀ W X) β = Xᵀ W y。但显式构造 Xᵀ W X 会放大数值误差,应避免。推荐直接调用带权重的求解接口:
arma::vec weighted_ls_solve(const arma::mat& X, const arma::vec& y, const arma::vec& w) {arma::mat Xw = arma::diagmat(arma::sqrt(w)) * X; arma::vec yw = arma::diagmat(arma::sqrt(w)) * y; return arma::solve(Xw, yw, arma::solve_opts::no_approx); }
-
arma::solve_opts::no_approx强制使用更稳的 QR(而非默认可能启用的快速近似) - 权重向量
w必须非负;若出现w(i) == 0,对应样本被完全忽略,sqrt(w)仍合法 - 若
X高度共线性,可改用arma::solve(Xw, yw, arma::solve_opts::svd)启用 SVD 降维
Logistic 回归中 IRLS 权重和工作响应怎么算
以 logistic 为例:链接函数 g(μ) = logit(μ) = log(μ/(1−μ)),方差函数 V(μ) = μ(1−μ)。第 k 轮迭代中:
- 当前预测均值:
mu = 1.0 / (1.0 + arma::exp(-X * beta_old)) - 工作响应:
z = X * beta_old + (y - mu) % (1.0 / (mu % (1.0 - mu))) - 权重:
w = mu % (1.0 - mu) - 然后调用
weighted_ls_solve(X, z, w)得到新beta_new
注意:% 是 Armadillo 的逐元乘法;所有中间量必须用 arma::vec/arma::mat,不能混用裸数组——否则 隐式转换 可能静默截断精度。
立即学习“C++ 免费学习笔记(深入)”;
常见崩溃点和绕过方式
实际调试时最常卡在三处:
-
std::bad_alloc:发生在构造超大arma::mat(如百万级样本 × 百维特征)时。解决方法 是改用分块计算或换arma::sp_mat(稀疏矩阵),但 GLM 权重通常稠密,慎用稀疏 -
arma::solve(): solution not found:说明当前Xw秩亏。不要急着加岭参数,先检查w是否全为零(比如初始beta过大导致mu全趋 0 或 1),或加入防溢出 clamp:mu = arma::clamp(mu, 1e-15, 1.0 - 1e-15) - 收敛震荡:IRLS 不保证全局收敛。建议限制最大迭代次数(通常 20–50 足够),并监控目标函数变化:
deviance = 2 * arma::sum(y % arma::log((y+1e-15)/(mu+1e-15)) + (1-y) % arma::log((1-y+1e-15)/(1-mu+1e-15)))
真正麻烦的是链接函数和分布不匹配——比如用 identity 链接配泊松响应却不做非负约束,这时候数值解出来也无意义。算法能跑通,不代表模型对。






























