edge_differential_privacy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "patent-broadcasting",
   "metadata": {},
   "source": [
    "# Edge Differential Privacy\n",
    "Reference implementations of techniques described in the paper [Smooth Sensitivity and Sampling in Private Data Analysis](https://cs-people.bu.edu/ads22/pubs/NRS07/NRS07-full-draft-v1.pdf) by Kobbi Nissim, Sofya Raskhodnikova, and Adam Smith."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "incomplete-solomon",
   "metadata": {},
   "source": [
    "## Library Imports\n",
    "We use [networkx](https://networkx.org) to perform the required graph computations and the implementation of the Cauchy mechanisms provided by [RelM](https://github.com/anusii/RelM) to release the differentially private query responses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "suspected-batman",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import networkx as nx\n",
    "from relm.mechanisms import CauchyMechanism"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "hybrid-magazine",
   "metadata": {},
   "source": [
    "<hr style=\"border:2px solid gray\"> </hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "expected-beginning",
   "metadata": {},
   "source": [
    "## Triangle Counts\n",
    "We construct a differentially private release mechanism for the number of triangles in a graph. This process is comprised of three steps:\n",
    "  1. Compute the exact query response,\n",
    "  2. Compute the smooth sensitivity of the query,\n",
    "  3. Add noise scaled according to the smooth sensitvity."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "committed-standard",
   "metadata": {},
   "source": [
    "### Compute the exact query response"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "introductory-isaac",
   "metadata": {},
   "source": [
    "#### Generate a random graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "major-trick",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 2**13\n",
    "p = 0.01\n",
    "g = nx.random_graphs.gnp_random_graph(n, p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "expired-granny",
   "metadata": {},
   "source": [
    "#### Compute the exact triangle count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "everyday-cooking",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Notice here that we divide by three becuase nx.triangles returns the\n",
    "# sum of the number of triangles each vertex is in.  Therefore, each triangle\n",
    "# gets counted three times, once for each of its vertices.\n",
    "triangle_count = np.array([sum(nx.triangles(g).values()) / 3.0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "macro-justice",
   "metadata": {},
   "source": [
    "### Compute the smooth sensitivity of the query"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fluid-operations",
   "metadata": {},
   "source": [
    "#### Compute the partial count matrices\n",
    "The algorithm that we use to compute the smooth sensitivity of the query is given in terms of the adjacency matrix $X = \\{x_{ij}\\}$ where $x_{ij} = \\mathbf{1}((i,j) \\in E)$. Let the matrix $A$ count the number of triangles that involve each potential edge.  That is, $A = \\{a_{ij}\\}$ where $a_{ij} = \\sum_{k \\in [n]} x_{ik} \\cdot x_{kj}$. Let the matrix $B$ count the number of half-built triangles that involve each potential edge. That is $B = \\{b_{ij}\\}$ where $b_{ij} = \\sum_{k \\in [n]} x_{ik} \\oplus x_{kj}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "sixth-winter",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the adjacency matrix\n",
    "X = nx.linalg.graphmatrix.adjacency_matrix(g).astype(np.int)\n",
    "\n",
    "# Compute A and B using matrix operations\n",
    "A = (X @ X).todense()\n",
    "B = X.sum(0) + X.sum(1) - 2 * A\n",
    "\n",
    "# Zero out the main diagonal because we are\n",
    "# interested only in indices (i,j) with i != j\n",
    "np.fill_diagonal(A, 0)\n",
    "np.fill_diagonal(B, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "loving-albuquerque",
   "metadata": {},
   "source": [
    "#### Compute the local sensitivity at distance $s$ for $0 \\leq s \\leq n$\n",
    "We recall from the paper that if $S_{f,\\epsilon}^*(G)$ is the smooth sensitivity of a query at an input $G$, $LS_f(H)$ is the local sensitivity of $f$ at an input $H$, and $$A^{(s)}(G) = \\max_{H: d(G,H) \\leq s} LS_f(H)$$\n",
    "is the sensitivity of $f$ at distance $s$, then we have\n",
    "$$S_{f,\\epsilon}^*(G) = \\max_{0 \\leq s \\leq n} e^{-\\epsilon s} A^{(s)}(G).$$\n",
    "\n",
    "Furthermore, Nissim et al show that if $f$ counts the number of triangles in a graph then we have\n",
    "$$A^{(s)}(G) = \\max_{0 \\leq i \\neq j \\leq n} c_{ij}(s) \\quad \\text{where} \\quad c_{ij} = \\min \\left(a_{ij} + \\left\\lfloor\\frac{s + \\min(s,b_{ij})}{2}\\right\\rfloor, n-2 \\right).$$\n",
    "\n",
    "They then describe an $O(n^2 \\log n)$ algorithm for computing $A^{(s)}(G)$ for $0 \\leq s \\leq n$ which works by efficiently identifying the pairs $(a_{ij}, b_{ij})$ needed to compute $c_{ij}(s)$ for all $s$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "complicated-sculpture",
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_maximal_pairs(A, B):\n",
    "    \"\"\"\n",
    "    Find maximal pairs A[i,j], B[i,j].\n",
    "    \n",
    "    Parameters:\n",
    "        A: A matrix describing the number of potential triangles involving each potential edge (i,j)\n",
    "        B: A matrix describing the number of half-built triangles involving each potential edge (i,j)\n",
    "        \n",
    "    Returns:\n",
    "        A generator of pairs (a, b) where for for each distinct value a in set(A)\n",
    "        b = max(B[i,j]) where the maximum is taken over all indices {i,j} where A[i,j] = a.\n",
    "    \"\"\"\n",
    "    W = np.array(A).flatten()\n",
    "    X = np.array(B).flatten()\n",
    "    idx = np.argsort(W)\n",
    "    Y = W[idx]\n",
    "    Z = X[idx]\n",
    "    delta = np.concatenate((np.zeros(1), Y[1:] != Y[:-1]))\n",
    "    new_val_idxs = np.concatenate((np.array([0]),\n",
    "                                   np.where(delta)[0],\n",
    "                                   np.array([len(Y)])))\n",
    "    S = np.empty((len(new_val_idxs)-1), dtype=np.int)\n",
    "    T = np.empty((len(new_val_idxs)-1), dtype=np.int)\n",
    "    for i in range(len(new_val_idxs)-1):\n",
    "        S[i] = Y[new_val_idxs[i]]\n",
    "        T[i] = np.max(Z[new_val_idxs[i]:new_val_idxs[i+1]])\n",
    "    return zip(S[::-1], T[::-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "japanese-blowing",
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_survivors(maximal_pairs, n):\n",
    "    \"\"\"\n",
    "    Find (i,j) used to compute local sensitivity at distance 0 <= s <= n.\n",
    "    \n",
    "    First identify all maximal pairs (A[i,j], B[i,j]).\n",
    "    Then find the maximal pairs that will produce the correct value for the\n",
    "    local sensitivity at distance s computations.\n",
    "    \n",
    "    Parameters:\n",
    "        maximal_pairs: A generator of maximal pairs as returned by find_maximal_pairs\n",
    "        n: The number of nodes in the underlying graph\n",
    "        \n",
    "    Returns:\n",
    "        A sorted list of survivors that can be used to compute the local sensitivity\n",
    "        at each distance 0 <= s <= n.\n",
    "    \"\"\"\n",
    "    prev_survivor = next(maximal_pairs)\n",
    "    break_points = {prev_survivor: n+1}\n",
    "    for survivor in maximal_pairs:\n",
    "        a0, b0 = prev_survivor\n",
    "        a1, b1 = survivor\n",
    "        intersection = 2*(a0 - a1) + b0\n",
    "        if b0 <= intersection <= b1:\n",
    "            break_points[prev_survivor] = intersection\n",
    "            break_points[survivor] = n+1\n",
    "            prev_survivor = survivor\n",
    "    survivors = sorted(break_points.items(), key=lambda _: _[1])\n",
    "    return survivors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "ultimate-nirvana",
   "metadata": {},
   "outputs": [],
   "source": [
    "def local_sensitivity_at_distance(s, break_list):\n",
    "    \"\"\"\n",
    "    Compute the local sensitivity at distance s.\n",
    "    \n",
    "    Parameters:\n",
    "        s: The distance at which to compute the local sensitivity.\n",
    "        break_list: A dictionary of candidate (key, value) pairs that\n",
    "                    are used to compute local sensitivities at various distances.\n",
    "                    \n",
    "    Returns:\n",
    "        The local sensitivity at distance s.\n",
    "    \"\"\"\n",
    "    a, b = next((k for k,v in break_list if s <= v))\n",
    "    return np.minimum(a + np.floor((s + np.minimum(s, b)) / 2.0), n - 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "crude-stanford",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the local sensitivity at distance s for 0 <= s <= n\n",
    "maximal_pairs = find_maximal_pairs(A, B)\n",
    "survivors = find_survivors(maximal_pairs, n)\n",
    "lsd = np.array([local_sensitivity_at_distance(s, survivors) for s in range(n+1)])\n",
    "\n",
    "# Compute the smooth sensitivity\n",
    "epsilon = 1.0\n",
    "smooth_scaling = np.exp(-epsilon * np.arange(n + 1))\n",
    "smooth_sensitivity = np.max(lsd * smooth_scaling)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "beautiful-convenience",
   "metadata": {},
   "source": [
    "### Add noise scaled according to the smooth sensitvity\n",
    "We create a differentially private release mechanism by adding noise from the Cauchy distribution scaled according to the smooth sensitivity. Because the Cauchy distributed random variables are real-valued, the differentially private query response will be real-valued despite the exact query response being integer-valued."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "outside-documentation",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a differentially private release mechanism\n",
    "mechanism = CauchyMechanism(epsilon=epsilon)\n",
    "\n",
    "# Compute the differentially private query response\n",
    "dp_triangle_count = mechanism.release(triangle_count, smooth_sensitivity)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "appointed-clarity",
   "metadata": {},
   "source": [
    "#### Display results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "korean-death",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact triangle count = 91765\n",
      "Differentially private triangle count = 91763.400858\n"
     ]
    }
   ],
   "source": [
    "print(\"Exact triangle count = %i\" % int(triangle_count))\n",
    "print(\"Differentially private triangle count = %f\" % dp_triangle_count)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "industrial-clearance",
   "metadata": {},
   "source": [
    "<hr style=\"border:2px solid gray\"> </hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "overhead-butter",
   "metadata": {},
   "source": [
    "## Cost of a Minimum Spanning Tree\n",
    "We construct a differentially private release mechanism for the cost of a minimum spanning tree. As before, this process is comprised of three steps:\n",
    "  1. Compute the exact query response,\n",
    "  2. Compute the smooth sensitivity of the query,\n",
    "  3. Add noise scaled according to the smooth sensitvity."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "finnish-identity",
   "metadata": {},
   "source": [
    "### Compute the exact query response"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "integral-replication",
   "metadata": {},
   "source": [
    "#### Generate a random graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "suspended-toronto",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 2 ** 8\n",
    "p = 0.1\n",
    "g = nx.random_graphs.gnp_random_graph(n=n, p=p)\n",
    "\n",
    "bound = 10.0  # An upper bound on the edge weights in the graph\n",
    "weights = {e: bound * np.random.randint(1, 11) / 10.0 for e in g.edges()}\n",
    "nx.set_edge_attributes(g, weights, \"weight\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "rising-combine",
   "metadata": {},
   "source": [
    "#### Compute the exact cost of a minimum spanning tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "difficult-iraqi",
   "metadata": {},
   "outputs": [],
   "source": [
    "mst = nx.minimum_spanning_tree(g)\n",
    "mst_cost = np.array([mst.size(weight=\"weight\")])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "endless-consensus",
   "metadata": {},
   "source": [
    "### Compute the smooth sensitivity of the query"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "iraqi-schedule",
   "metadata": {},
   "source": [
    "As per Nissim et al, let $G = (V, E)$ be a graph with vertex set $V$ and edge set $E$. For every $S \\subset V$ we can define a cut to be the partition $(S, V \\setminus S)$. We say an edge $(i,j)$ crosses the cut $S$ when $i \\in S$ and $j \\in V \\setminus S$. For a cut $S \\subset V$, let $w_t(S)$ denote the weight of the $t$-lightest edge in the cut. In the paper, the authors show that\n",
    "$$ A^{(k)}(G) = \\max \\left(\\max_{S \\subset V} w_{k+1}(S), \\max_{S \\subset V} (w_{k+2}(S) - w_1(S))\\right).$$\n",
    "Furthermore, the authors show that if $T$ is a minimum spanning tree of $G$ then\n",
    "$$ \\max_{S \\subset V} \\left(w_{k+2}(S) - w_1(S))\\right) = \\max_{e \\in T}\\left(\\max_{e-\\text{cuts} \\ S} w_{k+2}(S) - w_1(S)\\right)$$\n",
    "where an $(i,j)$-cut is a cut $S$ such that $(i,j)$ crosses $S$.\n",
    "\n",
    "The authors then describe an efficient method for computing $\\max_S w_k(S)$. We implement this algorithm with the function `find_max_cut_weight` below. We use the same algorithm to compute the max cut weight over $e$-cuts. This allows us to efficiently compute $A^{(k)}$ for all $k$ and thereby compute\n",
    "$$S_{f,\\epsilon}^*(G) = \\max_{0 \\leq k \\leq n} e^{-\\epsilon k} A^{(k)}(G).$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "editorial-services",
   "metadata": {},
   "outputs": [],
   "source": [
    "def retrieve_min_cut_cost(costs, key, g, w, **args):\n",
    "    \"\"\" Return the minimum cut cost for an unweighted graph related to g.\n",
    "    \n",
    "        If a key is in costs, then return costs[key]\n",
    "        Otherwise, generate the appropriate value for costs[key],\n",
    "        add this new (key, value pair) to costs and then\n",
    "        return costs[key].\n",
    "    \"\"\"\n",
    "    if key not in costs:\n",
    "        gw = nx.Graph()\n",
    "        gw.add_nodes_from(g)\n",
    "        gw.add_edges_from([e for e in g.edges() if g.edges[e][\"weight\"] <= w[key]])\n",
    "        costs[key] = nx.edge_connectivity(gw, **args)\n",
    "    return costs[key]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "formal-washer",
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_max_cut_weight(k, w, bound, costs, **kwargs):\n",
    "    \"\"\"\n",
    "    Compute max_{e-cuts S} w_k(S).\n",
    "    If kwargs is none then this is the same as max_{S \\subset V} w_k(S)\n",
    "    \n",
    "    Parameters:\n",
    "        k: The index of the weight of interest\n",
    "        w: A list of sorted edge-weights\n",
    "        bound: An upper bound for the edge-weights in w\n",
    "        costs: A dictionary of min_cut_costs\n",
    "        \n",
    "        kwargs (optional): Two nodes that define an edge e which\n",
    "                           must be separated by any cut under consideration.\n",
    "            s: a source node\n",
    "            t: a destination node\n",
    "    \n",
    "    Returns:\n",
    "        max_{e-cuts S} w_k(S)\n",
    "    \"\"\"\n",
    "    cost = retrieve_min_cut_cost(costs, len(w) - 1, g, w, **kwargs)\n",
    "    if cost <= k:\n",
    "        return bound\n",
    "\n",
    "    cost = retrieve_min_cut_cost(costs, 0, g, w, **kwargs)\n",
    "    if cost > k:\n",
    "        return w[0]\n",
    "\n",
    "    high = len(w) - 1\n",
    "    low = 0\n",
    "    while (high - low) > 1:\n",
    "        mid = (low + high) // 2\n",
    "        cost = retrieve_min_cut_cost(costs, mid, g, w, **kwargs)\n",
    "        if cost <= k:\n",
    "            low = mid\n",
    "        else:\n",
    "            high = mid\n",
    "    return w[high]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "inside-shanghai",
   "metadata": {},
   "source": [
    "#### Compute $\\max_{S \\subset V} w_{k+1}(S)$ for all $0 \\leq k \\leq n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "coastal-offer",
   "metadata": {},
   "outputs": [],
   "source": [
    "edge_weights = [g.edges[e][\"weight\"] for e in g.edges()]\n",
    "edge_weights = sorted(set(edge_weights))\n",
    "costs = dict()\n",
    "# Notice here that we use k instead of k+1 to correct for 0-up indexing in python\n",
    "lsd1 = np.array([find_max_cut_weight(k, edge_weights, bound, costs) for k in range(n + 1)])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "animal-lafayette",
   "metadata": {},
   "source": [
    "#### Compute $\\max_{e \\in T} \\left(\\max_{e-\\text{cuts} \\ S} \\left(w_{k+2} - w_1(S)\\right)\\right)$ for all $0 \\leq k \\leq n$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "unlikely-behalf",
   "metadata": {},
   "outputs": [],
   "source": [
    "mst = nx.minimum_spanning_tree(g)\n",
    "lsd2 = np.zeros(n + 1)\n",
    "for e in mst.edges():\n",
    "    costs = dict()\n",
    "    # Notice here that we use k+1 instead of k+2 to correct for 0-up indexing in python\n",
    "    max_cut_weights = [\n",
    "        find_max_cut_weight(k + 1, edge_weights, bound, costs, s=e[0], t=e[1])\n",
    "        for k in range(n + 1)\n",
    "    ]\n",
    "    lsd2_e = np.array(max_cut_weights) - mst.edges[e][\"weight\"]\n",
    "    lsd2 = np.maximum(lsd2, lsd2_e)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "juvenile-hunter",
   "metadata": {},
   "source": [
    "#### Compute smooth sensitivity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "dimensional-crest",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the local sensitivity at distance s for 0 <= s <= n\n",
    "lsd = np.maximum(lsd1, lsd2)\n",
    "\n",
    "# Compute the smooth sensitivity\n",
    "epsilon = 1.0\n",
    "smooth_scaling = np.exp(-epsilon * np.arange(n + 1))\n",
    "smooth_sensitivity = np.max(lsd * smooth_scaling)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "under-spirituality",
   "metadata": {},
   "source": [
    "### Add noise scaled according to the smooth sensitvity\n",
    "We create a differentially private release mechanism by adding noise from the Cauchy distribution scaled according to the smooth sensitivity. Because the Cauchy distributed random variables are real-valued, the differentially private query response will be real-valued despite the exact query response being integer-valued."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "standard-settle",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a differentiall private release mechanism\n",
    "mechanism = CauchyMechanism(epsilon=epsilon)\n",
    "\n",
    "# Compute the differentially private query response\n",
    "dp_mst_cost = mechanism.release(mst_cost, smooth_sensitivity)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "norman-hamburg",
   "metadata": {},
   "source": [
    "#### Display results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "printable-spokesman",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact MST cost = 279.000000\n",
      "Differentially private MST cost = 279.447695\n"
     ]
    }
   ],
   "source": [
    "print(\"Exact MST cost = %f\" % mst_cost)\n",
    "print(\"Differentially private MST cost = %f\" % dp_mst_cost)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "other-manufacturer",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}