Initial commit. Working implementation of FHMC.

Author: qiyangku

main_unimodal.m

@@ -0,0 +1,175 @@
+function main_unimodal(option,varargin)
+
+switch option
+    case 'plot_only'
+        plot_unimodal(varargin{1},[0 2]);
+        return
+    case 'analysis_only'
+        post_process_unimodal(varargin{1});
+        return
+    otherwise
+        gen_raw_data_unimodal;
+end
+end
+
+
+function gen_raw_data_unimodal(~)
+
+p.ntrials = 8*3; %for parallel processing,use multiples of 8
+p.target = 'unimodal';
+p.sigma = 1;
+p.dt = 1e-3; %1e-2 is no good for calculating Riesz derivative
+T = 1e3;
+
+%AIM: 1. Sample trajectory. 2. histogram 2. ACF. 3. MSE
+%For FHMC,HMC,
+
+a = [1.2 2];
+solver = {'Hamiltonian2','Langevin'}; %'Hamiltonian','Underdamped'};
+%for accurate results don't need the other two... too slow
+
+X = zeros(T/p.dt,p.ntrials,length(solver),length(a));
+
+disp('Generating raw data...')
+for i = 1:length(a) %group fractional/nonfractional together as this is less familiar to people
+    aa = a(i);
+    tic
+    for j = 1:length(solver)
+        p.methods.solver = solver{j};
+        temp = zeros(T/p.dt,p.ntrials);
+        parfor k = 1:p.ntrials
+            [temp(:,k),t] = fHMC(T,aa,p);
+        end
+        clc
+        fprintf('Solver: %u/%u\n',j,length(solver))
+        fprintf('alpha: %u/%u\n',i,length(a))
+        toc
+        X(:,:,j,i) = temp;
+    end
+    
+    
+    
+end
+
+n = floor(T/p.dt); %number of samples
+t = (0:n-1)*p.dt;
+
+save main_unimodal_raw_data.mat T a solver X t p
+end
+
+function post_process_unimodal(Z)
+n=50;
+binedge = linspace(-5,5,n+1);
+bin = binedge(1:end-1) + 0.5*(binedge(2)-binedge(1));
+disp('Evaluating histogram, ACF, and MSE')
+H = zeros(n,length(Z.solver),length(Z.a));
+
+m = 10/Z.p.dt;%1e4;
+tau = (0:m)*Z.p.dt;
+ACF = zeros(length(tau),length(Z.solver),length(Z.a));
+MSE = zeros(length(tau),length(Z.solver),length(Z.a));
+
+t = (0:5e4)*Z.p.dt;
+X = zeros(length(t),length(Z.solver),length(Z.a)); %sample X for visualization
+
+tic
+for i = 1:length(Z.a)
+    for j = 1:length(Z.solver)
+        x = Z.X(:,:,j,i); x = x(:);
+        H(:,j,i) = histcounts(x,binedge,'normalization','pdf');
+        temp_acf = 0;
+        temp_mse = 0;
+        for k = 1:Z.p.ntrials
+            temp_acf = temp_acf + autocorr(Z.X(:,k,j,i),m);
+            temp_mse = temp_mse + myMSE(Z.X(:,k,j,i),m);
+        end
+        ACF(:,j,i) = temp_acf/Z.p.ntrials;
+        MSE(:,j,i) = temp_mse/Z.p.ntrials;
+        X(:,j,i) = Z.X(1:length(t),1,j,i);
+        clc
+        toc
+        fprintf('%u / %u\n',i,length(Z.a))
+        fprintf('%u / %u\n',j,length(Z.solver))
+    end
+end
+save main_unimodal_post_process.mat H ACF MSE tau bin binedge t X
+end
+
+
+function plot_unimodal(Q,index)
+
+close all
+if any(index == 0) %sample trajectory
+    f0 = myfigure;
+    ttl = {'Fractional Hamiltonian Monte Carlo','Fractional Langevin Monte Carlo','Hamiltonian Monte Carlo','Langevin Monte Carlo'};
+    for k = 1:4
+        subplot(4,1,k)
+        [j,i]=ind2sub([2 2],k);
+        plot(Q.t(1:10:end),Q.X(1:10:end,j,i),'.','color',mycolor(1,k),'markersize',1);
+        title(ttl(k))
+        if k<4
+            set(gca,'xtick',[],'XColor','none');
+        else
+            xlabel('t')
+        end
+        ylabel('x_t')
+        %subplotmod;
+        box off
+        set(gca,'TickDir','out')
+        %hold on
+        %plot(xx,yy,'k--','linewidth',1.5)        
+        ylim([-4 4])
+        xlim([Q.t(1) Q.t(end)])
+    end
+    pos =get(f0,'Position');
+    set(gcf,'position',[pos(1) pos(2) pos(3) 600])
+    export_fig(f0,'figures/fig_main_unimodal_postprocess_sample_path.pdf','-pdf','-nocrop','-transparent','-painters');
+end
+
+if any(index == 1)    %ACF
+    f1 = myfigure;    
+    % plot(Q.tau,Q.ACF(:,1,1),'-','color',mycolor(1,2),'linewidth',1.5)
+    % hold on
+    % plot(Q.tau,Q.ACF(:,2,1),'-','color',mycolor(1,1),'linewidth',1.5)
+    % plot(Q.tau,Q.ACF(:,1,2),'--','color',mycolor(1,2) ,'linewidth',1.5)
+    % plot(Q.tau,Q.ACF(:,2,2),'--','color',mycolor(1,1),'linewidth',1.5)
+    
+    %color scheme default
+    plot(Q.tau,Q.ACF(:,1:2,1),Q.tau,Q.ACF(:,1:2,2),'--','linewidth',1.5)
+    hold on
+    plot(Q.tau([1 end]),[0 0],'color',[0.4 0.4 0.4],'HandleVisibility','off')
+    
+    xlim([0 8])
+    xlabel('Time Lag')
+    ylabel('Autocorrelation function')
+    legend('FHMC','FLMC','HMC','LMC','box','off')
+    subplotmod;
+    offsetAxes;
+    
+    export_fig(f1,'figures/fig_main_unimodal_postprocess_acf.pdf','-pdf','-nocrop','-transparent','-painters');
+end
+if any(index == 2)
+    f2 = myfigure;
+    gs = makedist('normal');
+    xx = linspace(-3,3,51);
+    yy = pdf(gs,xx);
+    
+    for k = 1:4
+        subplot(4,1,k)
+        [j,i]=ind2sub([2 2],k);
+        bar(Q.bin,Q.H(:,j,i),1,'FaceColor',mycolor(1,k),'EdgeColor',mycolor(1,k),'FaceAlpha',0.5);
+        %bar(bins,hc,1,'FaceColor',mycolor(2),'EdgeColor',mycolor(0));hold on
+        set(gca,'xtick',[])
+        set(gca,'ytick',[])
+        hold on
+        plot(xx,yy,'k--','linewidth',1.5)        
+        ylim([0 0.5])
+        xlim([-4 4])
+        set(gcf,'position',[100 100 200 600])
+    end
+    %export_fig(f2,'figures/fig_main_unimodal_postprocess_histogram.pdf','-pdf','-nocrop','-transparent','-painters');
+    print(gcf, '-dpdf', 'figures\fig_main_unimodal_postprocess_histogram.pdf'); 
+end
+
+
+end

