Proof Haar matrices satisfy JL lemma The Johnson-Lindenstrauss lemma says roughly that for any collection $S$ of $n$ points in $\mathbb{R}^d$, there exists a linear map $f:\mathbb{R}^d \rightarrow \mathbb{R}^k$ where $k = O(\log n/\epsilon^2)$ such that for all $x, y \in S$:
$$(1-\epsilon)\|x-y\|_2 \leq \|f(x)-f(y)\|_2 \leq (1+\epsilon)\|x-y\|_2$$

I've been trying to prove that the inequality above is satisfied if $f$ is picked using 
the Haar Transform as follows. Let D be a random $n\times n$ diagonal matrix with each diagonal element drawn uniformly and independently from $\{-1, 1\}$ (i.e. diagonal entries of $D$ are Rademacher random variables). Let $H$ be the standard $n\times n$ Haar matrix. Let finally $M$ be a random $k \times n$ binary matrix, such that each row $i$ has a single entry $M_{ij}$ equal to $1$ for $j$ picked uniformly at random from $\{1, \ldots, n\}$, and $M_{ij} = 0$ for all other $j$. In other words, $MH$ is equal to a random $k\times n$ matrix each of whose rows is a uniformly random row of $H$. Then the transform is defined as

$$
f(y) = MHDx
$$

When $f$ is picked as above, what is the probability that it satisfies the Johnson-Lindenstrauss condition for a fixed pair $x,y$ (when $f$ is given by a an appropriately scaled Gaussian matrix, the probability is $1 - \frac{1}{n^2}$)?
Your construction does not work in general for the value of $k$ given. 

Say $x = 0$ and $y= (1, 0, \ldots, 0)$ (or any other standard basis vector). Then $f(x) = 0$ and $HDy$ is a vector with $1 + \log_2 d$ nonzero entries. We have

$$
Pr[f(y) = 0] = \left(1 - \frac{1 + \log_2 d}{d}\right)^k \approx 1-O\left(\frac{k\log d}{d}\right),
$$
for $k \ll d$. Of course, in this case $\|f(y) - f(x)\|_2 = 0$ is a very bad approximation for $\|x - y\|_2 = 1$. So you almost certainly get infinite error.

This is why the Fast J-L Transform of Ailon and Chazelle uses a Walsh matrix rather than a Haar matrix: $WDx$, for $D$ picked as in your setup, is likely to have a lot of non-zero entries. This is related to the Uncertainty Principle. See also this CACM exposition.

More generally, Krahmer and Ward showed that if $M$ is a (possibly random) matrix that satisfies the restricted isometry property for submatrices of size at most $k$ (with sufficiently high probability), and $D$ is picked as above, then $f(x) = MDx$ satisfies the J-L condition for any $x, y$ with probability at least $1 - 2^{-\Omega(k)}$. This is in some sense a tight connection: a random matrix that satisfies the J-L condition for any $x, y$ with probability $1 - 2^{\Omega(k \log (n/k))}$, also satisfies the RIP property for submatrices of size up to $k$. 

For more on this subject check out the lecture notes from Jelani Nelson's course. In particular, lecture 17 contains the Krahmer-Ward theorem. Lectures 9-12 introduce the JL lemma with the sparse and fast variants.