Powell method - Multidimensional Optimization

powell.m

% Moved from: http://pastebin.com/YD8FGFL8
function out = powel(f, x0, e)
    [directions, ~] = size(x0);
    u = eye(directions);
    
    t(:, 1) = x0;
    xNew = x0;
    
    while true
        xOld = xNew
        
        % Pravim u(i)num za i = 1..broj pravaca
        for i = 1:directions
            teta = parabola(f, u(:, i), t(:, i), e);
            t(:, i + 1) = t(:, i) + teta * u(:, i);
        end
        
        % Pravim u(broj pravaca + 1)
        u(:, directions + 1) = t(:, directions + 1) - t(:, 1);
        teta = parabola(f, u(:, directions + 1), t(:, directions + 1), e);
        xNew = t(:, directions + 1) + teta * u(:, directions + 1);
        t(:, 1) = xNew;
        
        for i = 1:directions
            u(:, i) = u(:, i + 1);
        end
        
        if abs(xOld(1) - xNew(1)) < e
            out = xNew;
            break
        end
    end
end


function xopt = parabola(f, direction, add, e)
    x = [0 0.5 1];

    while true
        mata = [1 x(1) x(1).^2;
            1 x(2) x(2).^2;
            1 x(3) x(3).^2
        ];
        matb = [f(x(1)*direction + add); 
            f(x(2)*direction + add); 
            f(x(3)*direction + add)
        ];
    
        abc = linsolve(mata, matb);
        xopt = -abc(2)/(2*abc(3));
              
        x(4) = xopt;
        [vals indexes] = sort([f(direction*x(1) + add) 
            f(direction*x(2) + add) 
            f(direction*x(3) + add) 
            f(direction*x(4) + add)]);
        x = [x(indexes(1)) x(indexes(2)) x(indexes(3))];
        
        if abs(f(direction*x(1) + add) - f(direction*xopt +add)) < e
            break;
        end
    end
end