Posts

• Jan 17, 2018

  Dynamic-Graph Orientation

  tags:
  TheoryLunch Blog

  Computer science topics over lunch

  Dynamic-Graph Orientation

  Jan 17, 2018 • Eran Amar

  tags: $\newcommand{\lyxlock}{}$

  1 Orienting Acyclic Graphs

  Throughout this post, $G=\left(V,E\right)$ is an undirected acyclic graph on $n$ vertices and $m\le n-1$ edges. The problem we want to solve is called graph orientation:

  The input $G$ is given in the insertion-only streaming model. That is, we know $V$ in advance, and the set $E$ is being revealed one edge at a time. For each edge $\left(u,v\right)\in E$ that we get, we need to assign a direction (i.e. an orientation), meaning to decide that it either points to $u$ or to $v$ . When $u$ points to $v$ we denote it by $u\rightarrow v$ . Note that when we assign a direction to a new edge, we may as well update the orientation of previously seen edges. When the stream ends, we left with a directed version of $G$ , denote it $\hat{G}=\left(V,\hat{E}\right)$ .

  We are interested in two metrics: the first is the outgoing degree bound, which is defined by the term $$deg\left(\hat{G}\right):=\max_{v\in V}deg_{out}\left(v;\hat{G}\right)$$ where $deg_{out}(v;\hat{G}):=\left|\left\{ u\in V\mid v\rightarrow u\right\} \right|$ is the number of $v$ ’s neighbors which $v$ points to, in the oriented graph $\hat{G}$ .

  The second metric is a bound on the orientation cost for a single step of the stream. That is, a bound on the number of edge updates for every time a new edge is being revealed.

  Due to the fact that $G$ is acyclic, we know it is actually a collection of one or more trees, then there exists an orientation $\hat{G}$ in which $deg\left(\hat{G}\right)\le1$ (think that each node points to its parent in its corresponding tree).

  2 $\mathcal{O}\left(1\right)$ Orientation Cost

  Lets start with a very simple orientation rule: For every step of the stream, let $e_{i}=\left(u,v\right)$ be the edge given in the current step, and $E_{i}=\left\{ e_{1},..,e_{i}\right\}$ . Consider the endpoints of $e_{i}$ . If $u$ belongs to a larger connected component relative to the connected component of $v$ , under the graph induced by $E_{i}$ , then orient $e_{i}$ to be $v\rightarrow u$ , otherwise orient it $u\rightarrow v$ . Note that $u$ and $v$ cannot be in the same connected component in $G_{i}=\left(V,E_{i}\right)$ .

  Obviously, that update rule has constant orientation cost (in particular it never change the orientations of previous edges). However, the outgoing degree bound is $\mathcal{O}\left(\log n\right)$ . To prove that, consider some arbitrary vertex $v$ . We increase its outgoing degree by $1$ when we make it points to another larger connected component. So every time the outgoing degree of $v$ increases by $1$ , the size of its connected component is at least doubled. Since we have $n$ vertices, we cannot double that size more than $\log_{2}n$ times.

  3 $\mathcal{O}\left(1\right)$ Outgoing Degree Bound

  Now we will see an update rule that yields a constant degree bound, say $4$ :

  Upon receiving $e_{i}=\left(u,v\right)$ , orient it $v\rightarrow u$ . Then, if the outgoing degree of $v$ exceeds $4$ , we apply $makeSink\left(v\right)$ . This procedure will flip the orientation of any outgoing edge of $v$ . Next, check any neighbor $u$ that was affected by this change, and apply $makeSink\left(u\right)$ if it has more than $4$ outgoing edges. This process continues recursively until no orientation is changed.

  Interesting question that immediately pops to mind is, why this process of recursive updates eventually terminates?

  To see that infinite oscillation is impossible, consider some optimal orientation of $G$ , denote it $\hat{G}^{*}$ , in which $deg\left(\hat{G}^{*}\right)\le1$ . Suppose that $e_{i}=\left(v,u_{1}\right)$ was revealed to the stream, and let $u_{1},u_{2},...$ be the chain of vertices that we applied $makeSink\left(\cdot\right)$ to, after orienting it. Denote by $\phi_{k}$ the number of edges in $G_{i}=\left(V,E_{i}\right)$ that don’t agree with $\hat{G}^{*}$ , after applying $makeSink\left(\cdot\right)$ on $u_{1},..,u_{k}$ . So for example, $\phi_{0}$ is the number of conflicting edges right after we inserted $e_{i}$ but before we start the recursive fixing process, and $\phi_{1}$ is the number of conflicting edges after we flip all the outgoing edges of $u_{1}$ , and so on.

  We will show that $\phi_{0},\phi_{1},..$ is a strictly decreasing series, and since $\phi_{k}$ cannot be negative, that series must terminates. Fix some $\phi_{k}$ , which is the number of edges conflicting with $\hat{G}^{*}$ before applying $makeSink\left(u_{k+1}\right)$ . In that moment, $deg_{out}\left(u_{k+1}\right)>4$ , but in the optimal orientation it should be at most $1$ , which means that at least $4$ of $u_{k+1}$ ’s outgoing edges have wrong orientations. So applying $makeSink\left(u_{k+1}\right)$ adding up to one conflicting edge, but on the other hand, fixing the orientation of at least $4$ other edges, thus the total number of conflicting edges with $\hat{G}^{*}$ decreased by at least $3$ , i.e. $\phi_{k}\ge3+\phi_{k+1}$ .

  It is immediate that this algorithm yields an orientation with outgoing degree bound of $4$ . We now move to show that the amortize orientation cost is also constant. To show that we will use the accounting method, in which we assign positive and negative values to different operations that we perform during processing the stream. As long as those values are constants, if in every step of the stream the aggregated value is non-negative, then we may conclude that the amortize cost per operation is also constant. So suppose that the action of inserting a new edge has a value of $+2$ , and flipping orientation of an existing edge has a value of $-1$ . We first prove the following claim.

  Suppose we got a new edge $e_{i}$ , and let $u_{1},u_{2},...,u_{k}$ be the series of vertices that we applied on $makeSink$ during the recursive fixing process, after inserting $e_{i}$ to the graph. Then, all $u_{1},u_{2},...,u_{k}$ are distinct, and the value of $makeSink\left(u_{t}\right)$ is exactly $-5$ for any $t\in\left\{ 1,...,k\right\}$ .

  For the first part, if there is a vertex in $u_{1},u_{2},...,u_{k}$ appearing twice, it means that we found a cycle in $G$ which is impossible. For the second part, assume toward contradiction that we have $t\in\left\{ 1,..,k\right\}$ such that the outgoing degree of $u_{t}$ was at least 6. It must be due to $makeSink\left(\cdot\right)$ application on some vertices from $\left\{ u_{1},...,u_{t-1}\right\}$ , which means that $u_{t}$ has two neighbors in $\left\{ u_{1},...,u_{t-1}\right\}$ , and we know that those neighbors are connected, so we found a cycle in $G$ and that is a contradiction.

  It is left to show that at any point in time, the total value of previous operation is non negative. The following claim prove that.

  For every step, $i$ , of the stream, let $e_{i}$ be the edge revealed on that step. Let $u_{i,1},u_{i,2},...,u_{i,k_{i}}$ be the series of vertices that we applied on $makeSink$ during the recursive fixing process, with $\phi_{i,0,}...,\phi_{i,k_{i}}$ as defined before. Denote by $\psi_{i,t}$ the aggregated value of all the operations done until inserting $e_{i}$ (including) and after applying $makeSink$ to $\left\{ u_{1},..,u_{t}\right\}$ . Then, $\psi_{i,t}\ge2\cdot\phi_{i,t}$ for any $t$ .

  By induction on $i$ . For $i=1$ , note that $\phi_{1,0}=0$ and $\psi_{1,0}=3$ . For $i>1$ , we have $\psi_{i,0}=\psi_{i-1,k_{i-1}}+2$ and $\phi_{i,0}\le\phi_{i-1,k_{i-1}}+1$ , so using the induction hypothesis, $$\psi_{i,0}\ge2+2\cdot\phi_{i-1,k_{i-1}}\ge2\cdot\phi_{i,0}$$ for $t>0$ , assume by induction on $t$ that the claim holds, then $$\psi_{i,t}=\psi_{i,t-1}-5\ge2\cdot\phi_{i,t-1}-5\ge2\left(\phi_{i,t}+3\right)-5\ge2\phi_{i,t}$$

  4 Releasing the Assumptions

  This post dealt with a relative simple model, in which the stream only revealing new edges but never deleting a previously seen one. In addition, we assume that our graph is acyclic which is also a very strong assumption. There are works that shows more sophisticated streaming algorithm for fully dynamic graphs (i.e. with edges insertions and deletions), and also works on general graphs that are not necessarily acyclic. Further details can be found in this paper.

  $\newcommand{\lyxlock}{}$

  1 Orienting Acyclic Graphs

  Throughout this post, $G=\left(V,E\right)$ is an undirected acyclic graph on $n$ vertices and $m\le n-1$ edges. The problem we want to solve is called graph orientation:

  The input $G$ is given in the insertion-only streaming model. That is, we know $V$ in advance, and the set $E$ is being revealed one edge at a time. For each edge $\left(u,v\right)\in E$ that we get, we need to assign a direction (i.e. an orientation), meaning to decide that it either points to $u$ or to $v$ . When $u$ points to $v$ we denote it by $u\rightarrow v$ . Note that when we assign a direction to a new edge, we may as well update the orientation of previously seen edges. When the stream ends, we left with a directed version of $G$ , denote it $\hat{G}=\left(V,\hat{E}\right)$ .

  We are interested in two metrics: the first is the outgoing degree bound, which is defined by the term $$deg\left(\hat{G}\right):=\max_{v\in V}deg_{out}\left(v;\hat{G}\right)$$ where $deg_{out}(v;\hat{G}):=\left|\left\{ u\in V\mid v\rightarrow u\right\} \right|$ is the number of $v$ ’s neighbors which $v$ points to, in the oriented graph $\hat{G}$ .

  The second metric is a bound on the orientation cost for a single step of the stream. That is, a bound on the number of edge updates for every time a new edge is being revealed.

  Due to the fact that $G$ is acyclic, we know it is actually a collection of one or more trees, then there exists an orientation $\hat{G}$ in which $deg\left(\hat{G}\right)\le1$ (think that each node points to its parent in its corresponding tree).

  2 $\mathcal{O}\left(1\right)$ Orientation Cost

  Lets start with a very simple orientation rule: For every step of the stream, let $e_{i}=\left(u,v\right)$ be the edge given in the current step, and $E_{i}=\left\{ e_{1},..,e_{i}\right\}$ . Consider the endpoints of $e_{i}$ . If $u$ belongs to a larger connected component relative to the connected component of $v$ , under the graph induced by $E_{i}$ , then orient $e_{i}$ to be $v\rightarrow u$ , otherwise orient it $u\rightarrow v$ . Note that $u$ and $v$ cannot be in the same connected component in $G_{i}=\left(V,E_{i}\right)$ .

  Obviously, that update rule has constant orientation cost (in particular it never change the orientations of previous edges). However, the outgoing degree bound is $\mathcal{O}\left(\log n\right)$ . To prove that, consider some arbitrary vertex $v$ . We increase its outgoing degree by $1$ when we make it points to another larger connected component. So every time the outgoing degree of $v$ increases by $1$ , the size of its connected component is at least doubled. Since we have $n$ vertices, we cannot double that size more than $\log_{2}n$ times.

  3 $\mathcal{O}\left(1\right)$ Outgoing Degree Bound

  Now we will see an update rule that yields a constant degree bound, say $4$ :

  Upon receiving $e_{i}=\left(u,v\right)$ , orient it $v\rightarrow u$ . Then, if the outgoing degree of $v$ exceeds $4$ , we apply $makeSink\left(v\right)$ . This procedure will flip the orientation of any outgoing edge of $v$ . Next, check any neighbor $u$ that was affected by this change, and apply $makeSink\left(u\right)$ if it has more than $4$ outgoing edges. This process continues recursively until no orientation is changed.

  Interesting question that immediately pops to mind is, why this process of recursive updates eventually terminates?

  To see that infinite oscillation is impossible, consider some optimal orientation of $G$ , denote it $\hat{G}^{*}$ , in which $deg\left(\hat{G}^{*}\right)\le1$ . Suppose that $e_{i}=\left(v,u_{1}\right)$ was revealed to the stream, and let $u_{1},u_{2},...$ be the chain of vertices that we applied $makeSink\left(\cdot\right)$ to, after orienting it. Denote by $\phi_{k}$ the number of edges in $G_{i}=\left(V,E_{i}\right)$ that don’t agree with $\hat{G}^{*}$ , after applying $makeSink\left(\cdot\right)$ on $u_{1},..,u_{k}$ . So for example, $\phi_{0}$ is the number of conflicting edges right after we inserted $e_{i}$ but before we start the recursive fixing process, and $\phi_{1}$ is the number of conflicting edges after we flip all the outgoing edges of $u_{1}$ , and so on.

  We will show that $\phi_{0},\phi_{1},..$ is a strictly decreasing series, and since $\phi_{k}$ cannot be negative, that series must terminates. Fix some $\phi_{k}$ , which is the number of edges conflicting with $\hat{G}^{*}$ before applying $makeSink\left(u_{k+1}\right)$ . In that moment, $deg_{out}\left(u_{k+1}\right)>4$ , but in the optimal orientation it should be at most $1$ , which means that at least $4$ of $u_{k+1}$ ’s outgoing edges have wrong orientations. So applying $makeSink\left(u_{k+1}\right)$ adding up to one conflicting edge, but on the other hand, fixing the orientation of at least $4$ other edges, thus the total number of conflicting edges with $\hat{G}^{*}$ decreased by at least $3$ , i.e. $\phi_{k}\ge3+\phi_{k+1}$ .

  It is immediate that this algorithm yields an orientation with outgoing degree bound of $4$ . We now move to show that the amortize orientation cost is also constant. To show that we will use the accounting method, in which we assign positive and negative values to different operations that we perform during processing the stream. As long as those values are constants, if in every step of the stream the aggregated value is non-negative, then we may conclude that the amortize cost per operation is also constant. So suppose that the action of inserting a new edge has a value of $+2$ , and flipping orientation of an existing edge has a value of $-1$ . We first prove the following claim.

  Suppose we got a new edge $e_{i}$ , and let $u_{1},u_{2},...,u_{k}$ be the series of vertices that we applied on $makeSink$ during the recursive fixing process, after inserting $e_{i}$ to the graph. Then, all $u_{1},u_{2},...,u_{k}$ are distinct, and the value of $makeSink\left(u_{t}\right)$ is exactly $-5$ for any $t\in\left\{ 1,...,k\right\}$ .

  For the first part, if there is a vertex in $u_{1},u_{2},...,u_{k}$ appearing twice, it means that we found a cycle in $G$ which is impossible. For the second part, assume toward contradiction that we have $t\in\left\{ 1,..,k\right\}$ such that the outgoing degree of $u_{t}$ was at least 6. It must be due to $makeSink\left(\cdot\right)$ application on some vertices from $\left\{ u_{1},...,u_{t-1}\right\}$ , which means that $u_{t}$ has two neighbors in $\left\{ u_{1},...,u_{t-1}\right\}$ , and we know that those neighbors are connected, so we found a cycle in $G$ and that is a contradiction.

  It is left to show that at any point in time, the total value of previous operation is non negative. The following claim prove that.

  For every step, $i$ , of the stream, let $e_{i}$ be the edge revealed on that step. Let $u_{i,1},u_{i,2},...,u_{i,k_{i}}$ be the series of vertices that we applied on $makeSink$ during the recursive fixing process, with $\phi_{i,0,}...,\phi_{i,k_{i}}$ as defined before. Denote by $\psi_{i,t}$ the aggregated value of all the operations done until inserting $e_{i}$ (including) and after applying $makeSink$ to $\left\{ u_{1},..,u_{t}\right\}$ . Then, $\psi_{i,t}\ge2\cdot\phi_{i,t}$ for any $t$ .

  By induction on $i$ . For $i=1$ , note that $\phi_{1,0}=0$ and $\psi_{1,0}=3$ . For $i>1$ , we have $\psi_{i,0}=\psi_{i-1,k_{i-1}}+2$ and $\phi_{i,0}\le\phi_{i-1,k_{i-1}}+1$ , so using the induction hypothesis, $$\psi_{i,0}\ge2+2\cdot\phi_{i-1,k_{i-1}}\ge2\cdot\phi_{i,0}$$ for $t>0$ , assume by induction on $t$ that the claim holds, then $$\psi_{i,t}=\psi_{i,t-1}-5\ge2\cdot\phi_{i,t-1}-5\ge2\left(\phi_{i,t}+3\right)-5\ge2\phi_{i,t}$$

  4 Releasing the Assumptions

  This post dealt with a relative simple model, in which the stream only revealing new edges but never deleting a previously seen one. In addition, we assume that our graph is acyclic which is also a very strong assumption. There are works that shows more sophisticated streaming algorithm for fully dynamic graphs (i.e. with edges insertions and deletions), and also works on general graphs that are not necessarily acyclic. Further details can be found in this paper.
• Oct 8, 2017

  Last TheoryLunch at Weizmann - Cut Equivalent Trees

  tags:
  TheoryLunch Blog

  Computer science topics over lunch

  Last TheoryLunch at Weizmann - Cut Equivalent Trees

  Oct 8, 2017 • Eran Amar

  tags: $\newcommand{\lyxlock}{}$

  1 Last TheoryLunch Meeting at Weizmann Institute of Science

  For those of you who don’t know, this blog originated from my lunch meetings with my friends at Weizmann Institute of Science during my Master studies. As all good things comes to an end, also did my time at Weizmann, and last week our lunch group gathered for the last time. This post presents the topic discussed on that meeting.

  I want to thank my good friends, Yosef Pogrow and Ron Shiff, for the great time we had and for all the enriching lunch meetings (also those which didn’t include theory). Without you guys, this blog would not have been created, I wish you both best of luck ahead!

  This post, by the way, contains some results from the research of Yosef Pogrow (who also presented most of the meetings).

  2 Min-Cuts in a Graph

  This post deals with undirected, weighted graphs and cuts on them. We denote by $G=\left(V,E_{G},w_{G}\right)$ the graph equipped with a weight function on its edges. A cut in this graph is a partition of its vertices to $\left(S,V\backslash S\right)$ where $S\subseteq V$ . For simplicity, we will only mention one of the parts in the partition, usually $S$ . Note that it doesn’t matter whether we refer to the first part in the partition $S$ , or to the second part $V\backslash S$ , because both parts define the same cut.

  The value of the cut is defined as $$w_{G}\left(S\right):=\sum_{uv\in E_{G},u\in S,v\notin S}w_{G}\left(uv\right)$$

  A min $st$ -cut for a pair $s\ne t\in V$ is a cut defined by $S\subseteq V$ such that $$S=\arg\min_{S'\subseteq V}w_{G}\left(S'\right)$$ note that it is possible that several cuts will satisfy the definition above, but for our purposes they are equivalent (since they all have the same value).

  Before diving into the theorem about cut equivalent trees, we setup the following notation. For a given tree $T=\left(V,E\right)$ and any edge $uv\in E$ , removing that edge from the tree breaks it down into two connected components, one that contains $u$ , and another one that contains $v$ . Denote the component containing $u$ as $S_{T|u}\subseteq V$ , and then the other component is exactly $S_{T|v}:=V\backslash S_{T|u}$ . Since the only path connecting $u$ and $v$ was the removed edge, then $S_{T|u}$ is actually the min $uv$ -cut in $T$ (the same is true for $S_{T|v}$ since it defines the same cut).

  A simple fact about trees is, that for any tree $T$ and vertices $s\ne t$ inside it, the min $st$ -cut is given by $S_{T|u}$ where $uv\in E_{T}$ is the edge with minimum weight along the path that connects $s$ and $t$ . Since any $n$ vertices tree has exactly $n-1$ edges, there are at most $n-1$ unique values of min-cuts that the tree can yield.

  Cut-Equivalent Trees [Gomory and Hu, 1961]. For any weighted undirected graph $G=\left(V,E_{G},w_{G}\right)$ there is a weighted tree on the same set of vertices $T=\left(V,E_{T},w_{T}\right)$ such that:
  (a) for any $s\ne t\in V$ , the value of the min $st$ -cut in $G$ is equal the value of the min $st$ -cut in $T$ (with respect to their respective weight functions), and
  (b) for any edge $uv\in E_{T}$ , $$w_{T}\left(S_{T|u}\right)=w_{G}\left(S_{T|u}\right)$$ Such trees are known as Cut-Equivalent trees (also known as GH-trees).

  The proof is omitted from this post, and we move to discuss a nice implication of that theorem: it puts a limit on the number of different min-cut values that a graph can have. Since there are ${n \choose 2}$ different $st$ pairs in a graph there can be potentially ${n \choose 2}$ unique min-cut values, however, the theorem says that any value of a min-cut in $G$ is being represented by a (possibly different) min-cut in $T$ , and there are at most $n-1$ min-cut values in $T$ , then $G$ has at most $n-1$ unique min-cut values.

  3 General Cuts in a Graph

  Until now, we talked about minimum $st$ cuts. The theorem above state that for any graph $G$ , there is a cut-equivalent tree $T$ that can tell us much about all the minimum cuts in $G$ . It can tell us all the possible min-cut values that $G$ can have, moreover, for any min-cut value $x$ , there is a very easy way to find a cut in $G$ that have this exact value. All we need to do is to find an edge in the tree that have weight $x$ , suppose that edge is $uv\in E_{T}$ , and the desired cut in $G$ is given by $S_{T|u}$ .

  One interesting question to ask is, what cut-equivalent trees can tell us about “general” cuts? The answer is given in the following theorem.

  [Krauthgamer and Pogrow, 2017]. Let $G=\left(V,E_{G},w_{G}\right)$ be a weighted undirected graph with $n$ vertices, and let $T=\left(V,E_{T},w_{T}\right)$ be a cut-equivalent tree of $G$ . Then for any $S\subseteq V$ $$w_{G}\left(S\right)\le w_{T}\left(S\right)\le\left(n-1\right)\cdot w_{G}\left(S\right)$$

  We start with the left-hand side. Recall that the value of a cut $S$ is the sum of weights of all the “crossing edges”, meaning, edges which have endpoints in both parts of the partition. Then, combined with the second property of the tree, and the fact that $w_{T}\left(st\right)=w_{T}\left(S_{T|s}\right)$ for any edge $st\in E_{T}$ , we get $$\begin{aligned} w_{T}\left(S\right) & =\sum_{st\in E_{T},s\in S,t\notin S}w_{T}\left(st\right)\\ & =\sum_{st\in E_{T},s\in S,t\notin S}w_{G}\left(S_{T|s}\right) \end{aligned}$$ so it is suffice to show that any crossing edge in the sum of the original graph $w_{G}\left(S\right)$ , contributes also to a sum of $w_{G}\left(S_{T|s}\right)$ for some tree-edge $st\in E_{T}$ such that $s\in S,t\notin S$ . So fix some crossing edge $uv\in E_{G}$ ($u\in S,v\in V\backslash S$ ) from the original graph, and consider the unique path connecting $u$ and $v$ over the tree. Since we starting the path with a vertex from $S$ and ending with a vertex in $V\backslash S$ , there must be an edge $st\in E_{T}$ along the path such that $s\in S$ and $t\in V\backslash S$ . That edge $st$ separates $u$ from $v$ in the tree, so $u\in S_{T|s}$ and $v\notin S_{T|s}$ . Since $uv\in E_{G}$ then it contributes to $w_{G}\left(S_{T|s}\right)$ , which is what we wanted to prove.
  We proceed to show the second inequality. For each $st\in E_{T}$ in the sum of $w_{T}\left(S\right)$ , we know that $w_{T}\left(st\right)\le w_{G}\left(S\right)$ since the weight of the tree-edge $w_{T}\left(st\right)$ is equal to the min $st$ -cut in $G$ , thus in particular smaller or equal to the value of the cut $S$ (which also separates $s$ from $t$ but might not be of minimum weight). The inequality follows by summing over all the crossing edges in $S$ over $T$ , which is at most $n-1$ edges.

  4 The Inequalities are Tight

  We conclude this post with examples for cases in which one of the inequalities became an equality.

  For the first inequality, take any graph $G$ and a corresponding cut equivalent tree $T$ . Fix an arbitrary tree-edge $st\in E_{T}$ . The cut $S=S_{T|s}$ is the min $st$ -cut in $T$ and by the first property of the tree $w_{G}\left(S\right)=w_{T}\left(S\right)$ .

  For the second inequality, consider the complete graph on $n$ vertices with all its edges weighted $1$ . Its cut-equivalent tree must be a star in which every edge weighed $n-1$ . See example of an eight vertices star in the figure below (taken from Wikipedia).

  

  Suppose the vertex in the center of the star is $v$ . The cut $S=\left\{ v\right\}$ has weight of $n-1$ over $G$ , however, it has weight $\left(n-1\right)^{2}$ in the star $T$ which is exactly $\left(n-1\right)\cdot w_{G}\left(S\right)$ .

  $\newcommand{\lyxlock}{}$

  1 Last TheoryLunch Meeting at Weizmann Institute of Science

  For those of you who don’t know, this blog originated from my lunch meetings with my friends at Weizmann Institute of Science during my Master studies. As all good things comes to an end, also did my time at Weizmann, and last week our lunch group gathered for the last time. This post presents the topic discussed on that meeting.

  I want to thank my good friends, Yosef Pogrow and Ron Shiff, for the great time we had and for all the enriching lunch meetings (also those which didn’t include theory). Without you guys, this blog would not have been created, I wish you both best of luck ahead!

  This post, by the way, contains some results from the research of Yosef Pogrow (who also presented most of the meetings).

  2 Min-Cuts in a Graph

  This post deals with undirected, weighted graphs and cuts on them. We denote by $G=\left(V,E_{G},w_{G}\right)$ the graph equipped with a weight function on its edges. A cut in this graph is a partition of its vertices to $\left(S,V\backslash S\right)$ where $S\subseteq V$ . For simplicity, we will only mention one of the parts in the partition, usually $S$ . Note that it doesn’t matter whether we refer to the first part in the partition $S$ , or to the second part $V\backslash S$ , because both parts define the same cut.

  The value of the cut is defined as $$w_{G}\left(S\right):=\sum_{uv\in E_{G},u\in S,v\notin S}w_{G}\left(uv\right)$$

  A min $st$ -cut for a pair $s\ne t\in V$ is a cut defined by $S\subseteq V$ such that $$S=\arg\min_{S'\subseteq V}w_{G}\left(S'\right)$$ note that it is possible that several cuts will satisfy the definition above, but for our purposes they are equivalent (since they all have the same value).

  Before diving into the theorem about cut equivalent trees, we setup the following notation. For a given tree $T=\left(V,E\right)$ and any edge $uv\in E$ , removing that edge from the tree breaks it down into two connected components, one that contains $u$ , and another one that contains $v$ . Denote the component containing $u$ as $S_{T|u}\subseteq V$ , and then the other component is exactly $S_{T|v}:=V\backslash S_{T|u}$ . Since the only path connecting $u$ and $v$ was the removed edge, then $S_{T|u}$ is actually the min $uv$ -cut in $T$ (the same is true for $S_{T|v}$ since it defines the same cut).

  A simple fact about trees is, that for any tree $T$ and vertices $s\ne t$ inside it, the min $st$ -cut is given by $S_{T|u}$ where $uv\in E_{T}$ is the edge with minimum weight along the path that connects $s$ and $t$ . Since any $n$ vertices tree has exactly $n-1$ edges, there are at most $n-1$ unique values of min-cuts that the tree can yield.

  Cut-Equivalent Trees [Gomory and Hu, 1961]. For any weighted undirected graph $G=\left(V,E_{G},w_{G}\right)$ there is a weighted tree on the same set of vertices $T=\left(V,E_{T},w_{T}\right)$ such that:
  (a) for any $s\ne t\in V$ , the value of the min $st$ -cut in $G$ is equal the value of the min $st$ -cut in $T$ (with respect to their respective weight functions), and
  (b) for any edge $uv\in E_{T}$ , $$w_{T}\left(S_{T|u}\right)=w_{G}\left(S_{T|u}\right)$$ Such trees are known as Cut-Equivalent trees (also known as GH-trees).

  The proof is omitted from this post, and we move to discuss a nice implication of that theorem: it puts a limit on the number of different min-cut values that a graph can have. Since there are ${n \choose 2}$ different $st$ pairs in a graph there can be potentially ${n \choose 2}$ unique min-cut values, however, the theorem says that any value of a min-cut in $G$ is being represented by a (possibly different) min-cut in $T$ , and there are at most $n-1$ min-cut values in $T$ , then $G$ has at most $n-1$ unique min-cut values.

  3 General Cuts in a Graph

  Until now, we talked about minimum $st$ cuts. The theorem above state that for any graph $G$ , there is a cut-equivalent tree $T$ that can tell us much about all the minimum cuts in $G$ . It can tell us all the possible min-cut values that $G$ can have, moreover, for any min-cut value $x$ , there is a very easy way to find a cut in $G$ that have this exact value. All we need to do is to find an edge in the tree that have weight $x$ , suppose that edge is $uv\in E_{T}$ , and the desired cut in $G$ is given by $S_{T|u}$ .

  One interesting question to ask is, what cut-equivalent trees can tell us about “general” cuts? The answer is given in the following theorem.

  [Krauthgamer and Pogrow, 2017]. Let $G=\left(V,E_{G},w_{G}\right)$ be a weighted undirected graph with $n$ vertices, and let $T=\left(V,E_{T},w_{T}\right)$ be a cut-equivalent tree of $G$ . Then for any $S\subseteq V$ $$w_{G}\left(S\right)\le w_{T}\left(S\right)\le\left(n-1\right)\cdot w_{G}\left(S\right)$$

  We start with the left-hand side. Recall that the value of a cut $S$ is the sum of weights of all the “crossing edges”, meaning, edges which have endpoints in both parts of the partition. Then, combined with the second property of the tree, and the fact that $w_{T}\left(st\right)=w_{T}\left(S_{T|s}\right)$ for any edge $st\in E_{T}$ , we get $$\begin{aligned} w_{T}\left(S\right) & =\sum_{st\in E_{T},s\in S,t\notin S}w_{T}\left(st\right)\\ & =\sum_{st\in E_{T},s\in S,t\notin S}w_{G}\left(S_{T|s}\right) \end{aligned}$$ so it is suffice to show that any crossing edge in the sum of the original graph $w_{G}\left(S\right)$ , contributes also to a sum of $w_{G}\left(S_{T|s}\right)$ for some tree-edge $st\in E_{T}$ such that $s\in S,t\notin S$ . So fix some crossing edge $uv\in E_{G}$ ($u\in S,v\in V\backslash S$ ) from the original graph, and consider the unique path connecting $u$ and $v$ over the tree. Since we starting the path with a vertex from $S$ and ending with a vertex in $V\backslash S$ , there must be an edge $st\in E_{T}$ along the path such that $s\in S$ and $t\in V\backslash S$ . That edge $st$ separates $u$ from $v$ in the tree, so $u\in S_{T|s}$ and $v\notin S_{T|s}$ . Since $uv\in E_{G}$ then it contributes to $w_{G}\left(S_{T|s}\right)$ , which is what we wanted to prove.
  We proceed to show the second inequality. For each $st\in E_{T}$ in the sum of $w_{T}\left(S\right)$ , we know that $w_{T}\left(st\right)\le w_{G}\left(S\right)$ since the weight of the tree-edge $w_{T}\left(st\right)$ is equal to the min $st$ -cut in $G$ , thus in particular smaller or equal to the value of the cut $S$ (which also separates $s$ from $t$ but might not be of minimum weight). The inequality follows by summing over all the crossing edges in $S$ over $T$ , which is at most $n-1$ edges.

  4 The Inequalities are Tight

  We conclude this post with examples for cases in which one of the inequalities became an equality.

  For the first inequality, take any graph $G$ and a corresponding cut equivalent tree $T$ . Fix an arbitrary tree-edge $st\in E_{T}$ . The cut $S=S_{T|s}$ is the min $st$ -cut in $T$ and by the first property of the tree $w_{G}\left(S\right)=w_{T}\left(S\right)$ .

  For the second inequality, consider the complete graph on $n$ vertices with all its edges weighted $1$ . Its cut-equivalent tree must be a star in which every edge weighed $n-1$ . See example of an eight vertices star in the figure below (taken from Wikipedia).

  

  Suppose the vertex in the center of the star is $v$ . The cut $S=\left\{ v\right\}$ has weight of $n-1$ over $G$ , however, it has weight $\left(n-1\right)^{2}$ in the star $T$ which is exactly $\left(n-1\right)\cdot w_{G}\left(S\right)$ .
