bbukfopt.m

unscented kalman filter for optimization

function [x,e]=ukfopt(h,x,tol,P,Q,R)
%UKFOPT     Unconstrained optimization using the unscented Kalman filter
%
%       [x,e]=ukfopt(f,x,tol,P,Q,R) minimizes e=norm(f(x)) until e<tol. P, Q
%       and R are tuning parmeters relating to the Kalman filter
%       performance. P and Q should be n x n, where n is the number of 
%       decision variables, while R should be m x m, where m is the 
%       dimension of f. Normally, Q and R should be set to d*I e*I where 
%       both d and e are vary small positive scalars. P could be initially 
%       set to a*I where a is the estimated distence between the initial 
%       guess and the optimal valuve. This function can also be used to solve 
%       a set of nonlinear equation, f(x)=0.
%
% This is an example of what a nonlinear Kalman filter can do. It is
% strightforward to replace the unscented Kalman filter with the extended 
% Kalman filter to achieve this same functionality. The advantage of using 
% the unscented Kalman filter is that it is derivative-free. Therefore, it 
% is also suitable for non-analytical functions where no derivatives can 
% be obtained.
%
% Example 1: The Rosenbrock's function is solved within 100 iterations
%{
f=@(x)100*(x(2)-x(1)^2)^2+(1-x(1))^2;
x=ukfopt(f,[2;-1],1e-9,0.45*eye(2))
%}
% Example 2: Solver a set of nonlinear equations represented by MLP NN
%{
rand('state',0);
n=3;
nh=4;
W1=rand(nh,n);
b1=rand(nh,1);
W2=rand(n,nh);
b2=rand(n,1);
x0=zeros(n,1);
f=@(x)W2*tanh(W1*x+b1)+b2;
tol=1e-6;
P=1000*eye(n);
Q=1e-7*eye(n);
R=1e-6*eye(n);
x=ukfopt(f,x0,tol,P,Q,R);
%}
% Example 3: Training a MLP NN. You may have to try several time to get
% results. The best way is to use nnukf.
%{
rand('state',0)
randn('state',0)
N=100;        %training data length
x=randn(1,N); %training data x
y=sin(x);     %training data y=f(x)
nh=4;         %MLP NN hiden nodes
ns=2*nh+nh+1;
f=@(z)y-(z(2*nh+(1:nh))'*tanh(z(1:nh)*x+z(nh+1:2*nh,ones(1,N)))+z(end,ones(1,N)));
theta0=rand(ns,1);
theta=ukfopt(f,theta0,1e-3,0.5*eye(ns),1e-7*eye(ns),1e-6*eye(N));
W1=theta(1:nh);
b1=theta(nh+1:2*nh);
W2=theta(2*nh+(1:nh))';
b2=theta(ns);
M=200;         %Test data length
x1=randn(1,M);
y1=sin(x1);
z1=W2*tanh(W1*x1+b1(:,ones(1,M)))+b2(:,ones(1,M));
plot(x1,y1,'ob',x1,z1,'.r')
legend('Actual','NN model')
%}
%
% By Yi Cao at Cranfield University, 08 January 2008
%

n=numel(x);     %numer of decision variables
f=@(x)x;        %virtual state euqation to update the decision parameter
e=h(x);         %initial residual
m=numel(e);     %number of equations to be solved
if nargin<3,    %default values
    tol=1e-9;
end
if nargin<4
    P=eye(n);
end
if nargin<5
    Q=1e-6*eye(n);
end
if nargin<6
    R=1e-6*eye(m);
end
k=1;            %number of iterations
z=zeros(m,1);   %target vector
ne=norm(e);
while ne>tol
    [x,P]=ukf(f,x,P,h,z,Q,R);               %the unscented Kalman filter
    e=h(x);
    ne=norm(e);                                 %residual
    if mod(k,100)==1
        fprintf('k=%d e=%g\n',k,ne)    %display iterative information
    end
    k=k+1;                                  %iteration count
end
fprintf('k=%d e=%g\n',k,ne)            %final result