R语言CVXR包大数据集下拟合薄板样条报错问题咨询

问题现象 #

使用R语言CVXR包针对脑部扫描数据集拟合thin-plate spline（薄板样条）模型时，核心求解函数返回无法解读的报错。报错触发概率与数据集规模强相关：全量1567条观测下稳定触发，随机抽取数百条小样本运行时基本不会出现该问题。

可复现问题的完整示例代码 #

require(gamair)
require(CVXR)
require(npreg)
data(brain)
x = brain[, c(1, 2)]
x <- as.matrix(x)
z = brain$medFPQ
z <- as.vector(z)

m = 2
b.tp <- basis.tps(x, knots = x, m = m, rk = TRUE, intercept = TRUE)
pen.tp <- penalty.tps(x, m = m, rk = TRUE)
mstar <- choose(m+dim(x)[2]-1, dim(x)[2])
pen.tp <- rbind(matrix(0, ncol = dim(pen.tp)[2]+mstar, nrow = mstar), cbind( matrix(0, nrow = dim(pen.tp)[1], ncol = mstar ), pen.tp ) )
theta <- Variable(dim(b.tp)[2])
obj <- sum((z-b.tp%*%theta)^2) + 1e-01*quad_form(theta, pen.tp)
prob <- Problem(Minimize(obj))
result <- solve(prob, solver = "SCS")

运行返回的报错信息 #

Error in (function (cl, name, valueClass)  : 
  assignment of an object of class “complex” is not valid for @‘eigvals’ in an object of class “Constant”; is(value, "numeric") is not TRUE

报错原因 #

CVXR在处理quad_form传入的矩阵时，会自动计算矩阵特征值校验其是否符合凸优化要求的半正定属性。当样本量较大时，npreg包生成的薄板样条惩罚矩阵会因浮点计算累积误差出现两个问题：

• 矩阵上下三角存在极小的数值差，不是严格对称矩阵
• 存在幅值极小的负特征值，R在计算这类接近奇异的矩阵特征值时会产生带虚部的复数结果，触发CVXR内部的类型校验错误
  小样本下浮点误差幅值极低，不会触发上述问题，因此报错概率极低。

修复方法 #

在构造完惩罚矩阵pen.tp、定义目标函数之前，增加两行矩阵修正代码即可：

# 强制矩阵严格对称，消除浮点误差导致的上下三角不对称
pen.tp <- (pen.tp + t(pen.tp))/2
# 增加极小对角扰动，保证矩阵严格半正定，扰动幅值远小于惩罚系数，不影响拟合结果
diag(pen.tp) <- diag(pen.tp) + 1e-10

修正后全量1567条观测即可正常求解，无需更换数据集或调整模型参数。如果修正后仍存在兼容问题，可尝试将求解器从SCS更换为OSQP求解。

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