underdetermined-systems-sparsity.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "19eeaff9",
   "metadata": {
    "cell_marker": "r\"\"\""
   },
   "source": [
    "We are interested in solving systems of the form\n",
    "$Xb = y$, where $X$ is a matrix of size $n \\times p$ and $y$ is our $n$ observations.\n",
    "Our goal here is to recover $b$ when p > n.\n",
    "\n",
    "We will assume that $b$ is sparse meaning that $b$ has very few non-zero entries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e0f57fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import linprog\n",
    "import matplotlib.pyplot as plt\n",
    "from itertools import combinations\n",
    "from scipy.stats import norm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "809a4d92",
   "metadata": {},
   "source": [
    "# Question 1\n",
    "Generate a dataset.\n",
    "- Create a matrix $X$ of size $(n, p)$ where each entry of $X$ is sampled from $\\mathcal{N}(0, 1)$\n",
    "- Create a vector $b$ of size $p$ where only the $k$ first entries are non-zero\n",
    "- Generate the vector $y = Xb$ of size $n$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dba33732",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "def generate_dataset(k, n, p):\n",
    "    \"\"\"Generate a dataset\n",
    "    Generate a dataset where\n",
    "    X has random entries sampled from N(0, 1),\n",
    "    y such that Xb = y\n",
    "    b has its k first entry equal to 1 and other entries are 0.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    k : int\n",
    "        Number of non-zero components of b\n",
    "    n : int\n",
    "        Number of samples\n",
    "    p : int\n",
    "        Number of features\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    X : array of size (n, p)\n",
    "        Design matrix\n",
    "    y : array of size (n,)\n",
    "        Observations\n",
    "    b : array of size (p,)\n",
    "        Parameters\n",
    "    \"\"\"\n",
    "    raise NotImplementedError\n",
    "\n",
    "\n",
    "X, y, b =  generate_dataset(k=2, n=4, p=5)\n",
    "np.testing.assert_allclose(X.dot(b), y)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f8c8d50",
   "metadata": {},
   "source": [
    "# Question 2\n",
    "## Create a solver for linear regression\n",
    "Consider the problem of finding $b$ such that $Xb = y$\n",
    "- If the system has no solution return None\n",
    "- If the system has one or more solution, return a solution\n",
    "\n",
    "Help: you can remove redundant equations using the SVD of $X$.\n",
    "Taking $UDV$ an SVD of $X$, solving $Xb = y$ amount to solve $(DV)b = U^T y$\n",
    "in particular the zeroes in $D$ should yield zeroes in $U^T y$ (If this is not the case there is no solution and\n",
    "you should return None. Otherwise return a solution).\n",
    "In practice you won't encounter exact zeroes so it is enough to check that $U^Ty$ has absolute value lower than 1e-10 at coordinates where $D$ has absolute value lower than 1e-10.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0bdf704a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def solve_linear_regression(X, y):\n",
    "    \"\"\"Find a solution of linear regression\n",
    "\n",
    "    Find b such that Xb = y if it exists.\n",
    "    Otherwise return None\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    X : array of size (n, p)\n",
    "        Design matrix\n",
    "    y : array of size (n,)\n",
    "        Observations\n",
    "\n",
    "    Return\n",
    "    ------\n",
    "    solution: np array of size (p,) or None\n",
    "        A solution if there exists one or None if there isn't any\n",
    "    \"\"\"\n",
    "    raise NotImplementedError\n",
    "\n",
    "def test_solve_linear_regression():\n",
    "    # An overdetermined but redondant system\n",
    "    X = np.array([[1, 2], [3, 4], [2, 4], [6, 8]])\n",
    "    b = np.array([1, 2])\n",
    "    y = np.array([5, 11, 10, 22])\n",
    "    sol = solve_linear_regression(X, y)\n",
    "    np.testing.assert_allclose(b, sol)\n",
    "\n",
    "    # An underdetermined system\n",
    "    X = np.array([[2, 1, 2]])\n",
    "    b = np.array([0, 1, 2])\n",
    "    y = np.array([5])\n",
    "    sol = solve_linear_regression(X, y)\n",
    "    np.testing.assert_allclose(y, X.dot(sol))\n",
    "\n",
    "    # A system without a solution\n",
    "    X = np.array([[1, 2], [3, 4], [2, 4], [6, 8]])\n",
    "    b = np.array([1, 2])\n",
    "    y = np.array([5, 11, 11, 22])\n",
    "    sol = solve_linear_regression(X, y)\n",
    "    assert sol is None\n",
    "    \n",
    "test_solve_linear_regression()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aee0455a",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "# Question 3\n",
    "\n",
    "Consider a matrix $S$ with boolean entries (True or False).\n",
    "\n",
    "Build a solver for finding $b$ such that $Xb = y$ where $b$ has only non-zero coordinates at indexes $i$ such that $S[i] = True$ (if no such solution exists, return None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d00b4f55",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def search_in_support(S, X, y):\n",
    "    \"\"\"Find a solution in the support S\n",
    "\n",
    "    Solve Xb = y with constraints b_j = 0 for all j such \n",
    "    that S[j] = False\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    S : boolean array of size (p,)\n",
    "        Support of b\n",
    "    X : array of size (n, p)\n",
    "        Design matrix\n",
    "    y : array of size (n,)\n",
    "        Observations\n",
    "\n",
    "    Return\n",
    "    ------\n",
    "    solution: np array of size (p,) or None\n",
    "        A solution if there exists one or None if there isn't any\n",
    "    \"\"\"\n",
    "    raise NotImplementedError\n",
    "\n",
    "\n",
    "def test_search_in_support():\n",
    "    # An overdetermined but redondant system\n",
    "    X = np.array([[2, 1, 2], [4, 3, 4], [3, 2, 4], [1, 6, 8]])\n",
    "    b = np.array([0, 1, 2])\n",
    "    y = np.array([5, 11, 10, 22])\n",
    "    sol = search_in_support(np.array([False, True, True]), X, y)\n",
    "    np.testing.assert_allclose(b, sol)\n",
    "\n",
    "    # An underdetermined system\n",
    "    X = np.array([[2, 1, 2]])\n",
    "    b = np.array([0, 1, 2])\n",
    "    y = np.array([5])\n",
    "    sol = search_in_support(np.array([True, True, True]), X, y)\n",
    "    np.testing.assert_allclose(y, X.dot(sol))\n",
    "\n",
    "    # A system without a solution\n",
    "    X = np.array([[2, 1, 2], [4, 3, 4], [3, 2, 4], [1, 6, 8]])\n",
    "    b = np.array([0, 1, 2])\n",
    "    y = np.array([5, 11, 11, 22])\n",
    "    sol = search_in_support(np.array([False, True, True]), X, y)\n",
    "    assert sol is None\n",
    "\n",
    "test_search_in_support()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "245cdd6c",
   "metadata": {},
   "source": [
    "# Question 4\n",
    "\n",
    "- Build a solver that takes all possible support of a given size and tries them sequentially until it finds a solution (or nothing works).\n",
    "\n",
    "Example: the possible supports of size $2$ in dimension $3$ are\n",
    "$(True, True, False), (True, False, True), (True, True, False)$\n",
    "\n",
    "For all supports S of size `size_support`, build a solver for finding $b$ such that $Xb = y$ where $b$ has only non-zero coordinates at indexes $i$ such that $S[i] = True$, return the first solution found. If None of the support yield a solution, return None.\n",
    "\n",
    "- Use your solver to recover b\n",
    "\n",
    "    * Generate a dataset of with n=10, p=20,k=4\n",
    "    * Do an exhaustive search among supports of size 5\n",
    "    * Check that the solution found is the correct one"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d054abb",
   "metadata": {},
   "outputs": [],
   "source": [
    "def exhaustive_search(size_support, X, y):\n",
    "    \"\"\"For all supports S of size `size_support` try to solve Xb = y\n",
    "    If there exists a support that has at least one solution, return one solution.\n",
    "    If None of the support yield a solution, return None\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    size_support: int\n",
    "        The size of the support\n",
    "    X : array of size (n, p)\n",
    "        Design matrix\n",
    "    y : array of size (n,)\n",
    "        Observations\n",
    "\n",
    "    Return\n",
    "    ------\n",
    "    solution: np array of size (p,) or None\n",
    "        A solution if there exist one\n",
    "    \"\"\"\n",
    "    raise NotImplementedError\n",
    "\n",
    "\n",
    "def test_exhaustive_search():\n",
    "    # An overdetermined but redondant system\n",
    "    X = np.array([[2, 1, 2], [4, 3, 4], [3, 2, 4], [1, 6, 8]])\n",
    "    b = np.array([0, 1, 2])\n",
    "    y = np.array([5, 11, 10, 22])\n",
    "\n",
    "    sol = exhaustive_search(2, X, y)\n",
    "    np.testing.assert_allclose(b, sol)\n",
    "    \n",
    "    sol = exhaustive_search(1, X, y)\n",
    "    assert sol is None\n",
    "\n",
    "\n",
    "test_exhaustive_search()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07c586f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_sparse_b():\n",
    "    \"\"\"\n",
    "    - Generate a dataset of size n=10, p=20,k=4\n",
    "    - Do an exhaustive search among supports of size 5\n",
    "    - Check that the solution found is the correct one\n",
    "    \"\"\"\n",
    "    raise NotImplementedError\n",
    "\n",
    "\n",
    "find_sparse_b()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5729f12d",
   "metadata": {},
   "source": [
    "# Question 5\n",
    "\n",
    "- Use linprog to solve Basis pursuit:\n",
    "\n",
    "$\\min_b |b|$ with constraints $Xb = y$ where $| |$ is the l1 norm.\n",
    "\n",
    "Help: BP can be recast as a linear program with n equality constraints\n",
    "and 2p variables. Use variables for the positive and negative part of each entry of b.\n",
    "\n",
    "linprog solves \n",
    "\n",
    "```\n",
    "min_c c^T x\n",
    "under the constrains \n",
    "A_eq x = b_eq\n",
    "```\n",
    "\n",
    "The solution is given by\n",
    "`c_star = linprog(c=c, A_eq=A_eq, b_eq=b_eq)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d58b407f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def basis_pursuit(X, y):\n",
    "    \"\"\" Using linprog, solve\n",
    "    min_b |b| such that Xb = y where | | is the l1 norm.\n",
    "\n",
    "    Hint: BP can be recast as a linear program with n equality constraints\n",
    "    and 2p variables. Use variables for the positive and negative part\n",
    "    of each entry of b\n",
    "\n",
    "    Parameters\n",
    "    -----------\n",
    "\n",
    "    X : array of size (n, p)\n",
    "        Design matrix\n",
    "    y : array of size (n,)\n",
    "        Observations\n",
    "\n",
    "    Return\n",
    "    -------\n",
    "    sol: array of size p\n",
    "        The solution of basis pursuit\n",
    "    \"\"\"\n",
    "    raise NotImplementedError\n",
    "\n",
    "def test_basis_pursuit():\n",
    "    X = np.random.randn(10, 20)\n",
    "    b = np.zeros(20)\n",
    "    b[0] = 1\n",
    "    b[1] = 2\n",
    "    y = X.dot(b)\n",
    "    np.testing.assert_array_almost_equal(b, basis_pursuit(X, y))\n",
    "\n",
    "test_basis_pursuit()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4702a053",
   "metadata": {},
   "source": [
    "# Question 6\n",
    "\n",
    "Take p = 40\n",
    "Take k from 1 to p\n",
    "Take n from 1 to p\n",
    "Generate 2 datasets with the above parameters\n",
    "Solve with basis pursuit\n",
    "\n",
    "- Build the matrix `data` of size $p, p$\n",
    "where `data[i, j]` contains the number of times (0, 1 or 2) \n",
    "that basis pursuit succeeds when $k=i$ and $n=j$.\n",
    "\n",
    "Execute `test_basis_pursuit_success` to see the phase transition\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12e013b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_basis_pursuit_success():\n",
    "    \"\"\"\n",
    "    Take p = 40\n",
    "    Take k from 1 to p\n",
    "    Take n from 1 to p\n",
    "    Generate 2 datasets with the above parameters\n",
    "    Solve with basis pursuit\n",
    "    Store the result in a matrix of size (p, p) that counts\n",
    "    the number of times it succeeds for each value of n and k.\n",
    "\n",
    "    Return\n",
    "    -------\n",
    "    data : np array of size (p, p)\n",
    "        data[i, j] contains the number of success for n = i and k=j\n",
    "    \"\"\"\n",
    "\n",
    "    raise NotImplementedError\n",
    "\n",
    "\n",
    "data = compute_basis_pursuit_success()\n",
    "\n",
    "def test_basis_pursuit_success():\n",
    "    \"\"\"\n",
    "    Plot a heatmap of the result of basis pursuit\n",
    "    and check visually that it matches the theoretical results\n",
    "    \"\"\"\n",
    "\n",
    "    phi = norm.pdf\n",
    "    Phi = norm.cdf\n",
    "\n",
    "    F = lambda t: 2 * ((1 + t**2) * Phi(-t) - t * phi(t))\n",
    "    delta = lambda t: 4 * phi(t) / (2 * (t * (1 - 2 * Phi(-t)) + 2 * phi(t)))\n",
    "    s_max = lambda t: (delta(t) - F(t)) / ((1 + t**2) - F(t))\n",
    "\n",
    "    grid = np.linspace(0, 10, num=500)\n",
    "\n",
    "    plt. figure()\n",
    "    plt.imshow(data, extent=[0, 1, 1, 0])\n",
    "    plt.plot(s_max(grid), delta(grid))\n",
    "    plt.ylabel(\"sample size / dimension\")\n",
    "    plt.xlabel(\"sparsity / dimension\")\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "test_basis_pursuit_success()"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "formats": "py:percent,ipynb"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}