Extended Backslash autodiff

Extended Autodiff mldivide function to loop over gradient components to allow gradient propagation of batch likelihood objective function in VGP implementation.

MatlabGP/+BO/FUNBOmax.m

@@ -1,23 +1,37 @@
 function [alpha, dalpha] = FUNBOmax(Z,x)
 
+if nargout>1
+    x=AutoDiff(x);
+end
+
 ys = max(Z.Y);
 
 gamma = 1;
 
 %Calculate std at x
-[varf,dvarf] = Z.eval_var(x);
+[muf] = Z.eval_mu(x);
+[varf] = Z.eval_var(x);
+
+if nargout>1
+    dmuf = getderivs(muf);
+    muf = full(getvalue(muf));
 
-[muf,dmuf] = Z.eval_mu(x);
+    dvarf = getderivs(varf);
+    varf = full(getvalue(varf));
+end
 
 sigf = sqrt(abs(varf));
 dsigf = dvarf./(2*sigf+eps);
 
 %Calculate beta value at x
 beta = (ys - muf + gamma*sigf)/sigf;
-dbeta = -dmuf./sigf + (gamma - beta)*dsigf/sigf;
 
 %Calculate expected improvement over current best measurement
 alpha = -sigf*(beta*normcdf(beta)+normpdf(beta));
-dalpha = -dsigf*alpha/sigf - sigf*dbeta*normcdf(beta);
+
+if nargout>1
+    dbeta = (-dmuf./sigf + (gamma - beta)*dsigf/sigf);
+    dalpha = -1*(-dsigf*alpha/sigf - sigf*dbeta*normcdf(beta));
+end
 
 end

MatlabGP/+utils/contourf.m

@@ -0,0 +1,54 @@
+function contourf(kr,X0,dim1,dim2,varargin)
+
+%   dim1 - Vary this dimension
+%   dim2 - plot along this dimension
+%   f - plot this fidelity
+
+input=inputParser;
+input.KeepUnmatched=true;
+input.addOptional('new_fig',false,@islogical);  % Create a new figure
+input.addOptional('lb',kr.lb_x,@isnumeric);     % Lower bound of plot
+input.addOptional('ub',kr.ub_x,@isnumeric);     % Upper bound of plot
+input.addOptional('cmap','thermal');
+input.addOptional('LS','-');
+input.addOptional('nL',50);
+input.parse(varargin{:})
+in=input.Results;
+
+if in.new_fig
+    figure
+end
+
+a1=in.lb(dim1);
+b1=in.ub(dim1);
+
+
+a2=in.lb(dim2);
+b2=in.ub(dim2);
+
+n = in.nL;
+
+xx1 = linspace(a1,b1,n);
+xx2 = linspace(a2,b2,n);
+
+mb = X0;
+
+for i = 1:n
+    for j = 1:n
+
+        XX = mb;
+        XX(dim1) = xx1(i);
+        XX(dim2) = xx2(j);
+
+        Yj = kr.eval(XX);
+        YY(i,j) = Yj(1);
+
+    end
+end
+
+contourf(xx2,xx1,YY);
+hold on
+shading interp
+utils.cmocean(in.cmap);
+
+end

MatlabGP/AutoDiff/AutoDiff.m

@@ -572,10 +572,22 @@
             if isa(y, 'AutoDiff')
                 if isa(x, 'AutoDiff')
                     z.values = x.values \ y.values;