• Jul 18, 2017

  On the Hardness of Counting Problems - Part 2

  tags:
  TheoryLunch Blog

  Computer science topics over lunch

  On the Hardness of Counting Problems - Part 2

  Jul 18, 2017 • Eran Amar

  tags: $\newcommand{\lyxlock}{}$

  1 Recap

  The previous post discussed the $\mathcal{P},\mathcal{NP}$ and $\#\mathcal{P}$ families. We said that the hierarchy between decision problems from $\mathcal{P},\mathcal{NP}$ do not preserved when moving to their Counting versions. We defined $CNF$ and $DNF$ formulas (see section 2 in previous post), and the problems CNF-SAT that is $\mathcal{NP}$ -complete and DNF-SAT that is in $\mathcal{P}$ . We used the relation between their Counting versions to show that $\#DNF$ and $\#CNF$ are equally hard to be solved exactly. The relation is as follows: if $\phi$ is a $n$ variables $CNF$ formula, then $$\#CNF\left(\phi\right)+\#DNF\left(\neg\phi\right)=2^{n}$$ Furthermore, we showed that both problems are $\#\mathcal{P}$ -complete, demonstrating that even an “easy” decision problem like DNF-SAT can be very hard when considering its Counting counterpart.

  This post we will show that there is some hierarchy inside $\#\mathcal{P}$ when considering approximated solutions instead of exact solutions.

  2 Approximation Algorithm for #DNF-SAT

  2.1 Deterministic Approximation

  Let $\psi=\bigvee_{i=1}^{m}C_{i}$ be a $DNF$ formula, and $k^{*}=\#DNF\left(\psi\right)$ . We want to find a polynomial time algorithm that yields some $k$ which is an approximation for $k^{*}$ . Formally, for an approximation factor of $p>1$ , $k^{*}\le k\le p\cdot k^{*}$ . We will talk more about the approximation parameter later.

  First note that if a clause contains a variable and its negation then immediately we can conclude that there are no satisfying assignments (s.a.) for that clause as it contains a contradiction, so we will omit it from the formula (can be done in linear time). So from now on, assume that $\psi$ do not contain any contradiction. For the rest of this discussion, we say that a boolean variable appear in a clause $C$ if either it or its negation contained in $C$ .

  Let’s focus on an arbitrary clause $C_{i}\in\psi$ . We can easily count the number of s.a. it has; it is exactly $\#DNF\left(C_{i}\right)=2^{t}$ when $t\in\left\{ 0,1,..,n-1\right\}$ is the number of variables that do not appear in $C_{i}$ . That is true because the value of any variable in $C_{i}$ is determined by its appearance (if it has negation or not) so only variables which are not in $C_{i}$ can be either $True$ or $False$ and having two options.

  Moreover, we can efficiently sample at uniform any s.a. for that clause. Simply iterate over the $n$ variables, if a variable appears in $C_{i}$ then it have only one allowed value so assign it that value, and if it is not appear in $C_{i}$ then toss a coin and assign $True$ or $False$ accordingly.

  Trivial approximation algorithm will calculate the exact $\#DNF\left(C_{i}\right)$ for each $C_{i}$ in $\psi$ and returns their sum $k=\sum_{i=1}^{m}\#DNF\left(C_{i}\right)$ . However, we may have s.a. that were counted more than once (if they are satisfying more than a single clause). Thus the approximation factor of that algorithm is $m$ , because in the worst case all the clauses have the same set of satisfying assignments so we counted each s.a. $m$ times, yielding $k\le m\cdot k^{*}$ .

  We can improve that algorithm by using the inclusion-exclusion principle: instead of accumulating the number of s.a. for each clause, accumulate the number of s.a. that weren’t count so far.

  Unfortunately, finding which assignment was already counted forces us to enumerate all possible assignments, which obviously takes exponential time. Hence we will give a probabilistic estimation for that, instead. When working with estimations, we need to state the success probability for our estimation and it is also preferable that we will have a way to increase that probability (mostly at the cost of more computations and drawing more random samples).

  Next section we will describe how to get an estimation of $k$ such that $\left(1-p\right)k^{*}\le k\le\left(1+p\right)k^{*}$ for any $p$ , even a polynomially small one (i.e. $p=\frac{1}{poly\left(n\right)}$ ) with a very high probability.

  2.2 Probabilistic Approximation

  The overall procedure is still simple: iterate over each clause $C_{i}$ and estimate the fraction of s.a. for $C_{i}$ that were not counted before (as a number between zero and one), denote that estimation by $\tilde{k}_{i}$ , let $S_{i}$ be the total number of s.a. for $C_{i}$ calculated as mentioned in the deterministic approximation section. Then output $k=\sum_{i=1}^{m}\tilde{k}_{i}\cdot S_{i}$ .

  It left to show how to compute $\tilde{k}_{i}$ , so we do the following. Let $t=c\left(\frac{m}{p}\right)^{2}\ln m$ , when $c$ will be determined later. Sample $t$ independent times, an uniformly random s.a. for $C_{i}$ . For each assignment define a corresponding binary variable $X_{j}$ that is $1$ if it is not satisfying any previous clause, and $0$ otherwise. Note that evaluating all $\left\{ X_{j}\right\} _{j=1}^{t}$ takes $\mathcal{O}\left(poly\left(n,m\right)\cdot t\right)$ time. Then return $\tilde{k}_{i}=\frac{1}{t}\sum_{j=1}^{t}X_{j}$ . That completes our probabilistic approximation algorithm for $\#DNF$ .

  Analysis

  We need to analyze the quality of approximation of each $\tilde{k_{i}}$ , then show that with high probability $\sum_{i=1}^{m}\tilde{k}_{i}\cdot S_{i}$ is within the claimed approximation bounds. For the first part, we focus on a single estimation $\tilde{k}_{i}$ for the clause $C_{i}$ .

  Let $S_{i}$ be the number of all s.a. for $C_{i}$ and let $s_{i}$ be the number of s.a. for $C_{i}$ that aren’t satisfying any clause before $C_{i}$ . Out goal is to show that $\tilde{k}_{i}$ is, with high probability, close to $\frac{s_{i}}{S_{i}}$ . Observe that all the variables $\left\{ X_{j}\right\} _{j=1}^{t}$ for estimating $\tilde{k}_{i}$ were actually sampled i.i.d from a Bernoulli distribution with parameter $q_{i}=\frac{s_{i}}{S_{i}}$ (which is also its expectation). As a result $\mathbb{E}\left[\tilde{k}_{i}\right]=q_{i}$ , so apply Hoeffding inequality to get $$\begin{aligned} \mathbf{P}\left[\left|\tilde{k}_{i}-q_{i}\right|>\frac{p}{m}\right] & \le2\exp\left(\frac{-2\left(\frac{p}{m}\right)^{2}t^{2}}{t}\right)\\ & =2\exp\left(-2c\ln m\right)=2m^{-2c} \end{aligned}$$ which means that with high probability the estimator is “successful”, that is, $$q_{i}-\frac{p}{m}\le\tilde{k}_{i}\le q_{i}+\frac{p}{m}$$ and equivalently, $$s_{i}-S_{i}\frac{p}{m}\le\tilde{k}_{i}\cdot S_{i}\le s_{i}+S_{i}\frac{p}{m}$$ If all the estimators $\tilde{k}_{1},\tilde{k}_{2},..$ are successful, then $$k\le\sum_{i=1}^{m}\left(s_{i}+S_{i}\frac{p}{m}\right)=k^{*}+\sum_{i=1}^{m}S_{i}\frac{p}{m}\le k^{*}+\frac{p}{m}\left(mk^{*}\right)=\left(1+p\right)k^{*}$$ and also similarly $k\ge\left(1-p\right)k^{*}$ , which are exactly the approximation bounds that we want to achieve. To conclude the analysis, note that by applying the union bound, all the estimators will be successful with probability $\ge1-\sum_{i=1}^{m}2m^{-2c}=1-\frac{2}{m^{2c-1}}$ .

  Trade off between the parameters

  It is interesting to examine the relation between the different parameters, and how each affect the running time, success probability, and the quality of approximation.

  Consider the parameter $t$ , which is the number of samples the algorithm has to draw to estimate each $\tilde{k}_{i}$ . That $t$ depends linearly on the parameter $c$ which controls the “success probability” of all the estimators (simultaneously). On the one hand, large $c$ gives a success probability that is extremely close to $1$ . But on the other hand, it increase the number of samples $t$ per clause, and the total running time of the algorithm, which is $\mathcal{O}\left(poly\left(n,m\right)\cdot c\cdot p^{-2}\right)$ .

  Similarly, the better approximation we want for $k^{*}$ , the smaller $p$ that we should set (even smaller than $1$ ), then the number of samples and the running time increase inversely proportional to $p^{2}$ .

  3 What About Approximating #CNF-SAT?

  Previous section dealt with approximating, either in the deterministic or probabilistic manner, the $\#DNF$ problem. Unfortunately, as oppose to that we did with exact solutions, approximating $\#DNF$ do not help us approximating $\#CNF$ . We will prove that by proving a stronger claim.

  Let $P$ be any $\mathcal{NP}$ -complete problem, and denote $\#P$ its Counting versions. There is no polynomial-time approximation algorithm of $\#P$ to any factor, unless $\mathcal{P}=\mathcal{NP}$ .

  Assume that $\mathcal{P}\ne\mathcal{NP}$ . Suppose toward contradiction that we have some polynomial-time approximation algorithm for $\#P$ , and denote it by $Alg$ . Given an instance $x$ for the decision problem $P$ , simulate $Alg\left(x\right)$ and get an approximation for $\#P\left(x\right)$ . Then declare $No$ as the output of the decision problem if that approximation is zero, and $Yes$ otherwise. Note that this reduction runs in polynomial time. Moreover, the Counting versions have zero solutions if and only if $P\left(x\right)=No$ . Since any multiplicative approximation of zero is still zero, we can use the approximation algorithm for the counting problem to efficiently solve the $\mathcal{NP}$ -complete decision problem, contradicting that $\mathcal{P}\ne\mathcal{NP}$ .

  4 $DNF$ Formulas are Not Always Easy!

  During this and the previous post, we mentioned $DNF$ formulas as something which is easy to solve. We showed that $DNFSAT\in\mathcal{P}$ , and always compared it to its “difficult” variation - the $CNF$ version. However, a $DNF$ formula can also be “difficult”. It just depends what we are asking. For example, deciding if a $DNF$ formula $\psi$ is a tautology (i.e. evaluated to $True$ for any assignment) is an $\mathcal{NP}$ -complete problem. One can show that by noting that $\phi:=\neg\psi$ is a $CNF$ formula, and $\phi$ is a $No$ instance for CNF-SAT if and only if $\psi$ is a $Yes$ instance for the DNF-TAUTOLOGY problem. Thus, the problem are equally hard.

  Similarly, we can define an “easy” problem for $CNF$ formulas: decide whether a $CNF$ formula $\phi$ is a tautology can be done in polynomial time, by checking if its negation $\psi:=\neg\phi$ has a s.a. which is polynomial since $DNFSAT\in\mathcal{P}$ and $\psi$ is a $DNF$ formula.

  As a conclusion, the hardness of those decision problems depends not only on the form of the formula but also on what are we asking to decide.

  $\newcommand{\lyxlock}{}$

  1 Recap

  The previous post discussed the $\mathcal{P},\mathcal{NP}$ and $\#\mathcal{P}$ families. We said that the hierarchy between decision problems from $\mathcal{P},\mathcal{NP}$ do not preserved when moving to their Counting versions. We defined $CNF$ and $DNF$ formulas (see section 2 in previous post), and the problems CNF-SAT that is $\mathcal{NP}$ -complete and DNF-SAT that is in $\mathcal{P}$ . We used the relation between their Counting versions to show that $\#DNF$ and $\#CNF$ are equally hard to be solved exactly. The relation is as follows: if $\phi$ is a $n$ variables $CNF$ formula, then $$\#CNF\left(\phi\right)+\#DNF\left(\neg\phi\right)=2^{n}$$ Furthermore, we showed that both problems are $\#\mathcal{P}$ -complete, demonstrating that even an “easy” decision problem like DNF-SAT can be very hard when considering its Counting counterpart.

  This post we will show that there is some hierarchy inside $\#\mathcal{P}$ when considering approximated solutions instead of exact solutions.

  2 Approximation Algorithm for #DNF-SAT

  2.1 Deterministic Approximation

  Let $\psi=\bigvee_{i=1}^{m}C_{i}$ be a $DNF$ formula, and $k^{*}=\#DNF\left(\psi\right)$ . We want to find a polynomial time algorithm that yields some $k$ which is an approximation for $k^{*}$ . Formally, for an approximation factor of $p>1$ , $k^{*}\le k\le p\cdot k^{*}$ . We will talk more about the approximation parameter later.

  First note that if a clause contains a variable and its negation then immediately we can conclude that there are no satisfying assignments (s.a.) for that clause as it contains a contradiction, so we will omit it from the formula (can be done in linear time). So from now on, assume that $\psi$ do not contain any contradiction. For the rest of this discussion, we say that a boolean variable appear in a clause $C$ if either it or its negation contained in $C$ .

  Let’s focus on an arbitrary clause $C_{i}\in\psi$ . We can easily count the number of s.a. it has; it is exactly $\#DNF\left(C_{i}\right)=2^{t}$ when $t\in\left\{ 0,1,..,n-1\right\}$ is the number of variables that do not appear in $C_{i}$ . That is true because the value of any variable in $C_{i}$ is determined by its appearance (if it has negation or not) so only variables which are not in $C_{i}$ can be either $True$ or $False$ and having two options.

  Moreover, we can efficiently sample at uniform any s.a. for that clause. Simply iterate over the $n$ variables, if a variable appears in $C_{i}$ then it have only one allowed value so assign it that value, and if it is not appear in $C_{i}$ then toss a coin and assign $True$ or $False$ accordingly.

  Trivial approximation algorithm will calculate the exact $\#DNF\left(C_{i}\right)$ for each $C_{i}$ in $\psi$ and returns their sum $k=\sum_{i=1}^{m}\#DNF\left(C_{i}\right)$ . However, we may have s.a. that were counted more than once (if they are satisfying more than a single clause). Thus the approximation factor of that algorithm is $m$ , because in the worst case all the clauses have the same set of satisfying assignments so we counted each s.a. $m$ times, yielding $k\le m\cdot k^{*}$ .

  We can improve that algorithm by using the inclusion-exclusion principle: instead of accumulating the number of s.a. for each clause, accumulate the number of s.a. that weren’t count so far.

  Unfortunately, finding which assignment was already counted forces us to enumerate all possible assignments, which obviously takes exponential time. Hence we will give a probabilistic estimation for that, instead. When working with estimations, we need to state the success probability for our estimation and it is also preferable that we will have a way to increase that probability (mostly at the cost of more computations and drawing more random samples).

  Next section we will describe how to get an estimation of $k$ such that $\left(1-p\right)k^{*}\le k\le\left(1+p\right)k^{*}$ for any $p$ , even a polynomially small one (i.e. $p=\frac{1}{poly\left(n\right)}$ ) with a very high probability.

  2.2 Probabilistic Approximation

  The overall procedure is still simple: iterate over each clause $C_{i}$ and estimate the fraction of s.a. for $C_{i}$ that were not counted before (as a number between zero and one), denote that estimation by $\tilde{k}_{i}$ , let $S_{i}$ be the total number of s.a. for $C_{i}$ calculated as mentioned in the deterministic approximation section. Then output $k=\sum_{i=1}^{m}\tilde{k}_{i}\cdot S_{i}$ .

  It left to show how to compute $\tilde{k}_{i}$ , so we do the following. Let $t=c\left(\frac{m}{p}\right)^{2}\ln m$ , when $c$ will be determined later. Sample $t$ independent times, an uniformly random s.a. for $C_{i}$ . For each assignment define a corresponding binary variable $X_{j}$ that is $1$ if it is not satisfying any previous clause, and $0$ otherwise. Note that evaluating all $\left\{ X_{j}\right\} _{j=1}^{t}$ takes $\mathcal{O}\left(poly\left(n,m\right)\cdot t\right)$ time. Then return $\tilde{k}_{i}=\frac{1}{t}\sum_{j=1}^{t}X_{j}$ . That completes our probabilistic approximation algorithm for $\#DNF$ .

  Analysis

  We need to analyze the quality of approximation of each $\tilde{k_{i}}$ , then show that with high probability $\sum_{i=1}^{m}\tilde{k}_{i}\cdot S_{i}$ is within the claimed approximation bounds. For the first part, we focus on a single estimation $\tilde{k}_{i}$ for the clause $C_{i}$ .

  Let $S_{i}$ be the number of all s.a. for $C_{i}$ and let $s_{i}$ be the number of s.a. for $C_{i}$ that aren’t satisfying any clause before $C_{i}$ . Out goal is to show that $\tilde{k}_{i}$ is, with high probability, close to $\frac{s_{i}}{S_{i}}$ . Observe that all the variables $\left\{ X_{j}\right\} _{j=1}^{t}$ for estimating $\tilde{k}_{i}$ were actually sampled i.i.d from a Bernoulli distribution with parameter $q_{i}=\frac{s_{i}}{S_{i}}$ (which is also its expectation). As a result $\mathbb{E}\left[\tilde{k}_{i}\right]=q_{i}$ , so apply Hoeffding inequality to get $$\begin{aligned} \mathbf{P}\left[\left|\tilde{k}_{i}-q_{i}\right|>\frac{p}{m}\right] & \le2\exp\left(\frac{-2\left(\frac{p}{m}\right)^{2}t^{2}}{t}\right)\\ & =2\exp\left(-2c\ln m\right)=2m^{-2c} \end{aligned}$$ which means that with high probability the estimator is “successful”, that is, $$q_{i}-\frac{p}{m}\le\tilde{k}_{i}\le q_{i}+\frac{p}{m}$$ and equivalently, $$s_{i}-S_{i}\frac{p}{m}\le\tilde{k}_{i}\cdot S_{i}\le s_{i}+S_{i}\frac{p}{m}$$ If all the estimators $\tilde{k}_{1},\tilde{k}_{2},..$ are successful, then $$k\le\sum_{i=1}^{m}\left(s_{i}+S_{i}\frac{p}{m}\right)=k^{*}+\sum_{i=1}^{m}S_{i}\frac{p}{m}\le k^{*}+\frac{p}{m}\left(mk^{*}\right)=\left(1+p\right)k^{*}$$ and also similarly $k\ge\left(1-p\right)k^{*}$ , which are exactly the approximation bounds that we want to achieve. To conclude the analysis, note that by applying the union bound, all the estimators will be successful with probability $\ge1-\sum_{i=1}^{m}2m^{-2c}=1-\frac{2}{m^{2c-1}}$ .

  Trade off between the parameters

  It is interesting to examine the relation between the different parameters, and how each affect the running time, success probability, and the quality of approximation.

  Consider the parameter $t$ , which is the number of samples the algorithm has to draw to estimate each $\tilde{k}_{i}$ . That $t$ depends linearly on the parameter $c$ which controls the “success probability” of all the estimators (simultaneously). On the one hand, large $c$ gives a success probability that is extremely close to $1$ . But on the other hand, it increase the number of samples $t$ per clause, and the total running time of the algorithm, which is $\mathcal{O}\left(poly\left(n,m\right)\cdot c\cdot p^{-2}\right)$ .

  Similarly, the better approximation we want for $k^{*}$ , the smaller $p$ that we should set (even smaller than $1$ ), then the number of samples and the running time increase inversely proportional to $p^{2}$ .

  3 What About Approximating #CNF-SAT?

  Previous section dealt with approximating, either in the deterministic or probabilistic manner, the $\#DNF$ problem. Unfortunately, as oppose to that we did with exact solutions, approximating $\#DNF$ do not help us approximating $\#CNF$ . We will prove that by proving a stronger claim.

  Let $P$ be any $\mathcal{NP}$ -complete problem, and denote $\#P$ its Counting versions. There is no polynomial-time approximation algorithm of $\#P$ to any factor, unless $\mathcal{P}=\mathcal{NP}$ .

  Assume that $\mathcal{P}\ne\mathcal{NP}$ . Suppose toward contradiction that we have some polynomial-time approximation algorithm for $\#P$ , and denote it by $Alg$ . Given an instance $x$ for the decision problem $P$ , simulate $Alg\left(x\right)$ and get an approximation for $\#P\left(x\right)$ . Then declare $No$ as the output of the decision problem if that approximation is zero, and $Yes$ otherwise. Note that this reduction runs in polynomial time. Moreover, the Counting versions have zero solutions if and only if $P\left(x\right)=No$ . Since any multiplicative approximation of zero is still zero, we can use the approximation algorithm for the counting problem to efficiently solve the $\mathcal{NP}$ -complete decision problem, contradicting that $\mathcal{P}\ne\mathcal{NP}$ .

  4 $DNF$ Formulas are Not Always Easy!

  During this and the previous post, we mentioned $DNF$ formulas as something which is easy to solve. We showed that $DNFSAT\in\mathcal{P}$ , and always compared it to its “difficult” variation - the $CNF$ version. However, a $DNF$ formula can also be “difficult”. It just depends what we are asking. For example, deciding if a $DNF$ formula $\psi$ is a tautology (i.e. evaluated to $True$ for any assignment) is an $\mathcal{NP}$ -complete problem. One can show that by noting that $\phi:=\neg\psi$ is a $CNF$ formula, and $\phi$ is a $No$ instance for CNF-SAT if and only if $\psi$ is a $Yes$ instance for the DNF-TAUTOLOGY problem. Thus, the problem are equally hard.

  Similarly, we can define an “easy” problem for $CNF$ formulas: decide whether a $CNF$ formula $\phi$ is a tautology can be done in polynomial time, by checking if its negation $\psi:=\neg\phi$ has a s.a. which is polynomial since $DNFSAT\in\mathcal{P}$ and $\psi$ is a $DNF$ formula.

  As a conclusion, the hardness of those decision problems depends not only on the form of the formula but also on what are we asking to decide.
• Jul 6, 2017

  On the Hardness of Counting Problems - Part 1

  tags:
  TheoryLunch Blog

  Computer science topics over lunch

  On the Hardness of Counting Problems - Part 1

  Jul 6, 2017 • Eran Amar

  tags: $\newcommand{\lyxlock}{}$

  1 What is a Counting Problem?

  Every undergraduate student in Computer Science knows the $\mathcal{P}$ and $\mathcal{NP}$ families and the legendary question “is $\mathcal{P}$ equals $\mathcal{NP}$ ?”. However, there is another interesting family that is not always covered in undergraduate courses; that is the $\#\mathcal{P}$ (“sharp P”) family. In this post we will define and discuss this family and two concrete similar-but-different problems from it that demonstrate interesting properties about its hardness.

  Formally, $\mathcal{P}$ is the family of decision problems that can be solved in polynomial time. A known example is the Perfect Matching problem: given a graph, determine whether it has a perfect matching or not. A perfect matching is a set of edges that touches all the vertices in the graph and any pair of edges have disjoint endpoints.

  The family $\mathcal{NP}$ contains all the decision problems that has polynomial time algorithm over the non-deterministic computational model (we won’t get into that). We know that $\mathcal{P}\subseteq\mathcal{NP}$ , but we don’t know if the other direction also holds. A famous $\mathcal{NP}$ -complete example is the Hamiltonian Cycle problem: given a graph, determine whether it contains a cycle that visits every vertex exactly once. We say “complete” to note that the problem is harder at least as any other problem inside that family, that is, if we could solve it efficiently (i.e. in polynomial time) then we could solve any other problem in the family efficiently. When a problem $P_{1}$ is at least harder as other problem $P_{2}$ we denote it $P_{1}\ge P_{2}$ .

  Any problem in either $\mathcal{P}$ or $\mathcal{NP}$ can be rephrased slightly different. We can ask “how many perfect matching this graph contains?” rather than “is there a perfect matching in this graph?”. When asking how much instead of is there the problem become a Counting problem. Following that, the formal definition of $\#\mathcal{P}$ is: the family the Counting version of any problem in $\mathcal{NP}$ .

  Notice that for any problem $P\in\mathcal{P}\cup\mathcal{NP}$ we have $P\le\#P$ because having a count of the desired property can, in particular, distinguish between the cases that this count is $0$ or higher, that is, decide whether that property exists or not. Hence, we can safely state the general hierarchy, that $\mathcal{P}\le\mathcal{NP}\le\#\mathcal{P}$ .

  Let $P_{1},P_{2}$ be the Perfect Matching and Hamiltonian Cycle problems respectively, and let $\#P_{1},\#P_{2}$ be their Counting versions. It is interesting to see whether the hierarchy between problems from $\mathcal{P}$ and $\mathcal{NP}$ preserved for their Counting versions as well. For instance, for the problems above, we know that $P_{1}\le P_{2}$ since $\mathcal{P}\subseteq\mathcal{NP}$ and $P_{2}$ is $\mathcal{NP}$ -complete, is that implies also $\#P_{1}\le\#P_{2}$ ?

  The short answer is no. In the next section we will show that in a detailed discussion.

  2 Counting Makes Problems Much Harder

  We start with two basic definitions from the field of boolean logic.

  A Conjunctive Normal Form ($CNF$ ) formula $\phi$ other a set of $n$ boolean variables is defined as $$\phi=\bigwedge_{i=1}^{m}C_{i}$$ where each of its $m$ clauses is of the form $C_{i}=\bigvee_{j=1}^{\left|C_{i}\right|}\ell_{i,j}$ , that is, the clause is a disjunction of some number of literals. A literal $\ell_{i,j}$ is a variable or its negation. Note that since there are $n$ variables then $\left|C_{i}\right|\le2n$ .

  A Disjunctive Normal Form ($DNF$ ) formula $\psi$ other a set of $n$ boolean variables is defined analogously to the $CNF$ version, as $$\psi=\bigvee_{i=1}^{m}C_{i}$$ where each of its $m$ clauses is of the form $C_{i}=\bigwedge_{j=1}^{\left|C_{i}\right|}\ell_{i,j}$ .

  The CNF-SAT problem, gets as an input a $CNF$ boolean formula, and decides whether that formula has a satisfying assignment (s.a.) or not. That is, an assignment of $True/False$ for each of its variables such that the formula evaluates to $True$ . The DNF-SAT is defined the same but for $DNF$ formulas rather than $CNF$ ones.

  We will need the following facts in the coming discussion:

  1. The $\neg$ operator (negation) transforms any $CNF$ formula $\phi$ to a $DNF$ formula and vice versa.
  2. For any boolean formula $\phi$ over $n$ variables (regardless its form), there are exactly $2^{n}$ possible assignments each of which satisfy exactly one from $\left\{ \phi,\neg\phi\right\}$ .

  It is very easy to solve DNF-SAT efficiently, meaning, in polynomial time. Given a $DNF$ formula $\psi$ , iterate over its clauses until you find a satisfiable clause then declare that $\psi$ is satisfiable, however, if you failed to find such a clause then declare that $\psi$ is not satisfiable. Observe that any clause of the form $C_{i}=\bigwedge_{j=1}^{\left|C_{i}\right|}\ell_{i,j}$ is satisfiable if and only if it does not contains a variable and its negation; which we can check with a single pass over the literals of $C_{i}$ . As a result, $DNFSAT\in\mathcal{P}$ .

  In contrast to DNF-SAT, the CNF-SAT is know to be $\mathcal{NP}$ -complete, so we believe that it cannot be solved in polynomial time. Similar to the claim from the previous section, we conclude that CNF-SAT is “at least harder” as DNF-SAT, that is, $DNFSAT\le CNFSAT$ . However, as regard to the Counting versions of those problems, we can show the opposite inequality; that $\#CNFSAT\le\#DNFSAT$ .

  To prove that, we need to show that if we could solve $\#DNFSAT$ efficiently, then we could solve $\#CNFSAT$ efficiently. So suppose we have a magical algorithm that can solve $\#DNFSAT$ in polynomial time. Given a $CNF$ formula $\phi$ over $n$ boolean variables, construct $\psi=\neg\phi$ and feed it to the magical algorithm. Let its output be the number $t$ , then return $2^{n}-t$ .

  The correctness of that reduction follows from the two facts we stated earlier. Note that also all the procedures can be executed in polynomial time, so we solved $\#CNFSAT$ in polynomial time using an efficient solver for $\#DNFSAT$ , therefore, $\#CNFSAT\le\#DNFSAT$ .

  Actually, if we will look closely on that reduction, we will see that it should work perfectly for the other direction as well, that is, solving $\#DNFSAT$ by a an efficient algorithm for the $\#CNFSAT$ , so we may conclude that both problems are equally hard. This is an evidence that the hierarchy between decision problems do not preserved for their Counting versions. Putting it differently, we can say that Counting makes problems much harder, such that even an “easy” problems from $\mathcal{P}$ become harder as the entire $\mathcal{NP}$ when you consider its Counting version.

  3 $\#DNF$ is Also $\#\mathcal{P}$ -Complete

  Since $\#DNF$ and $\#CNF$ are equally hard, it is enough to show that $\#CNF$ is $\#\mathcal{P}$ -complete; for that we need a reduction from any problem in $\#\mathcal{P}$ to $\#CNF$ . Recall the Cook-Levin reduction which proves that CNF-SAT is $\mathcal{NP}$ -complete. That reduction can be applied to any $P\in\mathcal{NP}$ and produce a polynomial procedure that converts any input $x$ for the original problem $P$ , to a $CNF$ formula $\phi$ for the problem CNF-SAT, such that $\phi$ is satisfiable if and only if the decision problem $P$ says $Yes$ for the original input $x$ . That way, any efficient solver for CNF-SAT solves $P$ efficiently, so $P\le CNFSAT$ .

  An important property of the Cook-Levin reduction is that it keeps the number of s.a. for $\phi$ exactly as the number of witnesses for the original input $x$ of $P$ . That means, that we can use the same reduction to convert any Counting problem $\#P\in\#\mathcal{P}$ to $\#CNF$ . That is, the reduction will yields a procedure to convert any input $x$ for $\#P$ to a formula $\phi$ for $\#CNF$ such that the number of s.a. of $\phi$ is the solution for $\#P\left(x\right)$ . As a result, $\#CNF$ is $\#\mathcal{P}$ -complete, and then also $\#DNF$ .

  Following the conclusion from previous section, we can say further that even a simple problem from $\mathcal{P}$ , can have a Counting version that is not just “harder” but harder as any other problem in $\#\mathcal{P}$ .

  4 Is There a Structure inside $\#\mathcal{P}$ ?

  Going back to the $\#CNF$ and $\#DNF$ problems, we saw that they are related via the relation $$\#CNF\left(\phi\right)+\#DNF\left(\neg\phi\right)=2^{n}$$ for any $CNF$ formula $\phi$ over $n$ boolean variables. We exploited that relation to show that the problems are equally “hard” when talking about exact-solutions. However, an hierarchy emerged once we moved to consider approximated solutions for those problems. Therefore, there is some kind of structure for the family $\#\mathcal{P}$ under the appropriate notion of “solving” the problems.

  Next post we will discuss an approximation algorithm for $\#DNF$ and the hardness of approximating $\#CNF$ , showing that the additive relation from above it not enough for constructing an approximation algorithm for $\#CNF$ based on approximation for $\#DNF$ .

  $\newcommand{\lyxlock}{}$

  1 What is a Counting Problem?

  Every undergraduate student in Computer Science knows the $\mathcal{P}$ and $\mathcal{NP}$ families and the legendary question “is $\mathcal{P}$ equals $\mathcal{NP}$ ?”. However, there is another interesting family that is not always covered in undergraduate courses; that is the $\#\mathcal{P}$ (“sharp P”) family. In this post we will define and discuss this family and two concrete similar-but-different problems from it that demonstrate interesting properties about its hardness.

  Formally, $\mathcal{P}$ is the family of decision problems that can be solved in polynomial time. A known example is the Perfect Matching problem: given a graph, determine whether it has a perfect matching or not. A perfect matching is a set of edges that touches all the vertices in the graph and any pair of edges have disjoint endpoints.

  The family $\mathcal{NP}$ contains all the decision problems that has polynomial time algorithm over the non-deterministic computational model (we won’t get into that). We know that $\mathcal{P}\subseteq\mathcal{NP}$ , but we don’t know if the other direction also holds. A famous $\mathcal{NP}$ -complete example is the Hamiltonian Cycle problem: given a graph, determine whether it contains a cycle that visits every vertex exactly once. We say “complete” to note that the problem is harder at least as any other problem inside that family, that is, if we could solve it efficiently (i.e. in polynomial time) then we could solve any other problem in the family efficiently. When a problem $P_{1}$ is at least harder as other problem $P_{2}$ we denote it $P_{1}\ge P_{2}$ .

  Any problem in either $\mathcal{P}$ or $\mathcal{NP}$ can be rephrased slightly different. We can ask “how many perfect matching this graph contains?” rather than “is there a perfect matching in this graph?”. When asking how much instead of is there the problem become a Counting problem. Following that, the formal definition of $\#\mathcal{P}$ is: the family the Counting version of any problem in $\mathcal{NP}$ .

  Notice that for any problem $P\in\mathcal{P}\cup\mathcal{NP}$ we have $P\le\#P$ because having a count of the desired property can, in particular, distinguish between the cases that this count is $0$ or higher, that is, decide whether that property exists or not. Hence, we can safely state the general hierarchy, that $\mathcal{P}\le\mathcal{NP}\le\#\mathcal{P}$ .

  Let $P_{1},P_{2}$ be the Perfect Matching and Hamiltonian Cycle problems respectively, and let $\#P_{1},\#P_{2}$ be their Counting versions. It is interesting to see whether the hierarchy between problems from $\mathcal{P}$ and $\mathcal{NP}$ preserved for their Counting versions as well. For instance, for the problems above, we know that $P_{1}\le P_{2}$ since $\mathcal{P}\subseteq\mathcal{NP}$ and $P_{2}$ is $\mathcal{NP}$ -complete, is that implies also $\#P_{1}\le\#P_{2}$ ?

  The short answer is no. In the next section we will show that in a detailed discussion.

  2 Counting Makes Problems Much Harder

  We start with two basic definitions from the field of boolean logic.

  A Conjunctive Normal Form ($CNF$ ) formula $\phi$ other a set of $n$ boolean variables is defined as $$\phi=\bigwedge_{i=1}^{m}C_{i}$$ where each of its $m$ clauses is of the form $C_{i}=\bigvee_{j=1}^{\left|C_{i}\right|}\ell_{i,j}$ , that is, the clause is a disjunction of some number of literals. A literal $\ell_{i,j}$ is a variable or its negation. Note that since there are $n$ variables then $\left|C_{i}\right|\le2n$ .

  A Disjunctive Normal Form ($DNF$ ) formula $\psi$ other a set of $n$ boolean variables is defined analogously to the $CNF$ version, as $$\psi=\bigvee_{i=1}^{m}C_{i}$$ where each of its $m$ clauses is of the form $C_{i}=\bigwedge_{j=1}^{\left|C_{i}\right|}\ell_{i,j}$ .

  The CNF-SAT problem, gets as an input a $CNF$ boolean formula, and decides whether that formula has a satisfying assignment (s.a.) or not. That is, an assignment of $True/False$ for each of its variables such that the formula evaluates to $True$ . The DNF-SAT is defined the same but for $DNF$ formulas rather than $CNF$ ones.

  We will need the following facts in the coming discussion:

  1. The $\neg$ operator (negation) transforms any $CNF$ formula $\phi$ to a $DNF$ formula and vice versa.
  2. For any boolean formula $\phi$ over $n$ variables (regardless its form), there are exactly $2^{n}$ possible assignments each of which satisfy exactly one from $\left\{ \phi,\neg\phi\right\}$ .

  It is very easy to solve DNF-SAT efficiently, meaning, in polynomial time. Given a $DNF$ formula $\psi$ , iterate over its clauses until you find a satisfiable clause then declare that $\psi$ is satisfiable, however, if you failed to find such a clause then declare that $\psi$ is not satisfiable. Observe that any clause of the form $C_{i}=\bigwedge_{j=1}^{\left|C_{i}\right|}\ell_{i,j}$ is satisfiable if and only if it does not contains a variable and its negation; which we can check with a single pass over the literals of $C_{i}$ . As a result, $DNFSAT\in\mathcal{P}$ .

  In contrast to DNF-SAT, the CNF-SAT is know to be $\mathcal{NP}$ -complete, so we believe that it cannot be solved in polynomial time. Similar to the claim from the previous section, we conclude that CNF-SAT is “at least harder” as DNF-SAT, that is, $DNFSAT\le CNFSAT$ . However, as regard to the Counting versions of those problems, we can show the opposite inequality; that $\#CNFSAT\le\#DNFSAT$ .

  To prove that, we need to show that if we could solve $\#DNFSAT$ efficiently, then we could solve $\#CNFSAT$ efficiently. So suppose we have a magical algorithm that can solve $\#DNFSAT$ in polynomial time. Given a $CNF$ formula $\phi$ over $n$ boolean variables, construct $\psi=\neg\phi$ and feed it to the magical algorithm. Let its output be the number $t$ , then return $2^{n}-t$ .

  The correctness of that reduction follows from the two facts we stated earlier. Note that also all the procedures can be executed in polynomial time, so we solved $\#CNFSAT$ in polynomial time using an efficient solver for $\#DNFSAT$ , therefore, $\#CNFSAT\le\#DNFSAT$ .

  Actually, if we will look closely on that reduction, we will see that it should work perfectly for the other direction as well, that is, solving $\#DNFSAT$ by a an efficient algorithm for the $\#CNFSAT$ , so we may conclude that both problems are equally hard. This is an evidence that the hierarchy between decision problems do not preserved for their Counting versions. Putting it differently, we can say that Counting makes problems much harder, such that even an “easy” problems from $\mathcal{P}$ become harder as the entire $\mathcal{NP}$ when you consider its Counting version.

  3 $\#DNF$ is Also $\#\mathcal{P}$ -Complete

  Since $\#DNF$ and $\#CNF$ are equally hard, it is enough to show that $\#CNF$ is $\#\mathcal{P}$ -complete; for that we need a reduction from any problem in $\#\mathcal{P}$ to $\#CNF$ . Recall the Cook-Levin reduction which proves that CNF-SAT is $\mathcal{NP}$ -complete. That reduction can be applied to any $P\in\mathcal{NP}$ and produce a polynomial procedure that converts any input $x$ for the original problem $P$ , to a $CNF$ formula $\phi$ for the problem CNF-SAT, such that $\phi$ is satisfiable if and only if the decision problem $P$ says $Yes$ for the original input $x$ . That way, any efficient solver for CNF-SAT solves $P$ efficiently, so $P\le CNFSAT$ .

  An important property of the Cook-Levin reduction is that it keeps the number of s.a. for $\phi$ exactly as the number of witnesses for the original input $x$ of $P$ . That means, that we can use the same reduction to convert any Counting problem $\#P\in\#\mathcal{P}$ to $\#CNF$ . That is, the reduction will yields a procedure to convert any input $x$ for $\#P$ to a formula $\phi$ for $\#CNF$ such that the number of s.a. of $\phi$ is the solution for $\#P\left(x\right)$ . As a result, $\#CNF$ is $\#\mathcal{P}$ -complete, and then also $\#DNF$ .

  Following the conclusion from previous section, we can say further that even a simple problem from $\mathcal{P}$ , can have a Counting version that is not just “harder” but harder as any other problem in $\#\mathcal{P}$ .

  4 Is There a Structure inside $\#\mathcal{P}$ ?

  Going back to the $\#CNF$ and $\#DNF$ problems, we saw that they are related via the relation $$\#CNF\left(\phi\right)+\#DNF\left(\neg\phi\right)=2^{n}$$ for any $CNF$ formula $\phi$ over $n$ boolean variables. We exploited that relation to show that the problems are equally “hard” when talking about exact-solutions. However, an hierarchy emerged once we moved to consider approximated solutions for those problems. Therefore, there is some kind of structure for the family $\#\mathcal{P}$ under the appropriate notion of “solving” the problems.

  Next post we will discuss an approximation algorithm for $\#DNF$ and the hardness of approximating $\#CNF$ , showing that the additive relation from above it not enough for constructing an approximation algorithm for $\#CNF$ based on approximation for $\#DNF$ .
