Harmony Search Algorithm

Even with only a few parameters controlling a trading system, the optimization search space gets big, quickly. How do I find an optimal set of parameters via back-testing when my poor laptop takes an hour for 500 back-tests? The solution may be some form of genetic search algorithm, in this case the Harmony Search Algorithm.

I believe genetic optimization algorithms may provide a good solution. They work with discrete and continuous variables, they don’t care how “rough” the optimization surface is and they are resistant to being trapped by local maxima. A nice feature is the fact that, aside from the max and min values for each parameter, the algorithm looks where it pleases – no data snooping on my part!

My proposal is that I use a genetic search algorithm implemented in R to create each generation of parameter sets in a file. Read the file into TradingBlox Builder which steps through the parameter sets, writing results out to a file which will be passed back to the genetic algorithm. Repeat…

The process of setting up TBB to do stepped parameter tests outside of the built-in method will be covered in a another post, here I will outline my first experience with a simple genetic algorithm: “The Harmony Search Algorithm”

The Harmony Search Algorithm

Harmony is simple:

1. Randomly generate ‘harmony memory size’ (hms) parameter sets and arrange them as the rows of a table. Add a column for the results of each parameter set. This collection of information is called the Harmony Memory (HM)
2. Generate a new candidate parameter set based on some random variations on the existing sets in HM
3. Test the new candidate. If it is better than the worst performing parameter set in HM, then replace that one with this.
4. Repeat until some termination criteria is reached.

First thing to note is a BIG disadvantage of this algorithm: only ONE new candidate is generated with each cycle. Despite this, I pressed ahead as the big advantage is simplicity: I could code this up and test it out quickly. I like to learn on a simple system!

Some details on each of the steps above:

Creating HM

You have to decide how large to make HM. I guess that the larger your optimization search space the larger HM should be. If HM is too small the algorithm is more likely to get trapped at a local maxima. Larger HM leads to longer computing time.

The parameter values are just drawn independently using a uniform distribution between the user-specified minimum and maximum values.

Creating New Candidates

Start with the first parameter. With probability “harmony memory choosing rate” (hmcr) choose from amongst the existing values in HM. Otherwise, randomly generate a new value between the user-specified manimum and maximum values.

If the parameter was drawn from HM, then, with probability “pitch adjustment rate” (par), adjust the value previously drawn. Otherwise leave it as is. The adjustment depends on the whether the variable is discrete or continuous. If discrete adjust randomly up or down by one unit (“$delta $”). If continuous, adjust by the “fret width” (fw) multiplied by a random number between -1 and +1.

Here, I did something different. I think fw ought to bear some relation to the value of a parameter. I defined fw as a percent of the minimum to maximum range of the parameter (adjustment = fw x runif(-1, 1) x (min – max)). I could just as easily justify using fw as a percentage of the value of the selected parameter as well.

Repeat this process until you have constructed a complete set for testing.

Finally, there is the issue of accidentally duplicating parameter sets in HM. Not likely with multiple parameters with continuous values, but possible. I don’t want to test an already represented parameter set, so I repeat the parameter set. If I am unable to find a non-duplicate set I bail after so many attempts.

I put all this into a function in R called “new.candidate” which takes as input the current state of “harmony” and returns a new candidate parameter set.

Termination Criteria

I have to decide when to stop looking for better values. I use a simple decision process requiring at least a fixed fraction of the parameter set members of HM to agree within a fixed number of significant digits. I then report the current state of “harmony” along with the average values for the columns and the values that maximize the objective function.

I also implement a total cycle limit.

The Code

The following is the code for the heart of the algorithm: new.candidate.

new.candidate <- function(harmony){
  dimensions <- harmony$dimensions
  hms <- harmony$hms
  hmcr <- harmony$hmcr
  HM <- harmony$HM
  par <- harmony$par
  fw <- harmony$fw
  x.max <- harmony$xmax
  x.min <- harmony$xmin
  sig.dig <- harmony$sigDig
  x.new <- vector(mode="numeric", length=dimensions)
  dupe.test <- T
  loop.limit <- 50
  loop.count <- 0
  while(dupe.test){
    loop.count <- loop.count + 1
    if (loop.count >= loop.limit) stop("Can't find new candidate")
    for (i in 1:dimensions){
      if(runif(1, 0, 1) <= hmcr){
        x.new[i] <- HM[ceiling(runif(1, 0, hms)), i]
        if (runif(1, 0, 1) <= par){
          x.new[i] <- x.new[i] + 
            fw * (x.max[i] - x.min[i]) * runif(1, -1, 1)
          x.new[i] <- max(x.min[i], min(x.max[i], x.new[i]))
        }
      }else{
        x.new[i] <- runif(1, x.min[i], x.max[i])
      }
      x.new <- signif(x.new, sig.dig)
    }
    dupe.test <- (sum(duplicated(rbind(x.new, HM[ , 1:dimensions]))) != 0)
  }
  return(x.new)
}

Performance

The Search Optimization Torture Test

I came up with a function to create a surface with multiple local maxima and minima.  It is essentially two sine wave pulses, one in the x-direction [-1, 1] the other in the y-direction [-1, 1]. It is depicted in the 3-D (persp(x, y, z, theta=20, phi=30)) chart.

Harmony Search Parameters

I chose (using intuition) to set the Harmony Search parameters as follows:

• hms (harmony memory size) = 20
• hmcr (harmony memory change rate) = 75%
• par (pitch adjustment rate) = 25%
• fw (fret width) = 25%
• x.min (minimum values for each parameter) = c(-2, -2)
• x.max (maximum values for each parameter) = c(2, 2)
• max.count (maximum number of generations) = 5,000
• agree.dig (significant figures that qualify as “agreement”) = 3
• quorum (number of parameter sets that must agree for termination) = 50%

The “Successful Trials” plot shows the parameter sets that improved on existing sets in HM and thus displaced existing, inferior sets. You can clearly see clustering of successful solutions around the local maxima, with the majority homing in on the global maximum.

Pretty cool. My next post will explore the effect of adjusting the search parameters.