%% Reproducible file accompayining the paper
%% Linearly Embedding Sparse Vectors from \ell_2 to \ell_1 via Deterministic Dimension-Reducing Maps.
%% by Simon Foucart

close all; clear all; clc;

%
%% Some Low-Coherence Matrices (see Section 2)
%
%% The DeVore matrix
% To produce the "DeVore" matrix A_devo, we need an auxiliary matrix B containing the values (mod p) at {0,1,...,p-1} of all the polynomials of degree at most d. Its recursive construction is provided at end of the file.

p = 13;            % p is a prime number
d = 3;             % d is an integer < p
m = p^2;           % m is the number of rows
N = p^(d+1);       % N is the number of columns
B = value_matrix(p,d+1);
A_devo = zeros(m,N);
for j=0:p-1
    A_devo(j*p+1:(j+1)*p,:) = (B==j)/sqrt(p); 
end
Gram_devo = A_devo'*A_devo;
coh_devo = max(max(abs(Gram_devo - diag(diag(Gram_devo)))));
fprintf('The coherence of A_devo is %.3f, which indeed coincides with d/sqrt(m) = %.3f.', coh_devo, d/p)

%
%% The Weil matrix
% The "Weil" matrix A_weil capitalizes on the Weil bound. We also rely on the auxiliary matrix B for its swift construction.

p = 13;            % p is a prime number
d = 3;             % d is an integer < p
m = p;             % m is the number of rows
N = p^(d+1);       % N is the number of columns
B = value_matrix(p,d+1);
A_weil = exp( 1i*2*pi/p * diag(0:p-1) * B) / sqrt(p);
Gram_weil = A_weil'*A_weil;
coh_wei = max(max(abs(Gram_weil - diag(diag(Gram_weil)))));
fprintf('The coherence of A_weil is %.3f, which is indeed smaller than d/sqrt(m) = %.3f.', coh_wei, d/sqrt(p))

%
%% Explicit Isometric Embedding from \ell_2^p(\mathbb{C}) to \ell_4^{c p^2}(\mathbb{C}) (see Section 4)
% The matrix M_isom of the isometric embedding is a slight alteration of the "Golomb" matrix A_golo, whose construction is based on Golomb rulers (aka Sidon sets).

m = 6*p^2-6*p+1;
A_golo = zeros(m,p);
for k=0:p-1
    gk = 2*p*k + mod(k^2,p);
    A_golo(:,k+1) = exp(1i*2*pi*gk/m*(0:m-1));
end
M_isom = [ A_golo / (2*m)^(1/4); eye(p) / 2^(1/4) ];
x = randn(p,1);
fprintf('Verification on a random x: the 2-norm of x is %.3f and the 4-norm of M_isom x is also %.3f.',...
    norm(x,2), norm(M_isom*x,4))

%
%% Explicit Embedding of s-Sparse Vectors from \ell_2^N(\mathbb{C}) to \ell_1^{c s^4}(\mathbb{C}) (see Section 5)
% The matrix below is simply the (normalized) product of the Golomb matrix and the Weil matrix.
% Note that the parameter regime below does not meet the assumptions stated in Theorem 4.

A = A_golo*A_weil/m;
s = 4;
n_test = 100;
ratio = zeros(n_test,1);
for t=1:n_test
    x = zeros(N,1);
    supp = randsample(N,s); 
    x(supp) = randn(s,1);
    ratio(t) = norm(A*x,1)/norm(x,2);
end
fprintf('After %i tests on %i-sparse vectors x, the 1-norm of Ax divided by the 2-norm of x was always between %.3f and %.3f.', n_test, s, min(ratio), max(ratio));

%
%% Quasideterministic Embedding of s-Sparse Vectors from \ell_2^N(\mathbb{R}) to \ell_1^{c s^4}(\mathbb{R}) (see Section 3)
% We verify below the assumptions (a)-(b) from Theorem 1 for a realization of a Rademacher random matrix. The verification, despite being performed in "polynomial time", is very slow. Moreover, using 4 nested for loops would be quite impractical in MATLAB. For these reasons,  we relied on some not-so-intuitive workarounds.
% Note: To be convincing, the construction should be tested on matrices of much bigger sizes, but this would take a extremely long time. Below, the choice of m=800 and N=1400 (with 20 pieces) took a little less than 4 hours and the resulting cutoff sparsity s \approx ( m / \kappa^2)^{1/4} was only 2.

tic;
m = 800;              % m is the number of rows
N = 1400;             % N is the number of columns
A_rade = sign(randn(m,N));
% we verify assumptions (a)-(b) at once by calculating the coherence of a matrix 
% that contains the pointwise (Hadamard) products of pairs of columns of A;
% for efficient storage in memory, this matrix shall be manipulated in pieces
Idx = zeros(2,N*(N-1)/2);
for k=1:N-1
    for l=k+1:N
        Idx(:,(l-2)*(l-1)/2+k) = [k; l]; 
    end
end
n_pieces = 20;
NN = ceil(N*(N-1)/2/n_pieces);
AA = cell(n_pieces,1);
for i=1:n_pieces-1
    aux = zeros(m,NN);
    for j= 1:NN
        k = Idx(1,(i-1)*NN+j);
        l = Idx(2,(i-1)*NN+j);
        aux(:,j) = A_rade(:,k).*A_rade(:,l);
    end
    AA{i} = aux;
end
aux = zeros(m,N*(N-1)/2-(n_pieces-1)*NN);
for j=1:N*(N-1)/2-(n_pieces-1)*NN
    k = Idx(1,(n_pieces-1)*NN+j);
    l = Idx(2,(n_pieces-1)*NN+j);
    aux(:,j) = A_rade(:,k).*A_rade(:,l);
end
AA{n_pieces} = aux;
% now we calculate the coherence of AA piece by piece
kapp = zeros(n_pieces,n_pieces);
for i=1:n_pieces
    for j=i:n_pieces
        G = AA{i}'*AA{j};
        if i==j
            kapp(i,i) = max(max(abs(G - diag(diag(G))))) / sqrt(m);
        else
            kapp(i,j) = max(max(abs(G))) / sqrt(m);
        end
    end
end
t=toc;
fprintf('The value of kappa making (a)-(b) hold is %.3f.', max(max(kapp)))

%
%% Appendix: construction of the auxiliary matrix used for the DeVore and Weil matrices
% Given integers d \le p, the p \times p^d matrix M below has entries in \{0,1,\ldots,p-1\}. Its columns, indexed by polynomials f of degree <d, contain the values f(0),f(1),\ldots,f(p-1) modulo p.

function M = value_matrix(p,d)
M_1 = zeros(p,p);
for c=0:p-1
    M_1(:,c+1) = c*ones(p,1); 
end
if d==1
    M = M_1;
else
    M_aux = value_matrix(p,d-1);
    pp = p^(d-1);
    for c=0:p-1
        M(:,c*pp+1:(c+1)*pp) = mod( c+diag(0:p-1)*M_aux , p);
    end
end  
end