LAPACK下实现Moore-Penrose伪逆求解超定线性模型的最优例程选择

你遇到的是典型的超定线性最小二乘求解场景，同时N远大于I、O的特征非常明确，你原本采用的先构造XᵀX再调用DSYSV求解的方案可行，但存在更高效、数值稳定性更好的LAPACK例程组合，无需显式构造XᵀX即可完成计算。

最优例程选择 #

1. 最小二乘求解场景（你的核心需求） #

你不需要单独计算Moore-Penrose伪逆矩阵，直接调用LAPACK封装好的最小二乘驱动例程即可，性能比你现有方案高，精度更好：

• 如果你能保证X是满秩矩阵（即XᵀX为对称正定矩阵，无零特征值），首选DGELS例程：
  该例程基于QR分解实现超定最小二乘求解，内部自动完成所有运算，不需要你手动计算XᵀX和XᵀY。你只需要按接口要求传入X、Y的维度和数据，运算完成后Y的存储位置会直接返回你需要的系数矩阵A。相比你原有方案，省去了XᵀY的乘法开销，当O数值较大时优势非常明显，同时避免了显式构造XᵀX带来的数值误差放大问题。
• 如果X可能不满秩（即XᵀX为对称半正定矩阵，存在零特征值），需要得到Moore-Penrose伪逆对应的最小范数解，推荐两个例程：
  • DGELSD：基于分治策略的SVD求解例程，计算速度快，支持设置奇异值截断阈值，自动过滤小奇异值得到稳定的伪逆解，小规模I/O下性能表现优异。
  • DGELSS：同样基于SVD的最小二乘求解例程，和DGELSD的区别是SVD实现逻辑不同，大尺寸矩阵下DGELSD的性能更优。

2. 单独计算Moore-Penrose伪逆的实现方案 #

如果你需要单独得到X的伪逆矩阵，LAPACK没有提供直接的驱动例程，标准实现是基于SVD的例程组合：

1. 调用DGESVD或DGESDD计算X的奇异值分解：X = UΣVᵀ
2. 构造Σ的伪逆Σ⁺：所有大于截断阈值的奇异值取倒数，其余值置0
3. 计算伪逆矩阵：X⁺ = VΣ⁺Uᵀ

额外优化建议 #

如果你需要用同一个X矩阵多次求解不同Y对应的A矩阵，可以提前缓存X的分解结果，后续每次求解的开销可以降到O(I*O)：

• 满秩场景：先调用DGEQRF完成X的QR分解，缓存Q、R矩阵，后续每次求解仅需要调用DORMQR + DTRTRS即可得到A矩阵
• 不满秩场景：先调用DGESDD完成X的SVD分解，缓存U、S、V矩阵，后续直接用分解结果计算得到A矩阵即可

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