• Jul 2, 2017

  Approximation Algorithm for the Steiner Tree Problem

  tags:
  TheoryLunch Blog

  Computer science topics over lunch

  Approximation Algorithm for the Steiner Tree Problem

  Jul 2, 2017 • Eran Amar

  tags: $\newcommand{\lyxlock}{}$

  1 The Steiner Tree Problem

  An instance of the Steiner Tree problem is a connected graph $G$ on a set of $n$ vertices $V\left(G\right)$ and $m$ edges $E$ $\left(G\right)$ equipped with a positive weight function $w:E\left(G\right)\to\mathbb{R}_{>0}$ . Given a set of terminals $T\subseteq V\left(G\right)$ , a valid solution is s set $S\subseteq E\left(G\right)$ that induces a sub-graph in which any pair of terminals $t_{1},t_{2}\in T$ is connected. The weight of a solution $S$ is defined by $w\left(S\right)=\sum_{e\in S}w\left(e\right)$ , and we want to find a valid solution with minimum weight, usually denoted $S^{*}$ . Note that since the initial graph is connected there is always a valid solution.

  Finding an exact optimal solution is NP-Complete. However, there is polynomial time approximation algorithm that yields a factor 2 approximation. We say that a solution $S$ is 2-approximation for the problem if $w\left(S\right)\le2\cdot w\left(S^{*}\right)$ for some optimal solution $S^{*}$ .

  2 Adding Assumptions “for Free”

  Given a graph $G$ and weight-function $w$ , consider the following construction: Let $\tilde{G}$ be a clique on the same set of vertices $V\left(\tilde{G}\right)=V\left(G\right)$ , that is, $E\left(\tilde{G}\right)=\left\{ \left(u,v\right)\mid u\ne v\in V\left(G\right)\right\}$ . Define the weight function $\tilde{w}$ as follows, $\forall u,v\in V\left(\tilde{G}\right)$ let $\tilde{w}\left(u,v\right)$ be the weight of the shortest path from $u$ to $v$ in the original $G$ . When “shortest path” is defined to be the path with minimum sum of weights among all the simple paths that from $u$ to $v$ . Furthermore, it worth noting that by its definition, $\tilde{w}$ obeys the triangle inequality, meaning that for any triple $u,v,w\in V\left(\tilde{G}\right)$ , if $e_{1}=\left(u,v\right)$ , $e_{2}=\left(v,w\right)$ and $e_{3}=\left(u,w\right)$ then $\tilde{w}\left(e_{3}\right)\le\tilde{w}\left(e_{1}\right)+\tilde{w}\left(e_{2}\right)$ .

  Fix some $T\subseteq V\left(G\right)$ , we will show that the weight of any optimal solution for the original instance $\left(G,w\right)$ is exactly the same weight as of any optimal solution for $\left(\tilde{G},\tilde{w}\right)$ . We will do that by showing two inequalities: let $S^{*},\tilde{S}^{*}$ be optimal solutions for the original and the modified instances respectively. Then $w\left(S^{*}\right)\le w\left(\tilde{S}^{*}\right)$ and $w\left(\tilde{S}^{*}\right)\le w\left(S^{*}\right)$ .

  Obviously, the inequality $w\left(\tilde{S}^{*}\right)\le w\left(S^{*}\right)$ holds, because $\tilde{G}$ contains all the edges of $G$ and for any edge $e$ in $G$ its weight over the modified instance $\tilde{G}$ have $\tilde{w}\left(e\right)\le w\left(e\right)$ , by the definition of $\tilde{w}$ . Therefore, the optimal weight over $\tilde{G}$ can only be equal to or smaller than the one of $G$ .

  The other direction is less trivial. Take any $\tilde{S}^{*}$ and convert it to be a valid solution over the original $G$ as follows: for any $e\in\tilde{S}^{*}$ , if $e\in E\left(G\right)$ and $w\left(e\right)\le\tilde{w}\left(e\right)$ then add $e\in S^{*}$ . Otherwise, take all the edges along the shortest path connecting the endpoints of $e$ other the original graph $G$ and add them to $S^{*}$ . Note that the entire weight of that path cannot exceeds $\tilde{w}\left(e\right)$ by the definition of $\tilde{w}$ . Therefore we end up with a valid solution $S^{*}$ , that also have $w\left(S^{*}\right)\le w\left(\tilde{S}^{*}\right)$ . In conclusion, $w\left(S^{*}\right)=w\left(\tilde{S}^{*}\right)$ .

  The proof above gave us a concrete polynomial-time procedure to convert any valid solution $\tilde{S}$ over $\tilde{G}$ to a valid solution $S$ over $G$ such that $w\left(S\right)\le w\left(\tilde{S}\right)$ . From that we can conclude that it is suffice to find a 2-approximation for the “easier” instance $\left(\tilde{G},\tilde{w}\right)$ and then converts it to a 2-approximation solution for the original instance $\left(G,w\right)$ , note also that the construction of $\tilde{G},\tilde{w}$ from $G,w$ can be done in polynomial time (for example use Dijkstra's algorithm for computing the shortest paths).

  In other words, we can assume for free (i.e. without losing anything) that $G$ is a clique and $w$ have the triangle inequality, and design an approximation algorithm that deals only with such instances.

  3 Approximation Algorithm for the Easier Problem

  Let $G$ be a clique, $w$ be a weight function that obeys the triangle inequality, and $T$ is a set of terminals. Use any known polynomial time algorithm to get a Minimum Spanning Tree (MST) on $T$ (for instance, use Kruskal's algorithm), let the edges of that tree be $S\subseteq E\left(G\right)$ . Output $S$ as the approximated solution for the Steiner Tree problem.

  Putting it differently, we claim that $w\left(S\right)\le2w\left(S^{*}\right)$ , when $S^{*}$ is some optimal solution for the Steiner Tree problem for that “easier” instance.

  And now we will prove that.

  Let $S^{*}$ be some optimal solution for the Steiner Tree problem. It must be a tree because otherwise we have a redundant edge such that removing it will decrease the weight of the solution and keep it valid (the terminals will stay connected); in contradiction to the optimality of $S^{*}$ . Run a Depth First Search (DFS) at the root of $S^{*}$ , traverse all the nodes in the tree and return back to the root. Memorize all the nodes traversed that way (with repetitions) and denote that sequence $P$ . By the way we build $P$ , it is a cycle that touches all the vertices of $T$ (and maybe also others) and uses twice any edge of $S^{*}$ , once for the forward pass and once for the backward pass, so it has a length $k:=2\left|S^{*}\right|$ . Suppose that $P=\left(u_{1},...,u_{k}\right)$ , then we conclude that $$w\left(P\right)=\sum_{i=1}^{k-1}w\left(\left(u_{i},u_{i+1}\right)\right)=2\cdot w\left(S^{*}\right)$$

  Define the following order on $T$ : for any two terminals $t_{1},t_{2}\in T$ , we say that $t_{1}$ comes before $t_{2}$ if the first appearance of $t_{1}$ in $P$ comes before the first appearance of $t_{2}$ in $P$ . Let $\left(t_{1},...,t_{\left|T\right|}\right)$ be the sequence of ordered terminals, note that it is also defines a path because $G$ is a clique, that is, $\left(t_{i},t_{i+1}\right)\in E\left(G\right)$ for any $i$ . Denote that path $\tilde{P}$ .

  By the triangle inequality, for any pair $t_{i},t_{i+1}\in T$ the weighed subpath in $\tilde{P}$ , which is the weight of a single edge, is less or equal to the weighted length in $P$ (which may use more than one edge), thus $w\left(\tilde{P}\right)\le w\left(P\right)$ .

  Since $\tilde{P}$ touches all the terminals in $T$ exactly once, it has no cycles, hence it is also a spanning tree over the vertices $T$ and edges of the clique. Recall that the proposed solution $S$ is a minimum spanning tree on $T$ , therefore $w\left(S\right)\le w\left(\tilde{P}\right)$ . We conclude that $$w\left(S\right)\le w\left(\tilde{P}\right)\le w\left(P\right)=2\cdot w\left(S^{*}\right)$$ as required.

  $\newcommand{\lyxlock}{}$

  1 The Steiner Tree Problem

  An instance of the Steiner Tree problem is a connected graph $G$ on a set of $n$ vertices $V\left(G\right)$ and $m$ edges $E$ $\left(G\right)$ equipped with a positive weight function $w:E\left(G\right)\to\mathbb{R}_{>0}$ . Given a set of terminals $T\subseteq V\left(G\right)$ , a valid solution is s set $S\subseteq E\left(G\right)$ that induces a sub-graph in which any pair of terminals $t_{1},t_{2}\in T$ is connected. The weight of a solution $S$ is defined by $w\left(S\right)=\sum_{e\in S}w\left(e\right)$ , and we want to find a valid solution with minimum weight, usually denoted $S^{*}$ . Note that since the initial graph is connected there is always a valid solution.

  Finding an exact optimal solution is NP-Complete. However, there is polynomial time approximation algorithm that yields a factor 2 approximation. We say that a solution $S$ is 2-approximation for the problem if $w\left(S\right)\le2\cdot w\left(S^{*}\right)$ for some optimal solution $S^{*}$ .

  2 Adding Assumptions “for Free”

  Given a graph $G$ and weight-function $w$ , consider the following construction: Let $\tilde{G}$ be a clique on the same set of vertices $V\left(\tilde{G}\right)=V\left(G\right)$ , that is, $E\left(\tilde{G}\right)=\left\{ \left(u,v\right)\mid u\ne v\in V\left(G\right)\right\}$ . Define the weight function $\tilde{w}$ as follows, $\forall u,v\in V\left(\tilde{G}\right)$ let $\tilde{w}\left(u,v\right)$ be the weight of the shortest path from $u$ to $v$ in the original $G$ . When “shortest path” is defined to be the path with minimum sum of weights among all the simple paths that from $u$ to $v$ . Furthermore, it worth noting that by its definition, $\tilde{w}$ obeys the triangle inequality, meaning that for any triple $u,v,w\in V\left(\tilde{G}\right)$ , if $e_{1}=\left(u,v\right)$ , $e_{2}=\left(v,w\right)$ and $e_{3}=\left(u,w\right)$ then $\tilde{w}\left(e_{3}\right)\le\tilde{w}\left(e_{1}\right)+\tilde{w}\left(e_{2}\right)$ .

  Fix some $T\subseteq V\left(G\right)$ , we will show that the weight of any optimal solution for the original instance $\left(G,w\right)$ is exactly the same weight as of any optimal solution for $\left(\tilde{G},\tilde{w}\right)$ . We will do that by showing two inequalities: let $S^{*},\tilde{S}^{*}$ be optimal solutions for the original and the modified instances respectively. Then $w\left(S^{*}\right)\le w\left(\tilde{S}^{*}\right)$ and $w\left(\tilde{S}^{*}\right)\le w\left(S^{*}\right)$ .

  Obviously, the inequality $w\left(\tilde{S}^{*}\right)\le w\left(S^{*}\right)$ holds, because $\tilde{G}$ contains all the edges of $G$ and for any edge $e$ in $G$ its weight over the modified instance $\tilde{G}$ have $\tilde{w}\left(e\right)\le w\left(e\right)$ , by the definition of $\tilde{w}$ . Therefore, the optimal weight over $\tilde{G}$ can only be equal to or smaller than the one of $G$ .

  The other direction is less trivial. Take any $\tilde{S}^{*}$ and convert it to be a valid solution over the original $G$ as follows: for any $e\in\tilde{S}^{*}$ , if $e\in E\left(G\right)$ and $w\left(e\right)\le\tilde{w}\left(e\right)$ then add $e\in S^{*}$ . Otherwise, take all the edges along the shortest path connecting the endpoints of $e$ other the original graph $G$ and add them to $S^{*}$ . Note that the entire weight of that path cannot exceeds $\tilde{w}\left(e\right)$ by the definition of $\tilde{w}$ . Therefore we end up with a valid solution $S^{*}$ , that also have $w\left(S^{*}\right)\le w\left(\tilde{S}^{*}\right)$ . In conclusion, $w\left(S^{*}\right)=w\left(\tilde{S}^{*}\right)$ .

  The proof above gave us a concrete polynomial-time procedure to convert any valid solution $\tilde{S}$ over $\tilde{G}$ to a valid solution $S$ over $G$ such that $w\left(S\right)\le w\left(\tilde{S}\right)$ . From that we can conclude that it is suffice to find a 2-approximation for the “easier” instance $\left(\tilde{G},\tilde{w}\right)$ and then converts it to a 2-approximation solution for the original instance $\left(G,w\right)$ , note also that the construction of $\tilde{G},\tilde{w}$ from $G,w$ can be done in polynomial time (for example use Dijkstra's algorithm for computing the shortest paths).

  In other words, we can assume for free (i.e. without losing anything) that $G$ is a clique and $w$ have the triangle inequality, and design an approximation algorithm that deals only with such instances.

  3 Approximation Algorithm for the Easier Problem

  Let $G$ be a clique, $w$ be a weight function that obeys the triangle inequality, and $T$ is a set of terminals. Use any known polynomial time algorithm to get a Minimum Spanning Tree (MST) on $T$ (for instance, use Kruskal's algorithm), let the edges of that tree be $S\subseteq E\left(G\right)$ . Output $S$ as the approximated solution for the Steiner Tree problem.

  Putting it differently, we claim that $w\left(S\right)\le2w\left(S^{*}\right)$ , when $S^{*}$ is some optimal solution for the Steiner Tree problem for that “easier” instance.

  And now we will prove that.

  Let $S^{*}$ be some optimal solution for the Steiner Tree problem. It must be a tree because otherwise we have a redundant edge such that removing it will decrease the weight of the solution and keep it valid (the terminals will stay connected); in contradiction to the optimality of $S^{*}$ . Run a Depth First Search (DFS) at the root of $S^{*}$ , traverse all the nodes in the tree and return back to the root. Memorize all the nodes traversed that way (with repetitions) and denote that sequence $P$ . By the way we build $P$ , it is a cycle that touches all the vertices of $T$ (and maybe also others) and uses twice any edge of $S^{*}$ , once for the forward pass and once for the backward pass, so it has a length $k:=2\left|S^{*}\right|$ . Suppose that $P=\left(u_{1},...,u_{k}\right)$ , then we conclude that $$w\left(P\right)=\sum_{i=1}^{k-1}w\left(\left(u_{i},u_{i+1}\right)\right)=2\cdot w\left(S^{*}\right)$$

  Define the following order on $T$ : for any two terminals $t_{1},t_{2}\in T$ , we say that $t_{1}$ comes before $t_{2}$ if the first appearance of $t_{1}$ in $P$ comes before the first appearance of $t_{2}$ in $P$ . Let $\left(t_{1},...,t_{\left|T\right|}\right)$ be the sequence of ordered terminals, note that it is also defines a path because $G$ is a clique, that is, $\left(t_{i},t_{i+1}\right)\in E\left(G\right)$ for any $i$ . Denote that path $\tilde{P}$ .

  By the triangle inequality, for any pair $t_{i},t_{i+1}\in T$ the weighed subpath in $\tilde{P}$ , which is the weight of a single edge, is less or equal to the weighted length in $P$ (which may use more than one edge), thus $w\left(\tilde{P}\right)\le w\left(P\right)$ .

  Since $\tilde{P}$ touches all the terminals in $T$ exactly once, it has no cycles, hence it is also a spanning tree over the vertices $T$ and edges of the clique. Recall that the proposed solution $S$ is a minimum spanning tree on $T$ , therefore $w\left(S\right)\le w\left(\tilde{P}\right)$ . We conclude that $$w\left(S\right)\le w\left(\tilde{P}\right)\le w\left(P\right)=2\cdot w\left(S^{*}\right)$$ as required.
• Jun 15, 2017

  Probabilistic Approximation Algorithm for Set Cover

  tags:
  TheoryLunch Blog

  Computer science topics over lunch

  Probabilistic Approximation Algorithm for Set Cover

  Jun 15, 2017 • Eran Amar

  tags: $\newcommand{\lyxlock}{}$

  1 The Weighted Set Cover Problem

  In the Weighted Set Cover (WSC) problem, the input is a triple $\left(U,\mathcal{F},w\right)$ where $U$ is a universe of $n$ elements, $\mathcal{F}=\left\{ A_{1},A_{2},..,A_{m}\right\}$ is a family of $m$ sets such that each $A_{i}\subseteq U$ and $\bigcup_{A\in\mathcal{F}}A=U$ , and $w:\mathcal{F}\to\mathbb{R}_{>0}$ is a weight function that assigns positive weight to each set in $\mathcal{F}$ . A solution $S\subseteq\mathcal{F}$ for the WSC problem is valid if and only if $\bigcup_{A\in S}A=U$ , the weight of the solution is $w\left(S\right)=\sum_{A\in S}w\left(A\right)$ . Note that there is always a valid solution, and the goal in the WSC problem is to find such one with minimum weight.

  Weighted Set Cover is “hard” in many aspects. First, it is hard to efficiently find the exact solution, that is, the problem is NP-Complete so no time-efficient algorithm is known for it. Moreover, it is also hard to approximate the solution better than a $\log n$ factor. Formally, we said that a solution $S$ is a $\log n$ approximation for WSC if $w\left(S\right)\le w\left(S^{*}\right)\cdot\log n$ , when $S^{*}$ is an optimal solution with minimum weight. The WSC is hard to approximate in the sense that any $\left(1-\epsilon\right)\cdot\log n$ approximation for that problem is NP-Complete by itself, and that is for any $\epsilon>0$ .

  In this post we will describe a probabilistic $c\cdot\log n$ approximation algorithm for WSC which runs in $\mathcal{O}\left(poly\left(n,m\right)\right)$ time, and $c$ is some positive constant that controls the success probability of the algorithm.

  2 Equivalent Formulation

  We can reformulate the Weighted Set Cover problem as a system of linear equations as follows: For any set $A_{i}$ define an indicator variable $X_{i}\in\left\{ 0,1\right\}$ and let $w\left(X_{i}\right):=w\left(A_{i}\right)\cdot X_{i}$ , thus if $X_{i}=1$ then $w\left(X_{i}\right)=w\left(A_{i}\right)$ , and $w\left(X_{i}\right)=0$ otherwise. The objective is then to minimize $$\sum_{i=1}^{m}w\left(X_{i}\right)$$ subject to the $n$ following linear constrains $$\forall u\in U:\quad\sum_{i\text{ s.t. }u\in A_{i}}X_{i}\ge1$$

  Given a solution $l=\left(x_{1},..,x_{m}\right)\in\left\{ 0,1\right\} ^{m}$ to that system of equations, we can define the corresponding solution for WSC as $S=\left\{ A_{i}\in\mathcal{F}\mid x_{i}=1\right\}$ . Note that by this construction $S$ must be valid, moreover, the weights are the same: $w\left(l\right)=w\left(S\right)$ , so if $w\left(l\right)$ is minimal then also $w\left(S\right)$ . The proof of the other direction (from WSC to the system of equations) is similar.

  Because the linear-time mapping from any optimal solution of one problem to an optimal solution for the other, we can conclude that solving this system of equations is NP-Hard. However, if we relax the “integer requirement”, i.e. that $X_{i}\in\left\{ 0,1\right\}$ , then the problem becomes much easier. That is, when $X_{i}\in\left[0,1\right]$ for any $i$ , that system becomes a convex optimization problem, which we can solved with a polynomial-time algorithm. From now on, we will refer only to the convex version, and denote its (exact) solution as fractional solution.

  So suppose we have a fractional optimal solution for the linear system $l^{*}=\left(x_{1},..,x_{m}\right)\in\left[0,1\right]^{m}$ . Note that $w\left(l^{*}\right)\le w\left(S^{*}\right)$ for any optimal solution $S^{*}$ for WSC, because the convex formulation allows for more possible valid solutions so the minimum weight can only be smaller than of $w\left(S^{*}\right)$ .

  How do we get from $l^{*}$ a solution for the WSC? The nice trick here is probabilistic rounding. That is, we will treat each $x_{i}$ as the probability for selecting $A_{i}$ .

  We define the following procedure $Alg\left(U,\mathcal{F},w\right)$ as an helper function:

  1. Repeat for $t=1,...,\left(2\log n\right)$ times:
     1. Initialize $S_{t}=\emptyset$ .
     2. Solve the convex optimization version of the WSC to get an optimal fractional solution $l_{t}^{*}$ .
     3. For each $A_{i}\in\mathcal{F}$ , add it to $S_{t}$ with probability $x_{i}\in l_{t}^{*}$ .
     4. Let the “solution” for the current iteration be $S_{t}$ .
  2. Output $\bigcup_{t=1}^{2\log n}S_{t}$

  The final algorithm is just applying the $Alg\left(U,\mathcal{F},w\right)$ procedure $c$ independent times, and keeping the candidate solutions $S^{\left(1\right)},...,S^{\left(c\right)}$ , then returning the $S_{final}$ solution with minimum weight among those candidates.

  Note that all the procedures above have polynomial time (in $n$ and $m$ ), therefore the total running time of the algorithm is $\mathcal{O}\left(poly\left(n,m\right)\right)$ .

  3 The Solution is Approximately Optimal

  We want to show that with high probability $w\left(S_{final}\right)\le c\cdot\log n\cdot w\left(S^{*}\right)$ . It is easy to see that for each iteration of the sub-procedure $Alg$ , the intermediate set $S_{t}$ has expected weight equals to the weight of the optimal fractional solution. Formally, $$\mathbb{E}\left[w\left(S_{t}\right)\right]=\sum_{i=1}^{m}w\left(A_{i}\right)\mathbb{E}\left[X_{i}\right]=\sum_{i=1}^{m}w\left(A_{i}\right)x_{i}=w\left(l_{t}^{*}\right)$$ note that the optimal value of all the fractional solutions is the same so we will remove the subscript $t$ . Next, we continue to calculate the expected weight of the output solution of the sub-procedure, that is, the weight of the candidate $S^{\left(i\right)}$

  $$\begin{aligned} \mathbb{E}\left[w\left(S^{\left(i\right)}\right)\right] & =\mathbb{E}\left[w\left(\bigcup_{t=1}^{2\log n}S_{t}\right)\right]\\ & \le\sum_{t=1}^{2\log n}\mathbb{E}\left[w\left(S_{t}\right)\right]\\ & =2\log n\cdot w\left(l^{*}\right)\\ & \le2\log n\cdot w\left(S^{*}\right) \end{aligned}$$ It left to show that with high probability all the candidates $\left\{ S^{\left(i\right)}\right\} _{i}$ are within the weight bound of $c\cdot\log n\cdot w\left(S^{*}\right)$ , and because $S_{final}$ is being selected from those candidates, we conclude that it also has $\mathbb{E}\left[w\left(S_{final}\right)\right]\le c\cdot\log n\cdot w\left(S^{*}\right)$ . In order to get high probability for a single candidate, we can use Markov’s inequality as follows. Define a “failure” of a candidate $S^{\left(k\right)}$ as $w\left(S^{\left(k\right)}\right)\ge c\cdot\log n\cdot w\left(S^{*}\right)$ , then the probability of a single candidate to fail is bounded by $$\mathbb{P}\left[w\left(S^{\left(k\right)}\right)\ge c\cdot\log n\cdot w\left(S^{*}\right)\right]\le\frac{\mathbb{E}\left[w\left(S^{\left(k\right)}\right)\right]}{c\cdot\log n\cdot w\left(S^{*}\right)}\le\frac{2}{c}$$ So with probability of at least $1-\frac{2}{c}$ the candidate has $w\left(S^{\left(k\right)}\right)\le c\cdot\log n\cdot w\left(S^{*}\right)$ as required.

  When considering the final solution $S_{final}$ , it fails only when all the candidates fail and the probability for that is $\le\left(\frac{2}{c}\right)^{c}$ . So we have high constant probability to get a solution that is $c\cdot\log n$ approximation for the optimal solution.

  4 The Solution is Probably Valid

  What is the probability that the returned solution $S_{final}$ is not valid? that is, there is at least one element $u\in U$ that is not covered by $S_{final}$ . We analyze that by first considering an arbitrary element $u$ and an intermediate solution $S_{t}$ of the sub procedure $Alg$ . Suppose that $u$ belongs to $k$ sets, and for simplicity suppose they are $A_{1},...,A_{k}$ . We know that their corresponding variables have $\sum_{i=1}^{k}x_{i}\ge1$ by the requirements of the linear system. The probability to leave all these sets outside $S_{t}$ is $$\prod_{i=1}^{k}\left(1-x_{i}\right)\le\prod_{i=1}^{k}e^{-x_{i}}\le e^{-1}$$ that is, the probability that an intermediate solution $S_{t}$ will not cover $u$ is up to $e^{-1}$ .

  Next, we calculate the probability that $u$ is not covered by a $S_{final}$ where $S_{final}=S^{\left(i\right)}$ for some selected candidate $S^{\left(i\right)}$ . Then it must be that $u$ is not covered by any intermediate solution $S_{1},..,S_{2\log n}$ of the procedure $Alg$ , and that will happen with probability $\le\prod_{t=1}^{2\log n}e^{-1}=\frac{1}{n^{2}}$ . Applying the union-bound on all the $n$ elements, we conclude that with probability of at most $\frac{1}{n}$ , the selected candidate $S^{\left(i\right)}$ will not cover the entire universe. In other words, with at least $1-\frac{1}{n}$ probability, $S_{final}$ is a valid solution.

  $\newcommand{\lyxlock}{}$

  1 The Weighted Set Cover Problem

  In the Weighted Set Cover (WSC) problem, the input is a triple $\left(U,\mathcal{F},w\right)$ where $U$ is a universe of $n$ elements, $\mathcal{F}=\left\{ A_{1},A_{2},..,A_{m}\right\}$ is a family of $m$ sets such that each $A_{i}\subseteq U$ and $\bigcup_{A\in\mathcal{F}}A=U$ , and $w:\mathcal{F}\to\mathbb{R}_{>0}$ is a weight function that assigns positive weight to each set in $\mathcal{F}$ . A solution $S\subseteq\mathcal{F}$ for the WSC problem is valid if and only if $\bigcup_{A\in S}A=U$ , the weight of the solution is $w\left(S\right)=\sum_{A\in S}w\left(A\right)$ . Note that there is always a valid solution, and the goal in the WSC problem is to find such one with minimum weight.

  Weighted Set Cover is “hard” in many aspects. First, it is hard to efficiently find the exact solution, that is, the problem is NP-Complete so no time-efficient algorithm is known for it. Moreover, it is also hard to approximate the solution better than a $\log n$ factor. Formally, we said that a solution $S$ is a $\log n$ approximation for WSC if $w\left(S\right)\le w\left(S^{*}\right)\cdot\log n$ , when $S^{*}$ is an optimal solution with minimum weight. The WSC is hard to approximate in the sense that any $\left(1-\epsilon\right)\cdot\log n$ approximation for that problem is NP-Complete by itself, and that is for any $\epsilon>0$ .

  In this post we will describe a probabilistic $c\cdot\log n$ approximation algorithm for WSC which runs in $\mathcal{O}\left(poly\left(n,m\right)\right)$ time, and $c$ is some positive constant that controls the success probability of the algorithm.

  2 Equivalent Formulation

  We can reformulate the Weighted Set Cover problem as a system of linear equations as follows: For any set $A_{i}$ define an indicator variable $X_{i}\in\left\{ 0,1\right\}$ and let $w\left(X_{i}\right):=w\left(A_{i}\right)\cdot X_{i}$ , thus if $X_{i}=1$ then $w\left(X_{i}\right)=w\left(A_{i}\right)$ , and $w\left(X_{i}\right)=0$ otherwise. The objective is then to minimize $$\sum_{i=1}^{m}w\left(X_{i}\right)$$ subject to the $n$ following linear constrains $$\forall u\in U:\quad\sum_{i\text{ s.t. }u\in A_{i}}X_{i}\ge1$$

  Given a solution $l=\left(x_{1},..,x_{m}\right)\in\left\{ 0,1\right\} ^{m}$ to that system of equations, we can define the corresponding solution for WSC as $S=\left\{ A_{i}\in\mathcal{F}\mid x_{i}=1\right\}$ . Note that by this construction $S$ must be valid, moreover, the weights are the same: $w\left(l\right)=w\left(S\right)$ , so if $w\left(l\right)$ is minimal then also $w\left(S\right)$ . The proof of the other direction (from WSC to the system of equations) is similar.

  Because the linear-time mapping from any optimal solution of one problem to an optimal solution for the other, we can conclude that solving this system of equations is NP-Hard. However, if we relax the “integer requirement”, i.e. that $X_{i}\in\left\{ 0,1\right\}$ , then the problem becomes much easier. That is, when $X_{i}\in\left[0,1\right]$ for any $i$ , that system becomes a convex optimization problem, which we can solved with a polynomial-time algorithm. From now on, we will refer only to the convex version, and denote its (exact) solution as fractional solution.

  So suppose we have a fractional optimal solution for the linear system $l^{*}=\left(x_{1},..,x_{m}\right)\in\left[0,1\right]^{m}$ . Note that $w\left(l^{*}\right)\le w\left(S^{*}\right)$ for any optimal solution $S^{*}$ for WSC, because the convex formulation allows for more possible valid solutions so the minimum weight can only be smaller than of $w\left(S^{*}\right)$ .

  How do we get from $l^{*}$ a solution for the WSC? The nice trick here is probabilistic rounding. That is, we will treat each $x_{i}$ as the probability for selecting $A_{i}$ .

  We define the following procedure $Alg\left(U,\mathcal{F},w\right)$ as an helper function:

  1. Repeat for $t=1,...,\left(2\log n\right)$ times:
     1. Initialize $S_{t}=\emptyset$ .
     2. Solve the convex optimization version of the WSC to get an optimal fractional solution $l_{t}^{*}$ .
     3. For each $A_{i}\in\mathcal{F}$ , add it to $S_{t}$ with probability $x_{i}\in l_{t}^{*}$ .
     4. Let the “solution” for the current iteration be $S_{t}$ .
  2. Output $\bigcup_{t=1}^{2\log n}S_{t}$

  The final algorithm is just applying the $Alg\left(U,\mathcal{F},w\right)$ procedure $c$ independent times, and keeping the candidate solutions $S^{\left(1\right)},...,S^{\left(c\right)}$ , then returning the $S_{final}$ solution with minimum weight among those candidates.

  Note that all the procedures above have polynomial time (in $n$ and $m$ ), therefore the total running time of the algorithm is $\mathcal{O}\left(poly\left(n,m\right)\right)$ .

  3 The Solution is Approximately Optimal

  We want to show that with high probability $w\left(S_{final}\right)\le c\cdot\log n\cdot w\left(S^{*}\right)$ . It is easy to see that for each iteration of the sub-procedure $Alg$ , the intermediate set $S_{t}$ has expected weight equals to the weight of the optimal fractional solution. Formally, $$\mathbb{E}\left[w\left(S_{t}\right)\right]=\sum_{i=1}^{m}w\left(A_{i}\right)\mathbb{E}\left[X_{i}\right]=\sum_{i=1}^{m}w\left(A_{i}\right)x_{i}=w\left(l_{t}^{*}\right)$$ note that the optimal value of all the fractional solutions is the same so we will remove the subscript $t$ . Next, we continue to calculate the expected weight of the output solution of the sub-procedure, that is, the weight of the candidate $S^{\left(i\right)}$

  $$\begin{aligned} \mathbb{E}\left[w\left(S^{\left(i\right)}\right)\right] & =\mathbb{E}\left[w\left(\bigcup_{t=1}^{2\log n}S_{t}\right)\right]\\ & \le\sum_{t=1}^{2\log n}\mathbb{E}\left[w\left(S_{t}\right)\right]\\ & =2\log n\cdot w\left(l^{*}\right)\\ & \le2\log n\cdot w\left(S^{*}\right) \end{aligned}$$ It left to show that with high probability all the candidates $\left\{ S^{\left(i\right)}\right\} _{i}$ are within the weight bound of $c\cdot\log n\cdot w\left(S^{*}\right)$ , and because $S_{final}$ is being selected from those candidates, we conclude that it also has $\mathbb{E}\left[w\left(S_{final}\right)\right]\le c\cdot\log n\cdot w\left(S^{*}\right)$ . In order to get high probability for a single candidate, we can use Markov’s inequality as follows. Define a “failure” of a candidate $S^{\left(k\right)}$ as $w\left(S^{\left(k\right)}\right)\ge c\cdot\log n\cdot w\left(S^{*}\right)$ , then the probability of a single candidate to fail is bounded by $$\mathbb{P}\left[w\left(S^{\left(k\right)}\right)\ge c\cdot\log n\cdot w\left(S^{*}\right)\right]\le\frac{\mathbb{E}\left[w\left(S^{\left(k\right)}\right)\right]}{c\cdot\log n\cdot w\left(S^{*}\right)}\le\frac{2}{c}$$ So with probability of at least $1-\frac{2}{c}$ the candidate has $w\left(S^{\left(k\right)}\right)\le c\cdot\log n\cdot w\left(S^{*}\right)$ as required.

  When considering the final solution $S_{final}$ , it fails only when all the candidates fail and the probability for that is $\le\left(\frac{2}{c}\right)^{c}$ . So we have high constant probability to get a solution that is $c\cdot\log n$ approximation for the optimal solution.

  4 The Solution is Probably Valid

  What is the probability that the returned solution $S_{final}$ is not valid? that is, there is at least one element $u\in U$ that is not covered by $S_{final}$ . We analyze that by first considering an arbitrary element $u$ and an intermediate solution $S_{t}$ of the sub procedure $Alg$ . Suppose that $u$ belongs to $k$ sets, and for simplicity suppose they are $A_{1},...,A_{k}$ . We know that their corresponding variables have $\sum_{i=1}^{k}x_{i}\ge1$ by the requirements of the linear system. The probability to leave all these sets outside $S_{t}$ is $$\prod_{i=1}^{k}\left(1-x_{i}\right)\le\prod_{i=1}^{k}e^{-x_{i}}\le e^{-1}$$ that is, the probability that an intermediate solution $S_{t}$ will not cover $u$ is up to $e^{-1}$ .

  Next, we calculate the probability that $u$ is not covered by a $S_{final}$ where $S_{final}=S^{\left(i\right)}$ for some selected candidate $S^{\left(i\right)}$ . Then it must be that $u$ is not covered by any intermediate solution $S_{1},..,S_{2\log n}$ of the procedure $Alg$ , and that will happen with probability $\le\prod_{t=1}^{2\log n}e^{-1}=\frac{1}{n^{2}}$ . Applying the union-bound on all the $n$ elements, we conclude that with probability of at most $\frac{1}{n}$ , the selected candidate $S^{\left(i\right)}$ will not cover the entire universe. In other words, with at least $1-\frac{1}{n}$ probability, $S_{final}$ is a valid solution.
• Jun 11, 2017

  Visualizing high-dimensional data with t-SNE

  tags:
  TheoryLunch Blog

  Computer science topics over lunch

  Visualizing high-dimensional data with t-SNE

  Jun 11, 2017 • Aviv Netanyahu

  tags: $\newcommand{\lyxlock}{}$

  Guest author: thanks to Aviv Netanyahu for writing this post!

  1 What is visualization and why do we need it?

  Before working on a dataset we would like to understand it. After working on it, we want to communicate our conclusions with others easily. When dealing with high-dimensional data, it’s hard to get a ’feel’ for the data by just looking at numbers, and we will advance pretty slowly. For example, common visualization methods are histograms and scatter plots, which capture only one or two parameters at a time, meaning we would have to look at many histograms (one for each feature). However, by visualizing our data we may instantly be able to recognize structure. In order to do this we would like to find a representation of our data in $2$ D or $3$ D.

  2 Our goal

  Find an embedding of high-dimensional data to 2D or 3D (for simplicity, will assume $2$ D from now on).
  Formally, we want to learn $$\{\mathbf{x}_{1},...,\mathbf{x}_{N}\}\to\{\mathbf{y}_{1},...,\mathbf{y}_{N}\}$$ ($\{\mathbf{x}_{1},...,\mathbf{x}_{N}\}\subseteq\mathbb{R}^{n}$ is our data, $\{\mathbf{y}_{1},...,\mathbf{y}_{N}\}\subseteq\mathbb{R}^{2}$ is the embedding) such that the structure of the original data is preserved after the embedding.

  3 So why not PCA?

  PCA projects high-dimensional data onto a lower dimension determined by the directions where the data varies most in. The result is a linear projection which focuses on preserving distance between dissimilar data points. This is problematic:

  1. Linear sometimes isn’t good enough.
  2. We claim structure is mostly determined by similar data points.

  4 SNE (Hinton&Roweis)

  This algorithm for high-dimensional data visualization focuses on preserving high-dimensional ’neighbors’ in low-dimension.

  Assumption - ’neighbors’ in high dimension have small Euclidean distance. SNE converts the Euclidean distances to similarities using Gaussian distribution: $$p_{j|i}=\frac{\exp\left(-||\mathbf{x}_{i}-\mathbf{x}_{j}||^{2}/2\sigma_{i}^{2}\right)}{\sum\limits _{k\neq i}\exp\left(-||\mathbf{x}_{i}-\mathbf{x}_{k}||^{2}/2\sigma_{i}^{2}\right)}$$ $p_{j|i}$ is the probability that $j$ is the neighbor of $i$ in high-dimension. $$q_{j|i}=\frac{\exp\left(-||\mathbf{y}_{i}-\mathbf{y}_{j}||^{2}\right)}{\sum\limits _{k\neq i}\exp\left(-||\mathbf{y}_{i}-\mathbf{y}_{k}||^{2}\right)}$$ $q_{j|i}$ is the probability that $j$ is the neighbor of $i$ in low-dimension.
  Note - choosing $\sigma_{i}$ for $p_{j|i}$ is done by determining the perplexity, i.e. the effective number of neighbors for each data point. Then $\forall i$ , $\sigma_{i}$ is chosen such that $\mathbf{x}_{i}$ has according number of neighbors in the multivariate Gaussian with mean $\mathbf{x}_{i}$ and covariance matrix $\sigma_{i}\cdot I$ .
  
  By mapping each $\mathbf{x}_{i}\in\mathbb{R}^{n}$ to $\mathbf{y}_{i}\in\mathbb{R}^{2}$ such that $\forall i,j,\,p_{j|i}$ is close as possible to $q_{j|i}$ , we can preserve similarities. The algorithm finds such a mapping by gradient descent w.r.t. the following loss (cost) function: $$C=\sum_{i}KL(P_{i}||Q_{i})=\sum_{i}\sum_{j}p_{j|i}\log\frac{p_{j|i}}{q_{j|i}}$$ where $P_{i}=(p_{1|i},...,p_{N|i}),\,Q_{i}=(q_{1|i},...,q_{N|i})$ and the initialization of low-dimensional data is random around the origin.

  Notice the loss function, Kullback–Leibler (KL) divergence, is a natural way to measure similarity between two distributions.
  Also notice that this loss puts more emphasis on preserving the local structure of the data. KL-div is asymmetric, so different kinds of errors have a different affect:

  • Type 1 error: $p_{j|i}$ is large (i.e. two points were originally close) and $q_{j|i}$ small (i.e. these points were mapped too far), we pay a large penalty.
  • Type 2 error: $p_{j|i}$ is small (i.e. two points were originally far away) and $q_{j|i}$ large (i.e. these points were mapped too close), we pay a small penalty.

  This is good, as structure is preserved better by focusing on keeping close points together. This is easy to see in the following example: think of “unrolling” this swiss roll (like unrolling a carpet). Close points will stay within the same distance, far points won’t. This is true for many manifolds (spaces that locally resemble Euclidean space near each point) besides the swiss roll.

  

  4.1 The crowding problem

  The main problem with SNE - mapped points tend to “crowd” together, like in this example: images of 5 MNIST digits mapped with SNE to 2D and then marked by labels

  

  This happens because when reducing the dimensionality we need more “space” in which to represent our data-points. You can think of a sphere flattened to a circle where all distances are preserved (as much as possible) - the radius of the circle will have to be larger than the radius of the sphere!

  On one hand this causes type $1$ error (close points mapped too far), which forces far mapped points to come closer together. On the other, causes type $2$ error (far points mapped close) which forces them to move away from each other. Type $2$ error makes us pay a small price, however if many points cause error $2$ , the price will be higher, forcing the points to indeed move. What we just described causes the mapped points to correct themselves to the center such that they end up very close to each other. Consider the MNIST example above. If we did not know the labels we couldn’t make much sense of the mapping.

  5 T-SNE (van der Maaten&Hinton)

  5.1 Solving the crowding problem

  Instead of mapped points correcting themselves to minimize the loss, we can do this by increasing $q_{j|i}$ in advance. Recall, the loss is: $\sum\limits _{i}\sum\limits _{j}p_{j|i}\log\frac{p_{j|i}}{q_{j|i}}$ , thus larger $q_{j|i}$ minimizes it. We do this by defining: $$q_{ij}=\frac{(1+||\mathbf{y}_{i}-\mathbf{y}_{j}||)^{-1}}{\sum\limits _{k\neq l}(1+||\mathbf{y}_{k}-\mathbf{y}_{l}||)^{-1}}$$ using Cauchy distribution instead of Gaussian distribution. This increases $q_{j|i}$ since Cauchy distribution has a heavier tail than Gaussian distribution:

  

  (Hence the name of the algorithm - Cauchy distribution is a T-student distribution with $1$ degree of freedom).

  Note also that we changed the probability from conditional to joint probability (this was done for technical optimization reasons). For the same reason, $p_{j|i}$ was converted to $p_{ij}=\frac{p_{j|i}+p_{i|j}}{2N}$ . We can still use Gaussian distribution for the original data since we are trying to compensate for error in the mapped space.

  Now the performance on MNIST is much better:

  

  It’s obvious that global structure is preserved and crowding is resolved. Notice also that when zooming in on the mapping of $1$ digits, the local structure is preserved as well - digits with same orientation were mapped close to each other.

  6 Barnes-Hut approximation (van der Maaten)

  Computing the gradient for each mapped point each iteration is very costly, especially when the dataset is very large. At each iteration, $\forall i$ , we compute: $$\frac{\partial C}{\partial\mathbf{y}_{i}}=4\sum_{j\neq i}\frac{(p_{ij}-q_{ij})}{1+||\mathbf{y}_{i}-\mathbf{y}_{j}||^{2}}\cdot(\mathbf{y}_{i}-\mathbf{y}_{j})$$ We can think of $\mathbf{y}_{i}-\mathbf{y}_{j}$ as a spring between these two points (it might be too tight/loose depending on if the points were mapped to far/close accordingly). The fraction can be thought of as exertion/compression of the spring, in other words the correction we do in this iteration. We calculate this for $N$ mapped points, and derive $C$ by all data-points, so for each iteration we have $\mathcal{O}(N^{2})$ computations. This limits us to datasets with only few thousands of points! However, we can reduce the computation time to $\mathcal{O}(N\log N)$ by reducing the number of computations per iteration. The main idea - forces exerted by a group of points on a point that is relatively far away are all very similar. Therefore we can calculate $\mathbf{y}_{i}-\mathbf{y}_{j^{*}}$ where $\mathbf{y}_{j^{*}}$ is the average of the group of far away points, instead of $\mathbf{y}_{i}-\mathbf{y}_{j}$ for each $\mathbf{y}_{j}$ in the far away group.

  7 Some applications

  • Visualize word embeddings - http://nlp.yvespeirsman.be/blog/visualizing-word-embeddings-with-tsne/
  • Multiple maps word embeddings - http://homepage.tudelft.nl/19j49/multiplemaps/Multiple_maps_t-SNE/Multiple_maps_t-SNE.html
  • Visualization of genetic disease phenotype similarities by multiple maps t-SNE - https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4243097/

  8 When can t-SNE fail?

  When the data is noisy or on a highly varying manifold, the local linearity assumption we make (by converting Euclidean distances to similarities) does not hold. What can we do? - smooth data with PCA (in case of noise) or use Autoencoders (in case of highly varying manifold).

  $\newcommand{\lyxlock}{}$

  Guest author: thanks to Aviv Netanyahu for writing this post!

  1 What is visualization and why do we need it?

  Before working on a dataset we would like to understand it. After working on it, we want to communicate our conclusions with others easily. When dealing with high-dimensional data, it’s hard to get a ’feel’ for the data by just looking at numbers, and we will advance pretty slowly. For example, common visualization methods are histograms and scatter plots, which capture only one or two parameters at a time, meaning we would have to look at many histograms (one for each feature). However, by visualizing our data we may instantly be able to recognize structure. In order to do this we would like to find a representation of our data in $2$ D or $3$ D.

  2 Our goal

  Find an embedding of high-dimensional data to 2D or 3D (for simplicity, will assume $2$ D from now on).
  Formally, we want to learn $$\{\mathbf{x}_{1},...,\mathbf{x}_{N}\}\to\{\mathbf{y}_{1},...,\mathbf{y}_{N}\}$$ ($\{\mathbf{x}_{1},...,\mathbf{x}_{N}\}\subseteq\mathbb{R}^{n}$ is our data, $\{\mathbf{y}_{1},...,\mathbf{y}_{N}\}\subseteq\mathbb{R}^{2}$ is the embedding) such that the structure of the original data is preserved after the embedding.

  3 So why not PCA?

  PCA projects high-dimensional data onto a lower dimension determined by the directions where the data varies most in. The result is a linear projection which focuses on preserving distance between dissimilar data points. This is problematic:

  1. Linear sometimes isn’t good enough.
  2. We claim structure is mostly determined by similar data points.

  4 SNE (Hinton&Roweis)

  This algorithm for high-dimensional data visualization focuses on preserving high-dimensional ’neighbors’ in low-dimension.

  Assumption - ’neighbors’ in high dimension have small Euclidean distance. SNE converts the Euclidean distances to similarities using Gaussian distribution: $$p_{j|i}=\frac{\exp\left(-||\mathbf{x}_{i}-\mathbf{x}_{j}||^{2}/2\sigma_{i}^{2}\right)}{\sum\limits _{k\neq i}\exp\left(-||\mathbf{x}_{i}-\mathbf{x}_{k}||^{2}/2\sigma_{i}^{2}\right)}$$ $p_{j|i}$ is the probability that $j$ is the neighbor of $i$ in high-dimension. $$q_{j|i}=\frac{\exp\left(-||\mathbf{y}_{i}-\mathbf{y}_{j}||^{2}\right)}{\sum\limits _{k\neq i}\exp\left(-||\mathbf{y}_{i}-\mathbf{y}_{k}||^{2}\right)}$$ $q_{j|i}$ is the probability that $j$ is the neighbor of $i$ in low-dimension.
  Note - choosing $\sigma_{i}$ for $p_{j|i}$ is done by determining the perplexity, i.e. the effective number of neighbors for each data point. Then $\forall i$ , $\sigma_{i}$ is chosen such that $\mathbf{x}_{i}$ has according number of neighbors in the multivariate Gaussian with mean $\mathbf{x}_{i}$ and covariance matrix $\sigma_{i}\cdot I$ .
  
  By mapping each $\mathbf{x}_{i}\in\mathbb{R}^{n}$ to $\mathbf{y}_{i}\in\mathbb{R}^{2}$ such that $\forall i,j,\,p_{j|i}$ is close as possible to $q_{j|i}$ , we can preserve similarities. The algorithm finds such a mapping by gradient descent w.r.t. the following loss (cost) function: $$C=\sum_{i}KL(P_{i}||Q_{i})=\sum_{i}\sum_{j}p_{j|i}\log\frac{p_{j|i}}{q_{j|i}}$$ where $P_{i}=(p_{1|i},...,p_{N|i}),\,Q_{i}=(q_{1|i},...,q_{N|i})$ and the initialization of low-dimensional data is random around the origin.

  Notice the loss function, Kullback–Leibler (KL) divergence, is a natural way to measure similarity between two distributions.
  Also notice that this loss puts more emphasis on preserving the local structure of the data. KL-div is asymmetric, so different kinds of errors have a different affect:

  • Type 1 error: $p_{j|i}$ is large (i.e. two points were originally close) and $q_{j|i}$ small (i.e. these points were mapped too far), we pay a large penalty.
  • Type 2 error: $p_{j|i}$ is small (i.e. two points were originally far away) and $q_{j|i}$ large (i.e. these points were mapped too close), we pay a small penalty.

  This is good, as structure is preserved better by focusing on keeping close points together. This is easy to see in the following example: think of “unrolling” this swiss roll (like unrolling a carpet). Close points will stay within the same distance, far points won’t. This is true for many manifolds (spaces that locally resemble Euclidean space near each point) besides the swiss roll.

  

  4.1 The crowding problem

  The main problem with SNE - mapped points tend to “crowd” together, like in this example: images of 5 MNIST digits mapped with SNE to 2D and then marked by labels

  

  This happens because when reducing the dimensionality we need more “space” in which to represent our data-points. You can think of a sphere flattened to a circle where all distances are preserved (as much as possible) - the radius of the circle will have to be larger than the radius of the sphere!

  On one hand this causes type $1$ error (close points mapped too far), which forces far mapped points to come closer together. On the other, causes type $2$ error (far points mapped close) which forces them to move away from each other. Type $2$ error makes us pay a small price, however if many points cause error $2$ , the price will be higher, forcing the points to indeed move. What we just described causes the mapped points to correct themselves to the center such that they end up very close to each other. Consider the MNIST example above. If we did not know the labels we couldn’t make much sense of the mapping.

  5 T-SNE (van der Maaten&Hinton)

  5.1 Solving the crowding problem

  Instead of mapped points correcting themselves to minimize the loss, we can do this by increasing $q_{j|i}$ in advance. Recall, the loss is: $\sum\limits _{i}\sum\limits _{j}p_{j|i}\log\frac{p_{j|i}}{q_{j|i}}$ , thus larger $q_{j|i}$ minimizes it. We do this by defining: $$q_{ij}=\frac{(1+||\mathbf{y}_{i}-\mathbf{y}_{j}||)^{-1}}{\sum\limits _{k\neq l}(1+||\mathbf{y}_{k}-\mathbf{y}_{l}||)^{-1}}$$ using Cauchy distribution instead of Gaussian distribution. This increases $q_{j|i}$ since Cauchy distribution has a heavier tail than Gaussian distribution:

  

  (Hence the name of the algorithm - Cauchy distribution is a T-student distribution with $1$ degree of freedom).

  Note also that we changed the probability from conditional to joint probability (this was done for technical optimization reasons). For the same reason, $p_{j|i}$ was converted to $p_{ij}=\frac{p_{j|i}+p_{i|j}}{2N}$ . We can still use Gaussian distribution for the original data since we are trying to compensate for error in the mapped space.

  Now the performance on MNIST is much better:

  

  It’s obvious that global structure is preserved and crowding is resolved. Notice also that when zooming in on the mapping of $1$ digits, the local structure is preserved as well - digits with same orientation were mapped close to each other.

  6 Barnes-Hut approximation (van der Maaten)

  Computing the gradient for each mapped point each iteration is very costly, especially when the dataset is very large. At each iteration, $\forall i$ , we compute: $$\frac{\partial C}{\partial\mathbf{y}_{i}}=4\sum_{j\neq i}\frac{(p_{ij}-q_{ij})}{1+||\mathbf{y}_{i}-\mathbf{y}_{j}||^{2}}\cdot(\mathbf{y}_{i}-\mathbf{y}_{j})$$ We can think of $\mathbf{y}_{i}-\mathbf{y}_{j}$ as a spring between these two points (it might be too tight/loose depending on if the points were mapped to far/close accordingly). The fraction can be thought of as exertion/compression of the spring, in other words the correction we do in this iteration. We calculate this for $N$ mapped points, and derive $C$ by all data-points, so for each iteration we have $\mathcal{O}(N^{2})$ computations. This limits us to datasets with only few thousands of points! However, we can reduce the computation time to $\mathcal{O}(N\log N)$ by reducing the number of computations per iteration. The main idea - forces exerted by a group of points on a point that is relatively far away are all very similar. Therefore we can calculate $\mathbf{y}_{i}-\mathbf{y}_{j^{*}}$ where $\mathbf{y}_{j^{*}}$ is the average of the group of far away points, instead of $\mathbf{y}_{i}-\mathbf{y}_{j}$ for each $\mathbf{y}_{j}$ in the far away group.

  7 Some applications

  • Visualize word embeddings - http://nlp.yvespeirsman.be/blog/visualizing-word-embeddings-with-tsne/
  • Multiple maps word embeddings - http://homepage.tudelft.nl/19j49/multiplemaps/Multiple_maps_t-SNE/Multiple_maps_t-SNE.html
  • Visualization of genetic disease phenotype similarities by multiple maps t-SNE - https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4243097/

  8 When can t-SNE fail?

  When the data is noisy or on a highly varying manifold, the local linearity assumption we make (by converting Euclidean distances to similarities) does not hold. What can we do? - smooth data with PCA (in case of noise) or use Autoencoders (in case of highly varying manifold).
• May 26, 2017

  Greedy Algorithm for Multiway-Cut in Trees

  tags:
  TheoryLunch Blog

  Computer science topics over lunch

  Greedy Algorithm for Multiway-Cut in Trees

  May 26, 2017 • Eran Amar

  tags: $\newcommand{\lyxlock}{}$

  1 The Multiway Cut Problem

  The Multiway Cut (MC) problem is a generalization of the known Min-Cut problem. The input is a triple $\left(G,T,w\right)$ where $G$ is an undirected graph with a set $V\left(G\right)$ of $n$ vertices and a set $E\left(G\right)$ of $m$ edges, $T\subset V\left(G\right)$ is the set of terminals, and $w:E\left(G\right)\to\mathbb{R}_{>0}$ is a function that assign positive weight to each edge in the graph. A valid cut $C\subseteq E\left(G\right)$ is a set of edges such that any two vertices $t_{1},t_{2}\in T$ are disconnected over the sub-graph $H:=\left(V\left(G\right),E\left(G\right)\backslash C\right)$ . The weight of a given cut $C$ is the sum of wights of its edges, that is, $w\left(C\right)=\sum_{e\in C}w\left(e\right)$ . A solution for the MC problem is a valid cut of minimum weight.

  Let $k=\left|T\right|$ . Note that when $k=2$ the problem is reduced to the regular Min-Cut problem. For general graphs, any $k>2$ makes the problem NP-Hard. However, there are efficient algorithms for Multiway-Cut in trees. A tree is a connected acyclic graph. It is easy to see that if our graph is a tree, any two vertices in it must be connected by exactly one path. Moreover, $m=n-1$ . We will use these properties later when we prove that the algorithm gives an optimal solution.

  2 Efficient Algorithm for MC in Trees

  Let $\left(G,T,w\right)$ be the input to the algorithm, and $G$ is a tree. To simplify the analysis later on, we will work slightly differently; instead of trying to find a small-weight set $C\subseteq E\left(G\right)$ to discard from the graph, we will try to find a large-weight set $S\subseteq E\left(G\right)$ to keep in the graph. Note that obviously $S=V\left(G\right)\backslash C$ and $C=V\left(G\right)\backslash S$ . Therefore, finding $S$ with maximum weight is equivalent to finding $C$ with minimum weight.

  The idea is to start with a trivial (and probably not optimal) solution $S=\emptyset$ , that is, we don’t keep any edge in the graph such that all the terminals are separated, and we will improve $S$ iteratively.

  We start off with some definitions:

  A connected component (CC) of a vertex $u$ over a graph $G$ is the set of all vertices that can be reached from $u$ in $G$ . Note that if $\left|V\left(G\right)\right|=n$ then there are at most $n$ CCs in $G$ (each vertex has its own CC, which is the case when the graph has no edges).

  A vertex $t$ inside some CC is a master node if it is a terminal vertex. We then say that all the other vertices in that CC are dominated by $t$ .

  Now it is time to describe the algorithm:

  1. Start by sorting the $n-1$ edges of $E$ in descending order according to their weights. Let $e_{1},e_{2},..,e_{n-1}$ be the ordered sequence such that $\forall i,\:w\left(e_{i}\right)>w\left(e_{i+1}\right)$ .
  2. Initialize a trivial solution $S_{0}=\emptyset$ . Maintain a set $R_{0}=\left\{ \left\{ v\right\} \mid v\in V\left(G\right)\right\}$ of all the CCs in the resulting graph (each vertex is its own CC).
  3. For each $i=1,...,n-1$ do the following:
     1. Consider the endpoints of the edge $e_{i}=\left\{ v_{1},v_{2}\right\}$ and let $CC_{1},CC_{2}$ be the connected components of $v_{1},v_{2}$ respectively. If at least one of them do not has a master node, update $S_{i}=S_{i-1}\cup\left\{ e_{i}\right\}$ and let $R_{i}$ be the same as $R_{i-1}$ but with $CC_{1},CC_{2}$ merged into a new single connected component. Otherwise $S_{i}=S_{i-1}$ and $R_{i}=R_{i-1}$ .
  4. Return the set $S_{n-1}$ .

  3 The Solution is Correct and Optimal

  It is easy to see that $S_{n-1}$ is a valid solution; we started with all the terminals separated in their own CC, and in each iteration we never connect two CCs that have different masters (i.e. terminals). Therefore, in the returned set, all the connected components has at most one master, meaning that the terminals are separated.

  The fact that $S_{n-1}$ has maximum weight follows by the two following claims:

  1. Any optimal solution $S^{*}\subseteq E\left(G\right)$ has $\left|S^{*}\right|=n-k$ . Moreover, $\left|S_{n-1}\right|=n-k$ .
  2. For any iteration $i$ of the algorithm, there is an (possibly different) optimal solution $S^{*}$ such that $S_{i}\subseteq S^{*}$ .

  So we conclude that there exists some optimal $S^{*}$ such that $S_{n-1}=S^{*}$ , thus $S_{n-1}$ has maximum weight.

  We move to state and prove those claims formally.

  Let $\mathcal{X}$ be the set of all maximum-wight solutions for Multiway-Cut for the instance $\left(G,T,w\right)$ and $G$ is a tree. Then any $S^{*}\in\mathcal{X}$ have $\left|S^{*}\right|=n-k$ . Moreover, $\left|S_{n-1}\right|=n-k$ .

  Since we started from a graph that is a tree, the number of CCs in the final sub-graph is exactly $n-\left|S\right|$ , for any $S\subseteq E\left(G\right)$ . Let $S^{*}$ be some optimal solution. It left to show that the sub-graph $H^{*}=\left(V\left(G\right),S^{*}\right)$ has exactly $k$ connected components:
  
  Assume that there are less than $k$ , then we must have two terminals which live in the same CC - invalidating the solution $S^{*}$ . If we assume that there are more than $k$ CCs, then we must have a connected component without a master, so we can safely connect it to some other CC with some edge $e$ that is not in $S^{*}$ (and such edge exists since $G$ is connected). That way, we get a new valid solution with greater weight than the weight of $S^{*}$ which is a contradiction.
  
  Similarly, we prove that $S_{n-1}$ has $k$ CCs: If it has less than $k$ , there must be two terminals in the same CCs, and that contradicts the correctness of $S_{n-1}$ . If it has more than $k$ CCs after $n-1$ iterations, there must be a CC without a master, so we may connect it with some edge, let it be $e_{i}$ , to other CC without invalidating the solution. But that $e_{i}$ wasn’t picked when the algorithm considered it - meaning that it could cause a violation in the sub-graph of the $i$ th iteration, and since that sub-graph is contained in the sub-graph of the last iteration then we can’t add $e_{i}$ to our solution - a contradiction.

  Note that since any optimal solution have exactly $k$ CCs, we have that each CC must contains one terminal. So in the final sub-graph defined by the solution, any $u\in V\left(G\right)$ is dominated by exactly one master mode $t\in T$ .

  Let $\mathcal{X}$ be the set of all optimal solutions for the instance $\left(G,T,w\right)$ and $G$ in a tree. For any $i\in\left\{ 0,1,..,n-1\right\}$ let $\mathcal{X}_{i}=\left\{ S^{*}\in\mathcal{X}\mid S_{i}\subset S^{*}\right\}$ be the set of optimal solutions such that $S_{i}$ contained in them. Then, for any $i$ , $\mathcal{X}_{i}\ne\emptyset$ .

  For $i=0$ we have $S_{0}=\emptyset$ thus $\mathcal{X}_{0}=\mathcal{X}$ which is not empty. Assume toward contradiction that the claim is false, and let $i$ be the smallest such that $\mathcal{X}_{i}=\emptyset$ , so we know that $\mathcal{X}_{i-1}$ is not empty, and hence $S_{i-1}\ne S_{i}$ . In other words, during the $i$ th iteration, the algorithm added the edge $e_{i}=\left\{ v,u\right\}$ to $S_{i-1}$ , causing $\mathcal{X}_{i}$ to be empty.
  
  Fix some $S^{*}\in\mathcal{X}_{i-1}$ . Denote $H_{i-1}=\left(V\left(G\right),S_{i-1}\right)$ and $H^{*}=\left(V\left(G\right),S^{*}\right)$ the subgraphs of $G$ . From the update rule of the algorithm, at least one of $\left\{ u,v\right\}$ is not dominated by a master over $H_{i-1}$ , w.l.o.g assume it is $u$ . Let $r\in T$ be the master that dominates $u$ over $H^{*}$ , and consider the path connecting $r\rightsquigarrow u$ , note that any edge along that path belongs to $S^{*}$ . Since the path does not exists in $H_{n-1}$ and in particular in $H_{i-1}$ , there is at least one edge that is part of the path and not in $S_{i-1}$ , let $e_{j}$ be such one with minimum weight.
  
  Observe that it is impossible that $w\left(e_{j}\right)>w\left(e_{i}\right)$ because it means that the algorithm already rejected $e_{j}$ on the $j$ th iteration and since it was a violation for $H_{j}$ and $H_{j}\subset H_{i-1}$ then it is also a violation for $H_{i-1}$ therefore $e_{j}$ does not belongs to any optimal solution in $\mathcal{X}_{i-1}$ , in contradiction that $e_{j}\in S^{*}\in\mathcal{X}_{i-1}$ .
  
  So suppose $w\left(e_{j}\right)\le w\left(e_{i}\right)$ . Note that in the sub-graph $H^{*}$ induced by $S^{*}$ , $v$ must be dominated my some terminal $t\in T$ . Because $H^{*}$ is a sub-graph of a tree, if $t=r$ then we can close a cycle by adding $e_{i}$ to $S^{*}$ , thus $t\ne r$ . Consider the sub-graph $H'=\left(V\left(G\right),S^{*}\cup\left\{ e_{i}\right\} \right)$ . The edge $e_{i}$ connects now the terminals $t$ and $r$ , and that is the only violation in $H'$ . Since $H'$ is a sub-graph of a tree there is only one path connecting $t$ and $r$ and it must contains the path $r\rightsquigarrow u$ , so removing the edge $e_{j}$ will break the connection between $t,r$ leading to a valid solution $S'=\left(S^{*}\cup\left\{ e_{i}\right\} \right)\backslash\left\{ e_{j}\right\}$ . By the assumption that $w\left(e_{j}\right)\le w\left(e_{i}\right)$ the weight of $S'$ can only increase (in fact, the optimality of $S^{*}$ forces $w\left(e_{j}\right)=w\left(e_{i}\right)$ ) and we get an optimal solution that contains $S_{i}$ , therefore $\mathcal{X}_{i}\ne\emptyset$ , contradiction.

  4 Running Time and Similar Algorithms

  Note that $\left|E\left(G\right)\right|=n-1$ , so step 1 of the algorithm takes $\mathcal{O}\left(n\log n\right)$ . Step 2 is order of $n$ , times the cost of the operations inside the loop. There are data structures such as Union-Find that allow us to track the CCs of the current sub-graph in amortize small running time that is smaller than $\log n$ . Therefore, the total running time is bounded by $\mathcal{O}$ $\left(n\log n\right)$ .

  Note that the algorithm presented in section 2, is similar to Kruskal’s algorithm for finding the minimum spanning tree. In both cases, the greedy approach leads to an optimal solution.

  $\newcommand{\lyxlock}{}$

  1 The Multiway Cut Problem

  The Multiway Cut (MC) problem is a generalization of the known Min-Cut problem. The input is a triple $\left(G,T,w\right)$ where $G$ is an undirected graph with a set $V\left(G\right)$ of $n$ vertices and a set $E\left(G\right)$ of $m$ edges, $T\subset V\left(G\right)$ is the set of terminals, and $w:E\left(G\right)\to\mathbb{R}_{>0}$ is a function that assign positive weight to each edge in the graph. A valid cut $C\subseteq E\left(G\right)$ is a set of edges such that any two vertices $t_{1},t_{2}\in T$ are disconnected over the sub-graph $H:=\left(V\left(G\right),E\left(G\right)\backslash C\right)$ . The weight of a given cut $C$ is the sum of wights of its edges, that is, $w\left(C\right)=\sum_{e\in C}w\left(e\right)$ . A solution for the MC problem is a valid cut of minimum weight.

  Let $k=\left|T\right|$ . Note that when $k=2$ the problem is reduced to the regular Min-Cut problem. For general graphs, any $k>2$ makes the problem NP-Hard. However, there are efficient algorithms for Multiway-Cut in trees. A tree is a connected acyclic graph. It is easy to see that if our graph is a tree, any two vertices in it must be connected by exactly one path. Moreover, $m=n-1$ . We will use these properties later when we prove that the algorithm gives an optimal solution.

  2 Efficient Algorithm for MC in Trees

  Let $\left(G,T,w\right)$ be the input to the algorithm, and $G$ is a tree. To simplify the analysis later on, we will work slightly differently; instead of trying to find a small-weight set $C\subseteq E\left(G\right)$ to discard from the graph, we will try to find a large-weight set $S\subseteq E\left(G\right)$ to keep in the graph. Note that obviously $S=V\left(G\right)\backslash C$ and $C=V\left(G\right)\backslash S$ . Therefore, finding $S$ with maximum weight is equivalent to finding $C$ with minimum weight.

  The idea is to start with a trivial (and probably not optimal) solution $S=\emptyset$ , that is, we don’t keep any edge in the graph such that all the terminals are separated, and we will improve $S$ iteratively.

  We start off with some definitions:

  A connected component (CC) of a vertex $u$ over a graph $G$ is the set of all vertices that can be reached from $u$ in $G$ . Note that if $\left|V\left(G\right)\right|=n$ then there are at most $n$ CCs in $G$ (each vertex has its own CC, which is the case when the graph has no edges).

  A vertex $t$ inside some CC is a master node if it is a terminal vertex. We then say that all the other vertices in that CC are dominated by $t$ .

  Now it is time to describe the algorithm:

  1. Start by sorting the $n-1$ edges of $E$ in descending order according to their weights. Let $e_{1},e_{2},..,e_{n-1}$ be the ordered sequence such that $\forall i,\:w\left(e_{i}\right)>w\left(e_{i+1}\right)$ .
  2. Initialize a trivial solution $S_{0}=\emptyset$ . Maintain a set $R_{0}=\left\{ \left\{ v\right\} \mid v\in V\left(G\right)\right\}$ of all the CCs in the resulting graph (each vertex is its own CC).
  3. For each $i=1,...,n-1$ do the following:
     1. Consider the endpoints of the edge $e_{i}=\left\{ v_{1},v_{2}\right\}$ and let $CC_{1},CC_{2}$ be the connected components of $v_{1},v_{2}$ respectively. If at least one of them do not has a master node, update $S_{i}=S_{i-1}\cup\left\{ e_{i}\right\}$ and let $R_{i}$ be the same as $R_{i-1}$ but with $CC_{1},CC_{2}$ merged into a new single connected component. Otherwise $S_{i}=S_{i-1}$ and $R_{i}=R_{i-1}$ .
  4. Return the set $S_{n-1}$ .

  3 The Solution is Correct and Optimal

  It is easy to see that $S_{n-1}$ is a valid solution; we started with all the terminals separated in their own CC, and in each iteration we never connect two CCs that have different masters (i.e. terminals). Therefore, in the returned set, all the connected components has at most one master, meaning that the terminals are separated.

  The fact that $S_{n-1}$ has maximum weight follows by the two following claims:

  1. Any optimal solution $S^{*}\subseteq E\left(G\right)$ has $\left|S^{*}\right|=n-k$ . Moreover, $\left|S_{n-1}\right|=n-k$ .
  2. For any iteration $i$ of the algorithm, there is an (possibly different) optimal solution $S^{*}$ such that $S_{i}\subseteq S^{*}$ .

  So we conclude that there exists some optimal $S^{*}$ such that $S_{n-1}=S^{*}$ , thus $S_{n-1}$ has maximum weight.

  We move to state and prove those claims formally.

  Let $\mathcal{X}$ be the set of all maximum-wight solutions for Multiway-Cut for the instance $\left(G,T,w\right)$ and $G$ is a tree. Then any $S^{*}\in\mathcal{X}$ have $\left|S^{*}\right|=n-k$ . Moreover, $\left|S_{n-1}\right|=n-k$ .

  Since we started from a graph that is a tree, the number of CCs in the final sub-graph is exactly $n-\left|S\right|$ , for any $S\subseteq E\left(G\right)$ . Let $S^{*}$ be some optimal solution. It left to show that the sub-graph $H^{*}=\left(V\left(G\right),S^{*}\right)$ has exactly $k$ connected components:
  
  Assume that there are less than $k$ , then we must have two terminals which live in the same CC - invalidating the solution $S^{*}$ . If we assume that there are more than $k$ CCs, then we must have a connected component without a master, so we can safely connect it to some other CC with some edge $e$ that is not in $S^{*}$ (and such edge exists since $G$ is connected). That way, we get a new valid solution with greater weight than the weight of $S^{*}$ which is a contradiction.
  
  Similarly, we prove that $S_{n-1}$ has $k$ CCs: If it has less than $k$ , there must be two terminals in the same CCs, and that contradicts the correctness of $S_{n-1}$ . If it has more than $k$ CCs after $n-1$ iterations, there must be a CC without a master, so we may connect it with some edge, let it be $e_{i}$ , to other CC without invalidating the solution. But that $e_{i}$ wasn’t picked when the algorithm considered it - meaning that it could cause a violation in the sub-graph of the $i$ th iteration, and since that sub-graph is contained in the sub-graph of the last iteration then we can’t add $e_{i}$ to our solution - a contradiction.

  Note that since any optimal solution have exactly $k$ CCs, we have that each CC must contains one terminal. So in the final sub-graph defined by the solution, any $u\in V\left(G\right)$ is dominated by exactly one master mode $t\in T$ .

  Let $\mathcal{X}$ be the set of all optimal solutions for the instance $\left(G,T,w\right)$ and $G$ in a tree. For any $i\in\left\{ 0,1,..,n-1\right\}$ let $\mathcal{X}_{i}=\left\{ S^{*}\in\mathcal{X}\mid S_{i}\subset S^{*}\right\}$ be the set of optimal solutions such that $S_{i}$ contained in them. Then, for any $i$ , $\mathcal{X}_{i}\ne\emptyset$ .

  For $i=0$ we have $S_{0}=\emptyset$ thus $\mathcal{X}_{0}=\mathcal{X}$ which is not empty. Assume toward contradiction that the claim is false, and let $i$ be the smallest such that $\mathcal{X}_{i}=\emptyset$ , so we know that $\mathcal{X}_{i-1}$ is not empty, and hence $S_{i-1}\ne S_{i}$ . In other words, during the $i$ th iteration, the algorithm added the edge $e_{i}=\left\{ v,u\right\}$ to $S_{i-1}$ , causing $\mathcal{X}_{i}$ to be empty.
  
  Fix some $S^{*}\in\mathcal{X}_{i-1}$ . Denote $H_{i-1}=\left(V\left(G\right),S_{i-1}\right)$ and $H^{*}=\left(V\left(G\right),S^{*}\right)$ the subgraphs of $G$ . From the update rule of the algorithm, at least one of $\left\{ u,v\right\}$ is not dominated by a master over $H_{i-1}$ , w.l.o.g assume it is $u$ . Let $r\in T$ be the master that dominates $u$ over $H^{*}$ , and consider the path connecting $r\rightsquigarrow u$ , note that any edge along that path belongs to $S^{*}$ . Since the path does not exists in $H_{n-1}$ and in particular in $H_{i-1}$ , there is at least one edge that is part of the path and not in $S_{i-1}$ , let $e_{j}$ be such one with minimum weight.
  
  Observe that it is impossible that $w\left(e_{j}\right)>w\left(e_{i}\right)$ because it means that the algorithm already rejected $e_{j}$ on the $j$ th iteration and since it was a violation for $H_{j}$ and $H_{j}\subset H_{i-1}$ then it is also a violation for $H_{i-1}$ therefore $e_{j}$ does not belongs to any optimal solution in $\mathcal{X}_{i-1}$ , in contradiction that $e_{j}\in S^{*}\in\mathcal{X}_{i-1}$ .
  
  So suppose $w\left(e_{j}\right)\le w\left(e_{i}\right)$ . Note that in the sub-graph $H^{*}$ induced by $S^{*}$ , $v$ must be dominated my some terminal $t\in T$ . Because $H^{*}$ is a sub-graph of a tree, if $t=r$ then we can close a cycle by adding $e_{i}$ to $S^{*}$ , thus $t\ne r$ . Consider the sub-graph $H'=\left(V\left(G\right),S^{*}\cup\left\{ e_{i}\right\} \right)$ . The edge $e_{i}$ connects now the terminals $t$ and $r$ , and that is the only violation in $H'$ . Since $H'$ is a sub-graph of a tree there is only one path connecting $t$ and $r$ and it must contains the path $r\rightsquigarrow u$ , so removing the edge $e_{j}$ will break the connection between $t,r$ leading to a valid solution $S'=\left(S^{*}\cup\left\{ e_{i}\right\} \right)\backslash\left\{ e_{j}\right\}$ . By the assumption that $w\left(e_{j}\right)\le w\left(e_{i}\right)$ the weight of $S'$ can only increase (in fact, the optimality of $S^{*}$ forces $w\left(e_{j}\right)=w\left(e_{i}\right)$ ) and we get an optimal solution that contains $S_{i}$ , therefore $\mathcal{X}_{i}\ne\emptyset$ , contradiction.

  4 Running Time and Similar Algorithms

  Note that $\left|E\left(G\right)\right|=n-1$ , so step 1 of the algorithm takes $\mathcal{O}\left(n\log n\right)$ . Step 2 is order of $n$ , times the cost of the operations inside the loop. There are data structures such as Union-Find that allow us to track the CCs of the current sub-graph in amortize small running time that is smaller than $\log n$ . Therefore, the total running time is bounded by $\mathcal{O}$ $\left(n\log n\right)$ .

  Note that the algorithm presented in section 2, is similar to Kruskal’s algorithm for finding the minimum spanning tree. In both cases, the greedy approach leads to an optimal solution.
• May 18, 2017

  Recurrent Neural Network - Recent Advancements

  tags:
  TheoryLunch Blog

  Computer science topics over lunch

  Recurrent Neural Network - Recent Advancements

  May 18, 2017 • Eran Amar

  tags: $\newcommand{\lyxlock}{}$

  1 Last Post of the Series

  This post is the last in the series that discussed Recurrent Neural Networks (RNNs). The first post introduced the basic recurrent cells and their limitations, the second post presented the formulation of gated cells as LSTMs and GRUs, and the third post overviewed dropout variants for LSTMs. This post will present some of the recent advancements in the filed of RNNs (mostly related to LSTM units).

  We start by reviewing some of the notations that we used during the series.

  We defined $\mathcal{S}_{\sigma}$ to be the family of all affine transformations $\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}^{m}$ for any dimensions $n$ and $m$ , followed by an element-wise activation function $\sigma$ , that is,

  $$\mathcal{S}_{\sigma}=\left\{ \mathbf{x},\mathbf{h}\mapsto\sigma\left(W\mathbf{x}+U\mathbf{h}+\mathbf{b}\right)\mid\forall W,U,\mathbf{b}\right\}$$ and noted that with the $sigmoid\left(x\right)=\left(1+e^{-x}\right)^{-1}$ as the activation function, all the functions in $\mathcal{S}_{sig}$ are gate functions, that is, their output belongs to $\left[0,1\right]^{m}$ .

  Usually we do not specify $n$ and $m$ , we assume that whenever there is a matrix multiplication of element-wise operation the dimensions always match, as they are not important for the discussion.

  2 Variations of the Gate Functions for LSTMs

  2.1 Recap: LSTM Formulation - The Use of Gates

  LSTM layer is a recurrent layer that keeps a memory vector $\mathbf{c}^{\left(t\right)}$ for any time step $t$ . Upon receiving the current time-step input $\mathbf{x}^{\left(t\right)}$ and the output from previous time-step $\mathbf{h}^{\left(t-1\right)}$ , the layer computes gate activations $f^{\left(t\right)},i^{\left(t\right)},o^{\left(t\right)}$ by applying a the corresponding gate function $\mathbf{f},\mathbf{i},\mathbf{o}\in\mathcal{S}_{sig}$ on the tuple $\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)$ , note that these activations are vectors in $\left[0,1\right]^{m}$ .

  A new memory-vector is then being computed by $$\mathbf{c}^{\left(t\right)}=f^{\left(t\right)}\circ\mathbf{c}^{\left(t-1\right)}+i^{\left(t\right)}\circ\mathbf{s}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)$$ where $\circ$ is element-wise multiplication and $\mathbf{s}\in\mathcal{S}_{tanh}$ is the state-transition function. The layer then outputs $\mathbf{h}^{\left(t\right)}=o^{\left(t\right)}\circ tanh\left(\mathbf{c}^{\left(t\right)}\right)$ .

  The gate functions are of great importance for LSTM layer, they are what allow the layer to learn short and long term dependencies, and they also help avoiding the exploding and vanishing gradients problem.

  2.2 Simple Variations

  The more expressive power a gate function have, the better it can learn short \ long term dependencies. That is, if a gate function is very rich, it can account for subtle changes in the input $\mathbf{x}^{\left(t\right)}$ as well to subtle changes in the “past”, i.e. in $\mathbf{h}^{\left(t-1\right)}$ . One way to enrich a gate function is by making it depends on the previous memory vector $\mathbf{c}^{\left(t-1\right)}$ in addition to the regular tuple $\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)$ . That will redefine our gate functions family to be $$\mathcal{S}_{\sigma}=\left\{ \mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)},\mathbf{c}^{\left(t-1\right)}\mapsto\sigma\left(W\mathbf{x}^{\left(t\right)}+U\mathbf{h}^{\left(t-1\right)}+V\mathbf{c}^{\left(t-1\right)}+\mathbf{b}\right)\right\}$$ for any $W,U$ and diagonal $V$ . Such construction can be found in here, note that we can think of $V\mathbf{c}^{\left(t-1\right)}$ as a bias vector that is memory-dependent, thus there are formulation in which the vector $\mathbf{b}$ is omitted.

  2.3 Multiplicative LSTMs

  In this paper from October 2016, the authors took that idea one step further. They wanted to make a unique gate function for each possible input. Recall that any function in $\mathcal{S}_{\sigma}$ is of the form $\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\mapsto W\mathbf{x}^{\left(t\right)}+U\mathbf{h}^{\left(t-1\right)}+\mathbf{b}$ . For simplicity, consider the case where $\mathbf{b}=\mathbf{0}$ .

  We end up with a sum of two components: one that depends on the current input and one that depends on the past (i.e. $\mathbf{h}^{\left(t-1\right)}$ which encodes information about previous time steps), and the component with larger magnitude will dominate the transition. If $W\mathbf{x}^{\left(t\right)}$ is larger, the layer will not be sensitive enough to the past, and if $U\mathbf{h}^{\left(t-1\right)}$ is larger, then the layer will not be sensitive to subtle changes in the input.

  The author noted that since in most cases the input $\mathbf{x}^{\left(t\right)}$ is an 1-hot vector, then multiply by $W$ is just selecting a specific column. So we may think of any gate function $\mathbf{s}\in\mathcal{S}_{\sigma}$ as a fixed base affine transformation $\mathbf{h}\mapsto U\mathbf{h}$ combined with input-dependent bias vector $\mathbf{b}_{\mathbf{x}^{\left(t\right)}}=W\mathbf{x}^{\left(t\right)}$ . That is $$\mathbf{s}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)=U\mathbf{h}^{\left(t-1\right)}+\mathbf{b}_{\mathbf{x}^{\left(t\right)}}$$ that formula emphasis the additive effect of the current input-vector on the transition of $\mathbf{h}^{\left(t-1\right)}$ .

  The goal in multiplicative LSTM (mLSTM), is to have an unique affine transformation (in our case, unique matrix $U$ ) for each possible input. Obviously, if the number of possible inputs is very large, the number of the parameters will explode and it won’t be feasible to train the network. To overcome that, the authors of the paper suggested to learn shared intermediate matrices that will be used to construct an unique $U$ for each input, and because the factorization is shared, there are less parameters to learn.

  The factorization is defined as follow: the unique matrix $U_{\mathbf{x}^{\left(t\right)}}$ is constructed by $V_{1}diag\left(T\mathbf{x}^{\left(t\right)}\right)V_{2}$ where $V_{1},V_{2}$ and $T$ are intermediate matrices which are shared across all the possible inputs, and $diag\left(\mathbf{v}\right)$ is an operation that maps any vector $\mathbf{v}$ to a square diagonal matrix with the elements of $\mathbf{v}$ on its diagonal.

  Note that the target dimension of $T$ may be arbitrarily large (in the paper they chose it to be the same as $\mathbf{h}^{\left(t\right)}$ ).

  The difference between LSTM and mLSTM is only in the definition of the gate functions family, all the equations for updating the memory cell and the output are the same. In mLSTM the family $\mathcal{S}_{\sigma}$ is being replaced with $$\begin{aligned} \mathcal{S}'_{\sigma} & =\left\{ \mathbf{x},\mathbf{h}\mapsto\sigma\left(W\mathbf{x}+U_{\mathbf{x}}\mathbf{h}\right)\right\} \\ & =\left\{ \mathbf{x},\mathbf{h}\mapsto\sigma\left(\mathbf{b}_{\mathbf{x}}+V_{1}diag\left(T\mathbf{x}\right)V_{2}\mathbf{h}\right)\mid\forall W,V_{1},V_{2},T\right\} \end{aligned}$$ Note that we can reduce the number of parameters even further by forcing all the gate functions to use the same $V_{2}$ and $T$ matrices. That is, each gate will be parametrized only by $W$ and $V_{1}$ .

  Formally, define the following transformation $\tau:\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}^{m}$ such that $\tau\left(\mathbf{x},\mathbf{h}\right)=T\mathbf{x}\circ V_{2}\mathbf{h}$ and then we can define $\mathcal{S}'_{\sigma}$ for some fixed learned transformation $\tau$ as $$\begin{aligned} \mathcal{S}'_{\sigma,\tau} & =\left\{ \mathbf{x},\mathbf{h}\mapsto\sigma\left(W\mathbf{x}+V_{1}\tau\left(\mathbf{x},\mathbf{h}\right)\right)\mid\forall W,V_{1}\right\} \\ & =\left\{ \mathbf{x},\mathbf{h}\mapsto\mathbf{s}\left(\mathbf{x},\tau\left(\mathbf{x},\mathbf{h}\right)\right)\mid\mathbf{s}\in\mathcal{S}_{\sigma}\right\} \end{aligned}$$ in other words, mLSTM is an LSTM that apply its gates to the tuple $\left(\mathbf{x}^{\left(t\right)},\tau\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)\right)$ rather than $\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)$ , and $\tau$ is another learned transformation. If you look again at the exact formula of $\tau$ you will see the multiplicative effect of the input vector on the transformation of $\mathbf{h}$ . That way, the authors said, $\mathcal{S}'_{\sigma,\tau}$ can yield much richer family of input-dependent transformations for $\mathbf{h}^{\left(t\right)}$ , which can be sensitive to the past as well as to subtle changes in the current input.

  3 Recurrent Batch Normalization

  Batch normalization is an operator applied to a layer before going through the non-linearity function in order to “normalize” the values before activation. That operator has two hyper-parameters for tuning and two statistics that it accumulates internally. Formally, given a vector of pre-activation values $\mathbf{x}$ , batch-normalization is $$BN\left(\mathbf{x};\gamma,\beta\right)=\beta+\gamma\circ\frac{\mathbf{x}-\mathbb{\hat{E}\left[\mathbf{x}\right]}}{\sqrt{\hat{\mathbb{V}}\left[\mathbf{x}\right]+\epsilon}}$$ where $\hat{\mathbb{E}}$ and $\hat{\mathbb{V}}$ are the empirical mean and variance of the current batch respectively. $\gamma,\beta$ and $\epsilon$ are vectors, and the division in the operator is being computed element-wise. Note that at inference time, the population statistics are estimated by averaging the empirical statistics across all the batches.

  That paper from February 2017, applied batch normalization also to the recurrent connections in LSTM layers, and they showed empirically that it helped the network to converge faster. The usage is simple, each gate function $\mathbf{i},\mathbf{f},\mathbf{o}\in\mathcal{S}_{sig}$ and the transition function $\mathbf{s}\in\mathcal{S}_{tanh}$ computes $$\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\mapsto\sigma\left(BN\left(W\mathbf{h}^{\left(t-1\right)};\gamma_{h},\beta_{h}\right)+BN\left(W\mathbf{x}^{\left(t\right)};\gamma_{x},\beta_{x}\right)+\mathbf{b}\right)$$ when the hyper-parameters $\gamma$ and $\beta$ are shared across different gates. Then the output of the layer is $\mathbf{h}^{\left(t\right)}=o^{\left(t\right)}\circ tanh\left(BN\left(\mathbf{c}^{\left(t\right)};\gamma_{c},\beta_{c}\right)\right)$ . The authors suggested to set $\beta_{x}=\beta_{h}=\mathbf{0}$ because there is already a bias parameter $\mathbf{b}$ and to prevent redundancy. In addition, they said that sharing the internal $BN$ statistics across time degrades performance severely. Therefore, one should use “fresh” $BN$ operators for each time-step with their own internal statistics (but share the $\beta$ and $\gamma$ parameters).

  4 The Ever-growing Field

  RNNs is a rapidly growing field and this series only covers a small part of it. There are much more advanced models, some of them are only from few months ago. If you are interested in this field, you may find very interesting stuff here.

  $\newcommand{\lyxlock}{}$

  1 Last Post of the Series

  This post is the last in the series that discussed Recurrent Neural Networks (RNNs). The first post introduced the basic recurrent cells and their limitations, the second post presented the formulation of gated cells as LSTMs and GRUs, and the third post overviewed dropout variants for LSTMs. This post will present some of the recent advancements in the filed of RNNs (mostly related to LSTM units).

  We start by reviewing some of the notations that we used during the series.

  We defined $\mathcal{S}_{\sigma}$ to be the family of all affine transformations $\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}^{m}$ for any dimensions $n$ and $m$ , followed by an element-wise activation function $\sigma$ , that is,

  $$\mathcal{S}_{\sigma}=\left\{ \mathbf{x},\mathbf{h}\mapsto\sigma\left(W\mathbf{x}+U\mathbf{h}+\mathbf{b}\right)\mid\forall W,U,\mathbf{b}\right\}$$ and noted that with the $sigmoid\left(x\right)=\left(1+e^{-x}\right)^{-1}$ as the activation function, all the functions in $\mathcal{S}_{sig}$ are gate functions, that is, their output belongs to $\left[0,1\right]^{m}$ .

  Usually we do not specify $n$ and $m$ , we assume that whenever there is a matrix multiplication of element-wise operation the dimensions always match, as they are not important for the discussion.

  2 Variations of the Gate Functions for LSTMs

  2.1 Recap: LSTM Formulation - The Use of Gates

  LSTM layer is a recurrent layer that keeps a memory vector $\mathbf{c}^{\left(t\right)}$ for any time step $t$ . Upon receiving the current time-step input $\mathbf{x}^{\left(t\right)}$ and the output from previous time-step $\mathbf{h}^{\left(t-1\right)}$ , the layer computes gate activations $f^{\left(t\right)},i^{\left(t\right)},o^{\left(t\right)}$ by applying a the corresponding gate function $\mathbf{f},\mathbf{i},\mathbf{o}\in\mathcal{S}_{sig}$ on the tuple $\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)$ , note that these activations are vectors in $\left[0,1\right]^{m}$ .

  A new memory-vector is then being computed by $$\mathbf{c}^{\left(t\right)}=f^{\left(t\right)}\circ\mathbf{c}^{\left(t-1\right)}+i^{\left(t\right)}\circ\mathbf{s}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)$$ where $\circ$ is element-wise multiplication and $\mathbf{s}\in\mathcal{S}_{tanh}$ is the state-transition function. The layer then outputs $\mathbf{h}^{\left(t\right)}=o^{\left(t\right)}\circ tanh\left(\mathbf{c}^{\left(t\right)}\right)$ .

  The gate functions are of great importance for LSTM layer, they are what allow the layer to learn short and long term dependencies, and they also help avoiding the exploding and vanishing gradients problem.

  2.2 Simple Variations

  The more expressive power a gate function have, the better it can learn short \ long term dependencies. That is, if a gate function is very rich, it can account for subtle changes in the input $\mathbf{x}^{\left(t\right)}$ as well to subtle changes in the “past”, i.e. in $\mathbf{h}^{\left(t-1\right)}$ . One way to enrich a gate function is by making it depends on the previous memory vector $\mathbf{c}^{\left(t-1\right)}$ in addition to the regular tuple $\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)$ . That will redefine our gate functions family to be $$\mathcal{S}_{\sigma}=\left\{ \mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)},\mathbf{c}^{\left(t-1\right)}\mapsto\sigma\left(W\mathbf{x}^{\left(t\right)}+U\mathbf{h}^{\left(t-1\right)}+V\mathbf{c}^{\left(t-1\right)}+\mathbf{b}\right)\right\}$$ for any $W,U$ and diagonal $V$ . Such construction can be found in here, note that we can think of $V\mathbf{c}^{\left(t-1\right)}$ as a bias vector that is memory-dependent, thus there are formulation in which the vector $\mathbf{b}$ is omitted.

  2.3 Multiplicative LSTMs

  In this paper from October 2016, the authors took that idea one step further. They wanted to make a unique gate function for each possible input. Recall that any function in $\mathcal{S}_{\sigma}$ is of the form $\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\mapsto W\mathbf{x}^{\left(t\right)}+U\mathbf{h}^{\left(t-1\right)}+\mathbf{b}$ . For simplicity, consider the case where $\mathbf{b}=\mathbf{0}$ .

  We end up with a sum of two components: one that depends on the current input and one that depends on the past (i.e. $\mathbf{h}^{\left(t-1\right)}$ which encodes information about previous time steps), and the component with larger magnitude will dominate the transition. If $W\mathbf{x}^{\left(t\right)}$ is larger, the layer will not be sensitive enough to the past, and if $U\mathbf{h}^{\left(t-1\right)}$ is larger, then the layer will not be sensitive to subtle changes in the input.

  The author noted that since in most cases the input $\mathbf{x}^{\left(t\right)}$ is an 1-hot vector, then multiply by $W$ is just selecting a specific column. So we may think of any gate function $\mathbf{s}\in\mathcal{S}_{\sigma}$ as a fixed base affine transformation $\mathbf{h}\mapsto U\mathbf{h}$ combined with input-dependent bias vector $\mathbf{b}_{\mathbf{x}^{\left(t\right)}}=W\mathbf{x}^{\left(t\right)}$ . That is $$\mathbf{s}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)=U\mathbf{h}^{\left(t-1\right)}+\mathbf{b}_{\mathbf{x}^{\left(t\right)}}$$ that formula emphasis the additive effect of the current input-vector on the transition of $\mathbf{h}^{\left(t-1\right)}$ .

  The goal in multiplicative LSTM (mLSTM), is to have an unique affine transformation (in our case, unique matrix $U$ ) for each possible input. Obviously, if the number of possible inputs is very large, the number of the parameters will explode and it won’t be feasible to train the network. To overcome that, the authors of the paper suggested to learn shared intermediate matrices that will be used to construct an unique $U$ for each input, and because the factorization is shared, there are less parameters to learn.

  The factorization is defined as follow: the unique matrix $U_{\mathbf{x}^{\left(t\right)}}$ is constructed by $V_{1}diag\left(T\mathbf{x}^{\left(t\right)}\right)V_{2}$ where $V_{1},V_{2}$ and $T$ are intermediate matrices which are shared across all the possible inputs, and $diag\left(\mathbf{v}\right)$ is an operation that maps any vector $\mathbf{v}$ to a square diagonal matrix with the elements of $\mathbf{v}$ on its diagonal.

  Note that the target dimension of $T$ may be arbitrarily large (in the paper they chose it to be the same as $\mathbf{h}^{\left(t\right)}$ ).

  The difference between LSTM and mLSTM is only in the definition of the gate functions family, all the equations for updating the memory cell and the output are the same. In mLSTM the family $\mathcal{S}_{\sigma}$ is being replaced with $$\begin{aligned} \mathcal{S}'_{\sigma} & =\left\{ \mathbf{x},\mathbf{h}\mapsto\sigma\left(W\mathbf{x}+U_{\mathbf{x}}\mathbf{h}\right)\right\} \\ & =\left\{ \mathbf{x},\mathbf{h}\mapsto\sigma\left(\mathbf{b}_{\mathbf{x}}+V_{1}diag\left(T\mathbf{x}\right)V_{2}\mathbf{h}\right)\mid\forall W,V_{1},V_{2},T\right\} \end{aligned}$$ Note that we can reduce the number of parameters even further by forcing all the gate functions to use the same $V_{2}$ and $T$ matrices. That is, each gate will be parametrized only by $W$ and $V_{1}$ .

  Formally, define the following transformation $\tau:\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}^{m}$ such that $\tau\left(\mathbf{x},\mathbf{h}\right)=T\mathbf{x}\circ V_{2}\mathbf{h}$ and then we can define $\mathcal{S}'_{\sigma}$ for some fixed learned transformation $\tau$ as $$\begin{aligned} \mathcal{S}'_{\sigma,\tau} & =\left\{ \mathbf{x},\mathbf{h}\mapsto\sigma\left(W\mathbf{x}+V_{1}\tau\left(\mathbf{x},\mathbf{h}\right)\right)\mid\forall W,V_{1}\right\} \\ & =\left\{ \mathbf{x},\mathbf{h}\mapsto\mathbf{s}\left(\mathbf{x},\tau\left(\mathbf{x},\mathbf{h}\right)\right)\mid\mathbf{s}\in\mathcal{S}_{\sigma}\right\} \end{aligned}$$ in other words, mLSTM is an LSTM that apply its gates to the tuple $\left(\mathbf{x}^{\left(t\right)},\tau\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)\right)$ rather than $\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)$ , and $\tau$ is another learned transformation. If you look again at the exact formula of $\tau$ you will see the multiplicative effect of the input vector on the transformation of $\mathbf{h}$ . That way, the authors said, $\mathcal{S}'_{\sigma,\tau}$ can yield much richer family of input-dependent transformations for $\mathbf{h}^{\left(t\right)}$ , which can be sensitive to the past as well as to subtle changes in the current input.

  3 Recurrent Batch Normalization

  Batch normalization is an operator applied to a layer before going through the non-linearity function in order to “normalize” the values before activation. That operator has two hyper-parameters for tuning and two statistics that it accumulates internally. Formally, given a vector of pre-activation values $\mathbf{x}$ , batch-normalization is $$BN\left(\mathbf{x};\gamma,\beta\right)=\beta+\gamma\circ\frac{\mathbf{x}-\mathbb{\hat{E}\left[\mathbf{x}\right]}}{\sqrt{\hat{\mathbb{V}}\left[\mathbf{x}\right]+\epsilon}}$$ where $\hat{\mathbb{E}}$ and $\hat{\mathbb{V}}$ are the empirical mean and variance of the current batch respectively. $\gamma,\beta$ and $\epsilon$ are vectors, and the division in the operator is being computed element-wise. Note that at inference time, the population statistics are estimated by averaging the empirical statistics across all the batches.

  That paper from February 2017, applied batch normalization also to the recurrent connections in LSTM layers, and they showed empirically that it helped the network to converge faster. The usage is simple, each gate function $\mathbf{i},\mathbf{f},\mathbf{o}\in\mathcal{S}_{sig}$ and the transition function $\mathbf{s}\in\mathcal{S}_{tanh}$ computes $$\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\mapsto\sigma\left(BN\left(W\mathbf{h}^{\left(t-1\right)};\gamma_{h},\beta_{h}\right)+BN\left(W\mathbf{x}^{\left(t\right)};\gamma_{x},\beta_{x}\right)+\mathbf{b}\right)$$ when the hyper-parameters $\gamma$ and $\beta$ are shared across different gates. Then the output of the layer is $\mathbf{h}^{\left(t\right)}=o^{\left(t\right)}\circ tanh\left(BN\left(\mathbf{c}^{\left(t\right)};\gamma_{c},\beta_{c}\right)\right)$ . The authors suggested to set $\beta_{x}=\beta_{h}=\mathbf{0}$ because there is already a bias parameter $\mathbf{b}$ and to prevent redundancy. In addition, they said that sharing the internal $BN$ statistics across time degrades performance severely. Therefore, one should use “fresh” $BN$ operators for each time-step with their own internal statistics (but share the $\beta$ and $\gamma$ parameters).

  4 The Ever-growing Field

  RNNs is a rapidly growing field and this series only covers a small part of it. There are much more advanced models, some of them are only from few months ago. If you are interested in this field, you may find very interesting stuff here.
