Как выразить x1 * x2 в квадратике R

Question

У меня вопрос о том, как выразить x1x2 в целевой функции. Вот пример в Интернете.

## 
## min x1^2 +2x2^2 + 4x3^2 - x1 - x2 + 5x3 
## x1 + x3 <= 1 
## x1 >= 5 
## x2 <= 0 
## 

P = 2*diag (c (1, 2, 4)); 
d = c (-1, -1, 5); 
A = matrix (0, nrow=3, ncol=3); 
A[1,] = c(-1, 0, -1); 
A[2,] = c(1, 0, 0); 
A[3,] = c(0, -1, 0); 
b = c(-1, 5, 0); 

Из примера, цель finction является x1^2 +2x2^2 + 4x3^2 - x1 - x2 + 5x3

В R, это P = 2*diag (c (1, 2, 4)); d = c (-1, -1, 5);

Однако, если у меня есть целевая функция, как x1x2 или x1^2 Как ввести команды в R?

Большое спасибо.

Answer 1

solve.QP функция минимизирует следующее выражение:

min 1/2 * (t(x) %*% D %*% x) - t(d) %*% x 

, где х представляет собой вектор n переменных (например, x1, x2), ... D симметричной матрицы коэффициентов n x n и d вектора коэффициентов n ,

Для того, чтобы понять, как эти матричные продукты переводить выражение целевой функции, я создал следующую простую функцию для расширения цели 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 

нам нужно установить D2=1 and D4=1 (не только D2=2 поскольку матрицы D должны быть симметричными!), А все остальные коэффициенты = 0, следовательно:

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"