Code covered by the BSD License  

Highlights from
Unbounded Resolution for Function Approximation

image thumbnail

Unbounded Resolution for Function Approximation

by

 

A continuous function is optimized for varying outputs and increasing parameterization dimension

[Out_out,OutGrad]=OutputFunc_lwr_8(Pnt,Opt)
function [Out_out,OutGrad]=OutputFunc_lwr_8(Pnt,Opt)
global Memory_all
global Support_Flag Support_Flag_indy % goes to one when failed to find support

w=FindWeights(Opt.KernalWeight,Memory_all.Pop,Pnt,Opt.h);
temp_com=combnk(1:size(Memory_all.Pop,2),2);
Pop_now=[Opt.BasisScale_0*ones(sum(w~=0),1),...
    bsxfun(@minus,Memory_all.Pop(w~=0,:),Pnt)/(Opt.h*Opt.BasisScale_1),...
    ((bsxfun(@minus,Memory_all.Pop(w~=0,:),Pnt)/Opt.h).^2+Opt.BasisOff_2)/Opt.BasisScale_2,...
    ((bsxfun(@minus,Memory_all.Pop(w~=0,temp_com(:,1)),Pnt(temp_com(:,1))).*...
    bsxfun(@minus,Memory_all.Pop(w~=0,temp_com(:,2)),Pnt(temp_com(:,2))))/(Opt.h^2)+Opt.BasisOff_2)/Opt.BasisScale_2];
W=diag(w(w~=0));
temp_limit=0;
[dummy,S,V]=svd(Pop_now.'*(W.'*W)*Pop_now);
%         fprintf(['Initial eigenvalues',repmat(' %g',1,size(S,1)),'\n'],diag(S));
% while cond(Pop_now.'*(W.'*W)*Pop_now)>Opt.MinCond
while min(diag(S))<Opt.MinEig
    %     SupportNeeded=max(MinCond-diag(S),0)/(KernalSupportNew_mean*NewPointBatch*h);
    %         New_Pop=bsxfun(@plus,KernalSupportPDF_fun(rand(NewPointBatch,size(Population,2)))*...
    %             diag(SupportNeeded)*(V.'),Pnt);
%     New_Pop=bsxfun(@plus,Opt.KernalSupportPDF_fun(rand(Opt.NewPointBatch,length(Pnt)))*...
%         diag(min((diag(S)<=Opt.MinCond)./(diag(S)+~(diag(S)<=Opt.MinCond)),Opt.h))*(V.'),Pnt);
% reallly dumb gausian, see CostFunc_lwr_8 for more thoughts
%     New_Pop=bsxfun(@plus,randn(Opt.NewPointBatch*sum(diag(S)<=max(diag(S))/Opt.MinCond),length(Pnt))*Opt.h/2,Pnt);
%     New_Pop=PSO_ConditionNum(Pop_now.'*(W.'*W)*Pop_now,length(Pnt),Opt);
    New_Pop=bsxfun(@plus,Fmin_ConditionNum(Pop_now.'*(W.'*W)*Pop_now,length(Pnt),Opt),Pnt);
    New_Out=zeros(Opt.NewPointBatch,1);
    New_Cost=zeros(Opt.NewPointBatch,1);
    for jackie=1:size(New_Pop,1)
        New_Out(jackie)=Opt.OutputFunc_true(New_Pop(jackie,:));
        New_Cost(jackie)=Opt.CostFunc_true(New_Pop(jackie,:));
    end
    Memory_all.Pop=[Memory_all.Pop;New_Pop]; 
    Memory_all.Out=[Memory_all.Out;New_Out]; 
    Memory_all.Cost=[Memory_all.Cost;New_Cost]; 
    w=FindWeights(Opt.KernalWeight,Memory_all.Pop,Pnt,Opt.h);
Pop_now=[Opt.BasisScale_0*ones(sum(w~=0),1),...
    bsxfun(@minus,Memory_all.Pop(w~=0,:),Pnt)/(Opt.h*Opt.BasisScale_1),...
    ((bsxfun(@minus,Memory_all.Pop(w~=0,:),Pnt)/Opt.h).^2+Opt.BasisOff_2)/Opt.BasisScale_2,...
    ((bsxfun(@minus,Memory_all.Pop(w~=0,temp_com(:,1)),Pnt(temp_com(:,1))).*...
    bsxfun(@minus,Memory_all.Pop(w~=0,temp_com(:,2)),Pnt(temp_com(:,2))))/(Opt.h^2)+Opt.BasisOff_2)/Opt.BasisScale_2];
    W=diag(w(w~=0));
    temp_limit=temp_limit+1;
    [dummy,S,V]=svd(Pop_now.'*(W.'*W)*Pop_now);
%     fprintf(['%g pts, %g gen, Eigenval: ',repmat(' %g',1,size(S,1)),'\n'],size(Memory_all.Pop,1),temp_limit,diag(S));
    %         [Pop_now(:,2:end).'*(W.'*W)*Pop_now(:,2:end), S,V,diag(min((diag(S)<=MinCond)./(diag(S)+~(diag(S)<=MinCond)),h))]
    if temp_limit>200
        fprintf('Failed to find sufficient data to support estimation\n')
        try %#ok<TRYNC>
        Support_Flag(Support_Flag_indy)=1;
        end
        break
    end
end
Out_now=Memory_all.Out(w~=0,:);
Pop_now=[Opt.BasisScale_0*ones(sum(w~=0),1),...
    bsxfun(@minus,Memory_all.Pop(w~=0,:),Pnt)/(Opt.h*Opt.BasisScale_1),...
    ((bsxfun(@minus,Memory_all.Pop(w~=0,:),Pnt)/Opt.h).^2+Opt.BasisOff_2)/Opt.BasisScale_2,...
    ((bsxfun(@minus,Memory_all.Pop(w~=0,temp_com(:,1)),Pnt(temp_com(:,1))).*...
    bsxfun(@minus,Memory_all.Pop(w~=0,temp_com(:,2)),Pnt(temp_com(:,2))))/(Opt.h^2)+Opt.BasisOff_2)/Opt.BasisScale_2];
Gain_now=pinv(Pop_now.'*(W.'*W)*Pop_now)*Pop_now.'*(W.'*W);
Out_out=[Opt.BasisScale_0*1;0*Pnt(:);...
    Opt.BasisOff_2/Opt.BasisScale_2+zeros(size(temp_com,1)+length(Pnt),1)...
    ].'*Gain_now*Out_now;
OutGrad=[zeros(length(Pnt),1),eye(length(Pnt))/(Opt.h*Opt.BasisScale_1),...
    zeros(length(Pnt),size(temp_com,1)+length(Pnt))]*Gain_now*Out_now;

Contact us