• May 12, 2017

  Recurrent Neural Network - Dropout for LSTMs

  tags:
  TheoryLunch Blog

  Computer science topics over lunch

  Recurrent Neural Network - Dropout for LSTMs

  May 12, 2017 • Eran Amar

  tags: $\newcommand{\lyxlock}{}$

  1 Recap

  This is the third post in the series about Recurrent Neural Networks. The first post defined the basic notations and the formulation of a Recurrent Cell, the second post discussed its extensions such as LSTM and GRU. Recall that LSTM layer receives vector $\mathbf{x}^{\left(t\right)}$ - the input for the current time step $t$ , and $\mathbf{c}^{\left(t-1\right)},\mathbf{h}^{\left(t-1\right)}$ the memory-vector and the output-vector from previous time-step respectively, then the layer computes a new memory vector $$\mathbf{c}^{\left(t\right)}=\mathbf{f}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)\circ\mathbf{c}^{\left(t-1\right)}+\mathbf{i}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)\circ\mathbf{s}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)$$ where $\circ$ denotes element wise multiplication. Then the layer outputs $$\mathbf{h}^{\left(t\right)}=\mathbf{o}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)\circ tanh\left(\mathbf{c}^{\left(t\right)}\right)$$ where $\mathbf{f},\mathbf{i},\mathbf{o}\in\mathcal{S}_{sigmoid}$ are the forget, input and output gates respectively, and $\mathbf{s}\in\mathcal{S}_{tanh}$ is the state transition function. $\mathcal{S}_{\sigma}$ denotes a family of all affine transformations $\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}^{m}$ following by an element-wise activation function $\sigma$ , for any input-dimension $n$ and output-dimension $m$ .

  In this post, we will focus on the different variants of dropout regularization that are being used in LSTM networks.

  2 Dropout Regularization

  Dropout is a very popular regularization mechanism that is being attached to layers of a neural network in order to reduce overfitting and improve generalization. Applying dropout to a layer starts by fixing some $p\in\left(0,1\right)$ . Then, any time that layer produces an output during training, each one of its neurons is being zeroed-out independently with probability $p$ , and with probability $1-p$ it is being scaled by $\frac{1}{1-p}$ . Scaling the output is important because it keeps the expected output-value of the neuron unchanged. During test time, the dropout mechanism is being turned-off, that is, the output of each neuron is being passed as is.

  For fully-connected layers, we can formulate dropout as follows: Suppose that $\mathbf{h}$ should be the output of some layer that has dropout mechanism. Generate a mask $\mathbf{z}\in\left\{ 0,1\right\} ^{n}$ , by picking each coordinate i.i.d w.p. $p$ from $\left\{ 0,1\right\}$ , and output $\mathbf{h}\circ\mathbf{z}$ instead of $\mathbf{h}$ . Note that we hide here the output-scaling for simplicity, and that the mask $\mathbf{z}$ is being generated any time the layer required to produce an output.

  At first glance it seems trivial to apply dropout to LSTM layers. However, when we come to implement the mechanism we encounter some technical issues that should be addressed. For instance, we need to decide whether to apply the given mask only to the input $\mathbf{x}^{\left(t\right)}$ or also to the recurrent connection (i.e. to $\mathbf{h}^{\left(t-1\right)}$ ), and if we choose to apply on both, should it be the same mask or different masks? Should a mask be shared across time or be generated for each time-step?

  Several works tried to apply dropout to LSTM layers in several naive ways but without success. It seems that just dropping randomly some coordinates from the recurrent connections $\mathbf{h}^{\left(t-1\right)}$ impair the ability of the LSTM layer to learn long\short term dependencies and does not improve generalization. In the next section we will review some works that applied dropout to LSTMs in a way that successfully yields better generalization.

  3 Variants

  3.1 Mask Only Inputs; Regenerate Masks

  Wojciech Zaremba, Ilya Sutskever and Oriol Vinyals published in 2014 a paper that describe a successful dropout variation for LSTM layers. Their idea was that dropout should be applied only to the inputs of the layer and not to the recurrent connections. Moreover, a new mask should be generated for each time step.

  Formally, for each time-step $t$ , generate a mask $\mathbf{z}^{\left(t\right)}$ and compute $\mathbf{\hat{x}}^{\left(t\right)}:=\mathbf{x}^{\left(t\right)}\circ\mathbf{z}^{\left(t\right)}$ . Then continue to compute the LSTM layer as usual but use $\mathbf{\hat{x}}^{\left(t\right)}$ as the input to the layer rather than $\mathbf{x}^{\left(t\right)}$ .

  3.2 rnnDrop: Mask Only the Memory; Fixed Mask

  On 2015, Taesup Moon, Heeyoul Choi, Hoshik Lee and Inchul Song published a different dropout variation: rnnDrop.

  They suggested to generate a mask for each training sequence and fix it for all the time-steps in that sequence, that is, the mask is being shared across time. The mask is then being applied to the memory vector of the layer rather than the input. In their formulation, only the second formula of the LSTM changed: $$\mathbf{h}^{\left(t\right)}=\mathbf{o}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)\circ tanh\left(\mathbf{c}^{\left(t\right)}\circ\mathbf{z}\right)$$ where $\mathbf{z}$ if the fixed mask to the entire current training sequence.

  3.3 Mask Input and Hidden; Fixed Mask

  A relatively recent work of Yarin Gal and Zoubin Ghahramani from 2016 also use a make that is shared across time, however it is being applied to the inputs as well on the recurrent connections. This is one of the first successful dropout variants that actually apply the mask to the recurrent connection.

  Formally, for each training sequence generate two masks $\mathbf{z}_{1},\mathbf{z}_{2}$ and compute $\mathbf{\hat{x}}^{\left(t\right)}:=\mathbf{x}^{\left(t\right)}\circ\mathbf{z}_{1}$ and $\mathbf{\hat{h}}^{\left(t-1\right)}:=\mathbf{h}^{\left(t-1\right)}\circ\mathbf{z}_{2}$ . Then use $\mathbf{\hat{x}}^{\left(t\right)},\mathbf{\hat{h}}^{\left(t-1\right)}$ as the input to the “regular” LSTM layer.

  3.4 Mask Gates; Fixed Mask

  Another paper from 2016, of Stanislau Semeniuta, Aliaksei Severyn and Erhardt Barth demonstrate the mask being applied to some of the gates rather than the input\hidden vectors. For any time-step, generate $\mathbf{z}^{\left(t\right)}$ and then use it to mask the input gate $$\mathbf{c}^{\left(t\right)}=\mathbf{f}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)\circ\mathbf{c}^{\left(t-1\right)}+\mathbf{z}^{\left(t\right)}\circ\mathbf{i}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)\circ\mathbf{s}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)$$ and the second LSTM equation left unchanged.

  Small note, the authors also addressed in their paper some issues related to scaling the not-dropped coordinates, which won’t be covered here.

  3.5 Zoneout

  The very last dropout variation (to the day I wrote this post) is by David Krueger et al, from 2017. They suggested to treat the memory vector $\mathbf{c}^{\left(t\right)}$ and the output vector (a.k.a hidden vector) $\mathbf{h}^{\left(t\right)}$ as follows: each coordinate in each of the vectors either being updated as usual or preserved its value from previous time-step. As opposed to regular dropout, where “dropping a coordinate” meaning to make it zero, in zoneout it just keeps its previous value - acting as a random identity map that allows better gradients propagation through more time steps. Note that preserving a coordinate in the memory vector should not affect the computation of the hidden vector. Therefore we have to rewrite the formulas for the LSTM layer: given $\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}$ and $\mathbf{c}^{\left(t-1\right)}$ . Generate two masks $\mathbf{z}_{1}^{\left(t\right)},\mathbf{z}_{2}^{\left(t\right)}$ for the current time-step $t$ . Start by computing a candidate memory-vector with the regular formula, $$\mathbf{\hat{c}}^{\left(t\right)}=\mathbf{f}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)\circ\mathbf{c}^{\left(t-1\right)}+\mathbf{i}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)\circ\mathbf{s}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)$$ and then use only some of its coordinates to update the memory-vector, that is, $$\mathbf{c}^{\left(t\right)}=\mathbf{z}_{1}^{\left(t\right)}\circ\mathbf{c}^{\left(t-1\right)}+\left(1-\mathbf{z}_{1}^{\left(t\right)}\right)\circ\mathbf{\hat{c}}^{\left(t\right)}$$ Similarly, for the hidden-vector, compute a candidate $$\mathbf{\hat{h}}^{\left(t\right)}=\mathbf{o}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)\circ tanh\left(\mathbf{\hat{c}}^{\left(t\right)}\right)$$ and selectively update $$\mathbf{h}^{\left(t\right)}=\mathbf{z}_{2}^{\left(t\right)}\circ\mathbf{h}^{\left(t-1\right)}+\left(1-\mathbf{z}_{2}^{\left(t\right)}\right)\circ\mathbf{\hat{h}}^{\left(t\right)}$$ Observe again that the computation of the candidate hidden-vector is based on $\mathbf{\hat{c}}^{\left(t\right)}$ therefore unaffected by the mask $\mathbf{z}_{1}^{\left(t\right)}$ .

  4 More and More and More

  In the last year, dropout for LSTM networks gained more attention, and the list of different variants is growing quickly. This post tried to cover only some fraction of all the variants available out there. It is left for the curious reader to search the literature for more about this topic.

  Note that some of the dropout variations discussed above can be applied to basic RNN and GRU cells without much modifications - please refer to the papers themselves for more details.

  Next post, I will discuss recent advancement in the field of RNNs, such as Multiplicative LSTMs, Recurrent Batch Normalization and more.

  $\newcommand{\lyxlock}{}$

  1 Recap

  This is the third post in the series about Recurrent Neural Networks. The first post defined the basic notations and the formulation of a Recurrent Cell, the second post discussed its extensions such as LSTM and GRU. Recall that LSTM layer receives vector $\mathbf{x}^{\left(t\right)}$ - the input for the current time step $t$ , and $\mathbf{c}^{\left(t-1\right)},\mathbf{h}^{\left(t-1\right)}$ the memory-vector and the output-vector from previous time-step respectively, then the layer computes a new memory vector $$\mathbf{c}^{\left(t\right)}=\mathbf{f}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)\circ\mathbf{c}^{\left(t-1\right)}+\mathbf{i}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)\circ\mathbf{s}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)$$ where $\circ$ denotes element wise multiplication. Then the layer outputs $$\mathbf{h}^{\left(t\right)}=\mathbf{o}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)\circ tanh\left(\mathbf{c}^{\left(t\right)}\right)$$ where $\mathbf{f},\mathbf{i},\mathbf{o}\in\mathcal{S}_{sigmoid}$ are the forget, input and output gates respectively, and $\mathbf{s}\in\mathcal{S}_{tanh}$ is the state transition function. $\mathcal{S}_{\sigma}$ denotes a family of all affine transformations $\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}^{m}$ following by an element-wise activation function $\sigma$ , for any input-dimension $n$ and output-dimension $m$ .

  In this post, we will focus on the different variants of dropout regularization that are being used in LSTM networks.

  2 Dropout Regularization

  Dropout is a very popular regularization mechanism that is being attached to layers of a neural network in order to reduce overfitting and improve generalization. Applying dropout to a layer starts by fixing some $p\in\left(0,1\right)$ . Then, any time that layer produces an output during training, each one of its neurons is being zeroed-out independently with probability $p$ , and with probability $1-p$ it is being scaled by $\frac{1}{1-p}$ . Scaling the output is important because it keeps the expected output-value of the neuron unchanged. During test time, the dropout mechanism is being turned-off, that is, the output of each neuron is being passed as is.

  For fully-connected layers, we can formulate dropout as follows: Suppose that $\mathbf{h}$ should be the output of some layer that has dropout mechanism. Generate a mask $\mathbf{z}\in\left\{ 0,1\right\} ^{n}$ , by picking each coordinate i.i.d w.p. $p$ from $\left\{ 0,1\right\}$ , and output $\mathbf{h}\circ\mathbf{z}$ instead of $\mathbf{h}$ . Note that we hide here the output-scaling for simplicity, and that the mask $\mathbf{z}$ is being generated any time the layer required to produce an output.

  At first glance it seems trivial to apply dropout to LSTM layers. However, when we come to implement the mechanism we encounter some technical issues that should be addressed. For instance, we need to decide whether to apply the given mask only to the input $\mathbf{x}^{\left(t\right)}$ or also to the recurrent connection (i.e. to $\mathbf{h}^{\left(t-1\right)}$ ), and if we choose to apply on both, should it be the same mask or different masks? Should a mask be shared across time or be generated for each time-step?

  Several works tried to apply dropout to LSTM layers in several naive ways but without success. It seems that just dropping randomly some coordinates from the recurrent connections $\mathbf{h}^{\left(t-1\right)}$ impair the ability of the LSTM layer to learn long\short term dependencies and does not improve generalization. In the next section we will review some works that applied dropout to LSTMs in a way that successfully yields better generalization.

  3 Variants

  3.1 Mask Only Inputs; Regenerate Masks

  Wojciech Zaremba, Ilya Sutskever and Oriol Vinyals published in 2014 a paper that describe a successful dropout variation for LSTM layers. Their idea was that dropout should be applied only to the inputs of the layer and not to the recurrent connections. Moreover, a new mask should be generated for each time step.

  Formally, for each time-step $t$ , generate a mask $\mathbf{z}^{\left(t\right)}$ and compute $\mathbf{\hat{x}}^{\left(t\right)}:=\mathbf{x}^{\left(t\right)}\circ\mathbf{z}^{\left(t\right)}$ . Then continue to compute the LSTM layer as usual but use $\mathbf{\hat{x}}^{\left(t\right)}$ as the input to the layer rather than $\mathbf{x}^{\left(t\right)}$ .

  3.2 rnnDrop: Mask Only the Memory; Fixed Mask

  On 2015, Taesup Moon, Heeyoul Choi, Hoshik Lee and Inchul Song published a different dropout variation: rnnDrop.

  They suggested to generate a mask for each training sequence and fix it for all the time-steps in that sequence, that is, the mask is being shared across time. The mask is then being applied to the memory vector of the layer rather than the input. In their formulation, only the second formula of the LSTM changed: $$\mathbf{h}^{\left(t\right)}=\mathbf{o}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)\circ tanh\left(\mathbf{c}^{\left(t\right)}\circ\mathbf{z}\right)$$ where $\mathbf{z}$ if the fixed mask to the entire current training sequence.

  3.3 Mask Input and Hidden; Fixed Mask

  A relatively recent work of Yarin Gal and Zoubin Ghahramani from 2016 also use a make that is shared across time, however it is being applied to the inputs as well on the recurrent connections. This is one of the first successful dropout variants that actually apply the mask to the recurrent connection.

  Formally, for each training sequence generate two masks $\mathbf{z}_{1},\mathbf{z}_{2}$ and compute $\mathbf{\hat{x}}^{\left(t\right)}:=\mathbf{x}^{\left(t\right)}\circ\mathbf{z}_{1}$ and $\mathbf{\hat{h}}^{\left(t-1\right)}:=\mathbf{h}^{\left(t-1\right)}\circ\mathbf{z}_{2}$ . Then use $\mathbf{\hat{x}}^{\left(t\right)},\mathbf{\hat{h}}^{\left(t-1\right)}$ as the input to the “regular” LSTM layer.

  3.4 Mask Gates; Fixed Mask

  Another paper from 2016, of Stanislau Semeniuta, Aliaksei Severyn and Erhardt Barth demonstrate the mask being applied to some of the gates rather than the input\hidden vectors. For any time-step, generate $\mathbf{z}^{\left(t\right)}$ and then use it to mask the input gate $$\mathbf{c}^{\left(t\right)}=\mathbf{f}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)\circ\mathbf{c}^{\left(t-1\right)}+\mathbf{z}^{\left(t\right)}\circ\mathbf{i}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)\circ\mathbf{s}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)$$ and the second LSTM equation left unchanged.

  Small note, the authors also addressed in their paper some issues related to scaling the not-dropped coordinates, which won’t be covered here.

  3.5 Zoneout

  The very last dropout variation (to the day I wrote this post) is by David Krueger et al, from 2017. They suggested to treat the memory vector $\mathbf{c}^{\left(t\right)}$ and the output vector (a.k.a hidden vector) $\mathbf{h}^{\left(t\right)}$ as follows: each coordinate in each of the vectors either being updated as usual or preserved its value from previous time-step. As opposed to regular dropout, where “dropping a coordinate” meaning to make it zero, in zoneout it just keeps its previous value - acting as a random identity map that allows better gradients propagation through more time steps. Note that preserving a coordinate in the memory vector should not affect the computation of the hidden vector. Therefore we have to rewrite the formulas for the LSTM layer: given $\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}$ and $\mathbf{c}^{\left(t-1\right)}$ . Generate two masks $\mathbf{z}_{1}^{\left(t\right)},\mathbf{z}_{2}^{\left(t\right)}$ for the current time-step $t$ . Start by computing a candidate memory-vector with the regular formula, $$\mathbf{\hat{c}}^{\left(t\right)}=\mathbf{f}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)\circ\mathbf{c}^{\left(t-1\right)}+\mathbf{i}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)\circ\mathbf{s}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)$$ and then use only some of its coordinates to update the memory-vector, that is, $$\mathbf{c}^{\left(t\right)}=\mathbf{z}_{1}^{\left(t\right)}\circ\mathbf{c}^{\left(t-1\right)}+\left(1-\mathbf{z}_{1}^{\left(t\right)}\right)\circ\mathbf{\hat{c}}^{\left(t\right)}$$ Similarly, for the hidden-vector, compute a candidate $$\mathbf{\hat{h}}^{\left(t\right)}=\mathbf{o}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)\circ tanh\left(\mathbf{\hat{c}}^{\left(t\right)}\right)$$ and selectively update $$\mathbf{h}^{\left(t\right)}=\mathbf{z}_{2}^{\left(t\right)}\circ\mathbf{h}^{\left(t-1\right)}+\left(1-\mathbf{z}_{2}^{\left(t\right)}\right)\circ\mathbf{\hat{h}}^{\left(t\right)}$$ Observe again that the computation of the candidate hidden-vector is based on $\mathbf{\hat{c}}^{\left(t\right)}$ therefore unaffected by the mask $\mathbf{z}_{1}^{\left(t\right)}$ .

  4 More and More and More

  In the last year, dropout for LSTM networks gained more attention, and the list of different variants is growing quickly. This post tried to cover only some fraction of all the variants available out there. It is left for the curious reader to search the literature for more about this topic.

  Note that some of the dropout variations discussed above can be applied to basic RNN and GRU cells without much modifications - please refer to the papers themselves for more details.

  Next post, I will discuss recent advancement in the field of RNNs, such as Multiplicative LSTMs, Recurrent Batch Normalization and more.
