How to set parameters' sum to 1 in constrained optimization

Question

Here's the code (I'm sorry if it's so long, but it was the first example I had); I'm using the CVaR example from CreditMetrics package by A. Wittmann and DEoptim solver to optimize:

library(CreditMetrics)
library(DEoptim)

N <- 3
n <- 100000
r <- 0.003
ead <- rep(1/N,N)
rc <- c("AAA", "AA", "A", "BBB", "BB", "B", "CCC", "D")
lgd <- 0.99
rating <- c("BBB", "AA", "B")   
firmnames <- c("firm 1", "firm 2", "firm 3")
alpha <- 0.99

# correlation matrix
rho <- matrix(c(  1, 0.4, 0.6,
                  0.4,   1, 0.5,
                  0.6, 0.5,   1), 3, 3, dimnames = list(firmnames, firmnames),
              byrow = TRUE)

# one year empirical migration matrix from standard&poors website
rc <- c("AAA", "AA", "A", "BBB", "BB", "B", "CCC", "D")
M <- matrix(c(90.81,  8.33,  0.68,  0.06,  0.08,  0.02,  0.01,   0.01,
              0.70, 90.65,  7.79,  0.64,  0.06,  0.13,  0.02,   0.01,
              0.09,  2.27, 91.05,  5.52,  0.74,  0.26,  0.01,   0.06,
              0.02,  0.33,  5.95, 85.93,  5.30,  1.17,  1.12,   0.18,
              0.03,  0.14,  0.67,  7.73, 80.53,  8.84,  1.00,   1.06,
              0.01,  0.11,  0.24,  0.43,  6.48, 83.46,  4.07,   5.20,
              0.21,     0,  0.22,  1.30,  2.38, 11.24, 64.86,  19.79,
              0,     0,     0,     0,     0,     0,     0, 100
)/100, 8, 8, dimnames = list(rc, rc), byrow = TRUE)

cm.CVaR(M, lgd, ead, N, n, r, rho, alpha, rating)

y <- cm.cs(M, lgd)[which(names(cm.cs(M, lgd)) == rating)]

Now I write my function...

fun <- function(w) {
  # ... 
  - (t(w) %*% y - r) / cm.CVaR(M, lgd, ead = w, N, n, r, 
                           rho, alpha, rating)
}

...and I want to optimize it:

DEoptim(fn = fun, lower = rep(0, N), upper = rep(1, N), 
        control = DEoptim.control())

Can you tell me what do I have to insert in # ... to make sum(w) = 1 during optimization?

Below I show you optimization results according to flodel's tips:

# The first trick is to include B as large number to force the algorithm to put sum(w) = 1

fun <- function(w) {
  - (t(w) %*% y - r) / cm.CVaR(M, lgd, ead = w, N, n, r, rho, alpha, rating) + 
    abs(10000 * (sum(w) - 1))
}

DEoptim(fn = fun, lower = rep(0, N), upper = rep(1, N), 
        control = DEoptim.control())

$optim$bestval
[1] -0.05326055

$optim$bestmem
par1        par2        par3 
0.005046258 0.000201286 0.994752456

parsB <- c(0.005046258, 0.000201286, 0.994752456)

> fun(parsB)
            [,1]
[1,] -0.05326089

...and...

As you can see, the first trick works better in that he finds a results which is smaller than the second one. Unfortunately it seems he takes longer.

# The second trick needs you use w <- w / sum(w) in the function itself

fun <- function(w) {
  w <- w / sum(w)
  - (t(w) %*% y - r) / cm.CVaR(M, lgd, ead = w, N, n, r, rho, alpha, rating) #+ 
    #abs(10000 * (sum(w) - 1))
}

DEoptim(fn = fun, lower = rep(0, N), upper = rep(1, N), 
        control = DEoptim.control())

$optim$bestval
[1] -0.0532794

$optim$bestmem
par1         par2         par3 
1.306302e-15 2.586823e-15 9.307001e-01

parsC <- c(1.306302e-15, 2.586823e-15, 9.307001e-01)
parC <- parsC / sum(parsC)

> fun(parC)
           [,1]
[1,] -0.0532794

Any comment?

Should I increase the number of iterations because of a "too-stochastic" to-be-optimized-function?

Answer 1

Try:

w <- w / sum(w)

and if DEoptim gives you an optimal solution w* such that sum(w*) != 1 then w*/sum(w*) should be your optimal solution.

Another approach is to solve over all your variables but one. We know the value of the last variable must be 1 - sum(w) so in the body of the function, have:

w <- c(w, 1-sum(w))

and do the same to the optimal solution returned by DEoptim: w* <- c(w*, 1-sum(w*))

Both solutions require that you re-formulate your problem into an unconstrained (not counting for variable bounds) optimization so DEoptim can be used; which forces you to do a little extra work outside of DEoptim to recover the solution to the original problem.

In reply to your comment, if you want DEoptim to give you the correct answer right away (i.e. without the need for a post-transformation), you could also try to include a penalty cost to your objective function: for example add B * abs(sum(w)-1) where B is some arbitrary large number so sum(w) will be forced to 1.