add matlab code for volume estimation of the set of correlation matrices

README.md

@@ -15,7 +15,7 @@ Authors:
 ------------
 
 ####  Install Rcpp package  
- 
+
 * Install package-dependencies: `Rcpp`, `RcppEigen`, `BH`.  
 
 1. Then use devtools package to install `volesti` Rcpp package. In folder `/root/R-prog` Run:
@@ -116,7 +116,6 @@ The function prints the volume.
 `
 ./vol -file <filename> -sample
 `
-
 - The default settings are: `100` uniformly distributed points from the uniform distribution using billiard walk with walk length `1`.  
 - You can use the following flags: i) `-walk_length <walk_length>` to set the walk length of the random walk, ii) `-N <number_of_points>` to set the number of points to sample, iii) `-boltz` to sample from the Boltzmann distribution, iv) `-rdhr` to sample with random directions Hit and Run, `-cdhr` to sample with coordinate directions Hit and Run, `-hmc` to sample with Hamiltonian Monte Carlo with reflections, `-temperature <variance_of_boltzmann_distribution>` to set the variance of the Boltzmann distribution.  
 
@@ -130,10 +129,16 @@ The function prints the sampled points.
 ./vol -file <filename> -sdp
 `
 - The default settings are: `20` iterations are performed, with HMC sampling with walk length equal to `1`.  
+
 - You can use the following flags: i) `-N <number_of_iterations>` to set the number of iterations, ii) `-walk_length  <walk_length>` to set the walk length of the random walk iii) `-rdhr` to sample with random directions Hit and Run, `-hmc` to sample with Hamiltonian Monte Carlo with reflections.  
 
 - Example:  
-`./generate -n 10 -m 16`  
-`./vol -file sdp_prob_10_16.txt -sdp -N 30 -hmc -walk_length 3`  
-The function prints the values of the objective function of each iteration.
+  `./generate -n 10 -m 16`  
+  `./vol -file sdp_prob_10_16.txt -sdp -N 30 -hmc -walk_length 3`  
+  The function prints the values of the objective function of each iteration.
+
+
+
+## Volume estimation of correlation matrices (matlab code)
 
+In folder `matlab_code` we include the matlab code to estimate the volume of correlation matrices. To use the main function `volume` use the script `volume_example`. 

matlab_code/utils/get_billiard_L.m

@@ -0,0 +1,59 @@
+function L = get_billiard_L(m, N)
+
+    n = m^2 / 2 - m/2;
+
+    [upper, lower] = initialize_sampler(m);
+    x = zeros(n, 1);
+    A = eye(m);
+    B = zeros(m);
+
+    bc = ones(2 * n, 1);
+    L = 0;
+
+    for j = 1:N
+            v = get_direction(n);
+
+            %A(lower) = x;
+            %q = triu(A',1);
+            %A(upper) = q(upper);
+
+            B(lower) = v;
+            q = triu(B',1);
+            B(upper) = q(upper);
+
+            % compute the intersection of the line x + l*v with the
+            % boundary of the convex set that contains the PSD matrices
+            % with ones in the diagonal
+            eigenvalues = eig(B, -A);
+            l_max = 1 / max(eigenvalues);
+            l_min = 1 / min(eigenvalues);
+
+            % compute the intersection of the line x + l*v with the
+            % boundary of the hypercube [-1, 1]^n
+            lambdas = [v; -v] ./ (bc - [x; -x]);
+            l_max_temp = 1 / max(lambdas);
+            l_min_temp = 1 / min(lambdas);
+
+            % decide which boundary the line hits first
+            if (l_max_temp < l_max)
+                l_max = l_max_temp;
+            end
+            if (l_min_temp > l_min)
+                l_min = l_min_temp;
+            end
+
+            if (L < (l_max + l_min))
+                L = l_max + l_min;
+            end
+
+            % pick a uniformly distributed point from the segment
+            %lambda = l_min + rand * (l_max - l_min);
+
+            % update the current point of the random walk
+            %x = x + lambda * v;
+
+    end
+
+    L = 4 * L;
+
+end

matlab_code/utils/get_direction.m

@@ -0,0 +1,6 @@
+function [v] = get_direction(n)
+
+    v = randn(n,1);
+    v = v / norm(v);
+
+end

matlab_code/utils/get_gradient.m

@@ -0,0 +1,7 @@
+function [s] = get_gradient(q)
+
+    s = nchoosek(q, 2);
+    s = s(:,1) .* s(:,2);
+    s = s / norm(s);
+
+end

matlab_code/utils/initialize_sampler.m

@@ -0,0 +1,7 @@
+function [upper, lower] = initialize_sampler(m)
+
+    At = zeros(m, m);
+    upper  = (1:size(At,1)).' < (1:size(At,2));
+    lower = (1:size(At,1)).' > (1:size(At,2));
+
+end

matlab_code/utils/randsphere.m

@@ -0,0 +1,18 @@
+function X = randsphere(m,n,r)
+
+% This function returns an m by n array, X, in which 
+% each of the m rows has the n Cartesian coordinates 
+% of a random point uniformly-distributed over the 
+% interior of an n-dimensional hypersphere with 
+% radius r and center at the origin.  The function 
+% 'randn' is initially used to generate m sets of n 
+% random variables with independent multivariate 
+% normal distribution, with mean 0 and variance 1.
+% Then the incomplete gamma function, 'gammainc', 
+% is used to map these points radially to fit in the 
+% hypersphere of finite radius r with a uniform % spatial distribution.
+% Roger Stafford - 12/23/05
+
+X = randn(m,n);
+s2 = sum(X.^2,2);
+X = X.*repmat(r*(gammainc(s2/2,n/2).^(1/n))./sqrt(s2),1,n);

matlab_code/volume/billiard_walk_intersection.m

@@ -0,0 +1,127 @@
+function [points] = billiard_walk_intersection(m, J, L, R, N)
+
+    n = sum(sum(J>0));
+    points = zeros(n, N);
+
+    [upper, ~] = initialize_sampler_vol(m);
+    %x = zeros(n, 1);
+    A = eye(m);
+    B = zeros(m);
+
+    bc = ones(2 * n, 1);
+    pos = 1:(2*n);
+    pos = mod(pos, n);
+    pos(pos==0) = n;
+
+    for i = 1:N
+
+        for jj = 1:3
+
+        T = -log(rand) * L;
+        v = get_direction(n);
+
+        rho = 0;
+
+        while (true)
+
+            %A(lower) = x;
+            %q = triu(A',1);
+            %A(upper) = q(upper);
+
+            B(J) = v;
+            q = triu(B',1);
+            B(upper) = q(upper);
+
+            % compute the intersection of the line x + l*v with the
+            % boundary of the convex set that contains the PSD matrices
+            % with ones in the diagonal
+            [Q, eigenvalues] = eig(B, -A);
+            [max_eig, pos_max_eig] = max(diag(eigenvalues));
+            l_max = 1 / max_eig;
+
+            x = A(J);
+
+            % compute the intersection of the line x + l*v with the
+            % boundary of the hypercube [-1, 1]^n
+            lambdas = [v; -v] ./ (bc - [x; -x]);
+            [l_max_temp, pos_max] = max(lambdas);
+            l_max_temp = 1 / l_max_temp;
+
+            vrc = v'*x;
+            v2 = v'*v;
+            rc2 = x'*x;
+
+            disc_sqrt = sqrt(vrc^2 - v2 * (rc2 - R^2));
+
+            lR_max = max((-vrc + disc_sqrt)/v2, (-vrc - disc_sqrt)/v2);
+
+            [l_max, lmax_ind] = min([l_max l_max_temp lR_max]);
+
+            % decide which boundary the line hits first
+            if (lmax_ind == 2)
+                %l_max = l_max_temp;
+                % pick a uniformly distributed point from the segment
+                lambda = 0.995 * l_max;
+
+                % update the current point of the random walk
+                if (T <= l_max)
+                    %x = x + T * v;
+                    A = A + T * B;
+                    break;
+                end
+                A = A + lambda * B;
+                %x = x + lambda * v;
+
+                %reflevt the ray
+                %v = v - (2*AA(pos_max, :)*v)*AA(pos_max, :)'
+                p = pos(pos_max);
+                v(p) = -v(p);
+            elseif (lmax_ind == 1)
+                % pick a uniformly distributed point from the segment
+                lambda = 0.995 * l_max;
+                if (T <= l_max)
+                    %x = x + T * v;
+                    A = A + T* B;
+                    break;
+                end
+                s = get_gradient(Q(:, pos_max_eig));
+                % update the current point of the random walk
+                %x = x + lambda * v;
+                A = A + lambda * B;
+                %reflect the ray
+                v = v - (2*(v'*s))*s;
+            else
+                lambda = 0.995 * l_max;
+                if (T <= l_max)
+                    %x = x + T * v;
+                    A = A + T* B;
+                    break;
+                end
+                A = A + lambda * B;
+                x = x + lambda * v;
+
+                x = x / norm(x);
+                v = v - (2 * (v' * x)) * x;
+            end
+            rho = rho + 1;
+            T = T - lambda;
+        end
+        end
+        %A(lower) = x;
+        %q = triu(A',1);
+        %A(upper) = q(upper);
+        x = A(J);
+        points(:, i) = x;
+
+        %is_in = false;
+        %if (norm(x) <= R2)
+        %    is_in = true;
+        %end
+
+        %conv = update_window(ratio_parameters, is_in);
+
+    end
+
+    %ratio = ratio_parameters.count_in / ratio_parameters.tot_count;
+
+end

matlab_code/volume/check_convergence_for_ball.m

@@ -0,0 +1,46 @@
+function [conv, ratio, too_few] = check_convergence_for_ball(nu, X, r0, lb, ub, lastball)
+
+    conv = false;
+    too_few = false;
+
+    N = size(X, 2);
+
+    mm = N / nu;
+    countsIn = 0;
+    ratios = [];
+
+    for i=1:N
+
+        x = X(:,i);
+
+        if (norm(x) <= r0)
+            countsIn = countsIn + 1;
+        end
+
+        if (mod(i, mm) == 0)
+            ratios = [ratios countsIn/mm];
+            countsIn = 0;
+        end
+    end
+
+    ratio = mean(ratios);
+    t_a = tinv(0.90, nu);
+    rs = std(ratios);
+    T = rs * (t_a / sqrt(nu));
+
+    if (ratio > lb + T)
+        if (lastball)
+            conv = true;
+            return
+        end
+
+        if (ratio < ub + T)
+            conv = true;
+            return
+        end
+        return
+    end
+    too_few = true;
+
+end
+

matlab_code/volume/check_convergence_in_spectra.m

@@ -0,0 +1,70 @@
+function [conv, ratio, too_few] = check_convergence_in_spectra(m, J, X, nu, lb, ub)
+
+    too_few = false;
+    conv = false;
+
+    %n = sum(sum(J>0));
+
+    [upper, ~] = initialize_sampler_vol(m);
+    A = eye(m);
+
+    N = size(X, 2);
+    mm = N / nu;
+    countsIn = 0;
+    ratios = [];
+
+    for i=1:N
+
+        A(J) = X(:,i);
+        q = triu(A',1);
+        A(upper) = q(upper);
+
+        if (min(eig(A)) > 0)
+            countsIn = countsIn + 1;
+        end
+
+        if (mod(i, mm) == 0)
+            ratios = [ratios countsIn/mm];
+            countsIn = 0;
+            if (length(ratios) > 1)
+
+                ratio = mean(ratios);
+                t_a = tinv(0.995, nu);
+                rs = std(ratios);
+                T = rs * (t_a / sqrt(length(ratios)));
+
+
+                %[h, ~] = ttest(ratios, lb, 'Alpha', 0.005, 'Tail', 'right');
+                if (ratio + T < lb)
+                    too_few = true;
+                    conv = false;
+                    return
+                end
+
+                %[h, ~] = ttest(ratios, ub, 'Alpha', 0.005, 'Tail', 'left');
+                if (ratio - T > ub)
+                    conv = false;
+                    return
+                end
+            end
+        end
+    end
+
+    ratio = mean(ratios);
+    t_a = tinv(0.95,nu);
+    rs = std(ratios);
+    T = rs * (t_a / sqrt(nu));
+
+    if (ratio > lb + T)
+        if (ratio < ub - T)
+            conv = true;
+            return
+        end
+        conv = false;
+        return
+    end
+    too_few = true;
+    conv = false;
+
+end
+