matrix with elements $\pm1$. @avron actually consider $H$ a multiple of a unitary matrix, but the special case $H$ Hadamard makes their proposal most similar to ours. Then, using simple random sampling (with replacement), they choose $N$ columns from this matrix and multiply this $n\times N$ matrix by an $n\times
n$ diagonal matrix with diagonal entries made up of independent symmetric Bernoulli random variables. The columns of this resulting matrix are the $U_j$’s. We are also multiplying a subset of the columns of a Hadamard matrix by a random diagonal matrix, but we do not select the columns by simple random sampling from some arbitrary Hadamard matrix.

The extra structure we impose yields beneficial results in terms of the variance of the randomized trace approximation, as the following calculations show. Partitioning $M$ into an $m\times m$ array of $N\times N$ matrices with $k\ell$th block $M^b_{k\ell}$, we obtain the following: $$\label{UjMUj}
\frac{1}{N}\sum_{j=1}^N
U_j'MU_j = \frac{1}{N} \sum
_{k,\ell=1}^m\sum_{j=1}^N
Y_{jk}Y_{j\ell} \beta_j'X_k
M^b_{k\ell} X_\ell\beta_j.$$ Using $Y_{jk}^2=1$ and $X_k^2 = I$, we have $$\begin{aligned}
\frac{1}{N}\sum_{j=1}^N
Y_{jk}^2 \beta_j'X_k
M^b_{kk} X_k \beta_j & = &
\frac{1}{N}\operatorname{tr}\Biggl(X_k M^b_{kk}
X_k \sum_{j=1}^N
\beta_j\beta_j' \Biggr)
\\
& = & \operatorname{tr}\bigl(M^b_{kk}X_k^2
\bigr)
\\
& = & \operatorname{tr}\bigl(M^b_{kk}\bigr),\end{aligned}$$ which is not random. Thus, 