mains/subroutines/TrapTime.m

@@ -0,0 +1,23 @@
+function T = TrapTime(t)
+%INPUT: t = binary variable
+%OUTPUT: Time spent in one mode
+%t = x < 0;
+dt = diff(t);
+t_in = find(dt == 1);
+t_out = find(dt == -1);
+
+T = zeros(size(t_in));
+for j =1 :length(T)
+    if t_in(j)==1 %ignore influence from initial condition
+        T(j)=NaN;
+        continue
+    end
+    
+    jj = find(t_out>t_in(j),1);
+    if isempty(jj)  %ignore terminal condition
+        T(j) = NaN;
+    else
+        T(j) = t_out(jj)-t_in(j);
+    end
+end
+T(isnan(T))=[];

mains/subroutines/fHMC.m

@@ -0,0 +1,123 @@
+function [X,t] = fHMC(T,a,p)
+%Levy Monte Carlo
+%INPUT: T = time span (s), a = Levy characteristic exponent
+% p = other parameters
+%OUTPUT: X = samples, t = time
+%
+p.depth = 1/p.sigma/p.sigma;
+p.rmax = 5; %cutoff
+p.xmax = 5; %cutoff
+p.x0 = 0;
+p.r0 = 1;
+
+%p.methods.drift = 'approx'; %exact (numerically unstable)
+
+
+%
+dt = p.dt;%1e-3; %integration time step (s)
+dta = dt.^(1/a); %fractional integration step
+
+n = floor(T/dt); %number of samples
+t = (0:n-1)*dt;
+x = p.x0; %initial condition
+r = p.r0;
+X = zeros(n,1); %position
+%R = zeros(n,1); %momentum
+
+switch p.target
+    case 'unimodal'
+        U = @(x) p.depth*x.*x/2;
+        dU = @(x) p.depth*x;
+    case 'bimodal'
+        s1 = p.location;%2.5;
+        s2 = -s1;
+        w = 0.5;
+        sigma2 = p.sigma.*p.sigma;%0.32;%*(0.32/0.54);
+        U = @(x) -log( 1e-8 +w*exp(-(x-s1).*(x-s1)/2/sigma2) + (1-w)*exp(-(x-s2).*(x-s2)/2/sigma2))+0.5;
+        dx = 1e-3;
+        dU = @(x) (U(x+dx)-U(x-dx))/dx/2;
+end
+
+if a<2 %approximation of drift term
+    ba = @(x) frac_drift(U,dU,x,a); %approximate full Riesz D with partial Dx                
+else
+    ba = @(x) -dU(x); %approximation
+end
+
+c = gamma(a-1)*(gamma(a/2).^(-2)); %<--should be this!
+% if a < 2  %use this for a=1.2 as it seems to give better result???
+%     c = gamma(a+1)*(gamma(a/2+1).^(-2)); %a = a-2;
+% end
+%tic
+switch p.methods.solver
+    case 'Hamiltonian'
+        for i = 1:n
+            while true
+                %diffusion term
+                g = stblrnd(a,0,1,0); %Levy r.v. wt expo a and scale 1                
+                xnew = x + c*r*dt;
+                rnew = r -c*dU(x)*dt - c*r*dt + g*dta;
+                %accept criterion
+                if abs(rnew)<p.rmax %5
+                    x = xnew;
+                    r = rnew;
+                    break %accept
+                end
+            end
+            X(i) = x;
+        end
+    case 'Langevin'
+        for i = 1:n
+            while true
+                %diffusion term
+                g = stblrnd(a,0,1,0); %Levy r.v. wt expo a and scale 1
+                %NB. when a=2, g~sqrt(2)*gaussian r.v. as it should be
+                xnew = x + ba(x)*dt + g*dta;
+                %accept criterion
+                if abs(xnew)<p.xmax %5
+                    x = xnew;
+                    break %accept
+                end
+            end
+            X(i) = x;
+        end
+    case 'Hamiltonian2'
+        for i = 1:n
+            while true
+                %diffusion term
+                g = stblrnd(a,0,1,0); %Levy r.v. wt expo a and scale 1                
+                xnew = x + ba(x)*dt + r*dt + g*dta;
+                rnew = r + ba(x)*dt;
+                %accept criterion
+                if abs(xnew)<p.xmax %5
+                    x = xnew;
+                    r = rnew;
+                    break %accept
+                end
+            end
+            X(i) = x;
+        end
+    case 'Underdamped'
+        sas = makedist('Stable','alpha',1/a,'beta',0,'gam',1,'delta',0);
+        dv = 1e-3;
+        G = @(v) (pdf(sas,v+dv)-pdf(sas,v-dv))/dv/2; %rough approximation
+        for i = 1:n
+            while true
+                %diffusion term
+                g = stblrnd(a,0,1,0); %Levy r.v. wt expo a and scale 1
+                
+                rnew = r - dU(x)*dt - r*dt + g*dta;
+                xnew = x + G(rnew)*dt;
+                
+                %accept criterion
+                if abs(rnew)<p.rmax %5
+                    x = xnew;
+                    r = rnew;
+                    break %accept
+                end
+            end
+            X(i) = x;
+        end
+        
+end
+%toc

