Estimate RIP constant of matrix

rip_vec.m

n = 100;
m = 50;
A = 1/sqrt(m) * randn(m, n);

K = 10;
trials = 1e4;
rmin = ones(K, 1);
rmax = ones(K, 1);
for k = 1:K
  for i = 1:trials
    s = randperm(n);
    s = s(1:k);
    x = zeros(n, 1);
    x(s) = randn(k, 1);
    r = norm(A*x) / norm(x);
    rmin(k) = min(rmin(k), r);
    rmax(k) = max(rmax(k), r);
  end;
end;

alpha = rmax ./ rmin;
eps = (alpha - 1) ./ (alpha + 1);

bar(1:K, eps);