Moved GP class to linear solve over pseudoinv

MatlabGP/GP.m

@@ -54,10 +54,10 @@
 
             if nargout>1
 
-                %v = dot(ksf',(obj.K\ksf'))';
-                %sig = (abs(obj.kernel.scale - v));
+                v = dot(ksf',(obj.K\ksf'))';
+                sig = (abs(obj.kernel.scale - v));
 
-                sig = abs(obj.kernel.scale  + obj.kernel.signn - dot(ksf',obj.Kinv*ksf')');
+                %sig = abs(obj.kernel.scale  + obj.kernel.signn - dot(ksf',obj.Kinv*ksf')');
             end
         end
 
@@ -79,10 +79,13 @@
 
             ksf = obj.kernel.build(xs,xx);
 
-            %v = dot(ksf',(obj.K\ksf'));
-            %sig = abs(obj.kernel.scale - v);
+            % v = obj.L\ksf';
+            % sig = abs(obj.kernel.scale - sum(v.^2,1))';
 
-            sig = abs(obj.kernel.scale - dot(ksf',obj.Kinv*ksf')');
+            v = dot(ksf',(obj.K\ksf'));
+            sig = abs(obj.kernel.scale - v);
+
+            %sig = abs(obj.kernel.scale - dot(ksf',obj.Kinv*ksf')');
 
         end
 
@@ -131,7 +134,7 @@
 
         end
 
-        function [obj] = condition(obj,X,Y,lb,ub)
+        function [obj] = conditionOld(obj,X,Y,lb,ub)
 
             obj.X = X;
             obj.Y = Y;
@@ -147,7 +150,6 @@
             xx = (X - obj.lb_x)./(obj.ub_x - obj.lb_x);
 
             obj.kernel.scale = 1;
-
             [obj.K] = obj.kernel.build(xx,xx);
 
             obj.K = obj.K + diag(0*xx(:,1)+obj.kernel.signn);
@@ -177,6 +179,38 @@
 
         end
 
+        function [obj] = condition(obj,X,Y,lb,ub)
+
+            obj.X = X;
+            obj.Y = Y;
+
+            if nargin<4
+                obj.lb_x = min(X);
+                obj.ub_x = max(X);
+            else
+                obj.lb_x = lb;
+                obj.ub_x = ub;
+            end
+
+            xx = (X - obj.lb_x)./(obj.ub_x - obj.lb_x);
+
+            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));
+
+            res = obj.Y - obj.mean.eval(obj.X);
+
+            obj.alpha = obj.K\res;
+
+            sigp = sqrt(abs(res'*obj.alpha./(size(obj.Y,1))));
+
+            obj.kernel.scale = sigp^2;
+
+            obj.K = obj.kernel.scale*obj.K;
+            obj.alpha = obj.alpha/sigp^2;
+
+        end
+
         function [nll] = LL(obj,theta,regress,ntm)
 
             if regress

MatlabGP/MFGP.m

@@ -49,6 +49,8 @@
 
         function obj = condition(obj,GPs)
 
+            obj.GPs = GPs;
+
             nF = numel(obj.GPs);
 
             for i = nF:-1:2

MatlabGP/NLMFGP.m

@@ -1,8 +1,7 @@
 %{
-    Multi-Fidelity Gaussian Process
+    Non-Linear Multi-Fidelity Gaussian Process
     
-    An auto-regressive (AR(1)) model based on Kennedy and O'Hagan using the
-    recursive formulation of Le Gratiet.
+    An auto-regressive (AR(N)) model based on Perdikaris
 
     The model can be constructed using data from nF fidelities.
 

MatlabGP/VGP.m

@@ -64,8 +64,11 @@
             y = obj.mean.eval(x) + ksu*obj.alpha;
 
             if nargout>1
+
+                v1 = obj.Kuu\ksu';
+                v2 = obj.M\ksu';
 
-                sigs = -dot(ksu',(obj.Kuuinv)*ksu') + dot(ksu',obj.Minv*ksu');
+                sigs = dot(ksu',v2)-dot(ksu',v1);
 
                 sig = obj.kernel.scale + obj.kernel.signn + sigs';
             end
@@ -102,7 +105,10 @@
 
             ksu = obj.kernel.build(xs,xu);
 
-            sigs = -dot(ksu',(obj.Kuuinv)*ksu') + dot(ksu',obj.Minv*ksu');
+            v1 = obj.Kuu\ksu';
+            v2 = obj.M\ksu';
+
+            sigs = dot(ksu',v2)-dot(ksu',v1);
 
             sig = obj.kernel.scale + obj.kernel.signn + sigs';
 
@@ -205,16 +211,15 @@
 
             obj.kernel.scale = std(Y)/2;
 
-            obj.Kuu = obj.kernel.build(xu,xu);
-            obj.Kuuinv = pinv(obj.Kuu,1*10^(-7));
+            obj.Kuu = obj.kernel.build(xu,xu) + (1e-6)*eye(size(xu,1));
+            %obj.Kuuinv = pinv(obj.Kuu,1*10^(-7));
 
             obj.Kuf = obj.kernel.build(xu,xf);
 
             obj.B = obj.Kuf*obj.Kuf'/obj.kernel.signn;
             obj.M = obj.Kuu + obj.B;
-            obj.Minv = pinv(obj.M,1*10^(-7));
 
-            obj.alpha = obj.Minv*obj.Kuf*(Y - obj.mean.eval(X))/obj.kernel.signn;
+            obj.alpha = obj.M\(obj.Kuf*(Y - obj.mean.eval(X))/obj.kernel.signn);
 
         end
 
@@ -397,7 +402,6 @@
 
                 obj.B = obj.B + k2s*k2s'/obj.kernel.signn;
                 obj.M = obj.Kuu + obj.B;
-                obj.Minv = pinv(obj.M);
 
                 obj.alpha = obj.alpha + k2s*y/obj.kernel.signn;
             end