A simple way to optimize

SeanS6

You can try out this simple optimization algorithm if you like:

// Continuous Gray code optimization
//http://citeseerx.ist.psu.edu/viewdoc/download;jsessionid=14BEBFED079E983F49AD75C3CA4AFF82?doi=10.1.1.100.6090&rep=rep1&type=pdf
// The most simple way to evolve solutions to engineering or other problems.  
// And effective too.
void setup() {
  int iterations=100000;
  int precision=1000;
  int vecLen=50;
  float range=5.12f;

  float[] parent=new float[vecLen];
  float[] child=new float[vecLen];
  for (int i=0; i<parent.length; i++) {
    parent[i]=random(-range, range);  // random from -5.12 to 5.12
  }
  float parentCost=deYong(parent);
  for (int i=0; i<iterations; i++) {
    for (int j=0; j<parent.length; j++) {
      float mutation=exp(-precision*random(0f, 1f));
        if (random(-1f, 1f)<0) {
        mutation=-mutation;
      }
      mutation*=2*range;
      float c=parent[j]+mutation;
      if (c>range) c=parent[j];
      if (c<-range) c=parent[j];
      child[j]=c;
    }
    float childCost=deYong(child);
    if (childCost<parentCost) {
      parentCost=childCost;
      for (int j=0; j<parent.length; j++) {
        parent[j]=child[j];
      }
    }
  }
  println("Parent Cost:"+parentCost);
  for(int i=0;i<parent.length;i++){
    println("Vector: "+parent[i]);
  }
}

float deYong(float[] x){
  float result=0f;
  for(int i=0;i<x.length;i++){
    result+=x[i]*x[i];
  }
  return result;
}

SeanS6

Result:
Parent Cost:1.0998695E-24
Vector: -4.446074E-15
Vector: -1.3206545E-14
Vector: 2.5333451E-14
Vector: 1.3092901E-13
Vector: 4.1820844E-16
Vector: -1.5196274E-13
Vector: 2.5448588E-14
Vector: -5.298662E-15
Vector: -2.8565362E-15
Vector: 2.7322212E-13
Vector: 6.90185E-16
Vector: 2.670862E-16
Vector: -6.143601E-13
Vector: 4.357555E-13
Vector: 3.0938024E-14
Vector: -5.0206584E-13
Vector: 1.3262619E-14
Vector: -9.462891E-14
Vector: 2.920175E-15
Vector: -5.9401404E-15
Vector: 2.101738E-14
Vector: -2.2861597E-15
Vector: 3.104379E-14
Vector: -1.4589519E-15
Vector: 5.691257E-17
Vector: -1.7216068E-14
Vector: -3.290261E-14
Vector: 5.133789E-14
Vector: 5.4137467E-14
Vector: 2.9214705E-14
Vector: 1.9756596E-13
Vector: 2.437591E-15
Vector: 2.0602704E-15
Vector: -2.8701068E-14
Vector: -1.2330929E-15
Vector: -1.6886144E-14
Vector: 1.6468369E-14
Vector: 3.1314783E-15
Vector: 3.2038975E-15
Vector: -9.365162E-16
Vector: -1.5792505E-15
Vector: 8.157097E-15
Vector: -9.7503885E-15
Vector: 7.7525444E-15
Vector: 1.015922E-14
Vector: -4.8009404E-15
Vector: -2.1792574E-16
Vector: 2.0833282E-14
Vector: 3.2158198E-13
Vector: 9.1788065E-15

SeanS6

