R 函数 solve.QP 错误 "constraints are inconsistent, no solution!"
R function solve.QP error "constraints are inconsistent, no solution!"
我正在尝试 运行 使用具有以下参数的 solve.QP 函数(来自 quadprog 包)进行优化
R = matrix( c( 2.231113e-05,-4.816095e-05,-5.115287e-05, 0,2.989584e-05,4.212173e-06,0,0, 5.504990e-05), ncol=3, byrow=T)
b = c(-1,0,rep(0,ncol(R)))
C = cbind(rep(1,ncol(R)), diag(ncol(R)))
C = cbind(-rep(1,ncol(R)),C)
d = as.matrix(c(57621264,78057622,171342351),ncol=1)
H = solve.QP(Dmat = R, factorized = TRUE, dvec = d, Amat = C, bvec = b)
但是我收到的错误是
Error in solve.QP(Dmat = R, factorized = TRUE, dvec = d, Amat = C, bvec = b) :
constraints are inconsistent, no solution!
但是,当我为 R 使用不同的矩阵时
R2 = matrix( c( 0.05365071,-0.06364421,-0.04102565, 0, 0.08423283,-0.04048879,0,0,0.09659707), ncol=3, byrow=T)
solve.QP 通话
H = solve.QP(Dmat = R2, factorized = TRUE, dvec = d, Amat = C, bvec = b)
不会造成任何问题。
我现在的问题是为什么在前一个案例中会出现问题。
非常感谢任何帮助!
可能有两个问题。首先,您还没有输入参数 meq,它告诉求解器您的约束方程中有多少是等式约束,有多少是不等式,因此它默认为 meq=0,这使得它们都是等式约束,因此您已经过度确定了您的解决方案。看看你的问题,我可能猜想至少最后三个约束是不平等的;即解向量的分量都应该> 0。接下来,第二个方程似乎有点混乱。如果是不等式约束,则与后三者是多余的。如果是等式约束,则与第一个冲突。也许它不应该存在或应该以某种方式改变。
------------编辑---------------------------- --------------
我对您的第一个 post 回复有点太快了,现在明白您的所有约束都是不等式约束,因此您可以使用 meq 的默认值。在我看来,第二个约束与其余约束是多余的,但这似乎不会造成任何问题,因此目前并不重要。您的编辑还帮助我更好地理解了您的问题,并且我同意您使用 R 矩阵的示例应该可以在给定的约束条件下解决。可能是 R 中矩阵元素的大小可能导致 solve.QP 的解决方案问题,所以我尝试了一个更通用的非线性约束优化器 constrOptim。这为您的 R 和 R2 矩阵提供了解决方案,其中 R2 解决方案非常接近 R2 的 solve.QP 解决方案。
R2 = matrix( c( 0.05365071,-0.06364421,-0.04102565, 0, 0.08423283,-0.04048879,0,0,0.09659707), ncol=3, byrow=T)
R = matrix( c( 2.231113e-05,-4.816095e-05,-5.115287e-05, 0,2.989584e-05,4.212173e-06,0,0, 5.504990e-05), ncol=3, byrow=T)
d = as.matrix(c(57621264,78057622,171342351),ncol=1)
start <- rep(1/(ncol(R)+1), ncol(R))
min_fn <- function(b, dvec, Dmat) -t(dvec)%*%b +t(b)%*%Dmat%*%b/2
grad_min_fn <- function(b, dvec, Dmat) -dvec + Dmat%*%b
b = c(-1., 0, rep(0,ncol(R)))
C = cbind(rep(1,ncol(R)), diag(ncol(R)))
C = cbind(-rep(1,ncol(R)),C)
D <- t(solve(R))%*%solve(R)
constrOptim(theta=start, f=min_fn, grad=grad_min_fn, ui=t(C), ci=b, control=list(reltol=10*.Machine$double.eps),
dvec=d, Dmat=D )
solve.QP solution for R2
$solution
[1] -1.025463e-10 0.000000e+00 1.000000e+00
constrOptim solution for R2
$par
[1] 2.479855e-15 1.178762e-14 1.000000e+00
以及给出 R
的解决方案
$par
[1] 9.272547e-17 5.958225e-14 1.040137e-01
constrOptim 比 solve.QP 提供了更多有关其解决方案路径的信息,因此可能对数字敏感问题更有帮助。
我不知道 constrOptim 是否可以代替 solve.QP 为您工作,但至少它可以帮助您调查问题解决方案的属性。
你的矩阵 dvec 中的元素非常大。您的矩阵 R 是 Dmat
的 cholesky 分解中上三角矩阵的逆矩阵,即 Dmat = M^T M,其中 M = R^{-1}。所以以未分解的形式:
M <- solve(R)
Dmat <- t(M)%*%M
Dmat
unfactorized 很大,与 dvec
:
的规模大致相同
# [,1] [,2] [,3]
#[1,] 2008893283 3236243201 1619059085
#[2,] 3236243201 6332319710 2522625866
#[3,] 1619059085 2522625866 1641403882
因此您的问题可能是由于某种溢出错误引起的。要解决此问题,您可以缩放 Dmat
(未分解)和 dvec
:
sc <- norm(Dmat,"2")
solve.QP(Dmat = Dmat/sc, dvec=d/sc, Amat=C, bvec=b, meq=0, factorized=FALSE )
解决方案是
# $solution
# [1] -1.220832e-17 0.000000e+00 1.043877e-01
这与 constrOptim
在另一个答案中找到的解决方案非常匹配。
缩放 Dmat
和 dvec
很好,因为 -d^T b + 1/2 b^T D b
的约束最小值与任何常量 sc 的 sc*(-d^T b + 1/2 b^T D b)
的约束最小值相同。
编辑:
要以分解形式解决问题,您可以尝试以下缩放比例:
nn2 = sqrt(norm(d,"2"))
H = solve.QP(Dmat = R*nn2, dvec = d/(nn2^2), Amat = C, bvec = b, factorized=TRUE)
#$solution
#[1] 0.0000000 0.0000000 0.1043877
我正在尝试 运行 使用具有以下参数的 solve.QP 函数(来自 quadprog 包)进行优化
R = matrix( c( 2.231113e-05,-4.816095e-05,-5.115287e-05, 0,2.989584e-05,4.212173e-06,0,0, 5.504990e-05), ncol=3, byrow=T)
b = c(-1,0,rep(0,ncol(R)))
C = cbind(rep(1,ncol(R)), diag(ncol(R)))
C = cbind(-rep(1,ncol(R)),C)
d = as.matrix(c(57621264,78057622,171342351),ncol=1)
H = solve.QP(Dmat = R, factorized = TRUE, dvec = d, Amat = C, bvec = b)
但是我收到的错误是
Error in solve.QP(Dmat = R, factorized = TRUE, dvec = d, Amat = C, bvec = b) :
constraints are inconsistent, no solution!
但是,当我为 R 使用不同的矩阵时
R2 = matrix( c( 0.05365071,-0.06364421,-0.04102565, 0, 0.08423283,-0.04048879,0,0,0.09659707), ncol=3, byrow=T)
solve.QP 通话
H = solve.QP(Dmat = R2, factorized = TRUE, dvec = d, Amat = C, bvec = b)
不会造成任何问题。 我现在的问题是为什么在前一个案例中会出现问题。
非常感谢任何帮助!
可能有两个问题。首先,您还没有输入参数 meq,它告诉求解器您的约束方程中有多少是等式约束,有多少是不等式,因此它默认为 meq=0,这使得它们都是等式约束,因此您已经过度确定了您的解决方案。看看你的问题,我可能猜想至少最后三个约束是不平等的;即解向量的分量都应该> 0。接下来,第二个方程似乎有点混乱。如果是不等式约束,则与后三者是多余的。如果是等式约束,则与第一个冲突。也许它不应该存在或应该以某种方式改变。
------------编辑---------------------------- --------------
我对您的第一个 post 回复有点太快了,现在明白您的所有约束都是不等式约束,因此您可以使用 meq 的默认值。在我看来,第二个约束与其余约束是多余的,但这似乎不会造成任何问题,因此目前并不重要。您的编辑还帮助我更好地理解了您的问题,并且我同意您使用 R 矩阵的示例应该可以在给定的约束条件下解决。可能是 R 中矩阵元素的大小可能导致 solve.QP 的解决方案问题,所以我尝试了一个更通用的非线性约束优化器 constrOptim。这为您的 R 和 R2 矩阵提供了解决方案,其中 R2 解决方案非常接近 R2 的 solve.QP 解决方案。
R2 = matrix( c( 0.05365071,-0.06364421,-0.04102565, 0, 0.08423283,-0.04048879,0,0,0.09659707), ncol=3, byrow=T)
R = matrix( c( 2.231113e-05,-4.816095e-05,-5.115287e-05, 0,2.989584e-05,4.212173e-06,0,0, 5.504990e-05), ncol=3, byrow=T)
d = as.matrix(c(57621264,78057622,171342351),ncol=1)
start <- rep(1/(ncol(R)+1), ncol(R))
min_fn <- function(b, dvec, Dmat) -t(dvec)%*%b +t(b)%*%Dmat%*%b/2
grad_min_fn <- function(b, dvec, Dmat) -dvec + Dmat%*%b
b = c(-1., 0, rep(0,ncol(R)))
C = cbind(rep(1,ncol(R)), diag(ncol(R)))
C = cbind(-rep(1,ncol(R)),C)
D <- t(solve(R))%*%solve(R)
constrOptim(theta=start, f=min_fn, grad=grad_min_fn, ui=t(C), ci=b, control=list(reltol=10*.Machine$double.eps),
dvec=d, Dmat=D )
solve.QP solution for R2
$solution
[1] -1.025463e-10 0.000000e+00 1.000000e+00
constrOptim solution for R2
$par
[1] 2.479855e-15 1.178762e-14 1.000000e+00
以及给出 R
的解决方案 $par
[1] 9.272547e-17 5.958225e-14 1.040137e-01
constrOptim 比 solve.QP 提供了更多有关其解决方案路径的信息,因此可能对数字敏感问题更有帮助。
我不知道 constrOptim 是否可以代替 solve.QP 为您工作,但至少它可以帮助您调查问题解决方案的属性。
你的矩阵 dvec 中的元素非常大。您的矩阵 R 是 Dmat
的 cholesky 分解中上三角矩阵的逆矩阵,即 Dmat = M^T M,其中 M = R^{-1}。所以以未分解的形式:
M <- solve(R)
Dmat <- t(M)%*%M
Dmat
unfactorized 很大,与 dvec
:
# [,1] [,2] [,3]
#[1,] 2008893283 3236243201 1619059085
#[2,] 3236243201 6332319710 2522625866
#[3,] 1619059085 2522625866 1641403882
因此您的问题可能是由于某种溢出错误引起的。要解决此问题,您可以缩放 Dmat
(未分解)和 dvec
:
sc <- norm(Dmat,"2")
solve.QP(Dmat = Dmat/sc, dvec=d/sc, Amat=C, bvec=b, meq=0, factorized=FALSE )
解决方案是
# $solution
# [1] -1.220832e-17 0.000000e+00 1.043877e-01
这与 constrOptim
在另一个答案中找到的解决方案非常匹配。
缩放 Dmat
和 dvec
很好,因为 -d^T b + 1/2 b^T D b
的约束最小值与任何常量 sc 的 sc*(-d^T b + 1/2 b^T D b)
的约束最小值相同。
编辑: 要以分解形式解决问题,您可以尝试以下缩放比例:
nn2 = sqrt(norm(d,"2"))
H = solve.QP(Dmat = R*nn2, dvec = d/(nn2^2), Amat = C, bvec = b, factorized=TRUE)
#$solution
#[1] 0.0000000 0.0000000 0.1043877