• Apr 30, 2017

  Recurrent Neural Network - LSTM and GRU

  tags:
  TheoryLunch Blog

  Computer science topics over lunch

  Recurrent Neural Network - LSTM and GRU

  Apr 30, 2017 • Eran Amar

  tags: $\newcommand{\lyxlock}{}$

  1 Recap

  The first post in the series discussed the basic structure of recurrent cells and their limitations. We defined two families of functions, the first is $\mathcal{F}_{\sigma}$ which contains all the affine transformations $\mathbb{R}^{n}\to\mathbb{R}^{m}$ for any $n$ and $m,$ followed by an element-wise activation function $\sigma.$ And another family $\mathcal{S}_{\sigma}$ which is some kind of extension of $\mathcal{F}_{\sigma}$ in the sense that the input space is $\mathbb{R}^{n}\times\mathbb{R}^{n}$ rather than $\mathbb{R}^{n}.$ Formally, the definition of $\mathcal{S}_{\sigma}$ is $$\mathcal{S}_{\sigma}=\left\{ \mathbf{x},\mathbf{h}\mapsto\sigma\left(W\mathbf{x}+U\mathbf{h}+\mathbf{b}\right)\mid\forall W,U,\mathbf{b}\right\}$$ In this post, which is the second in the series, we will focus on two activation functions: $tanh\left(x\right)=\frac{e^{x}-e^{-x}}{e^{x}+e^{-x}}$ and $sig\left(x\right)=\left(1+e^{-x}\right)^{-1}$ . Next section describes what are gate functions, then we move to review two common improvements to the basic recurrent cell: the LSTM and GRU cells.

  2 Gate Functions

  Recall that recurrent layer has two steps. The first is updating its inner state based on both the current input and the previous state vector, then producing an output by applying some other function to the new state. So the input to the layer at each time step $t$ is actually the tuple $\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)$ .

  A gate is a function that takes such a tuple and produces a vector of values between zero and one. Note that any function in $\mathcal{S}_{sig}$ and in $\mathcal{S}_{tanh}$ is actually a gate because of the properties of the $tanh$ and $sig$ functions. We add gates into a recurrent neuron in order to control how much information will flow over the “recurrent connection”. That is, suppose that $\mathbf{h}^{\left(t\right)}$ is the layer’s current state that should be used at the next time step, a gated layer will use the vector $\hat{\mathbf{h}}^{\left(t\right)}=\mathbf{h}^{\left(t\right)}\circ\mathbf{s}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right),$ instead of using $\mathbf{h}^{\left(t\right)}$ (when the symbol $\circ$ denotes element-wise multiplication). Note that any coordinate in $\hat{\mathbf{h}}^{\left(t\right)}$ is a “moderated” version of the corresponding coordinate in $\mathbf{h}^{\left(t\right)}$ because each entry in the gate’s output is in $\left(0,1\right).$ We will see a concrete example in the next section.

  3 LSTM Cell

  Long Short Term Memory Cell (LSTM cell) is an improved version of the recurrent neuron that was proposed on 1997. It went through minor modifications until the version presented below (which is from 2013). LSTM cells solve both limitations of the basic recurrent neuron; it prevents the exploding and vanishing gradients problem, and it can remember as well as it can forget.

  The main addition to the recurrent layer structure is the use of gates and a memory vector for each time step, denoted by $\mathbf{c}^{\left(t\right)}$ . The LSTM layer gets as inputs the tuple $\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)$ and the previous memory vector $\mathbf{c}^{\left(t-1\right)},$ then outputs $\mathbf{h}^{\left(t\right)}$ and an updated memory vector $\mathbf{c}^{\left(t\right)}.$

  Here is how an LSTM layer computes its outputs: it have the following gates $\mathbf{f},\mathbf{i},\mathbf{o}\in\mathcal{S}_{sig}$ named the forget, input and output gate respectively, and a state-transition function $\mathbf{s}\in\mathcal{S}_{tanh}$ . It first updates the memory vector $$\mathbf{c}^{\left(t\right)}=\mathbf{f}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)\circ\mathbf{c}^{\left(t-1\right)}+\mathbf{i}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)\circ\mathbf{s}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)$$ and then computes $$\mathbf{h}^{\left(t\right)}=\mathbf{o}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)\circ tanh\left(\mathbf{c}^{\left(t\right)}\right)$$ Those equations can be explained as follows: the first equation is an element-wise summation of two terms. The first term is the previous memory vector $\mathbf{c}^{\left(t-1\right)}$ moderated by the forget gate. That is, the layer uses the current input $\mathbf{x}^{\left(t\right)}$ and previous output $\mathbf{h}^{\left(t-1\right)}$ to determine how much to shrink each coordinate of the previous memory vector. The second term is the candidate for the new state, i.e. $\mathbf{s}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right),$ moderated by the input gate. Note that all the gates operate on the same input tuple $\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right).$

  The input and forget gates control the long-short term dependencies (i.e. the recurrent connection), and allow the LSTM layer to adaptively control the balance between new information that come from the state-transition function and the history information that comes from the memory vector, hence the names for the gates: input and forget.

  Another difference is that LSTM layer controls how much of its inner-memory to expose by using the output gate. That is being formulated in the second equation.

  The addition of gates is what preventing the exploding and vanishing gradients problem. It makes the LSTM layer able to learn long and short term dependencies at the cost of increasing the number of parameters that are needed to be trained, and that makes the network harder to optimize.

  4 GRU

  Another gated cell proposed on 2014 is the Gated Recurrent Unit (GRU). It has similar advantages as the LSTM cell, but fewer parameters to train because the memory vector was removed and one of the gates.

  GRU has two gate functions $\mathbf{z},\mathbf{r}\in\mathcal{S}_{sig}$ named the update and reset gate respectively, and a state transition $\mathbf{s}\in\mathcal{S}_{tanh}.$ The input to GRU layer is only the tuple $\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)$ and the output is $\mathbf{h}^{\left(t\right)}$ computed as follow: first the layer computes its gates for the current time step, denoted by $$\mathbf{z}^{\left(t\right)}:=\mathbf{z}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)$$ for the update gate and by $$\mathbf{r}^{\left(t\right)}:=\mathbf{r}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)$$ for the reset gate. Then the output of the layer is $$\mathbf{h}^{\left(t\right)}=\left(1-\mathbf{z}^{\left(t\right)}\right)\circ\mathbf{h}^{\left(t-1\right)}+\mathbf{z}^{\left(t\right)}\circ\mathbf{s}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\circ\mathbf{r}^{\left(t\right)}\right)$$ with all the arithmetic operations being done element wise.

  The term $\mathbf{s}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\circ\mathbf{r}^{\left(t\right)}\right)$ is a candidate for the next state. Note that the state-transition function accepts the previous state moderated by the reset gate, allowing it to forget past states (therefore the name of the gate: reset). Then the output of the layer is a linear interpolation between the previous state and the candidate state, controlled by the update gate.

  As oppose to the LSTM cell, the GRU doesn’t have an output gate to control how much of its inner state to expose therefore the entire state is being exposed at each time step.

  Both LSTM and GRU are very common in many recurrent network architectures and achieve great results in many tasks. LSTMs and GRUs can learn dependencies of various lengths which make the network very expressive. However, a too expressive network can leads sometimes to overfitting, to prevent that it is common to use some type of regularization, such as Dropout.

  Next post I will discuss the dropout variations that are specific for RNNs.

  $\newcommand{\lyxlock}{}$

  1 Recap

  The first post in the series discussed the basic structure of recurrent cells and their limitations. We defined two families of functions, the first is $\mathcal{F}_{\sigma}$ which contains all the affine transformations $\mathbb{R}^{n}\to\mathbb{R}^{m}$ for any $n$ and $m,$ followed by an element-wise activation function $\sigma.$ And another family $\mathcal{S}_{\sigma}$ which is some kind of extension of $\mathcal{F}_{\sigma}$ in the sense that the input space is $\mathbb{R}^{n}\times\mathbb{R}^{n}$ rather than $\mathbb{R}^{n}.$ Formally, the definition of $\mathcal{S}_{\sigma}$ is $$\mathcal{S}_{\sigma}=\left\{ \mathbf{x},\mathbf{h}\mapsto\sigma\left(W\mathbf{x}+U\mathbf{h}+\mathbf{b}\right)\mid\forall W,U,\mathbf{b}\right\}$$ In this post, which is the second in the series, we will focus on two activation functions: $tanh\left(x\right)=\frac{e^{x}-e^{-x}}{e^{x}+e^{-x}}$ and $sig\left(x\right)=\left(1+e^{-x}\right)^{-1}$ . Next section describes what are gate functions, then we move to review two common improvements to the basic recurrent cell: the LSTM and GRU cells.

  2 Gate Functions

  Recall that recurrent layer has two steps. The first is updating its inner state based on both the current input and the previous state vector, then producing an output by applying some other function to the new state. So the input to the layer at each time step $t$ is actually the tuple $\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)$ .

  A gate is a function that takes such a tuple and produces a vector of values between zero and one. Note that any function in $\mathcal{S}_{sig}$ and in $\mathcal{S}_{tanh}$ is actually a gate because of the properties of the $tanh$ and $sig$ functions. We add gates into a recurrent neuron in order to control how much information will flow over the “recurrent connection”. That is, suppose that $\mathbf{h}^{\left(t\right)}$ is the layer’s current state that should be used at the next time step, a gated layer will use the vector $\hat{\mathbf{h}}^{\left(t\right)}=\mathbf{h}^{\left(t\right)}\circ\mathbf{s}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right),$ instead of using $\mathbf{h}^{\left(t\right)}$ (when the symbol $\circ$ denotes element-wise multiplication). Note that any coordinate in $\hat{\mathbf{h}}^{\left(t\right)}$ is a “moderated” version of the corresponding coordinate in $\mathbf{h}^{\left(t\right)}$ because each entry in the gate’s output is in $\left(0,1\right).$ We will see a concrete example in the next section.

  3 LSTM Cell

  Long Short Term Memory Cell (LSTM cell) is an improved version of the recurrent neuron that was proposed on 1997. It went through minor modifications until the version presented below (which is from 2013). LSTM cells solve both limitations of the basic recurrent neuron; it prevents the exploding and vanishing gradients problem, and it can remember as well as it can forget.

  The main addition to the recurrent layer structure is the use of gates and a memory vector for each time step, denoted by $\mathbf{c}^{\left(t\right)}$ . The LSTM layer gets as inputs the tuple $\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)$ and the previous memory vector $\mathbf{c}^{\left(t-1\right)},$ then outputs $\mathbf{h}^{\left(t\right)}$ and an updated memory vector $\mathbf{c}^{\left(t\right)}.$

  Here is how an LSTM layer computes its outputs: it have the following gates $\mathbf{f},\mathbf{i},\mathbf{o}\in\mathcal{S}_{sig}$ named the forget, input and output gate respectively, and a state-transition function $\mathbf{s}\in\mathcal{S}_{tanh}$ . It first updates the memory vector $$\mathbf{c}^{\left(t\right)}=\mathbf{f}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)\circ\mathbf{c}^{\left(t-1\right)}+\mathbf{i}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)\circ\mathbf{s}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)$$ and then computes $$\mathbf{h}^{\left(t\right)}=\mathbf{o}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)\circ tanh\left(\mathbf{c}^{\left(t\right)}\right)$$ Those equations can be explained as follows: the first equation is an element-wise summation of two terms. The first term is the previous memory vector $\mathbf{c}^{\left(t-1\right)}$ moderated by the forget gate. That is, the layer uses the current input $\mathbf{x}^{\left(t\right)}$ and previous output $\mathbf{h}^{\left(t-1\right)}$ to determine how much to shrink each coordinate of the previous memory vector. The second term is the candidate for the new state, i.e. $\mathbf{s}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right),$ moderated by the input gate. Note that all the gates operate on the same input tuple $\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right).$

  The input and forget gates control the long-short term dependencies (i.e. the recurrent connection), and allow the LSTM layer to adaptively control the balance between new information that come from the state-transition function and the history information that comes from the memory vector, hence the names for the gates: input and forget.

  Another difference is that LSTM layer controls how much of its inner-memory to expose by using the output gate. That is being formulated in the second equation.

  The addition of gates is what preventing the exploding and vanishing gradients problem. It makes the LSTM layer able to learn long and short term dependencies at the cost of increasing the number of parameters that are needed to be trained, and that makes the network harder to optimize.

  4 GRU

  Another gated cell proposed on 2014 is the Gated Recurrent Unit (GRU). It has similar advantages as the LSTM cell, but fewer parameters to train because the memory vector was removed and one of the gates.

  GRU has two gate functions $\mathbf{z},\mathbf{r}\in\mathcal{S}_{sig}$ named the update and reset gate respectively, and a state transition $\mathbf{s}\in\mathcal{S}_{tanh}.$ The input to GRU layer is only the tuple $\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)$ and the output is $\mathbf{h}^{\left(t\right)}$ computed as follow: first the layer computes its gates for the current time step, denoted by $$\mathbf{z}^{\left(t\right)}:=\mathbf{z}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)$$ for the update gate and by $$\mathbf{r}^{\left(t\right)}:=\mathbf{r}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\right)$$ for the reset gate. Then the output of the layer is $$\mathbf{h}^{\left(t\right)}=\left(1-\mathbf{z}^{\left(t\right)}\right)\circ\mathbf{h}^{\left(t-1\right)}+\mathbf{z}^{\left(t\right)}\circ\mathbf{s}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\circ\mathbf{r}^{\left(t\right)}\right)$$ with all the arithmetic operations being done element wise.

  The term $\mathbf{s}\left(\mathbf{x}^{\left(t\right)},\mathbf{h}^{\left(t-1\right)}\circ\mathbf{r}^{\left(t\right)}\right)$ is a candidate for the next state. Note that the state-transition function accepts the previous state moderated by the reset gate, allowing it to forget past states (therefore the name of the gate: reset). Then the output of the layer is a linear interpolation between the previous state and the candidate state, controlled by the update gate.

  As oppose to the LSTM cell, the GRU doesn’t have an output gate to control how much of its inner state to expose therefore the entire state is being exposed at each time step.

  Both LSTM and GRU are very common in many recurrent network architectures and achieve great results in many tasks. LSTMs and GRUs can learn dependencies of various lengths which make the network very expressive. However, a too expressive network can leads sometimes to overfitting, to prevent that it is common to use some type of regularization, such as Dropout.

  Next post I will discuss the dropout variations that are specific for RNNs.