More of the same:

    Xor128 rnd=new Xor128();
    float[] parent=new float[2];
    float[] child=new float[2];
    float parentCost;
    int precision=100;
    int x1, x2, y1, y2;

    void setup() {
      background(0);
      size(512, 512);
      frameRate(60);
      for (int i=0; i<512; i++) {
        for (int j=0; j<512; j++) {
          parent[0]=i/256f-1f;
          parent[1]=j/256f-1f;
          stroke(boothFunction(parent)*2f);
          point(i, j);
        }
      }
      parent[0]=rnd.nextFloatSym();  // random starting position
      parent[1]=rnd.nextFloatSym();
      parentCost=boothFunction(parent); // cost at that position
      x1=(int)((parent[0]+1f)*256f);
      y1=(int)((parent[1]+1f)*256f);
      stroke(0f, 255f, 0f);
      ellipse(x1, y1, 3, 3);
    }

    void draw() {
      child[0]=rnd.mutateXSym(parent[0], precision);
      child[1]=rnd.mutateXSym(parent[1], precision);
      float childCost=boothFunction(child);
      if (childCost<parentCost) {
        parentCost=childCost;
        parent[0]=child[0];
        parent[1]=child[1];
        x2=(int)((parent[0]+1f)*256f);
        y2=(int)((parent[1]+1f)*256f);
        stroke(0f, 255f, 0f);
        line(x1, y1, x2, y2);
        x1=x2;
        y1=y2;
        println("Cost:"+parentCost+" X:"+10f*parent[0]+" Y:"+10f*parent[1]);
      }
    }

    // Min. 0 at 1,3 over x,y domain -10 to 10
    float boothFunction(float[] vec) {
      float x=10f*vec[0];  // get correct domain
      float y=10f*vec[1];
      float a=x+2f*y-7f;
      float b=2f*x+y-5f;
      return a*a+b*b;
    }

    //http://xoroshiro.di.unimi.it/
    class Xor128 {

      private long s0;
      private long s1;

      Xor128() {
        setSeed(System.nanoTime());
      }

      public long nextLong() {
        final long s0 = this.s0;
        long s1 = this.s1;
        final long result = s0 + s1;
        s1 ^= s0;
        this.s0 = Long.rotateLeft(s0, 55) ^ s1 ^ s1 << 14;
        this.s1 = Long.rotateLeft(s1, 36);
        return result;
      }

      public float nextFloat() {
        return (nextLong()&0x7FFFFFFFFFFFFFFFL)*1.0842021e-19f;
      }

      public float nextFloatSym() {
        return nextLong()*1.0842021e-19f;
      }

      public boolean nextBoolean() {
        return nextLong() < 0;
      }

      // Mutation between -1 and 1
      public float mutate(long precision) {
        long ra=nextLong();
        int e=126-(int)(((ra>>>32)*precision)>>>32);
        if (e<0) return 0f;
        return Float.intBitsToFloat((e<<23)|((int)ra&0x807fffff));
      }

      // For parameters x>=0, x<=1
      public float mutateX(float x, long precision) {
        float mx=x+mutate(precision);
        if (mx>1f) return x;
        if (mx<0f) return x;
        return mx;
      }

      // For parameters x>=-1, x<=1
      public float mutateXSym(float x, long precision) {
        float mx=x+2f*mutate(precision);
        if (mx>1f) return x;
        if (mx<-1f) return x;
        return mx;
      }

      public void setSeed(long seed ) {
        s0 = seed;
        s1 = seed+0x9E3779B97F4A7C15L;
        nextLong();
      }
    }

Chrisir

Thanks!

Is it a genetic algorithm?

SeanS6

It is related to genetic algorithms: It's called an evolutionary strategy. However it doesn't use the Gaussian distribution for mutation as most of them do: https://blog.openai.com/evolution-strategies/ To use the mutation operations you choose to have all the parameters to the function you want to optimized to be either between 0 and 1 and use the mutateX function, or to be between -1 and 1 and use the mutateXSym function. The only control parameter to choose is the precision. Maybe pick 30 for simple problems, 100 for more complex problems etc.

SeanS6

For more complex problems you can add a grandparent. Basically you are adding local (or not so local) restarts to the algorithm. Helping move it out of local minimum it could get trapped in.

Xor128 rnd=new Xor128();
float[] grandparent=new float[2];
float[] parent=new float[2];
float[] child=new float[2];
float grandparentCost;
int precision=100;
int iter=5000;
int x1, x2, y1, y2;

