experiments.py

# %%
from cvxopt import spmatrix, matrix, solvers
import numpy as np
from math import sqrt, log
import numpy as np
import numpy as np
from sklearn.linear_model import Lasso


# %%
# setup
def sensing_matrix(n, N, gaussian=False):
    """
    Sensing matrix of shape (n, N). By default, entries are independent Rademacher with unit variance.
    """
    if gaussian:
        return 1 / sqrt(n) * np.random.randn(n, N)
    else:
        return 1 / sqrt(n) * (2 * (np.random.randn(n, N) > 0) - 1)


def subsets(N, k, d):
    """
    Dataset of k-sparse subsets of shape (N, d).
    """
    return np.array([
        np.bincount(np.random.choice(N, k, replace=False), minlength=N) for _ in range(d)
    ]).T


# %%
# ended up not using this one
def lp_experiment(N, k, n, gaussian=False):
    f = sensing_matrix(n, N, gaussian) # (n, N)
    x = subsets(N, k, 1).reshape((-1)) # (N)
    y = f @ x # (n)

    # Set up an LP
    #   minimize <c, x>
    #   s.t. Gx - h <= 0  (N)
    #        f x = y      (n)

    c = matrix(1., (N, 1))
    G = spmatrix(-1., range(N), range(N))
    h = matrix(0., (N, 1))

    sol = solvers.lp(c, G, h, matrix(f), matrix(y))
    decoded = np.array(sol['x']).reshape((-1)) > 0.5

    return np.sum(decoded != x).item()


# %%
def lasso_experiment(N, k, n, gaussian=False):
    f = sensing_matrix(n, N, gaussian) # (n, N)
    x = subsets(N, k, 1).reshape((-1)) # (N)
    y = f @ x # (n)

    lasso_cv = Lasso(fit_intercept=False, max_iter=10000, positive=True, alpha=1e-5)
    lasso_cv.fit(f, y)
    x_reconstructed = lasso_cv.coef_

    decoded = (x_reconstructed > 0.5).astype(int)
    return np.sum(decoded != x).item()


# %%
def linear_experiment(N, k, n, gaussian=False):
    f = sensing_matrix(n, N, gaussian) # (n, N)
    x = subsets(N, k, 1).reshape((-1)) # (N)
    y = f @ x # (n)

    correlation = f.T @ y
    decoded = correlation > 0.5
    return np.sum(decoded != x).item()