• Apr 20, 2017

  Recurrent Neural Network - Introduction

  tags:
  TheoryLunch Blog

  Computer science topics over lunch

  Recurrent Neural Network - Introduction

  Apr 20, 2017 • Eran Amar

  tags: $\newcommand{\lyxlock}{}$

  1 The First Post in the Series

  Recurrent Neural Network (RNN) is a neural network architecture that leads to state of the art results in many tasks, such as machine translation, speech recognition, image captioning and many more. Because it is a rapidly evolving field, it is difficult to find all the information (both basic and advanced material) in one place, so I decided to write a series of posts that will discuss Recurrent Neural Networks from their basic formulation to their most recent variations.

  This post, as the first in the series, will present the concept of a recurrent neuron and its limitations. Next posts will discuss gated neurons like LSTM and GRU, variants of dropout regularization for RNNs and some advanced stuff that I found interesting (bi-directional RNN, Multiplicative LSTM and etc).

  The next section is dedicated to review the notations and basic terminology that will be used during the entire series. I will assume that the reader is familiar with the basics of neural network (fully connected architecture, back-propagation algorithm, etc). If you are not familiar with one of those concepts, I strongly recommend to review them online before going further. Those notes provides concise introduction to the field. There is also an excellent video-lecture (in Hebrew) by Prof. Shai Shalev-Shwartz from 2016. If you really want to dive deeper, you can find broad introduction in Chapter 6 in Deep Learning book (legally free electronic version here).

  2 Notations and Basic Terminology

  Artificial Neural Network (ANN) is a mathematical hierarchal model inspired from neuroscience. The basic computation unit is a neuron, which is a function from $\mathbb{R}^{d}$ to $\mathbb{R}$ , $$g:\mathbf{x}\mapsto\sigma\left(\left\langle \mathbf{x},\mathbf{w}\right\rangle +b\right)$$ where $\mathbf{x},\mathbf{w}\in\mathbb{R}^{d}$ are the input and weight vectors respectively, and $b\in\mathbb{R}$ is the bias. The function $\sigma$ is a non-linear activation function $\mathbb{R}\to\mathbb{R}$ , for example $\sigma\left(z\right)=ReLU\left(z\right):=\max\left(z,0\right)$ or $\sigma\left(z\right)=tanh\left(z\right):=\frac{e^{z}-e^{-z}}{e^{z}+e^{-z}}$ . Note that a neuron can ignore some of its input coordinates by setting their corresponding entries in its weight vector to zero.

  Neurons are organized within layers. In the very basic ANN architecture, called a fully connected (FC) neural network, a layer is just a collection of neurons that will be computed in parallel (meaning that they are not connected to each other). For any layer $i$ , each of its neurons accepts as inputs all neurons from layer $i-1$ .

  A layer in FC architecture can be formulate mathematically as a function from $\mathbb{R}^{d}$ to $\mathbb{R}^{n}$ such that $$\mathbf{f}:\mathbf{x}\mapsto\left(\begin{array}{c} g_{1}\left(\mathbf{x}\right)\\ \vdots\\ g_{n}\left(\mathbf{x}\right) \end{array}\right)$$ where each $g_{i}$ is a function of a neuron as defined above, with its own weight vector and a bias scalar. The formulation can be simplified using matrix notations: $$\mathbf{f}:\mathbf{x}\mapsto\sigma\left(W\mathbf{x}+\mathbf{b}\right)$$ where $W\in\mathbb{R}^{n\times d}$ contains in the i-th row the weight vector of $g_{i}$ , and $\mathbf{b}_{i}$ is the bias of $g_{i}$ . The function $\sigma$ is the same as before but since its input is now a vector, it apply element-wise on each coordinate of the input, producing an output vector of the same length.

  We denote the family of all possible FC functions based on some fixed activation function as $$\mathcal{F}_{\sigma}=\left\{ \mathbf{x}\mapsto\sigma\left(W\mathbf{x}+\mathbf{b}\right)\mid\forall W,\mathbf{b}\right\}$$ that family contains all possible layers for any input\output dimension. For abbreviation, I will not mention the dimensions of layers or the specific activation function $\sigma$ when they are not important for the discussion.

  Applying several layers sequentially yield a feedforward FC neural network. That is, given $\mathbf{x}\in\mathbb{R}^{d}$ , a FC network with $t$ layers is $$\mathcal{N}\left(\mathbf{x}\right)=\mathbf{f}_{t}\left(...\mathbf{f}_{2}\left(\mathbf{f}_{1}\left(\mathbf{x}\right)\right)\right)$$ where $\forall i:\:\mathbf{f}_{i}\in\mathcal{F}_{\sigma}$ for some fixed $\sigma$ .

  The number of layers in the network is its depth, and the size of the layers is its width (assuming all the layers are of the same dimension). Training the network is optimizing all the weight matrices and bias vectors such that the network will approximate some target (usually unknown) function. We can think of the optimization process as a process of iteratively updating the choice of layer functions $\mathbf{f}_{1},\mathbf{f}_{2},..$ from the family $\mathcal{F}_{\sigma}$ .

  The FC architecture is simple, yet yields a very powerful model. When wide enough, even a single hidden layer network can approximate to any arbitrary precision almost any continuous function (under some mild assumptions). However, for some tasks it might not be optimal, and different architectures are being used, instead of the FC one.

  3 The Need for Context

  There are tasks that require some notion of context in order to make a decent prediction. For example, tasks that involve time-steps can consider the history until the current time step as a context (a.k.a time-series tasks). For concrete example, consider predicting the next word in a given sequence. Each word has as a context all words preceding it. Note that there are more ways to define the notion of a context, and it mostly depends on the task in hand. For didactic purposes we will focus on the example of next-word prediction.

  How can we build a model for next-word prediction using FC neural network?

  The obvious answer is to fix some window size $d$ , and build a network that gets the last $d$ words in the sentence as features, in order to predict the next word. One immediate limitation of this approach is that the window size is the same for all the examples in the training set, and it is not simple to find a suitable window size. On the one hand, choosing small window may cause the model to miss some important long term relations between words that are more than $d$ steps apart (for sentences longer than $d$ ). On the other hand, choosing very large $d$ increases the number of parameters in the model, which make the optimization harder, and increase the computation time per example, which affect both training-time and prediction-time.

  Another limitation is that FC layers regards their inputs as ordered information. That is, a FC layer assigns a specific weight for each coordinate in its input vector, causing the order of the features to be meaningful. But in some tasks the order of the features may vary, for example, the subject of a sentence doesn’t limit to specific location inside an English sentence. If we want FC layer to account for that, we have to get a lot of examples (and they need to cover all different cases).

  What we would like to have here is called parameter sharing. We will come back to that term later in this post.

  4 Recurrent Cell - Learning to Remember

  We want our network to be able to remember information about previous invocations (i.e. the past inputs). In order to do that, our network have to maintain a state. So we need to extend our neuron formulation to include another variable, $h\in\mathbb{R}$ , that will be its state and be kept during consecutive invocations.That state variable is usually initialized to be zero so that it doesn’t affect the first computation.

  We will use an upper-script to denote the time step, so $h^{\left(t\right)}$ is the state of a neuron at time $t$ . Since the neurons need to pass information to themselves (i.e. their state from previous invocation) they become recurrent neurons, a.k.a recurrent cells.

  Recurrent neuron uses internally two functions. The first function updates the inner state based on current input and previous state, which is the transition from time $t-1$ to time $t$ , $$s:\mathbf{x},h^{\left(t-1\right)}\mapsto tanh\left(\left\langle \mathbf{x},\mathbf{w}\right\rangle +b+h^{\left(t-1\right)}\cdot u\right)$$ $u\in\mathbb{R}$ is a new parameter denoting the weight for the previous state.

  And the second function is for computing the output of the neuron at current time $t$ , based on its new state $$g^{\left(t\right)}:h^{\left(t\right)}\mapsto\sigma\left(h^{\left(t\right)}\cdot v+c\right)$$ with $v,c\in\mathbb{R}$ as the weight and bias of the computation respectively.

  We can move again to matrix notations to denote several parallel recurrent cells that forms a recurrent layer. The first function updates the vector of states for the entire layer $$\mathbf{s}:\mathbf{x},\mathbf{h}^{\left(t-1\right)}\mapsto tanh\left(W\mathbf{x}+\mathbf{b}+U\mathbf{h}^{\left(t-1\right)}\right)$$ and the second function for computing the output of the layer based on the current state $$\mathbf{f}:\mathbf{h}^{\left(t\right)}\mapsto\sigma\left(V\mathbf{h}^{\left(t\right)}+\mathbf{c}\right)$$ When we come to compare that to the FC formulation, we see that actually the function $\mathbf{f}$ itself hasn’t changed, but now it accepts as input the state of the layer $\mathbf{h}$ , rather than the input to the layer $\mathbf{x}$ . Another change is that before the layer computes its output, it first update its hidden state.

  Note that we can extend the definition of the family $\mathcal{F}_{tanh}$ to contains functions like the state-transforming function $\mathbf{s}$ , by making each $\mathbf{f}\in\mathcal{F}_{tanh}$ accepts two vectors as input instead of one (and we also need another weight matrix for the second input). So denote the family of all hidden-state transition functions as $$\mathcal{S}_{\sigma}=\left\{ \mathbf{x},\mathbf{h}^{\left(t-1\right)}\mapsto\sigma\left(W\mathbf{x}+\mathbf{b}+U\mathbf{h}^{\left(t-1\right)}\right)\mid\forall W,U,\mathbf{b}\right\}$$ usually the activation function is $\sigma=tanh$ . We will use that family later on, when discuss gated units as LSTM and GRU.

  5 Comparing the Architectures - an Example

  Consider again the next-word prediction problem. Suppose we have an ordered dictionary that contains the entire vocabulary we will ever need. We will represent any word by an integer indicating its index in the dictionary. In this section we will construct two single-neuron networks, once using the FC architecture and once using the Recurrent architecture.

  FC network for next-word prediction: Let the window size be $d=3$ . The input to the network can be represent as a vector $\mathbf{x}\in\mathbb{N}^{3}$ and the output as $y\in\mathbb{N}$ . The FC network amounts to $$\mathcal{N}_{FC}\left(\mathbf{x}\right)=f\left(x\right)$$ for some $f\in\mathcal{F}_{\sigma}$ that maps from $\mathbb{R}^{3}$ to $\mathbb{R}$ .

  Recurrent neural network (RNN), however, leads to different representation: the input and the output both integers $x,y\in\mathbb{N},$ so the state transition function $s\in\mathcal{S}$ and the output function $f\in\mathcal{F}_{\sigma}$ are both from $\mathbb{R}$ to $\mathbb{R}$ . The idea of the recurrent architecture is to apply the same recurrent layer several times on different inputs (each input from different time-step) in order to get a prediction.

  Here is a detailed description of how to do that: we take a series of three consecutive words from the beginning of a sentence: $x^{\left(1\right)},x^{\left(2\right)},x^{\left(3\right)}$ and a target next word: $y$ . First, feed the network with $x^{\left(1\right)}$ , it updates the inner state from $h^{\left(0\right)}=0$ to some hidden state $h^{\left(1\right)}$ and outputs $\hat{y}^{\left(1\right)}$ . We ignore that output and feed the next word $x^{\left(2\right)}$ , the network updates its state to $h^{\left(2\right)}$ and produces $\hat{y}^{\left(2\right)}$ that is also ignored. Lastly, we feed $x^{\left(3\right)}$ to the network, the inner state is being updated to $h^{\left(3\right)}$ and another output is produced $\hat{y}^{\left(3\right)}$ which is now considered as the prediction of the network. Note that in each invocation of the network, the output depends on both the current input and the previous state, therefore $\hat{y}^{\left(3\right)}$ is the result of the entire history $x^{\left(1\right)},x^{\left(2\right)},x^{\left(3\right)}$ and $h^{\left(0\right)}$ (which was initialized to zero).

  At first glance it looks like there is no meaningful differences between the architectures. We have an implicit sliding window over the words in the sentence for the recurrent neuron rather than explicit sliding window as with the FC neuron. Moreover, the output $y$ of FC neuron also depends on all previous time steps because the input $\mathbf{x}$ is just the values of $\left(x^{\left(1\right)},x^{\left(2\right)},x^{\left(3\right)}\right)$ .

  However, note that in the FC architecture, each word in the input vector has its own weight and bias whereas in recurrent architecture we use the same weights over and over again (by applying the same function). In other words, in recurrent layers we share the parameters across time. The formula below demonstrate how we applied the same function $s\in\mathcal{S}$ to updates the hidden-state across different time steps $$h^{\left(t\right)}=s\left(x^{\left(t\right)},h^{\left(t-1\right)}\right)\Longrightarrow h^{\left(3\right)}=s\left(x^{\left(3\right)},s\left(x^{\left(2\right)},s\left(x^{\left(1\right)},0\right)\right)\right)$$ then making a prediction using some $f\in\mathcal{F}_{\sigma}$ on the current state $$\hat{y}^{\left(3\right)}=f\left(h^{\left(3\right)}\right)$$ Parameter sharing has a very powerful effect on the learnability of the network both because it reduces the number of parameters that needed to be optimized and it make the network robust to small changes in the context (recall the example about the order of features in the input). In addition, the recurrent layer practically has a dynamic window size. Going back to our next-word prediction network, we can continue to apply the network on consecutive inputs, teaching it dependencies between words that are further than 3 time steps.

  Note that in our example, we ignored the “intermediate” predictions of the network: $\hat{y}^{\left(1\right)}$ and $\hat{y}^{\left(2\right)}$ , but we could keep them and feed them as inputs to another recurrent layer making our RNN network deeper. We will discuss it in a later post.

  6 Limitations of Recurrent Cells

  Unfortunately, RNNs are not perfect and maintaining memory is not enough, sometimes the networks needs to be able to forget, for instance, when predicting next-word in a sentence the network needs to reset its state once it reach the end of the sentence.

  Furthermore, parameter sharing has a serious drawback. Recall that applying the network for multiple time steps leads to a composition of the same $\mathbf{s}\in\mathcal{S}$ upon itself $t$ times, something like this $$\mathbf{h}^{\left(t\right)}=\mathbf{s}\left(\mathbf{x}^{\left(t\right)},\mathbf{s}\left(\mathbf{x}^{\left(t-1\right)},\mathbf{s}\left(\mathbf{x}^{\left(t-2\right)},...\right)\right)\right)$$ which may cause the gradient during back-propagation to either explode or vanish. This problem is known as the vanishing and exploding gradients. In order to get an intuition about that, we make the following simplifications: assume that the activation function inside $\mathbf{s}$ is the identity instead of $tanh$ , and $\mathbf{x}^{\left(t\right)}=\mathbf{0}$ for all $t$ . Suppose that $\mathbf{h}^{\left(0\right)}$ is given to us and it is not all zeros (that is, it encodes some information about the past).

  We have that $\mathbf{h}^{\left(t\right)}=U\mathbf{h}^{\left(t-1\right)}$ so applying the recurrent neuron for $t$ time steps amounts to $\mathbf{h}^{\left(t\right)}=U^{t}\mathbf{h}^{\left(0\right)}$ . What happen to $U^{t}$ as $t$ increases? it easier to see the answer when considering the scalar case. If $U$ was a real positive number, then either $U^{t}\to0$ when $U<1$ , or $U^{t}\to\infty$ when $U>1$ . The same goes with matrices, suppose $U$ has an eigenvalue decomposition such that $U=V^{\top}\Sigma V$ for orthogonal $V$ , so $U^{t}=V^{\top}\Sigma^{t}V$ causing the small entries on the diagonal of $\Sigma$ to vanish, and the large entries to explode. Since the gradients in our network are being scaled according to $\Sigma^{t}$ , they are either explode or vanished. The problem with very large gradients is that they cause the network to diverge (large updates to the parameters), on the other hand, very small gradients cause infinitesimally small updates to the parameters which prevent the network from learning. That is describe much more formally and rigorously here.

  Next post I will discuss two popular improvements to the recurrent neurons that helps avoid the limitations above: the LSTM and GRU cells.

  Acknowledgments: I want to thank Ron Shiff for his useful comments and corrections to this post.

  $\newcommand{\lyxlock}{}$

  1 The First Post in the Series

  Recurrent Neural Network (RNN) is a neural network architecture that leads to state of the art results in many tasks, such as machine translation, speech recognition, image captioning and many more. Because it is a rapidly evolving field, it is difficult to find all the information (both basic and advanced material) in one place, so I decided to write a series of posts that will discuss Recurrent Neural Networks from their basic formulation to their most recent variations.

  This post, as the first in the series, will present the concept of a recurrent neuron and its limitations. Next posts will discuss gated neurons like LSTM and GRU, variants of dropout regularization for RNNs and some advanced stuff that I found interesting (bi-directional RNN, Multiplicative LSTM and etc).

  The next section is dedicated to review the notations and basic terminology that will be used during the entire series. I will assume that the reader is familiar with the basics of neural network (fully connected architecture, back-propagation algorithm, etc). If you are not familiar with one of those concepts, I strongly recommend to review them online before going further. Those notes provides concise introduction to the field. There is also an excellent video-lecture (in Hebrew) by Prof. Shai Shalev-Shwartz from 2016. If you really want to dive deeper, you can find broad introduction in Chapter 6 in Deep Learning book (legally free electronic version here).

  2 Notations and Basic Terminology

  Artificial Neural Network (ANN) is a mathematical hierarchal model inspired from neuroscience. The basic computation unit is a neuron, which is a function from $\mathbb{R}^{d}$ to $\mathbb{R}$ , $$g:\mathbf{x}\mapsto\sigma\left(\left\langle \mathbf{x},\mathbf{w}\right\rangle +b\right)$$ where $\mathbf{x},\mathbf{w}\in\mathbb{R}^{d}$ are the input and weight vectors respectively, and $b\in\mathbb{R}$ is the bias. The function $\sigma$ is a non-linear activation function $\mathbb{R}\to\mathbb{R}$ , for example $\sigma\left(z\right)=ReLU\left(z\right):=\max\left(z,0\right)$ or $\sigma\left(z\right)=tanh\left(z\right):=\frac{e^{z}-e^{-z}}{e^{z}+e^{-z}}$ . Note that a neuron can ignore some of its input coordinates by setting their corresponding entries in its weight vector to zero.

  Neurons are organized within layers. In the very basic ANN architecture, called a fully connected (FC) neural network, a layer is just a collection of neurons that will be computed in parallel (meaning that they are not connected to each other). For any layer $i$ , each of its neurons accepts as inputs all neurons from layer $i-1$ .

  A layer in FC architecture can be formulate mathematically as a function from $\mathbb{R}^{d}$ to $\mathbb{R}^{n}$ such that $$\mathbf{f}:\mathbf{x}\mapsto\left(\begin{array}{c} g_{1}\left(\mathbf{x}\right)\\ \vdots\\ g_{n}\left(\mathbf{x}\right) \end{array}\right)$$ where each $g_{i}$ is a function of a neuron as defined above, with its own weight vector and a bias scalar. The formulation can be simplified using matrix notations: $$\mathbf{f}:\mathbf{x}\mapsto\sigma\left(W\mathbf{x}+\mathbf{b}\right)$$ where $W\in\mathbb{R}^{n\times d}$ contains in the i-th row the weight vector of $g_{i}$ , and $\mathbf{b}_{i}$ is the bias of $g_{i}$ . The function $\sigma$ is the same as before but since its input is now a vector, it apply element-wise on each coordinate of the input, producing an output vector of the same length.

  We denote the family of all possible FC functions based on some fixed activation function as $$\mathcal{F}_{\sigma}=\left\{ \mathbf{x}\mapsto\sigma\left(W\mathbf{x}+\mathbf{b}\right)\mid\forall W,\mathbf{b}\right\}$$ that family contains all possible layers for any input\output dimension. For abbreviation, I will not mention the dimensions of layers or the specific activation function $\sigma$ when they are not important for the discussion.

  Applying several layers sequentially yield a feedforward FC neural network. That is, given $\mathbf{x}\in\mathbb{R}^{d}$ , a FC network with $t$ layers is $$\mathcal{N}\left(\mathbf{x}\right)=\mathbf{f}_{t}\left(...\mathbf{f}_{2}\left(\mathbf{f}_{1}\left(\mathbf{x}\right)\right)\right)$$ where $\forall i:\:\mathbf{f}_{i}\in\mathcal{F}_{\sigma}$ for some fixed $\sigma$ .

  The number of layers in the network is its depth, and the size of the layers is its width (assuming all the layers are of the same dimension). Training the network is optimizing all the weight matrices and bias vectors such that the network will approximate some target (usually unknown) function. We can think of the optimization process as a process of iteratively updating the choice of layer functions $\mathbf{f}_{1},\mathbf{f}_{2},..$ from the family $\mathcal{F}_{\sigma}$ .

  The FC architecture is simple, yet yields a very powerful model. When wide enough, even a single hidden layer network can approximate to any arbitrary precision almost any continuous function (under some mild assumptions). However, for some tasks it might not be optimal, and different architectures are being used, instead of the FC one.

  3 The Need for Context

  There are tasks that require some notion of context in order to make a decent prediction. For example, tasks that involve time-steps can consider the history until the current time step as a context (a.k.a time-series tasks). For concrete example, consider predicting the next word in a given sequence. Each word has as a context all words preceding it. Note that there are more ways to define the notion of a context, and it mostly depends on the task in hand. For didactic purposes we will focus on the example of next-word prediction.

  How can we build a model for next-word prediction using FC neural network?

  The obvious answer is to fix some window size $d$ , and build a network that gets the last $d$ words in the sentence as features, in order to predict the next word. One immediate limitation of this approach is that the window size is the same for all the examples in the training set, and it is not simple to find a suitable window size. On the one hand, choosing small window may cause the model to miss some important long term relations between words that are more than $d$ steps apart (for sentences longer than $d$ ). On the other hand, choosing very large $d$ increases the number of parameters in the model, which make the optimization harder, and increase the computation time per example, which affect both training-time and prediction-time.

  Another limitation is that FC layers regards their inputs as ordered information. That is, a FC layer assigns a specific weight for each coordinate in its input vector, causing the order of the features to be meaningful. But in some tasks the order of the features may vary, for example, the subject of a sentence doesn’t limit to specific location inside an English sentence. If we want FC layer to account for that, we have to get a lot of examples (and they need to cover all different cases).

  What we would like to have here is called parameter sharing. We will come back to that term later in this post.

  4 Recurrent Cell - Learning to Remember

  We want our network to be able to remember information about previous invocations (i.e. the past inputs). In order to do that, our network have to maintain a state. So we need to extend our neuron formulation to include another variable, $h\in\mathbb{R}$ , that will be its state and be kept during consecutive invocations.That state variable is usually initialized to be zero so that it doesn’t affect the first computation.

  We will use an upper-script to denote the time step, so $h^{\left(t\right)}$ is the state of a neuron at time $t$ . Since the neurons need to pass information to themselves (i.e. their state from previous invocation) they become recurrent neurons, a.k.a recurrent cells.

  Recurrent neuron uses internally two functions. The first function updates the inner state based on current input and previous state, which is the transition from time $t-1$ to time $t$ , $$s:\mathbf{x},h^{\left(t-1\right)}\mapsto tanh\left(\left\langle \mathbf{x},\mathbf{w}\right\rangle +b+h^{\left(t-1\right)}\cdot u\right)$$ $u\in\mathbb{R}$ is a new parameter denoting the weight for the previous state.

  And the second function is for computing the output of the neuron at current time $t$ , based on its new state $$g^{\left(t\right)}:h^{\left(t\right)}\mapsto\sigma\left(h^{\left(t\right)}\cdot v+c\right)$$ with $v,c\in\mathbb{R}$ as the weight and bias of the computation respectively.

  We can move again to matrix notations to denote several parallel recurrent cells that forms a recurrent layer. The first function updates the vector of states for the entire layer $$\mathbf{s}:\mathbf{x},\mathbf{h}^{\left(t-1\right)}\mapsto tanh\left(W\mathbf{x}+\mathbf{b}+U\mathbf{h}^{\left(t-1\right)}\right)$$ and the second function for computing the output of the layer based on the current state $$\mathbf{f}:\mathbf{h}^{\left(t\right)}\mapsto\sigma\left(V\mathbf{h}^{\left(t\right)}+\mathbf{c}\right)$$ When we come to compare that to the FC formulation, we see that actually the function $\mathbf{f}$ itself hasn’t changed, but now it accepts as input the state of the layer $\mathbf{h}$ , rather than the input to the layer $\mathbf{x}$ . Another change is that before the layer computes its output, it first update its hidden state.

  Note that we can extend the definition of the family $\mathcal{F}_{tanh}$ to contains functions like the state-transforming function $\mathbf{s}$ , by making each $\mathbf{f}\in\mathcal{F}_{tanh}$ accepts two vectors as input instead of one (and we also need another weight matrix for the second input). So denote the family of all hidden-state transition functions as $$\mathcal{S}_{\sigma}=\left\{ \mathbf{x},\mathbf{h}^{\left(t-1\right)}\mapsto\sigma\left(W\mathbf{x}+\mathbf{b}+U\mathbf{h}^{\left(t-1\right)}\right)\mid\forall W,U,\mathbf{b}\right\}$$ usually the activation function is $\sigma=tanh$ . We will use that family later on, when discuss gated units as LSTM and GRU.

  5 Comparing the Architectures - an Example

  Consider again the next-word prediction problem. Suppose we have an ordered dictionary that contains the entire vocabulary we will ever need. We will represent any word by an integer indicating its index in the dictionary. In this section we will construct two single-neuron networks, once using the FC architecture and once using the Recurrent architecture.

  FC network for next-word prediction: Let the window size be $d=3$ . The input to the network can be represent as a vector $\mathbf{x}\in\mathbb{N}^{3}$ and the output as $y\in\mathbb{N}$ . The FC network amounts to $$\mathcal{N}_{FC}\left(\mathbf{x}\right)=f\left(x\right)$$ for some $f\in\mathcal{F}_{\sigma}$ that maps from $\mathbb{R}^{3}$ to $\mathbb{R}$ .

  Recurrent neural network (RNN), however, leads to different representation: the input and the output both integers $x,y\in\mathbb{N},$ so the state transition function $s\in\mathcal{S}$ and the output function $f\in\mathcal{F}_{\sigma}$ are both from $\mathbb{R}$ to $\mathbb{R}$ . The idea of the recurrent architecture is to apply the same recurrent layer several times on different inputs (each input from different time-step) in order to get a prediction.

  Here is a detailed description of how to do that: we take a series of three consecutive words from the beginning of a sentence: $x^{\left(1\right)},x^{\left(2\right)},x^{\left(3\right)}$ and a target next word: $y$ . First, feed the network with $x^{\left(1\right)}$ , it updates the inner state from $h^{\left(0\right)}=0$ to some hidden state $h^{\left(1\right)}$ and outputs $\hat{y}^{\left(1\right)}$ . We ignore that output and feed the next word $x^{\left(2\right)}$ , the network updates its state to $h^{\left(2\right)}$ and produces $\hat{y}^{\left(2\right)}$ that is also ignored. Lastly, we feed $x^{\left(3\right)}$ to the network, the inner state is being updated to $h^{\left(3\right)}$ and another output is produced $\hat{y}^{\left(3\right)}$ which is now considered as the prediction of the network. Note that in each invocation of the network, the output depends on both the current input and the previous state, therefore $\hat{y}^{\left(3\right)}$ is the result of the entire history $x^{\left(1\right)},x^{\left(2\right)},x^{\left(3\right)}$ and $h^{\left(0\right)}$ (which was initialized to zero).

  At first glance it looks like there is no meaningful differences between the architectures. We have an implicit sliding window over the words in the sentence for the recurrent neuron rather than explicit sliding window as with the FC neuron. Moreover, the output $y$ of FC neuron also depends on all previous time steps because the input $\mathbf{x}$ is just the values of $\left(x^{\left(1\right)},x^{\left(2\right)},x^{\left(3\right)}\right)$ .

  However, note that in the FC architecture, each word in the input vector has its own weight and bias whereas in recurrent architecture we use the same weights over and over again (by applying the same function). In other words, in recurrent layers we share the parameters across time. The formula below demonstrate how we applied the same function $s\in\mathcal{S}$ to updates the hidden-state across different time steps $$h^{\left(t\right)}=s\left(x^{\left(t\right)},h^{\left(t-1\right)}\right)\Longrightarrow h^{\left(3\right)}=s\left(x^{\left(3\right)},s\left(x^{\left(2\right)},s\left(x^{\left(1\right)},0\right)\right)\right)$$ then making a prediction using some $f\in\mathcal{F}_{\sigma}$ on the current state $$\hat{y}^{\left(3\right)}=f\left(h^{\left(3\right)}\right)$$ Parameter sharing has a very powerful effect on the learnability of the network both because it reduces the number of parameters that needed to be optimized and it make the network robust to small changes in the context (recall the example about the order of features in the input). In addition, the recurrent layer practically has a dynamic window size. Going back to our next-word prediction network, we can continue to apply the network on consecutive inputs, teaching it dependencies between words that are further than 3 time steps.

  Note that in our example, we ignored the “intermediate” predictions of the network: $\hat{y}^{\left(1\right)}$ and $\hat{y}^{\left(2\right)}$ , but we could keep them and feed them as inputs to another recurrent layer making our RNN network deeper. We will discuss it in a later post.

  6 Limitations of Recurrent Cells

  Unfortunately, RNNs are not perfect and maintaining memory is not enough, sometimes the networks needs to be able to forget, for instance, when predicting next-word in a sentence the network needs to reset its state once it reach the end of the sentence.

  Furthermore, parameter sharing has a serious drawback. Recall that applying the network for multiple time steps leads to a composition of the same $\mathbf{s}\in\mathcal{S}$ upon itself $t$ times, something like this $$\mathbf{h}^{\left(t\right)}=\mathbf{s}\left(\mathbf{x}^{\left(t\right)},\mathbf{s}\left(\mathbf{x}^{\left(t-1\right)},\mathbf{s}\left(\mathbf{x}^{\left(t-2\right)},...\right)\right)\right)$$ which may cause the gradient during back-propagation to either explode or vanish. This problem is known as the vanishing and exploding gradients. In order to get an intuition about that, we make the following simplifications: assume that the activation function inside $\mathbf{s}$ is the identity instead of $tanh$ , and $\mathbf{x}^{\left(t\right)}=\mathbf{0}$ for all $t$ . Suppose that $\mathbf{h}^{\left(0\right)}$ is given to us and it is not all zeros (that is, it encodes some information about the past).

  We have that $\mathbf{h}^{\left(t\right)}=U\mathbf{h}^{\left(t-1\right)}$ so applying the recurrent neuron for $t$ time steps amounts to $\mathbf{h}^{\left(t\right)}=U^{t}\mathbf{h}^{\left(0\right)}$ . What happen to $U^{t}$ as $t$ increases? it easier to see the answer when considering the scalar case. If $U$ was a real positive number, then either $U^{t}\to0$ when $U<1$ , or $U^{t}\to\infty$ when $U>1$ . The same goes with matrices, suppose $U$ has an eigenvalue decomposition such that $U=V^{\top}\Sigma V$ for orthogonal $V$ , so $U^{t}=V^{\top}\Sigma^{t}V$ causing the small entries on the diagonal of $\Sigma$ to vanish, and the large entries to explode. Since the gradients in our network are being scaled according to $\Sigma^{t}$ , they are either explode or vanished. The problem with very large gradients is that they cause the network to diverge (large updates to the parameters), on the other hand, very small gradients cause infinitesimally small updates to the parameters which prevent the network from learning. That is describe much more formally and rigorously here.

  Next post I will discuss two popular improvements to the recurrent neurons that helps avoid the limitations above: the LSTM and GRU cells.

  Acknowledgments: I want to thank Ron Shiff for his useful comments and corrections to this post.
