nm.js

//Implementation of Nelder-Mead Simplex Linear Optimizer
//	TODO: Robust Unit Test of 2D Function Optimizations
//	TODO: Extend to support functions beyond the 2D Space

function Simplex(vertices) {
    this.vertices = vertices;
    this.centroid = null;
    this.reflect_point = null; //Reflection point, updated on each iteration
    this.reflect_cost = null;
    this.expand_point = null;
    this.expand_cost = null;
    this.contract_point = null;
    this.contract_cost = null;
}

//sort the vertices of Simplex by their objective value as defined by objFunc
Simplex.prototype.sortByCost = function (objFunc) {
    this.vertices.sort(function (a, b) {
        var a_cost = objFunc(a), b_cost = objFunc(b);
            
        if (a_cost < b_cost) {
            return -1;
        } else if (a_cost > b_cost) {
            return 1;
        } else {
            return 0;
        }
    });
};

//find the centroid of the simplex (ignoring the vertex with the worst objective value)
Simplex.prototype.updateCentroid = function (objFunc) {
    this.sortByCost(objFunc); //vertices must be in order of best..worst

    var centroid_n = this.vertices.length - 1, centroid_sum = 0, i;
    for (i = 0; i < centroid_n; i += 1) {
        centroid_sum += this.vertices[i];
    }
    
    this.centroid = centroid_sum / centroid_n;
};

Simplex.prototype.updateReflectPoint = function (objFunc) {
    var worst_point = this.vertices[this.vertices.length - 1];
    this.reflect_point = this.centroid + (this.centroid - worst_point); // 1*(this.centroid - worst_point), 1 removed to make jslint happy
    this.reflect_cost = objFunc(this.reflect_point);
};

Simplex.prototype.updateExpandPoint = function (objFunc) {
    var worst_point = this.vertices[this.vertices.length - 1];
    this.expand_point = this.centroid + 2 * (this.centroid - worst_point);
    this.expand_cost = objFunc(this.expand_point);
};

Simplex.prototype.updateContractPoint = function (objFunc) {
    var worst_point = this.vertices[this.vertices.length - 1];
    this.contract_point = this.centroid + 0.5 * (this.centroid - worst_point);
    this.contract_cost = objFunc(this.contract_point);
};

//assumes sortByCost has been called prior!
Simplex.prototype.getVertexCost = function (objFunc, option) {
    if (option === 'worst') {
        return objFunc(this.vertices[this.vertices.length - 1]);
    } else if (option === 'secondWorst') {
        return objFunc(this.vertices[this.vertices.length - 2]);
    } else if (option === 'best') {
        return objFunc(this.vertices[0]);
    }
};

Simplex.prototype.reflect = function () {    
    this.vertices[this.vertices.length - 1] = this.reflect_point; //replace the worst vertex with the reflect vertex
};

Simplex.prototype.expand = function () {
    this.vertices[this.vertices.length - 1] = this.expand_point; //replace the worst vertex with the expand vertex
};

Simplex.prototype.contract = function () {    
    this.vertices[this.vertices.length - 1] = this.contract_point; //replace the worst vertex with the contract vertex
};

Simplex.prototype.reduce = function () {    
    var best_x = this.vertices[0],  a;
    for (a = 1; a < this.vertices.length; a += 1) {
        this.vertices[a] = best_x + 0.5 * (this.vertices[a] - best_x); //0.1 + 0.5(0.1-0.1)
    }
};

function NM(objFunc, x0, numIters) {

	//This is our Simplex object that will mutate based on the behavior of the objective function objFunc
    var S = new Simplex([x0, x0 + 1, x0 + 2]), itr, x;

    for (itr = 0; itr < numIters; itr += 1) {
        
        S.updateCentroid(objFunc); //needs to know which objFunc to hand to sortByCost
        S.updateReflectPoint(objFunc);

        x = S.vertices[0];
        document.write('Iteration: ' + itr + ' x = ' + x + ' cost = ' + objFunc(x) + '<br/>');

        if (S.reflect_cost < S.getVertexCost(objFunc, 'secondWorst') && S.reflect_cost > S.getVertexCost(objFunc, 'best')) {
            S.reflect();
        } else if (S.reflect_cost < S.getVertexCost(objFunc, 'best')) { //new point is better than previous best: expand

            S.updateExpandPoint(objFunc);
           
            if (S.expand_cost < S.reflect_cost) {
                S.expand();
            } else {           
                S.reflect();
            }
        } else { //new point was worse than all current points: contract

            S.updateContractPoint(objFunc);

            if (S.contract_cost < S.getVertexCost(objFunc, 'worst')) {
                S.contract();
            } else {                
                S.reduce();            
            }
        }
    }
}

//function that we are currently trying to minimize: 5(x^4) + 6x + 8
function parabola(x) {
    return 5 * Math.pow(x, 4) + 6 * x + 8;
}

//objective function that Nelder Mead will seek to minimize by mutating the simplex
function parabolicCost(x) {
    var residual = parabola(x); 
    return Math.pow(residual, 2);
}

//TODO: Robust Unit Testing
NM(parabolicCost, 100, 200);
document.write('<br>This particular function is minimized when x = -0.669433 according to <a href="http://www.wolframalpha.com/input/?i=5%28x^4%29+%2B+6x+%2B+8">WolframAlpha</a><br/>');