void setup() {
  background(0);
  size(512, 512);
  frameRate(600000);
  for (int i=0; i<512; i++) {
    for (int j=0; j<512; j++) {
      grandparent[0]=i/256f-1f;
      grandparent[1]=j/256f-1f;
      stroke(127+trefethenFunction(grandparent)*100f);
      point(i, j);
    }
  }
  grandparent[0]=rnd.nextFloatSym();  // random starting position
  grandparent[1]=rnd.nextFloatSym();
  grandparentCost=trefethenFunction(grandparent); // cost at that position
  x1=(int)((grandparent[0]+1f)*256f);
  y1=(int)((grandparent[1]+1f)*256f);
  stroke(0f, 255f, 0f);
  ellipse(x1, y1, 3, 3);
}

void draw() {
  parent[0]=rnd.mutateXSym(grandparent[0], precision);
  parent[1]=rnd.mutateXSym(grandparent[1], precision);
  float parentCost=trefethenFunction(parent);
  for (int i=0; i<iter; i++) {
    child[0]=rnd.mutateXSym(parent[0], precision);
    child[1]=rnd.mutateXSym(parent[1], precision);
    float childCost=trefethenFunction(child);
    if (childCost<parentCost) {
      parentCost=childCost;
      parent[0]=child[0];
      parent[1]=child[1];
    }
  }
  if (parentCost<grandparentCost) {
    grandparentCost=parentCost;
    grandparent[0]=parent[0];
    grandparent[1]=parent[1];
    x2=(int)((grandparent[0]+1f)*256f);
    y2=(int)((grandparent[1]+1f)*256f);
    stroke(0f, 255f, 0f);
    line(x1, y1, x2, y2);
    x1=x2;
    y1=y2;
    println("Cost:"+grandparentCost+" X:"+10f*grandparent[0]+" Y:"+10f*grandparent[1]);
  }
}

// The global minimum is located at x ∗ = f (−0.024403, 0.210612),
// f(x*) = −3.30686865.
float trefethenFunction(float[] vec) {
  float x=10f*vec[0];  // get correct domain
  float y=10f*vec[1];
  float a=exp(sin(50f*x))+sin(60*exp(y))+sin(70f*sin(x))+sin(sin(80f*y))-sin(10f*(x+y))+.25f*(x*x+y*y);
  return a;
}

//http://xoroshiro.di.unimi.it/
class Xor128 {

  private long s0;
  private long s1;

  Xor128() {
    setSeed(System.nanoTime());
  }

  public long nextLong() {
    final long s0 = this.s0;
    long s1 = this.s1;
    final long result = s0 + s1;
    s1 ^= s0;
    this.s0 = Long.rotateLeft(s0, 55) ^ s1 ^ s1 << 14;
    this.s1 = Long.rotateLeft(s1, 36);
    return result;
  }

  public float nextFloat() {
    return (nextLong()&0x7FFFFFFFFFFFFFFFL)*1.0842021e-19f;
  }

  public float nextFloatSym() {
    return nextLong()*1.0842021e-19f;
  }

  public boolean nextBoolean() {
    return nextLong() < 0;
  }

  // Mutation between -1 and 1
  public float mutate(long precision) {
    long ra=nextLong();
    int e=126-(int)(((ra>>>32)*precision)>>>32);
    if (e<0) return 0f;
    return Float.intBitsToFloat((e<<23)|((int)ra&0x807fffff));
  }

  // For parameters x>=0, x<=1
  public float mutateX(float x, long precision) {
    float mx=x+mutate(precision);
    if (mx>1f) return x;
    if (mx<0f) return x;
    return mx;
  }

  // For parameters x>=-1, x<=1
  public float mutateXSym(float x, long precision) {
    float mx=x+2f*mutate(precision);
    if (mx>1f) return x;
    if (mx<-1f) return x;
    return mx;
  }

  public void setSeed(long seed ) {
    s0 = seed;
    s1 = seed+0x9E3779B97F4A7C15L;
    nextLong();
  }
}