Максимизируйте квадратичный объект с линейными ограничениями в R: Rsolnp или Auglag

Question

Я пытаюсь найти решение следующей задачи оптимизации, используя либо auglag, либо Rsolnp.

Max t(w1 - w2) * Kf * Sf * t(Kf) * (w1 - w2) 
subject to Kc * w1 = Kc * w2 
and sum(w1) = 1 and sum(w2) = 1 and w1,w2 >= 0 
Sc and Sf are variance covariance matrices at the coarse and fine level respectively. 
Kc and Kf are exposure matrices as the coarse and fine level respectively. 
Nc and Nf are nodes at which exposure nodes at the coarse and fine level. 

Это эффективно пытается найти WTS двух портфелей w1 и w2, что бы максимизировать ТэВ более тонкого уровня экспозиции, при условии сумме WTS = 1 и все WTS> 0. Существует еще одно ограничение равенства слишком (Это эффективно означает, что экспозиция на грубом уровне идентична для двух портфелей). Rsolnp не может максимизировать и возвращает решение, в котором целевая функция равна 0, auglag полностью взрывается и не отвечает ограничениям с рядом предупреждений.

Может кто-нибудь, пожалуйста, помогите мне понять, где я иду не так?

seqFineNodes <- c(1, 2, 3, 4, 5, 6) 

Nc <- c(2, 3, 5) 

Kc <- matrix(c(0.2481316799436,0.495478766935844,0,0,0,0,0,0,0.743360061619584,0.497321712603124,0,0,0,0,0,0.497321712603124,1.23913608908603,1.48240730986596), nrow=length(seqFineNodes), ncol=length(Nc)) 
dimnames(Kc) <- list(as.character(seqFineNodes), as.character(Nc)) 

Sc <- matrix(c(619.806079280659,627.832850585004,549.805085990891,627.832850585004,668.726833059322,624.524848194842,549.805085990891,624.524848194842,696.498483673357), nrow=length(Nc), ncol=length(Nc)) 
dimnames(Sc) <- list(as.character(Nc), as.character(Nc)) 

Nf <- c(2, 3, 4, 5) 

Kf <- matrix(c(0.2481316799436,0.495478766935844,0,0,0,0,0,0,0.743360061619584,0,0,0,0,0,0,0.994643425206249,0,0,0,0,0,0,1.23913608908603,1.48240730986596), nrow=length(seqFineNodes), ncol=length(Nf)) 
dimnames(Kf) <- list(as.character(seqFineNodes), as.character(Nf)) 

Sf <- matrix(c(619.806079280659,627.832850585004,602.504944834256,549.805085990891,627.832850585004,668.726833059322,666.196728425214,624.524848194842,602.504944834256,666.196728425214,696.688027074344,681.064062606848,549.805085990891,624.524848194842,681.064062606848,696.498483673357), nrow=length(Nf), ncol=length(Nf)) 
dimnames(Sf) <- list(as.character(Nf), as.character(Nf)) 

KRD_fine <- Kf 
KRD_coarse <- Kc 
VC_fine <- Sf 
VC_coarse <- Sc 
countw <- length(seqFineNodes) 


t1 <- diag(x = 1, nrow = countw, ncol = countw) 
t2 <- diag(x = -1, nrow = countw, ncol = countw) 
tr <- cbind(t1,t2) 

D_fine <- t(tr) %*% KRD_fine %*% VC_fine %*% t(KRD_fine) %*% tr 
#round(eigen(Dmat)$values, 4) 
D_fine <- as.matrix(nearPD(D_fine)$mat) 
#round(eigen(Dmat)$values, 4) 

eq_coarse_krd_A <- t(KRD_coarse) %*% tr 
eq_coarse_krd_b <- rep(0, nrow(VC_coarse)) 

# Equality constraints 
eq_A1 <- c(rep(1, countw), rep(0,countw)) 
eq_A2 <- c(rep(0, countw), rep(1,countw)) 
eq_b <- c(1 , 1) 

# Constraint wts greater than zero 
ineq_A <- diag(x = 1, nrow = 2 * countw, ncol = 2 * countw) 
ineq_b <- rep(0, 2 * countw) 

# Combine constraints 
heq <- rbind(eq_coarse_krd_A, eq_A1, eq_A2) 
beq <- c(eq_coarse_krd_b, eq_b) 

hin <- ineq_A 

theta <- c(1, rep(0, countw - 1), 1, rep(0, countw - 1)) 

krdsol <- solnp(par = theta, 
       fun = function(x) -c(t(x) %*% D_fine %*% x), 
       ineqfun = function(x) c(hin %*% x), 
       ineqLB = rep(0, 2 * countw), 
       ineqUB = rep(1, 2 * countw), 
       eqfun = function(x) c(heq %*% x), 
       eqB = beq) 


krdFine <- auglag(par = theta, 
        fn = function(x) c(t(x) %*% D_fine %*% x), 
        hin = function(x) c(hin %*% x), 
        heq = function(x) c(heq %*% x) - beq, 
        control.outer = list(method = "nlminb"), 
        control.optim=list(fnscale=-1)) 

Answer 1

Я решил решить вашу проблему: solnp. ?solnp говорит fun, ineqfun и eqfun возвращение vector но ваш возврат matrix. Поэтому я добавил c(...) им.

library(Rsolnp) 

krdsol <- solnp(par = theta, 
       fun = function(x) c(-t(x) %*% D_fine %*% x), 
       ineqfun = function(x) c(hin %*% x), 
       ineqLB = rep(0, 2 * countw), 
       ineqUB = rep(1, 2 * countw), 
       eqfun = function(x) c(heq %*% x), 
       eqB = beq) 

[Изменено]

Элементы auglag(control.optim=list(...)) принимает в качестве аргументов перечислены в ?nlminb() (и увидеть ?auglag())