• Apr 6, 2017

  Memory Efficient Algorithm for Finding Cycles in Graphs

  tags:
  TheoryLunch Blog

  Computer science topics over lunch

  Memory Efficient Algorithm for Finding Cycles in Graphs

  Apr 6, 2017 • Eran Amar

  tags: $\newcommand{\lyxlock}{}$

  1 Cycle Detection in Graphs

  Somewhere in 2015, I encountered a nice algorithmic question during an interview for an algotreading company. The question was how one can decide if there is a cycle in an undirected connected graph? There are two immediate answers here, the first is running Breadth First Search (BFS) on arbitrary node and memorizing the nodes that were already visited. The search ends either when the algorithm complete scanning the whole graph or when it reaches some node twice, outputting “No” or “Yes” respectively. Another possible solution is by running Depth First Search (DFS) and keeping a variable to count the length of the longest path explored. The algorithm declare a cycle if at any point in time the length of the path reaches the number of vertices in the graph. For the rest of this post, denote that number by $n$ .

  Those two solutions are simple, and more important for this discussion, have the same space complexity. That is, both consume around the same amount of working memory (in bits) ignoring the memory needed for the input and output of the algorithm. Even though DFS have implementation that requires less memory than BFS, ultimately both requires $\Omega\left(n\log n\right)$ bits. One can (informally) justify the aforementioned space complexity as follow: both algorithms have to keep track of their search (a queue in BFS and a stack is DFS) which potentially can be as large as the entire graph thus $\Omega\left(n\right)$ , and each node need $\Omega\left(\log n\right)$ bits for encoding so the total space complexity is $\Omega\left(n\log n\right)$ bits.

  During the rest of this post we will focus on space complexity rather than time complexity. We will present an algorithm that can determine if there is a cycle in $\mathcal{O}\left(\log^{2}n\right)$ bits [base on this paper] and then discuss the importance of related problem: the s-t connectivity.

  2 The s-t Connectivity Problem

  2.1 Solving s-t connectivity as a Sub-step

  We will start from solving a slightly different problem then uses its solution as a procedure in our cycle-detection algorithm.

  Introducing the undirected s-t connectivity problem: The input is an undirected graph $G=\left(V,E\right)$ of size $n$ and two vertices $s,t\in V$ . The output should be $1$ if there is path in $G$ that connects $s$ and $t$ , and $0$ otherwise. The directed s-t connectivity is about the same but the input graph is directed (i.e its edges has directions) and the output should be $1$ whether there is a directed path from $s$ to $t$ .

  Before we present the algorithm, we first define the following predicate: for any $u,v\in V$ and $k\in\mathbb{N}$ define $conn\left(u,v;k\right)=1$ if exists a path from $u$ to $v$ of at most $k$ steps (i.e. edges), otherwise $conn\left(u,v;k\right)=0$ .

  The next observations are key points needed to solve the s-t connectivity efficiently (both for directed and undirected case).

  1. For any $u,v\in V$ and any $k\in\mathbb{N}$ , we have that $conn\left(u,v;2^{k}\right)=1$ if and only if $\exists w\in V$ such that $conn\left(u,w;2^{k-1}\right)=1$ and $conn\left(w,v;2^{k-1}\right)=1$
  2. Let $u,v\in V$ , they are connected if and only if $conn\left(u,v;2^{\lceil\log n\rceil}\right)=1$

  From here the algorithm is straightforward: denote the algorithm by $Alg\left(u,v;2^{k}\right)$ that upon receiving $u,v$ from the graph $G=\left(V,E\right)$ and $k\in\mathbb{N}$ :

  1. If $k=0$ , it outputs $\mathbf{1}_{\left[u=v\right]}$ , and if $k=1$ it outputs $\mathbf{1}_{\left[u=v\vee\left(u,v\right)\in E\right]}$
  2. For each $w\in V$ : if $Alg\left(u,w;2^{k-1}\right)=1$ and $Alg\left(w,v;2^{k-1}\right)=1$ then output $1$
  3. Output $0$

  The solution for s-t connectivity can be achieved by outputting “Yes” if and only if $Alg\left(s,t,2^{\lceil\log n\rceil}\right)=1$ .

  Obviously, the depth of the recursion is $\mathcal{O}\left(\log n\right)$ and in each level of the recursion the algorithm only needs to keep track of $w$ and $k$ , which is $\mathcal{O}\left(\log n\right)$ bits, therefore the total memory complexity of the algorithm is $\mathcal{O}\left(\log^{2}n\right)$ bits.

  2.2 Extending the algorithm to detect cycles

  How solving the s-t connectivity problem helps us with cycles detection? The next observation answer that: let $G=\left(V,E\right)$ be the graph in question and $u,v\in V$ two endpoints of some edge, they must stay connected after removing that edge from $G$ if and only if both participating in a cycle contained in $G$ .

  For the directed case the algorithm is as follow: iterate over any edge $\left(u,v\right)\in E$ , delete it temporarily from $G$ and check if $Alg\left(v,u;2^{\lceil\log n\rceil}\right)=1$ on the resulting graph, that is, if $v$ is still connected to $u$ even after you delete the edge from $u$ to $v$ . If $Alg\left(v,u;2^{\lceil\log n\rceil}\right)=1$ then declare a cycle, else, restore back the edge and continue to the next one. Note that this solution also works for undirected graphs. The memory complexity is $\mathcal{O}\left(\log^{2}n\right)+\log n=\mathcal{O}\left(\log^{2}n\right)$ and we done.

  2.3 Is it the best we can do?

  Actually, for undirected graphs there is an algorithm with $\mathcal{O}\left(\log n\right)$ bits, improving the square on the log in the construction from previous section. The curious reader can read more about it here. For the directed case it is an open question whether we can improve the $\log^{2}n$ or not.

  2.4 Beyond Cycles in Graphs

  The directed s-t connectivity problem is an important theoretical problem because it captures the “hardness” of any problem that requires logarithmic amount of memory space. The reduction to show that is actually very easy, suppose we have an algorithm that solves some problem in $\mathcal{O}\left(s\right)$ space where $s$ is logarithmic in $n$ (and $n$ is the size of the input), generate the following directed graph: create a vertex $v$ for any state (memory configuration) of the algorithm. Since its memory is $\mathcal{O}\left(s\right)$ , there are at most $2^{\mathcal{O}\left(s\right)}=\mathcal{O}\left(poly\left(n\right)\right)$ states (thus vertices). Create an edge from $u$ to $v$ iff $v$ is a valid possible successor state of $u$ . Add a vertex $t$ and connects to it all the states that leads to outputting “Yes”. Denote by $s$ the starting-state of the algorithm and now solving the original problem is equivalent to finding a directed path from $s$ to $t$ in the states-graph.

  This reduction was presented by Walter Savitch in 1970, and it also works for Nondeterministic algorithms. The conclusion is that any nondeterministic problem that can be solved in $\mathcal{O}\left(s\right)$ space (where $s$ is logarithmic in the size of the input) can be solved deterministically in $\mathcal{O}\left(s^{2}\right)$ . Combining that conclusion with the algorithm for directed s-t connectivity yields that any problem with nondeterministic $\mathcal{O}\left(\log n\right)$ space complexity can be solved deterministically (via reduction to directed s-t connectivity) with $\mathcal{O}\left(\log^{2}n\right)$ space complexity.

  $\newcommand{\lyxlock}{}$

  1 Cycle Detection in Graphs

  Somewhere in 2015, I encountered a nice algorithmic question during an interview for an algotreading company. The question was how one can decide if there is a cycle in an undirected connected graph? There are two immediate answers here, the first is running Breadth First Search (BFS) on arbitrary node and memorizing the nodes that were already visited. The search ends either when the algorithm complete scanning the whole graph or when it reaches some node twice, outputting “No” or “Yes” respectively. Another possible solution is by running Depth First Search (DFS) and keeping a variable to count the length of the longest path explored. The algorithm declare a cycle if at any point in time the length of the path reaches the number of vertices in the graph. For the rest of this post, denote that number by $n$ .

  Those two solutions are simple, and more important for this discussion, have the same space complexity. That is, both consume around the same amount of working memory (in bits) ignoring the memory needed for the input and output of the algorithm. Even though DFS have implementation that requires less memory than BFS, ultimately both requires $\Omega\left(n\log n\right)$ bits. One can (informally) justify the aforementioned space complexity as follow: both algorithms have to keep track of their search (a queue in BFS and a stack is DFS) which potentially can be as large as the entire graph thus $\Omega\left(n\right)$ , and each node need $\Omega\left(\log n\right)$ bits for encoding so the total space complexity is $\Omega\left(n\log n\right)$ bits.

  During the rest of this post we will focus on space complexity rather than time complexity. We will present an algorithm that can determine if there is a cycle in $\mathcal{O}\left(\log^{2}n\right)$ bits [base on this paper] and then discuss the importance of related problem: the s-t connectivity.

  2 The s-t Connectivity Problem

  2.1 Solving s-t connectivity as a Sub-step

  We will start from solving a slightly different problem then uses its solution as a procedure in our cycle-detection algorithm.

  Introducing the undirected s-t connectivity problem: The input is an undirected graph $G=\left(V,E\right)$ of size $n$ and two vertices $s,t\in V$ . The output should be $1$ if there is path in $G$ that connects $s$ and $t$ , and $0$ otherwise. The directed s-t connectivity is about the same but the input graph is directed (i.e its edges has directions) and the output should be $1$ whether there is a directed path from $s$ to $t$ .

  Before we present the algorithm, we first define the following predicate: for any $u,v\in V$ and $k\in\mathbb{N}$ define $conn\left(u,v;k\right)=1$ if exists a path from $u$ to $v$ of at most $k$ steps (i.e. edges), otherwise $conn\left(u,v;k\right)=0$ .

  The next observations are key points needed to solve the s-t connectivity efficiently (both for directed and undirected case).

  1. For any $u,v\in V$ and any $k\in\mathbb{N}$ , we have that $conn\left(u,v;2^{k}\right)=1$ if and only if $\exists w\in V$ such that $conn\left(u,w;2^{k-1}\right)=1$ and $conn\left(w,v;2^{k-1}\right)=1$
  2. Let $u,v\in V$ , they are connected if and only if $conn\left(u,v;2^{\lceil\log n\rceil}\right)=1$

  From here the algorithm is straightforward: denote the algorithm by $Alg\left(u,v;2^{k}\right)$ that upon receiving $u,v$ from the graph $G=\left(V,E\right)$ and $k\in\mathbb{N}$ :

  1. If $k=0$ , it outputs $\mathbf{1}_{\left[u=v\right]}$ , and if $k=1$ it outputs $\mathbf{1}_{\left[u=v\vee\left(u,v\right)\in E\right]}$
  2. For each $w\in V$ : if $Alg\left(u,w;2^{k-1}\right)=1$ and $Alg\left(w,v;2^{k-1}\right)=1$ then output $1$
  3. Output $0$

  The solution for s-t connectivity can be achieved by outputting “Yes” if and only if $Alg\left(s,t,2^{\lceil\log n\rceil}\right)=1$ .

  Obviously, the depth of the recursion is $\mathcal{O}\left(\log n\right)$ and in each level of the recursion the algorithm only needs to keep track of $w$ and $k$ , which is $\mathcal{O}\left(\log n\right)$ bits, therefore the total memory complexity of the algorithm is $\mathcal{O}\left(\log^{2}n\right)$ bits.

  2.2 Extending the algorithm to detect cycles

  How solving the s-t connectivity problem helps us with cycles detection? The next observation answer that: let $G=\left(V,E\right)$ be the graph in question and $u,v\in V$ two endpoints of some edge, they must stay connected after removing that edge from $G$ if and only if both participating in a cycle contained in $G$ .

  For the directed case the algorithm is as follow: iterate over any edge $\left(u,v\right)\in E$ , delete it temporarily from $G$ and check if $Alg\left(v,u;2^{\lceil\log n\rceil}\right)=1$ on the resulting graph, that is, if $v$ is still connected to $u$ even after you delete the edge from $u$ to $v$ . If $Alg\left(v,u;2^{\lceil\log n\rceil}\right)=1$ then declare a cycle, else, restore back the edge and continue to the next one. Note that this solution also works for undirected graphs. The memory complexity is $\mathcal{O}\left(\log^{2}n\right)+\log n=\mathcal{O}\left(\log^{2}n\right)$ and we done.

  2.3 Is it the best we can do?

  Actually, for undirected graphs there is an algorithm with $\mathcal{O}\left(\log n\right)$ bits, improving the square on the log in the construction from previous section. The curious reader can read more about it here. For the directed case it is an open question whether we can improve the $\log^{2}n$ or not.

  2.4 Beyond Cycles in Graphs

  The directed s-t connectivity problem is an important theoretical problem because it captures the “hardness” of any problem that requires logarithmic amount of memory space. The reduction to show that is actually very easy, suppose we have an algorithm that solves some problem in $\mathcal{O}\left(s\right)$ space where $s$ is logarithmic in $n$ (and $n$ is the size of the input), generate the following directed graph: create a vertex $v$ for any state (memory configuration) of the algorithm. Since its memory is $\mathcal{O}\left(s\right)$ , there are at most $2^{\mathcal{O}\left(s\right)}=\mathcal{O}\left(poly\left(n\right)\right)$ states (thus vertices). Create an edge from $u$ to $v$ iff $v$ is a valid possible successor state of $u$ . Add a vertex $t$ and connects to it all the states that leads to outputting “Yes”. Denote by $s$ the starting-state of the algorithm and now solving the original problem is equivalent to finding a directed path from $s$ to $t$ in the states-graph.

  This reduction was presented by Walter Savitch in 1970, and it also works for Nondeterministic algorithms. The conclusion is that any nondeterministic problem that can be solved in $\mathcal{O}\left(s\right)$ space (where $s$ is logarithmic in the size of the input) can be solved deterministically in $\mathcal{O}\left(s^{2}\right)$ . Combining that conclusion with the algorithm for directed s-t connectivity yields that any problem with nondeterministic $\mathcal{O}\left(\log n\right)$ space complexity can be solved deterministically (via reduction to directed s-t connectivity) with $\mathcal{O}\left(\log^{2}n\right)$ space complexity.
