Using Householder reflectors to speed up random matrix experiments is not a new idea. It is well-known [1]}, [2]} that a Haar-distributed random orthogonal matrix can be factorized as a product of Householder reflectors. This leads to an efficient way of generating a random orthogonal matrix with \(\mathcal {O}(n^2)\)  operations (rather than the \(\mathcal {O}(n^3)\)  cost associated with a full QR decomposition on a Gaussian matrix). Householder reflectors have also been applied to reduce a Gaussian matrix to a particularly simple random bidiagonal form [3]}, [4]}. This clever factorization leads to an \(\mathcal {O}(n^2)\)  algorithm for simulating the spectrum densities of Gaussian and Wishart matrices. (Recall that a standard eigenvalue decomposition on a dense matrix requires \(\mathcal {O}(n^3)\)  work in practice.) The proposed HD algorithm differs from the previous work in that it is a truly matrix-free construction. With the progress of the dynamics, it gradually builds a recursive set of (random) Householder reflectors based on the current iteration vector \(x_t\)  and the history of the iterations up to this point. This adaptive, “on-the-fly” construction is essential for us to capture the correlation structures generated by the dynamics without fixing the matrix in advance.