-                    if size(y, 2) > 1
-                        error('not yet implemented')
+                    % if size(y, 2) > 1
+                    %     error('not yet implemented')
+                    % end
+
+                    for i = 1:size(y.derivatives,2)
+
+                        Xi = reshape(kron(z.values', speye(size(x, 1))) * x.derivatives(:,i),size(y));
+                        Yi = reshape(y.derivatives(:,i),size(y));
+
+                        Zi = x.values\(Yi - Xi);
+
+                        z.derivatives(:,i) = Zi(:);
+
                     end
-                    z.derivatives = x.values \ (y.derivatives - kron(z.values', speye(size(x, 1))) * x.derivatives);
+
+                    %z.derivatives = x.values \ (y.derivatives - kron(z.values', speye(size(x, 1))) * x.derivatives);
                     z = AutoDiff(z);
 
                 else

MatlabGP/GP.m

@@ -196,7 +196,7 @@
 
             obj.kernel.scale = 1;
             [obj.K] = obj.kernel.build(xx,xx);
-            obj.K = obj.K + diag(0*xx(:,1)+obj.kernel.signn) + (1e-6)*eye(size(xx,1));
+            obj.K = obj.K + diag(0*xx(:,1)+obj.kernel.signn) + (1e-14)*eye(size(xx,1));
 
             res = obj.Y - obj.mean.eval(obj.X);
 
@@ -225,7 +225,7 @@
 
             detk = det(obj.K + diag(0*obj.K(:,1) + obj.kernel.signn));
 
-            if isnan(detk)
+            if isnan(detk)||isinf(detk)
                 detk = eps;
             end
 
@@ -250,7 +250,7 @@
 
             [obj] = obj.condition(obj.X,obj.Y);
 
-            detk = det(obj.K + diag(0*obj.K(:,1) + obj.kernel.signn));
+            detk = det(obj.K/obj.kernel.scale + diag(0*obj.K(:,1) + obj.kernel.signn));
 
             loss_nll = -0.5*log(sqrt(obj.kernel.scale)) - 0.5*log(abs(detk)+eps) + 0.01*sum(log(eps+gampdf(abs(theta(ntm+1:end)),1.1,0.5)));
 
@@ -315,15 +315,15 @@
 
             func = @(x) obj.LL(x,regress,ntm);
 
-            [xxt,LL] = optim.AdaptiveGridSampling(func,tlb,tub,10,20,4);
+            [xxt,LL1] = optim.AdaptiveGridSampling(func,tlb,tub,10,20,4);
 
             % xxt = tlb + (tub - tlb).*lhsdesign(200*length(tlb),length(tlb));
             % 
             % for ii = 1:size(xxt,1)
             %     LL(ii) = func(xxt(ii,:));
             % end
 
-            LL = exp(1 + LL - max(LL));
+            LL = exp(1 + LL1 - max(LL1));
 
             theta = sum(xxt.*LL(:))/sum(LL);
 

MatlabGP/MFGP.m

@@ -69,6 +69,14 @@
             end
         end
 
+        function obj = train2(obj)
+            nF = numel(obj.GPs);
+
+            for i = nF:-1:2
+                obj.Zd{i} = obj.Zd{i}.train2();
+            end
+        end
+
         function obj = resolve(obj,x,y,f)
 
             nF = numel(obj.GPs);

MatlabGP/examples/TestAppliedAviationProblem_clean.m

@@ -75,11 +75,14 @@
 %%
 for jj = 1:100
 
-    %[xn,Rn] = BO.argmax(@BO.MFSFDelta,MF);
-    [xn,Rn] = BO.argmax(@BO.maxVAR,MF);
+    [xn,Rn] = BO.argmax(@BO.MFSFDelta,MF);
+    %[xn,Rn] = BO.argmax(@BO.maxVAR,MF);
 
-    siggn(1) = abs(MF.expectedReward(xn,1))/(C(1));
-    siggn(2) = abs(MF.expectedReward(xn,2))/(C(2));
+    %siggn(1) = abs(MF.expectedReward(xn,1))/(C(1));
+    %siggn(2) = abs(MF.expectedReward(xn,2))/(C(2));
+
+    siggn(1) = abs(Z{1}.eval_var(xn))/(C(1));
+    siggn(2) = abs(Z{2}.eval_var(xn))/(C(2));
 
     [~,nu] = dec.action();
 

MatlabGP/examples/TestBackslashGradient.m

@@ -0,0 +1,30 @@
+
+clear
+clc
+
+xmesh = lhsdesign(10,2);
+
+xs = lhsdesign(1000,2);
+
+b = kernels.EQ(1,[1 1]);
+
+tt = 0.001:0.01:2;
+
+for i = 1:200
+
+theta = AutoDiff([tt(i) 1]);
+
+b = b.setHPs(theta);
+
+KXX = b.build(xmesh,xmesh);
+
+KXs = b.build(xmesh,xs);
+
+v = KXX\KXs;
+
+LL = sum(sum(v*v',1));
+
+lv(i) = LL.values;
+ld(:,i) = LL.derivatives;
+
+end

MatlabGP/examples/TestForresterProblem_clean.m

@@ -12,13 +12,13 @@
 xx = lb + (ub - lb).*lhsdesign(500,1);
 yy = testFuncs.Forrester(xx,1);
 
-x1 = [0;lb + (ub - lb).*lhsdesign(1,1);1];
+x1 = [0;1];
 y1 = testFuncs.Forrester(x1,1);
 
-x2 = [0;lb + (ub - lb).*lhsdesign(5,D);1];%20
+x2 = [0;lb + (ub - lb).*lhsdesign(1,D);1];%20
 y2 = testFuncs.Forrester(x2,2);
 
-x3 = [0;lb + (ub - lb).*lhsdesign(6,D);1];%100
+x3 = [0;lb + (ub - lb).*lhsdesign(1,D);1];%100
 y3 = testFuncs.Forrester(x3,3);
 
 x{1} = x1;
@@ -49,8 +49,8 @@
 
 %%
 tic
-MF = NLMFGP(Z,ma,a);%
-MF = MF.condition();
+MF = MFGP(Z,ma,a);%
+MF = MF.condition(Z);
 MF = MF.train();
 toc
 
@@ -80,19 +80,19 @@
 
 %%
 
-C = [50 30 1];%20
+C = [500 30 1];%20
 
 for jj = 1:200
 
-    [xn,Rn] = BO.argmax(@BO.MFSFDelta,MF);
+    [xn,Rn] = BO.argmax(@BO.maxVAR,MF);
 
-    siggn(1) = abs(MF.expectedReward(xn,1))/(C(1));
-    siggn(2) = abs(MF.expectedReward(xn,2))/(C(2));
-    siggn(3) = abs(MF.expectedReward(xn,3))/(C(3));
+    % siggn(1) = abs(MF.expectedReward(xn,1))/(C(1));
+    % siggn(2) = abs(MF.expectedReward(xn,2))/(C(2));
+    % siggn(3) = abs(MF.expectedReward(xn,3))/(C(3));
 
-    % siggn(1) = abs(Z{1}.eval_var(xn))/(C(1));
-    % siggn(2) = abs(Z{2}.eval_var(xn))/(C(2));
-    % siggn(3) = abs(Z{3}.eval_var(xn))/(C(3));
+    siggn(1) = abs(Z{1}.eval_var(xn))/(C(1));
+    siggn(2) = abs(Z{2}.eval_var(xn))/(C(2));
+    siggn(3) = abs(Z{3}.eval_var(xn))/(C(3));
 
     [~,nu] = dec.action();
 
@@ -127,8 +127,8 @@
 
     yh1 = MF.eval(xn);
 
-    MF.GPs = Z;
-    MF = MF.condition();
+    %MF.GPs = Z;
+    MF = MF.condition(Z);
 
     yh2 = MF.eval(xn);
 
@@ -150,27 +150,32 @@
 
     figure(3)
     clf(3)
+    utils.plotLineOut(MF,0,1)
     hold on
-    plot(cost,R2z)
-    plot(cost,R2MF)
-
-    figure(4)
-    clf(4)
-    hold on
-    plot(cost,RMAEz)
-    plot(cost,RMAEMF)
-    plot(cost,maxeMF)
-
-    figure(5)
-    clf(5)
-    hold on
-    plot(cost,Rie)
-    plot(cost,Ri)
-
-    figure(6)
-    clf(6)
-    hold on
-    dec.plotDists
+    plot(xn,yh2,'x')
+    % figure(3)
+    % clf(3)
+    % hold on
+    % plot(cost,R2z)
+    % plot(cost,R2MF)
+    % 
+    % figure(4)
+    % clf(4)
+    % hold on
+    % plot(cost,RMAEz)
+    % plot(cost,RMAEMF)
+    % plot(cost,maxeMF)
+    % 
+    % figure(5)
+    % clf(5)
+    % hold on
+    % plot(cost,Rie)
+    % plot(cost,Ri)
+    % 
+    % figure(6)
+    % clf(6)
+    % hold on
+    % dec.plotDists
 
     drawnow
 

MatlabGP/examples/TestMultiOutput.m

@@ -49,9 +49,9 @@
 ytest2 = Z1.eval(xtest2);
 ytest3 = Z1.eval(xtest3);
 
-sigtest1 = Z1.eval_var(xtest1);
-sigtest2 = Z1.eval_var(xtest2);
-sigtest3 = Z1.eval_var(xtest3);
+sigtest1 = Z1.eval_var(xtest1)';
+sigtest2 = Z1.eval_var(xtest2)';
+sigtest3 = Z1.eval_var(xtest3)';
 
 figure
 subplot(1,3,1)