solve.QP
関数は、次式を最小化:
min 1/2 * (t(x) %*% D %*% x) - t(d) %*% x
xは
n
変数(例えば、X1、X2 ...)、
n x n
係数の
D
対称行列と
n
係数の
d
ベクトルのベクトルである
。これらの行列積は、私が入力与えられたquadprogの目標を拡張するには、次の簡単な関数を作成しました、目的関数式に変換する方法を理解するために
:
# Given matrix D and vector d that will be passed to solve.QP,
# it returns the objective function expression as string
expandQPObj <- function(D, d, var.name = 'x') {
# helper function that expands a matrix product. NB it does not compute the result,
# it just returns the expressions of each element of the resulting matrix
expandMatrixProd <- function(m1, m2) {
m1 <- as.matrix(m1)
m2 <- as.matrix(m2)
if (ncol(m1) != nrow(m2)) { stop("incompatible dimensions") }
n <- nrow(m1)
m <- ncol(m2)
res <- matrix('', nrow = n, ncol = m)
for (i in 1:n) {
for (j in 1:m) {
a <- m1[i, ]
b <- m2[, j]
a <- ifelse(grepl('[+*]', a), paste0('(', a, ')'), a)
b <- ifelse(grepl('[+*]', b), paste0('(', b, ')'), b)
res[i, j] <- gsub('+-','-',paste(a, b, sep = '*', collapse = '+'),fixed = TRUE)
}
}
return(res)
}
D <- as.matrix(D)
d <- as.vector(d)
n <- length(d)
if (!all(dim(D) == n) || n == 0) {
stop('Dimensions problem: D should be an nxn matrix and d a vector of length n (n>0)')
}
xvec <- paste(var.name, 1:n, sep = '')
quadComp <- as.vector(Reduce(list(t(xvec), D, xvec),f=expandMatrixProd))
linearComp <- as.vector(strMatrixMult(t(d), xvec))
return(paste0('1/2*(', quadComp, ') - (', linearComp, ')'))
}
リテラル係数でそれを呼び出すことによって、我々得るために、そう
Dliteral <- matrix(paste0('D',1:9),nrow=3,byrow = T)
# [,1] [,2] [,3]
#[1,] "D1" "D2" "D3"
#[2,] "D4" "D5" "D6"
#[3,] "D7" "D8" "D9"
dliteral <- paste0('d',1:3)
#[1] "d1" "d2" "d3"
expandQPObj(Dliteral,dliteral)
#"1/2*((x1*D1+x2*D4+x3*D7)*x1+(x1*D2+x2*D5+x3*D8)*x2+(x1*D3+x2*D6+x3*D9)*x3) - d1*x1+d2*x2+d3*x3"
:簡単に私たちの目的の結果を得るために、これらの係数を設定する方法を理解することができます
min x1*x2
(マトリックス
D
は対称でなければならないためだけでなく、
D2=2
!)、従って他のすべての係数= 0、我々は
D2=1 and D4=1
を設定する必要が
:得るために、同様
D = rbind(c(0,1,0),c(1,0,0),c(0,0,0))
# [,1] [,2] [,3]
#[1,] 0 1 0
#[2,] 1 0 0
#[3,] 0 0 0
d <- rep(0,3)
#[1] 0 0 0
expandQPObj(D,d)
[1] "1/2*((x1*0+x2*1+x3*0)*x1+(x1*1+x2*0+x3*0)*x2+(x1*0+x2*0+x3*0)*x3) - 0*x1+0*x2+0*x3"
:
min x1^2 (== min x1*x1)
D1=2
と他のすべての係数= 0を設定する必要があります。
D = rbind(c(2,0,0),c(0,0,0),c(0,0,0))
# [,1] [,2] [,3]
#[1,] 2 0 0
#[2,] 0 0 0
#[3,] 0 0 0
d <- rep(0,3)
#[1] 0 0 0
expandQPObj(D,d)
# "1/2*((x1*2+x2*0+x3*0)*x1+(x1*0+x2*0+x3*0)*x2+(x1*0+x2*0+x3*0)*x3) - 0*x1+0*x2+0*x3"
BTW、私はなぜP = 2 * diag(c(1,2,4))この関数 "2" diagを掛ける必要がある別の質問を思いますか? – GGININDERRRRRRRRRRRRRRRRRRRRRR
目的関数の非線形部分に1/2を掛けるため、 "2"が必要です。詳細はsolve.QPのドキュメントと私の答えを参照してください。 – digEmAll