[x,e]=ekfopt(h,x,tol,P,Q,R)

function [x,e]=ekfopt(h,x,tol,P,Q,R)
%EKFOPT     Unconstrained optimization using the extended Kalman filter
%
%       [x,e]=ekfopt(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 extended Kalman filter with the unscented 
% 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 200 iterations
% f=@(x)100*(x(2)-x(1)^2)^2+(1-x(1))^2;
% x=ekfopt(f,[2;-1],1e-9,0.5*eye(2),eye(2),1e-9)
%
% 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-7*eye(n);
% x=ekfopt(f,x0,tol,P,Q,R);
%
%
% By Yi Cao at Cranfield University, 10 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
[m1,m2]=size(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(m1,m2); %target vector
ne=norm(e);
while ne>tol
    [x,P]=ekf(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