大型数据集CVXR中的一个错误

Question

我一直试图使用流行的CVXR包在R中将薄板样条拟合到脑扫描数据集，但不幸的是，该包的主要功能返回了一个我无法破解的错误。下面是我一直在研究的例子。

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

我一直在想这是什么原因，因为我没有找到任何相关的信息。然而，我注意到这个错误不太可能出现在较小的数据集中。例如，在1567个可用的观测中，我们只使用了几百个随机抽样的观测结果。

如果有人有更多有关如何解决这个问题的资料，我可否向他们求助？谢谢。

Answer 1

我无法重现错误(GNU/Linux上的R4.2.0)。但要大胆猜测:显然计算了一个复杂的值，例如，当一个平方根是从一个负数计算出来的时候。这样的负数可能只是“数字噪声”(即由舍入误差引起的，实际上为零)。这样接近于零的负数可能是需要满秩矩阵的计算结果，但确实得到了秩不足矩阵。

我的会议信息：

sessionInfo()
R version 4.2.0 (2022-04-22)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 21.10

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-openmp/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-openmp/libopenblasp-r0.3.13.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] CVXR_1.0-10  npreg_1.0-8  gamair_1.0-2

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.8.3     lattice_0.20-45  codetools_0.2-18 gmp_0.6-5       
 [5] slam_0.1-50      grid_4.2.0       R6_2.5.1         Rmpfr_0.8-7     
 [9] Matrix_1.4-1     tools_4.2.0      bit64_4.0.5      bit_4.0.4       
[13] compiler_4.2.0   scs_3.0-0        cccp_0.2-7       Rglpk_0.6-4     

Answer 2

在我看来，我没有看到你的问题是随机的，所以结果应该是@JohnK得到的。我还可以确认CVXR解决了这个问题。你为什么要用SCS作为解算器？对于默认的求解器(OSQP)，速度是原来的三倍，并且返回最优状态，顺便说一句，应该始终检查！

> system.time(result <- solve(prob))
   user  system elapsed 
 49.118  24.082  42.783 
> cat(sprintf("Results---solver: %s, status: %s, value: %f\n", result$solver, result$status, result$value))
Results---solver: OSQP, status: optimal, value: 1172.806162
> system.time(result <- solve(prob, solver = "SCS"))
   user  system elapsed 
130.451  26.029 122.103 
> cat(sprintf("Results---solver: %s, status: %s, value: %f\n", result$solver, result$status, result$value))
Results---solver: SCS, status: optimal_inaccurate, value: 1172.570113

这是我的课程：

> sessionInfo()
R version 4.2.0 (2022-04-22)
Platform: x86_64-apple-darwin21.3.0 (64-bit)
Running under: macOS Monterey 12.4

Matrix products: default
BLAS:   /usr/local/Cellar/openblas/0.3.20/lib/libopenblasp-r0.3.20.dylib
LAPACK: /usr/local/Cellar/r/4.2.0/lib/R/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] npreg_1.0-8  CVXR_1.0-10  gamair_1.0-2

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.8.3     bit_4.0.4        lattice_0.20-45  R6_2.5.1        
 [5] tools_4.2.0      grid_4.2.0       osqp_0.6.0.5     scs_3.0-0       
 [9] bit64_4.0.5      assertthat_0.2.1 cccp_0.2-7       Matrix_1.4-1    
[13] gmp_0.6-5        Rglpk_0.6-4      codetools_0.2-18 slam_0.1-50     
[17] Rcplex_0.3-5     gurobi_9.5-1     compiler_4.2.0   Rmpfr_0.8-7     
[21] rcbc_0.1.0.9001  Rmosek_9.3.2    
>