mains/subroutines/frac_drift.m

@@ -0,0 +1,23 @@
+function b = frac_drift(U,dU,x,a)
+% Method from the paper <Fractional Langevin Monte Carlo...>
+% Stable implementation of drift term to avoid 0/0 error for large x
+% INPUT: x = scalar or column vector; a = noise expo; type = well type
+% U,dU are function handles
+
+a = a - 2; %Levy exponent converted to that used for Riesz derivative
+
+h = 0.1; %discretization
+L = 8; %spatial range
+K = floor(L/h+1);
+k = -K:K; %row vector
+g = (-1).^k.*gamma(a+1)./gamma(a/2 -k+1)./gamma(a/2 +k +1);
+
+l = U(x)-U(x-k*h);
+l2 = max(l,[],2);
+f = -dU(x-k*h).*exp(l-l2);
+%f(f<0)=0;
+b = sum(g.*f,2)/(h^a).*exp(l2);
+
+%clipping
+b(b>1e3) = 500;%1e3; %limit the height of the well
+b(b<-1e3) = -500;%1e3; % this is set such that abs(b*dt) = 1

mains/subroutines/myMSE.m

@@ -0,0 +1,18 @@
+function MSE = myMSE(s,dT)
+%calculate MSE of a sampler from true mean
+
+%dT = 1e4; %ms
+dT=dT+1; %to be consistent with matlab autocorr
+%if ~exist('MSE','var')
+MSE = 0;
+k = 1;
+for i=1:1e3:length(s)-dT
+    %ss  =  angle(cumsum(  exp(1i*s(i + (1:dT)-1)) ));
+    %dss = wrapToPi(ss - true_mean);
+    dss = cumsum(s(i + (1:dT)-1))./(1:dT)'; %cum mean
+    MSE = MSE + dss.*dss;
+    k = k+1;
+end
+MSE = MSE/k;%(length(s)-dT);
+%end
+%t = 1:dT;