• Mar 30, 2017

  Derandomization for Pairwise Independent Seed

  tags:
  TheoryLunch Blog

  Computer science topics over lunch

  Derandomization for Pairwise Independent Seed

  Mar 30, 2017 • Eran Amar

  tags: $\newcommand{\lyxlock}{}$

  1 Randomized Algorithm Settings

  Randomized algorithm is an algorithm that in addition to its inputs also uses internally $t$ coin-tosses (i.e. uniformly random bits) in order to decide on its output, so different invocations of the algorithm with same inputs may yield different outputs that depend on the results of these coins. One can think about the internal coin-tosses as a binary string $r$ of length $t$ that is being fixed before the algorithm starts and being consumed sequentially during the algorithm’s execution. Such string is often called seed, and the randomness within the algorithm now comes from choosing a different string at the beginning of each execution.

  Note that tossing $t$ fair coins, is exactly like drawing uniformly at random a string $r$ from the set $\left\{ 0,1\right\} ^{t}$ (i.e. each string has a probability of $2^{-t}$ ).

  Suppose we have a randomized algorithm $A$ that solves a problem of interest. We have a guarantee that for any input $x$ of length $n$ , algorithm $A$ runs in $poly\left(n\right)$ , and outputs either a valid solution $S$ if it succeeds or output $failed$ . The success probability of the algorithm is $p$ (which may be even infinitesimally small).

  The first question to ask is, can we convert the randomized algorithm into a deterministic one?

  2 Brute Force Derandomization

  The answer is yes, and there is a simple way to do that: Upon receiving an input $x$ , iterate over all possible strings $r\in\left\{ 0,1\right\} ^{t}$ and run $A_{r}\left(x\right)$ . That is, simulate $A$ on $x$ with $r$ as the random seed. Return the first solution that $A$ outputs (i.e. anything that is not $failed$ ).

  The idea behind this method is that if there is a non-zero probability for $A$ to output a valid solution for the input $x$ , then there must be a random seed that leads to that output when executing algorithm $A$ on $x$ . So all we need to do is iterating over all possible seeds of length $t$ , and that leads to a deterministic algorithm.

  You probably already see the problem here. If $t=\mathcal{O}\left(n\right)$ then the running time of the deterministic procedure is $2^{\mathcal{O}\left(n\right)}\cdot poly\left(n\right)$ .

  That is, of course, terrible running time.

  On the other hand, if $t$ is a constant or even up to $\mathcal{O}\left(\log n\right)$ , then the resulting running time will stay polynomial, and everyone is happy. In the next section we will see that with an additional assumption, we can “decrease” the seed length from $n$ to $\log n$ bits.

  3 Pairwise Independent Seed

  We start off with a definition.

  k-wise Independent Family. Let $\mathcal{X}$ be a set of random variables. We say that $\mathcal{X}$ is a k-wise independent family if for any subset of $k$ variables $X_{1},..,X_{k}\in\mathcal{X}$ , and values $y_{1},...,y_{k}$ $$\mathbf{P}\left[\bigwedge_{i\in I}X_{i}=y_{i}\right]=\prod_{i\in I}\mathbf{P}\left[X_{i}=y_{i}\right]\qquad\forall I\subseteq\left[k\right]$$ that is, any subset of $k$ random variables from $\mathcal{X}$ is fully independent. The family $\mathcal{X}$ is called pairwise independent family when $k=2$ .

  What this definition have with our settings? Observe that up until now we implicitly assumed that the algorithm needs its coins to be fully independent. That is, if we think of $s$ as a series of Bernoulli RVs then we actually require $s$ to be t-wise independent family, which is a much stronger property compare to the pairwise independence property. If we relax that requirement and assume that the algorithm only needs its coins to be pairwise independent, then there is a way (described in next section) to construct $n$ pairwise independent bits using only $\mathcal{O}\left(\log n\right)$ fully independent bits. That construction can be “embedded” into our algorithm such that we ultimately end up with a procedure that only requires $\mathcal{O}\left(\log n\right)$ fully independent seed. Applying the derandomization method from previous section will result with a deterministic algorithm that will still run in $poly\left(n\right)$ time.

  Of course, the aforementioned derandomization method is valid only if the assumption of pairwise independent seed doesn’t break the correctness analysis of the randomized algorithm.

  4 Stretching the Seed

  Suppose we have $\approx\log n$ fair fully independent coins, we want to construct a method to generate from them $\approx n$ fair pairwise independent bits. The following claim shows such construction.

  Suppose we have $t:=\log\left(n\right)+1$ fully independent Bernoulli random variables, each with probability $\frac{1}{2}$ . Fix them in some arbitrary order to get a (random) vector $\mathbf{s}$ of length $t$ . Let $\mathcal{S}=\left\{ \left\langle \mathbf{s},\mathbf{x}\right\rangle _{mod2}\mid\mathbf{x}\in\left\{ 0,1\right\} ^{t},\:\:x\ne\mathbf{0}\right\}$ be the set of all possible sums modulo 2 of those RVs and note that $\left|\mathcal{S}\right|=2^{t}-1=\Theta\left(n\right)$ , then $\mathcal{S}$ is a family of pairwise independent bits.

  To prove that we start with a very basic lemma.

  Let $X,Y\sim Bernoulli\left(0.5\right)$ be two independent random variables, let $Z=X+Y$ then $Z\sim Bernoulli\left(0.5\right)$ . Moreover, $Z$ is independent from $X$ and from $Y$ .

  For any $b\in\left\{ 0,1\right\}$ $$\begin{aligned} \mathbf{P}\left[X+Y=b\right] & =\mathbf{P}\left[X=1,Y=1-b\right]+\mathbf{P}\left[X=0,Y=b\right]\\ & =\mathbf{P}\left[X=1\right]\mathbf{P}\left[Y=1-b\right]+\mathbf{P}\left[X=0\right]\mathbf{P}\left[Y=b\right]\\ & =\frac{1}{2}\left(\mathbf{P}\left[Y=1-b\right]+\mathbf{P}\left[Y=1-b\right]\right)\\ & =\frac{1}{2} \end{aligned}$$ thus $Z\sim Bernoulli\left(0.5\right)$ . Next w.l.o.g we prove that $Z$ independent from $X$ . Take any $b_{1},b_{2}\in\left\{ 0,1\right\}$ . In case $b_{1}=b_{2}$ then $$\begin{aligned} \mathbf{P}\left[X=b_{1},Z=b_{2}\right] & =\mathbf{P}\left[X=b_{1},X+Y=b_{1}\right]\\ & =\mathbf{P}\left[X=b_{1},Y=0\right]\\ & =\mathbf{P}\left[X=b_{1}\right]\mathbf{P}\left[Y=0\right]\\ & =\mathbf{P}\left[X=b_{1}\right]\cdot\frac{1}{2}\\ & =\mathbf{P}\left[X=b_{1}\right]\mathbf{P}\left[Z=b_{1}\right] \end{aligned}$$ so $X$ and $Z$ are independent. And in case that $b_{2}=1-b_{1}$ then similarly we get $$\begin{aligned} \mathbf{P}\left[X=b_{1},Z=b_{2}\right] & =\mathbf{P}\left[X=b_{1},X+Y=1-b_{1}\right]\\ & =\mathbf{P}\left[X=b_{1},Y=1\right]\\ & =\mathbf{P}\left[X=b_{1}\right]\mathbf{P}\left[Z=b_{1}\right] \end{aligned}$$ which conclude the proof.

  Now we prove the claim from above.

  Take two different variables $X_{1}=\left\langle \mathbf{s},\mathbf{x}_{1}\right\rangle \in\mathcal{S}$ and $X_{2}=\left\langle \mathbf{s},\mathbf{x}_{2}\right\rangle \in\mathcal{S}$ , and let $b_{1},b_{2}\in\left\{ 0,1\right\}$ , we want to show that $X_{1}$ and $X_{2}$ are independent. Obviously, we only care about the set of indices which are non-zero in the vectors $\mathbf{x}_{1}$ and $\mathbf{x}_{2}$ . Let these sets be $I_{1},I_{2}$ respectively. Denote $M_{0}=\sum_{i\in I_{1}\cap I_{2}}\mathbf{s}_{i}$ , $M_{1}=\sum_{i\in I_{1}\backslash I_{2}}\mathbf{s}_{i}$ and $M_{2}=\sum_{i\in I_{2}\backslash I_{1}}\mathbf{s}_{i}$ where summation over empty set is defined as the constant $0$ (all the operations in the proof are modulo 2). Note that any pair of variables from $\left\{ M_{0},M_{1},M_{2}\right\}$ are independent because they defined on mutually disjoint set of indices. Moreover, by the lemma from above, each of them if it is not the constant zero then is distributed as $Bernoulli\left(0.5\right)$ . When plugging-in the new notations we get $$\begin{aligned} \mathbf{P}\left[X_{1}=b_{1},\:X_{2}=b_{2}\right] & =\mathbf{P}\left[M_{0}+M_{1}=b_{1},\:M_{0}+M_{2}=b_{2}\right] \end{aligned}$$ Because $I_{1}\ne I_{2}$ there must be at most one variable from $\left\{ M_{0},M_{1},M_{2}\right\}$ which may be constant. Therefore, no matter which $M_{i}$ it is (if any), either $X_{1}$ or $X_{2}$ is a sum of two non-constant variables. W.l.o.g let it be $X_{1}$ , then it has the $M_{1}$ component which does not appear in $X_{2}$ and can flip the result of the sum with probability $0.5$ regardless of what happen anywhere else, therefore $X_{1}$ is independent of $X_{2}$ , and that completes the proof.

  Last one remark. In our settings the randomized algorithm either return a solution or says “failed”, if we have a polynomial procedure to validate a solution we can generalize the settings to the case when the algorithm always output a solution, and then we verify it to determine if it is “succeeds” case of “failed” case.

  $\newcommand{\lyxlock}{}$

  1 Randomized Algorithm Settings

  Randomized algorithm is an algorithm that in addition to its inputs also uses internally $t$ coin-tosses (i.e. uniformly random bits) in order to decide on its output, so different invocations of the algorithm with same inputs may yield different outputs that depend on the results of these coins. One can think about the internal coin-tosses as a binary string $r$ of length $t$ that is being fixed before the algorithm starts and being consumed sequentially during the algorithm’s execution. Such string is often called seed, and the randomness within the algorithm now comes from choosing a different string at the beginning of each execution.

  Note that tossing $t$ fair coins, is exactly like drawing uniformly at random a string $r$ from the set $\left\{ 0,1\right\} ^{t}$ (i.e. each string has a probability of $2^{-t}$ ).

  Suppose we have a randomized algorithm $A$ that solves a problem of interest. We have a guarantee that for any input $x$ of length $n$ , algorithm $A$ runs in $poly\left(n\right)$ , and outputs either a valid solution $S$ if it succeeds or output $failed$ . The success probability of the algorithm is $p$ (which may be even infinitesimally small).

  The first question to ask is, can we convert the randomized algorithm into a deterministic one?

  2 Brute Force Derandomization

  The answer is yes, and there is a simple way to do that: Upon receiving an input $x$ , iterate over all possible strings $r\in\left\{ 0,1\right\} ^{t}$ and run $A_{r}\left(x\right)$ . That is, simulate $A$ on $x$ with $r$ as the random seed. Return the first solution that $A$ outputs (i.e. anything that is not $failed$ ).

  The idea behind this method is that if there is a non-zero probability for $A$ to output a valid solution for the input $x$ , then there must be a random seed that leads to that output when executing algorithm $A$ on $x$ . So all we need to do is iterating over all possible seeds of length $t$ , and that leads to a deterministic algorithm.

  You probably already see the problem here. If $t=\mathcal{O}\left(n\right)$ then the running time of the deterministic procedure is $2^{\mathcal{O}\left(n\right)}\cdot poly\left(n\right)$ .

  That is, of course, terrible running time.

  On the other hand, if $t$ is a constant or even up to $\mathcal{O}\left(\log n\right)$ , then the resulting running time will stay polynomial, and everyone is happy. In the next section we will see that with an additional assumption, we can “decrease” the seed length from $n$ to $\log n$ bits.

  3 Pairwise Independent Seed

  We start off with a definition.

  k-wise Independent Family. Let $\mathcal{X}$ be a set of random variables. We say that $\mathcal{X}$ is a k-wise independent family if for any subset of $k$ variables $X_{1},..,X_{k}\in\mathcal{X}$ , and values $y_{1},...,y_{k}$ $$\mathbf{P}\left[\bigwedge_{i\in I}X_{i}=y_{i}\right]=\prod_{i\in I}\mathbf{P}\left[X_{i}=y_{i}\right]\qquad\forall I\subseteq\left[k\right]$$ that is, any subset of $k$ random variables from $\mathcal{X}$ is fully independent. The family $\mathcal{X}$ is called pairwise independent family when $k=2$ .

  What this definition have with our settings? Observe that up until now we implicitly assumed that the algorithm needs its coins to be fully independent. That is, if we think of $s$ as a series of Bernoulli RVs then we actually require $s$ to be t-wise independent family, which is a much stronger property compare to the pairwise independence property. If we relax that requirement and assume that the algorithm only needs its coins to be pairwise independent, then there is a way (described in next section) to construct $n$ pairwise independent bits using only $\mathcal{O}\left(\log n\right)$ fully independent bits. That construction can be “embedded” into our algorithm such that we ultimately end up with a procedure that only requires $\mathcal{O}\left(\log n\right)$ fully independent seed. Applying the derandomization method from previous section will result with a deterministic algorithm that will still run in $poly\left(n\right)$ time.

  Of course, the aforementioned derandomization method is valid only if the assumption of pairwise independent seed doesn’t break the correctness analysis of the randomized algorithm.

  4 Stretching the Seed

  Suppose we have $\approx\log n$ fair fully independent coins, we want to construct a method to generate from them $\approx n$ fair pairwise independent bits. The following claim shows such construction.

  Suppose we have $t:=\log\left(n\right)+1$ fully independent Bernoulli random variables, each with probability $\frac{1}{2}$ . Fix them in some arbitrary order to get a (random) vector $\mathbf{s}$ of length $t$ . Let $\mathcal{S}=\left\{ \left\langle \mathbf{s},\mathbf{x}\right\rangle _{mod2}\mid\mathbf{x}\in\left\{ 0,1\right\} ^{t},\:\:x\ne\mathbf{0}\right\}$ be the set of all possible sums modulo 2 of those RVs and note that $\left|\mathcal{S}\right|=2^{t}-1=\Theta\left(n\right)$ , then $\mathcal{S}$ is a family of pairwise independent bits.

  To prove that we start with a very basic lemma.

  Let $X,Y\sim Bernoulli\left(0.5\right)$ be two independent random variables, let $Z=X+Y$ then $Z\sim Bernoulli\left(0.5\right)$ . Moreover, $Z$ is independent from $X$ and from $Y$ .

  For any $b\in\left\{ 0,1\right\}$ $$\begin{aligned} \mathbf{P}\left[X+Y=b\right] & =\mathbf{P}\left[X=1,Y=1-b\right]+\mathbf{P}\left[X=0,Y=b\right]\\ & =\mathbf{P}\left[X=1\right]\mathbf{P}\left[Y=1-b\right]+\mathbf{P}\left[X=0\right]\mathbf{P}\left[Y=b\right]\\ & =\frac{1}{2}\left(\mathbf{P}\left[Y=1-b\right]+\mathbf{P}\left[Y=1-b\right]\right)\\ & =\frac{1}{2} \end{aligned}$$ thus $Z\sim Bernoulli\left(0.5\right)$ . Next w.l.o.g we prove that $Z$ independent from $X$ . Take any $b_{1},b_{2}\in\left\{ 0,1\right\}$ . In case $b_{1}=b_{2}$ then $$\begin{aligned} \mathbf{P}\left[X=b_{1},Z=b_{2}\right] & =\mathbf{P}\left[X=b_{1},X+Y=b_{1}\right]\\ & =\mathbf{P}\left[X=b_{1},Y=0\right]\\ & =\mathbf{P}\left[X=b_{1}\right]\mathbf{P}\left[Y=0\right]\\ & =\mathbf{P}\left[X=b_{1}\right]\cdot\frac{1}{2}\\ & =\mathbf{P}\left[X=b_{1}\right]\mathbf{P}\left[Z=b_{1}\right] \end{aligned}$$ so $X$ and $Z$ are independent. And in case that $b_{2}=1-b_{1}$ then similarly we get $$\begin{aligned} \mathbf{P}\left[X=b_{1},Z=b_{2}\right] & =\mathbf{P}\left[X=b_{1},X+Y=1-b_{1}\right]\\ & =\mathbf{P}\left[X=b_{1},Y=1\right]\\ & =\mathbf{P}\left[X=b_{1}\right]\mathbf{P}\left[Z=b_{1}\right] \end{aligned}$$ which conclude the proof.

  Now we prove the claim from above.

  Take two different variables $X_{1}=\left\langle \mathbf{s},\mathbf{x}_{1}\right\rangle \in\mathcal{S}$ and $X_{2}=\left\langle \mathbf{s},\mathbf{x}_{2}\right\rangle \in\mathcal{S}$ , and let $b_{1},b_{2}\in\left\{ 0,1\right\}$ , we want to show that $X_{1}$ and $X_{2}$ are independent. Obviously, we only care about the set of indices which are non-zero in the vectors $\mathbf{x}_{1}$ and $\mathbf{x}_{2}$ . Let these sets be $I_{1},I_{2}$ respectively. Denote $M_{0}=\sum_{i\in I_{1}\cap I_{2}}\mathbf{s}_{i}$ , $M_{1}=\sum_{i\in I_{1}\backslash I_{2}}\mathbf{s}_{i}$ and $M_{2}=\sum_{i\in I_{2}\backslash I_{1}}\mathbf{s}_{i}$ where summation over empty set is defined as the constant $0$ (all the operations in the proof are modulo 2). Note that any pair of variables from $\left\{ M_{0},M_{1},M_{2}\right\}$ are independent because they defined on mutually disjoint set of indices. Moreover, by the lemma from above, each of them if it is not the constant zero then is distributed as $Bernoulli\left(0.5\right)$ . When plugging-in the new notations we get $$\begin{aligned} \mathbf{P}\left[X_{1}=b_{1},\:X_{2}=b_{2}\right] & =\mathbf{P}\left[M_{0}+M_{1}=b_{1},\:M_{0}+M_{2}=b_{2}\right] \end{aligned}$$ Because $I_{1}\ne I_{2}$ there must be at most one variable from $\left\{ M_{0},M_{1},M_{2}\right\}$ which may be constant. Therefore, no matter which $M_{i}$ it is (if any), either $X_{1}$ or $X_{2}$ is a sum of two non-constant variables. W.l.o.g let it be $X_{1}$ , then it has the $M_{1}$ component which does not appear in $X_{2}$ and can flip the result of the sum with probability $0.5$ regardless of what happen anywhere else, therefore $X_{1}$ is independent of $X_{2}$ , and that completes the proof.

  Last one remark. In our settings the randomized algorithm either return a solution or says “failed”, if we have a polynomial procedure to validate a solution we can generalize the settings to the case when the algorithm always output a solution, and then we verify it to determine if it is “succeeds” case of “failed” case.
• Mar 24, 2017

  Comparing Chernoff - Hoeffding bounds

  tags:
  TheoryLunch Blog

  Computer science topics over lunch

  Comparing Chernoff - Hoeffding bounds

  Mar 24, 2017 • Eran Amar

  tags: $\newcommand{\lyxlock}{}$

  1 The Definitions and Motivation

  Hoeffding and Chernoff bounds (a.k.a “inequalities”) are very common concentration measures that are being used in many fields in computer science. A concentration measure is a way to bound the probability for the event in which the sum of random variables is “far” from the sum of their means.

  For example: suppose we have $n$ (fair) coins to toss and we want to know what is the probability that the total number of heads (i.e sum of indicators) will be “far” from half of the tosses (which is the sum of their means), say that either $\#heads>\frac{3}{4}n$ or $\#heads<\frac{1}{4}n$ . Chernoff’s and Hoeffding’s bounds can be used to bound the probability for that to happen.

  Both bounds are similar in their settings and their results, and sometimes it is not clear which bound is better for a given setting; maybe even both bounds are equivalent. In this post we will try to get an intuition which bound is stronger (and when).

  Note that each bounds have several common variations and we will discuss only the ones that are stated below. Let’s start with the first bound, a multiplicative version of Chernoff’s bound.

  Chernoff bound. Let $\left\{ X_{i}\right\} _{i=1}^{n}$ be a collection of independent random variables ranging in $\left[0,1\right]$ , denote $X:=\sum_{i=1}^{n}X_{i}$ and let $\mu:=\mathbb{E}\left[X\right]=\sum_{i=1}^{n}\mathbb{E}\left[X_{i}\right]$ , then for any $\epsilon\in\left(0,1\right)$ $$\mathbf{P}\left[\left|X-\mu\right|>\epsilon\mu\right]\le2\exp\left(-\epsilon^{2}\frac{\mu}{3}\right)$$

  Hoeffding bound. Let $\left\{ X_{i}\right\} _{i=1}^{n}$ be independent random variables ranging in $\left[a,b\right]$ where $a<b$ , denote $X:=\sum_{i=1}^{n}X_{i}$ and let $\mu:=\mathbb{E}\left[X\right]=\sum_{i=1}^{n}\mathbb{E}\left[X_{i}\right]$ , then for any $t$ : $$\mathbf{P}\left[\left|X-\mu\right|>t\right]\le2\exp\left(\frac{-t^{2}}{n\left(b-a\right)^{2}}\right)$$

  2 Writing Hoeffding’s bound as a Chernoff’s bound

  In first glance it seems that we cannot compare between the definitions, mostly because this small technical issue: Hoeffding’s bound allows the random variables to be in any close interval $\left[a,b\right]$ rather than $\left[0,1\right]$ . The solution for that is to scale and shift the variables to make them takes values in $\left[0,1\right]$ . We will start with exactly that, and then transform Hoeffding’s inequality into “Chernoff’s”.

  So for any $i\in\left[n\right]$ , denote $Y_{i}=\frac{X_{i}}{b-a}-a$ and let $Y=\sum_{i=1}^{n}Y_{i}=\frac{X}{b-a}-na$ . Observe that $\mathbb{E}\left[Y\right]=\frac{\mu}{b-a}-na$ . Clearly, now the range of the variables $\left\{ Y_{i}\right\} _{i=1}^{n}$ is $\left[0,1\right]$ so we can write Hoeffding’s bound in terms of $Y_{i}$ ’s

  $$\begin{aligned} \mathbf{P}\left[\left|X-\mu\right|>t\right] & =\mathbf{P}\left[\left|\frac{X}{b-a}-na-\frac{\mu}{b-a}+na\right|>\frac{t}{b-a}\right]\\ & =\mathbf{P}\left[\left|Y-\mathbb{E}\left[Y\right]\right|>\frac{t}{b-a}\cdot\frac{\mathbb{E}\left[Y\right]}{\mathbb{E}\left[Y\right]}\right]\\ & \le2\exp\left(\frac{-t^{2}}{\left(b-a\right)^{2}\mathbb{E}\left[Y\right]^{2}}\cdot\frac{\mathbb{E}\left[Y\right]}{3}\right)\\ & =2\exp\left(\frac{-t^{2}}{3\left(b-a\right)^{2}}\left(\frac{\mu}{b-a}-na\right)^{-1}\right)\\ & =2\exp\left(\frac{-t^{2}}{3\mu\left(b-a\right)-3a\cdot n\left(b-a\right)^{2}}\right) \end{aligned}$$

  Where the inequality in the middle of the process is where we apply Chernoff’s bound with $\epsilon=\frac{t}{\left(b-a\right)\mathbb{E}\left[Y\right]}$ , note that there is an implicit assumption here that $\epsilon<1$ .

  So we managed to formulate Hoeffding’s bound in Chernoff’s settings, ending with a very messy formula. In order to explore the behavior of the bounds further, we will simplify the analysis by adding another assumption: that $a=0$ .

  3 Which bound to used?

  Now we assuming that $a=0$ . The restriction from previous section that $\epsilon<1$ is now holds whenever $t<\mu$ , which make sense because Chernoff’s bound, as we defined it, only concerned with deviations that are up to $\mu$ .

  Next, the term we got for Hoeffding is simplified further into $$\mathbf{P}\left[\left|X-\mu\right|>t\right]\le2\exp\left(\frac{-t^{2}}{3\mu b}\right)$$ Now easy to compare between the bounds. With Hoeffding we can achieve bound of $2\exp\left(\frac{-t^{2}}{nb^{2}}\right)$ and with Chernoff $2\exp\left(\frac{-t^{2}}{3\mu b}\right)$ , it just left to compare between $nb$ and $3\mu$ . If $nb<3\mu$ then Chernoff is stronger and vise versa. We conclude that none of the bounds is always preferred upon the other (which is not so surprising), and the answer to the question “which one to use” depends on the parameters of the distribution in hand.

  4 Chernoff bound is tight!

  One surprising fact about Chernoff’s bound is that in some cases it is tight (i.e. gives a lower bound). The theorem below is taken from these notes (see the link for the full proof).

  Let $\left\{ X_{i}\right\} _{i=1}^{n}$ be i.i.d random variables ranging in $\left\{ 0,1\right\}$ with $\mathbf{P}\left[X_{i}=1\right]=p$ , denote $X:=\sum_{i=1}^{n}X_{i}$ and let $\mu:=\mathbb{E}\left[X\right]=np$ .
  If $p\le\frac{1}{4}$ , then for any $t>0$ $$\mathbf{P}\left[X-\mu>t\right]\ge\frac{1}{4}\exp\left(-t^{2}\frac{2}{\mu}\right)$$ If $p<\frac{1}{2}$ , then for any $t\in\left[0,n\left(1-2p\right)\right]$ $$\mathbf{P}\left[X-\mu>t\right]\ge\frac{1}{4}\exp\left(-t^{2}\frac{2}{\mu}\right)$$

  Note the differences between that regime and the one from the definition at the beginning of the post: in this theorem the variables are i.i.d (rather than only independent) and they are discrete.

  $\newcommand{\lyxlock}{}$

  1 The Definitions and Motivation

  Hoeffding and Chernoff bounds (a.k.a “inequalities”) are very common concentration measures that are being used in many fields in computer science. A concentration measure is a way to bound the probability for the event in which the sum of random variables is “far” from the sum of their means.

  For example: suppose we have $n$ (fair) coins to toss and we want to know what is the probability that the total number of heads (i.e sum of indicators) will be “far” from half of the tosses (which is the sum of their means), say that either $\#heads>\frac{3}{4}n$ or $\#heads<\frac{1}{4}n$ . Chernoff’s and Hoeffding’s bounds can be used to bound the probability for that to happen.

  Both bounds are similar in their settings and their results, and sometimes it is not clear which bound is better for a given setting; maybe even both bounds are equivalent. In this post we will try to get an intuition which bound is stronger (and when).

  Note that each bounds have several common variations and we will discuss only the ones that are stated below. Let’s start with the first bound, a multiplicative version of Chernoff’s bound.

  Chernoff bound. Let $\left\{ X_{i}\right\} _{i=1}^{n}$ be a collection of independent random variables ranging in $\left[0,1\right]$ , denote $X:=\sum_{i=1}^{n}X_{i}$ and let $\mu:=\mathbb{E}\left[X\right]=\sum_{i=1}^{n}\mathbb{E}\left[X_{i}\right]$ , then for any $\epsilon\in\left(0,1\right)$ $$\mathbf{P}\left[\left|X-\mu\right|>\epsilon\mu\right]\le2\exp\left(-\epsilon^{2}\frac{\mu}{3}\right)$$

  Hoeffding bound. Let $\left\{ X_{i}\right\} _{i=1}^{n}$ be independent random variables ranging in $\left[a,b\right]$ where $a<b$ , denote $X:=\sum_{i=1}^{n}X_{i}$ and let $\mu:=\mathbb{E}\left[X\right]=\sum_{i=1}^{n}\mathbb{E}\left[X_{i}\right]$ , then for any $t$ : $$\mathbf{P}\left[\left|X-\mu\right|>t\right]\le2\exp\left(\frac{-t^{2}}{n\left(b-a\right)^{2}}\right)$$

  2 Writing Hoeffding’s bound as a Chernoff’s bound

  In first glance it seems that we cannot compare between the definitions, mostly because this small technical issue: Hoeffding’s bound allows the random variables to be in any close interval $\left[a,b\right]$ rather than $\left[0,1\right]$ . The solution for that is to scale and shift the variables to make them takes values in $\left[0,1\right]$ . We will start with exactly that, and then transform Hoeffding’s inequality into “Chernoff’s”.

  So for any $i\in\left[n\right]$ , denote $Y_{i}=\frac{X_{i}}{b-a}-a$ and let $Y=\sum_{i=1}^{n}Y_{i}=\frac{X}{b-a}-na$ . Observe that $\mathbb{E}\left[Y\right]=\frac{\mu}{b-a}-na$ . Clearly, now the range of the variables $\left\{ Y_{i}\right\} _{i=1}^{n}$ is $\left[0,1\right]$ so we can write Hoeffding’s bound in terms of $Y_{i}$ ’s

  $$\begin{aligned} \mathbf{P}\left[\left|X-\mu\right|>t\right] & =\mathbf{P}\left[\left|\frac{X}{b-a}-na-\frac{\mu}{b-a}+na\right|>\frac{t}{b-a}\right]\\ & =\mathbf{P}\left[\left|Y-\mathbb{E}\left[Y\right]\right|>\frac{t}{b-a}\cdot\frac{\mathbb{E}\left[Y\right]}{\mathbb{E}\left[Y\right]}\right]\\ & \le2\exp\left(\frac{-t^{2}}{\left(b-a\right)^{2}\mathbb{E}\left[Y\right]^{2}}\cdot\frac{\mathbb{E}\left[Y\right]}{3}\right)\\ & =2\exp\left(\frac{-t^{2}}{3\left(b-a\right)^{2}}\left(\frac{\mu}{b-a}-na\right)^{-1}\right)\\ & =2\exp\left(\frac{-t^{2}}{3\mu\left(b-a\right)-3a\cdot n\left(b-a\right)^{2}}\right) \end{aligned}$$

  Where the inequality in the middle of the process is where we apply Chernoff’s bound with $\epsilon=\frac{t}{\left(b-a\right)\mathbb{E}\left[Y\right]}$ , note that there is an implicit assumption here that $\epsilon<1$ .

  So we managed to formulate Hoeffding’s bound in Chernoff’s settings, ending with a very messy formula. In order to explore the behavior of the bounds further, we will simplify the analysis by adding another assumption: that $a=0$ .

  3 Which bound to used?

  Now we assuming that $a=0$ . The restriction from previous section that $\epsilon<1$ is now holds whenever $t<\mu$ , which make sense because Chernoff’s bound, as we defined it, only concerned with deviations that are up to $\mu$ .

  Next, the term we got for Hoeffding is simplified further into $$\mathbf{P}\left[\left|X-\mu\right|>t\right]\le2\exp\left(\frac{-t^{2}}{3\mu b}\right)$$ Now easy to compare between the bounds. With Hoeffding we can achieve bound of $2\exp\left(\frac{-t^{2}}{nb^{2}}\right)$ and with Chernoff $2\exp\left(\frac{-t^{2}}{3\mu b}\right)$ , it just left to compare between $nb$ and $3\mu$ . If $nb<3\mu$ then Chernoff is stronger and vise versa. We conclude that none of the bounds is always preferred upon the other (which is not so surprising), and the answer to the question “which one to use” depends on the parameters of the distribution in hand.

  4 Chernoff bound is tight!

  One surprising fact about Chernoff’s bound is that in some cases it is tight (i.e. gives a lower bound). The theorem below is taken from these notes (see the link for the full proof).

  Let $\left\{ X_{i}\right\} _{i=1}^{n}$ be i.i.d random variables ranging in $\left\{ 0,1\right\}$ with $\mathbf{P}\left[X_{i}=1\right]=p$ , denote $X:=\sum_{i=1}^{n}X_{i}$ and let $\mu:=\mathbb{E}\left[X\right]=np$ .
  If $p\le\frac{1}{4}$ , then for any $t>0$ $$\mathbf{P}\left[X-\mu>t\right]\ge\frac{1}{4}\exp\left(-t^{2}\frac{2}{\mu}\right)$$ If $p<\frac{1}{2}$ , then for any $t\in\left[0,n\left(1-2p\right)\right]$ $$\mathbf{P}\left[X-\mu>t\right]\ge\frac{1}{4}\exp\left(-t^{2}\frac{2}{\mu}\right)$$

  Note the differences between that regime and the one from the definition at the beginning of the post: in this theorem the variables are i.i.d (rather than only independent) and they are discrete.
• Mar 16, 2017

  Introduction to Hypergraphs

  tags:
  TheoryLunch Blog

  Computer science topics over lunch

  Introduction to Hypergraphs

  Mar 16, 2017 • Eran Amar

  tags: $\newcommand{\lyxlock}{}$

  In this post, I will discuss basic definitions and theory related to hypergraphs and cut sparsifiers. I will start from the basic definitions of Graphs, then generalize these definitions to Hypergraphs. Later, I will discuss some related theory about Cut-Sparsifiers.

  1 Basic Definitions in Graphs

  Let’s start from the very beginning - graphs. Graphs are fundamental objects in computer science, and they are used to model relations between objects. Formally, a graph $G=\left(V,E,w\right)$ is a tuple of two sets, $V$ and $E$ and a function $w$ . $V$ is the set of vertices (i.e. the elements), usually denoted as $\left[n\right]:=\left\{ 1,2,..,n\right\}$ , and $E$ is a set of edges between those vertices (i.e. indicating which elements belongs to the relation the graph represents). An edge is just a set of two vertices, that is $e=\left\{ u,v\right\} \in E\subseteq\left\{ A\mid A\subseteq V,\:\:\left|A\right|=2\right\}$ . The function $w:\:E\rightarrow\mathbb{R}$ assigns a weight to every edge in the graph (we will focus on non-negative weights). For instance, we can model friendships in Facebook as a graph, the vertices will be all the users in Facebook, and there will be an edge between any pair of users if they are “friends” of each other.

  A cut in a graph $G$ is a partition of $V$ into two sets, $\left(S,V\backslash S\right)$ . We say that an edge $e=\left\{ u,v\right\}$ cross the cut defined by $S$ if and only if $S\cap e\ne\emptyset$ and $\left(V\backslash S\right)\cap e\ne\emptyset$ , that is, the edge “touches” both parts of the graph. Now we can talk about the weight of the cut, that is, if we have a cut that is defined by $S$ and a weight function $w$ then, $$w_{G}\left(S\right)=\sum_{\begin{array}{c} e\in E\\ e\cap S\notin\{\emptyset,e\} \end{array}}w(e)$$ in words, the weight of a cut is the sum of weights of all edges that cross it.

  2 Generalizing to Hyper-Graphs

  The model of a graph, as defined above, can’t be used to describe non binary relationships. For example, say we want to model Facebook groups via a graph; the only way to do that is by including in the set of vertices also any group in Facebook, and then connecting each user to any of its groups with an edge. The main issue with such construction is that the vertices represents both users and facebook-groups. In order to avoid such ambiguity we can use a richer model, we can use Hypergraphs.

  Hypergraphs are generalization of graphs in the sense that edges may be of arbitrary size. Meaning that now $E\subseteq2^{V}\backslash\emptyset=\left\{ A\mid A\subseteq V,\:\:\left|A\right|>0\right\}$ . Going back to our example, we can model Facebook groups by the graph $\left(V,E\right)$ when $V$ is the set of all users, and any group in Facebook will be an edge $e\in2^{V}$ (that notation stands for the Power Set of $V$ ) such that $e$ contains the users belongs to that group.

  Here is another example, in purely mathematical form, let $G=\left(\left[5\right],\left\{ \left\{ 4\right\} ,\left\{ 1,2,5\right\} \right\} \right)$ , which is a valid hypergraph over five vertices $\left\{ 1,2,..,5\right\}$ and with two edges: $\left\{ 4\right\}$ and $\left\{ 1,2,5\right\}$ .

  We also require the domain of the weight function to be $2^{V}\backslash\emptyset$ (rather than $E$ in the case of regular graphs). That requirement allows us to apply the same definitions of cuts and cut weight as they were phrased for regular graphs.

  Usually we are interested in family of hypergraphs with limited size, that is, hypergraphs where each edge is of size at most $r$ . Such hypergraphs are said to be r-uniform.

  3 Cuts Sparsification

  We move to define Cut-Sparsifiers: Let $\epsilon\in\left(0,1\right)$ . Given a hypergraph $H=\left(V,E\right)$ and a weight function $w$ , we say that a hypergraph $K=\left(V,E_{\epsilon},w\right)$ is $\epsilon$ -cut-sparsifier of $H$ if $$\forall S\subset V\qquad(1-\epsilon)\cdot w_{H}(S)\le w_{K}(S)\le(1+\epsilon)\cdot w_{H}(S)$$ and the set $E_{\epsilon}$ may be any set in $2^{V}\backslash\emptyset$ . The subscript “$H$ ” in the function $w_{H}\left(\cdot\right)$ indicates that the weights are over the hypergraph $H$ (and similarly with $w_{K}\left(\cdot\right)$ and $K$ ).

  It is not part of the definition, but the goal is to find cut-sparsifier that shrink the number of edges, that is, with the smallest $\left|E_{\epsilon}\right|$ possible. Whereas in regular graphs it is understood that we want to minimize the number of edges, in hypergraphs we should also consider the size of the edges in $E_{\epsilon}$ . That is, the quality of the sparsifier will be measured also by the size of the edges in the resulting graph. For that we define the total edge size of a hypergraph to be $\sum_{e\in E_{\epsilon}}\left|e\right|$ .

  4 Upper & Lower Bounds for Sparsification

  For regular graphs, we already know how to reduce the number of edges to be $\mathcal{O}\left(n/\epsilon^{2}\right)$ for any $\epsilon\in\left(0,1\right)$ . That is surprising, because no matter how many edges there are in the original graph, which can be up to $n^{2}$ , there is an algorithms that can reduce the number of edges to $\mathcal{O}\left(n/\epsilon^{2}\right)$ , while maintaining approximately the same weights for all possible cuts in the graph. Moreover, that algorithm find the cut-sparsifier in polynomial time.

  Later results show that this upper bound is also a lower bound, that is, there are graphs that cannot be reduced into less that $\mathcal{O}\left(n/\epsilon^{2}\right)$ edges.

  However, when working with hypergraphs it is not quite clear if one can reduce the total edges size of the graph to be even order of $n^{2}$ . Think about it, potentially the number of edges in a hypergraph can be $2^{n}$ , so $n^{2}$ can be considered quite small for such hypergraphs.

  An interesting result from 2015, showed lower bound of $O(\epsilon^{-2}n\cdot r)$ edges for r-uniform hypergraphs (assuming that $r>\log\left(n\right)$ ), which, when converting to total edges size is actually $O(\epsilon^{-2}n\cdot r^{2})$ .

  The question that one can ask now is, can we do better? Is there an algorithm that can construct a cut-sparsifier for hypergraphs with smaller total edges size?

  That is, actually, an open question which is studied nowdays.

  $\newcommand{\lyxlock}{}$

  In this post, I will discuss basic definitions and theory related to hypergraphs and cut sparsifiers. I will start from the basic definitions of Graphs, then generalize these definitions to Hypergraphs. Later, I will discuss some related theory about Cut-Sparsifiers.

  1 Basic Definitions in Graphs

  Let’s start from the very beginning - graphs. Graphs are fundamental objects in computer science, and they are used to model relations between objects. Formally, a graph $G=\left(V,E,w\right)$ is a tuple of two sets, $V$ and $E$ and a function $w$ . $V$ is the set of vertices (i.e. the elements), usually denoted as $\left[n\right]:=\left\{ 1,2,..,n\right\}$ , and $E$ is a set of edges between those vertices (i.e. indicating which elements belongs to the relation the graph represents). An edge is just a set of two vertices, that is $e=\left\{ u,v\right\} \in E\subseteq\left\{ A\mid A\subseteq V,\:\:\left|A\right|=2\right\}$ . The function $w:\:E\rightarrow\mathbb{R}$ assigns a weight to every edge in the graph (we will focus on non-negative weights). For instance, we can model friendships in Facebook as a graph, the vertices will be all the users in Facebook, and there will be an edge between any pair of users if they are “friends” of each other.

  A cut in a graph $G$ is a partition of $V$ into two sets, $\left(S,V\backslash S\right)$ . We say that an edge $e=\left\{ u,v\right\}$ cross the cut defined by $S$ if and only if $S\cap e\ne\emptyset$ and $\left(V\backslash S\right)\cap e\ne\emptyset$ , that is, the edge “touches” both parts of the graph. Now we can talk about the weight of the cut, that is, if we have a cut that is defined by $S$ and a weight function $w$ then, $$w_{G}\left(S\right)=\sum_{\begin{array}{c} e\in E\\ e\cap S\notin\{\emptyset,e\} \end{array}}w(e)$$ in words, the weight of a cut is the sum of weights of all edges that cross it.

  2 Generalizing to Hyper-Graphs

  The model of a graph, as defined above, can’t be used to describe non binary relationships. For example, say we want to model Facebook groups via a graph; the only way to do that is by including in the set of vertices also any group in Facebook, and then connecting each user to any of its groups with an edge. The main issue with such construction is that the vertices represents both users and facebook-groups. In order to avoid such ambiguity we can use a richer model, we can use Hypergraphs.

  Hypergraphs are generalization of graphs in the sense that edges may be of arbitrary size. Meaning that now $E\subseteq2^{V}\backslash\emptyset=\left\{ A\mid A\subseteq V,\:\:\left|A\right|>0\right\}$ . Going back to our example, we can model Facebook groups by the graph $\left(V,E\right)$ when $V$ is the set of all users, and any group in Facebook will be an edge $e\in2^{V}$ (that notation stands for the Power Set of $V$ ) such that $e$ contains the users belongs to that group.

  Here is another example, in purely mathematical form, let $G=\left(\left[5\right],\left\{ \left\{ 4\right\} ,\left\{ 1,2,5\right\} \right\} \right)$ , which is a valid hypergraph over five vertices $\left\{ 1,2,..,5\right\}$ and with two edges: $\left\{ 4\right\}$ and $\left\{ 1,2,5\right\}$ .

  We also require the domain of the weight function to be $2^{V}\backslash\emptyset$ (rather than $E$ in the case of regular graphs). That requirement allows us to apply the same definitions of cuts and cut weight as they were phrased for regular graphs.

  Usually we are interested in family of hypergraphs with limited size, that is, hypergraphs where each edge is of size at most $r$ . Such hypergraphs are said to be r-uniform.

  3 Cuts Sparsification

  We move to define Cut-Sparsifiers: Let $\epsilon\in\left(0,1\right)$ . Given a hypergraph $H=\left(V,E\right)$ and a weight function $w$ , we say that a hypergraph $K=\left(V,E_{\epsilon},w\right)$ is $\epsilon$ -cut-sparsifier of $H$ if $$\forall S\subset V\qquad(1-\epsilon)\cdot w_{H}(S)\le w_{K}(S)\le(1+\epsilon)\cdot w_{H}(S)$$ and the set $E_{\epsilon}$ may be any set in $2^{V}\backslash\emptyset$ . The subscript “$H$ ” in the function $w_{H}\left(\cdot\right)$ indicates that the weights are over the hypergraph $H$ (and similarly with $w_{K}\left(\cdot\right)$ and $K$ ).

  It is not part of the definition, but the goal is to find cut-sparsifier that shrink the number of edges, that is, with the smallest $\left|E_{\epsilon}\right|$ possible. Whereas in regular graphs it is understood that we want to minimize the number of edges, in hypergraphs we should also consider the size of the edges in $E_{\epsilon}$ . That is, the quality of the sparsifier will be measured also by the size of the edges in the resulting graph. For that we define the total edge size of a hypergraph to be $\sum_{e\in E_{\epsilon}}\left|e\right|$ .

  4 Upper & Lower Bounds for Sparsification

  For regular graphs, we already know how to reduce the number of edges to be $\mathcal{O}\left(n/\epsilon^{2}\right)$ for any $\epsilon\in\left(0,1\right)$ . That is surprising, because no matter how many edges there are in the original graph, which can be up to $n^{2}$ , there is an algorithms that can reduce the number of edges to $\mathcal{O}\left(n/\epsilon^{2}\right)$ , while maintaining approximately the same weights for all possible cuts in the graph. Moreover, that algorithm find the cut-sparsifier in polynomial time.

  Later results show that this upper bound is also a lower bound, that is, there are graphs that cannot be reduced into less that $\mathcal{O}\left(n/\epsilon^{2}\right)$ edges.

  However, when working with hypergraphs it is not quite clear if one can reduce the total edges size of the graph to be even order of $n^{2}$ . Think about it, potentially the number of edges in a hypergraph can be $2^{n}$ , so $n^{2}$ can be considered quite small for such hypergraphs.

  An interesting result from 2015, showed lower bound of $O(\epsilon^{-2}n\cdot r)$ edges for r-uniform hypergraphs (assuming that $r>\log\left(n\right)$ ), which, when converting to total edges size is actually $O(\epsilon^{-2}n\cdot r^{2})$ .

  The question that one can ask now is, can we do better? Is there an algorithm that can construct a cut-sparsifier for hypergraphs with smaller total edges size?

  That is, actually, an open question which is studied nowdays.

Tags