We thank the reviewer for his/her response. We answer here to his/her question.
Note that in the literature, different methods to generate orthogonal matrices were proposed - see e.g. [1,2,3,4,5,6,7,8,9].
***Direction matrix as Householder reflection.*** An efficient method to generate an orthogonal matrix is the following: At iteration $k \in \mathbb{N}$, we generate the direction matrix $G_k$ as a single random Householder reflection i.e.
\begin{equation}
G_k := I - 2 v_k v_k^\intercal
\end{equation}
where $I$ is the identity $\mathbb{R}^{d \times d}$ and $v_k$ is a vector uniformly sampled from the sphere (i.e $v_k \in \mathbb{S}^{d- 1}$). Note that the cost of this method consists of two parts:
- Generation of $v_k$: which consists in generating a Gaussian vector and in normalizing it.
- The outer product $v_k v_k^\intercal$ for which modern implementations exploit parallelization/vectorization.
The identity matrix can be generated and stored offline. Note that since it is very sparse, it can be stored using a sparse format (e.g. COO[10]). In this way, we can save resources in high-dimensional settings. The computation of the gradient approximation $g_k$ at iteration $k \in \mathbb{N}$ is
\begin{equation*}
g_k(x_k) := \frac{d}{\ell} \sum\limits_{i = 1}^\ell \frac{f(x_k + h_k G_k e_i) - f(x_k - h_k G_k e_i) }{2h_k}G_k e_i.
\end{equation*}
Thus, it "uses" only the first $\ell$ columns of $G_k$. Therefore, considering $\ell$ constant (as we proposed), we can store offline a (truncated) identity $I_{d, \ell}$ and compute the outer product truncating $v_k^\intercal$. We report the time cost of generating a set of directions using this procedure and random (Gaussian and Spherical) directions in the case $\ell = d$ (i.e. the most time-expensive setting). Mean and standard deviation are indicated using 500 repetitions.
***
d | Random Gaussian | Random Spherical | Householder
***
2 | 9.27e-7 $\pm$ 7.96e-7 | 5.49e-6 $\pm$ 2.05e-6 | 9.32e-6 $\pm$ 3.43e-6
4 | 1.30e-6 $\pm$ 7.21e-7 | 6.56e-6 $\pm$ 2.63e-6 | 1.12e-5 $\pm$ 5.79e-6
8 | 2.18e-6 $\pm$ 6.06e-7 | 8.01e-6 $\pm$ 5.32e-6 | 1.11e-5 $\pm$ 5.20e-6
16 | 5.69e-6 $\pm$ 1.61e-6 | 1.15e-5 $\pm$ 4.10e-6 | 1.18e-5 $\pm$ 7.20e-6
32 | 1.78e-5 $\pm$ 6.42e-6 | 2.49e-5 $\pm$ 1.33e-5 | 1.16e-5 $\pm$ 7.25e-6
64 | 6.58e-5 $\pm$ 7.03e-6 | 7.74e-5 $\pm$ 1.95e-5 | 1.62e-5 $\pm$ 3.79e-6
128 | 2.73e-4 $\pm$ 2.37e-5 | 2.98e-4 $\pm$ 2.45e-5 | 3.32e-5 $\pm$ 4.02e-6
256 | 1.26e-3 $\pm$ 2.79e-5 | 1.36e-3 $\pm$ 2.90e-5 | 1.20e-4 $\pm$ 1.04e-4
512 | 5.50e-3 $\pm$ 1.63e-4 | 5.91e-3 $\pm$ 1.22e-4 | 1.22e-3 $\pm$ 3.82e-4
1024 | 2.16e-2 $\pm$ 6.92e-4 | 2.41e-2 $\pm$ 7.35e-4 | 4.83e-3 $\pm$ 2.26e-3
2048 | 8.92e-2 $\pm$ 8.19e-2 | 1.04e-1 $\pm$ 1.03e-1 | 2.40e-2 $\pm$ 3.87e-2
***
The resources of the machine used to perform this experiment are described in Appendix C. Note that, our procedure is more expensive than random directions only in small dimensional settings (i.e. for $d \leq 16$). However, for higher dimensional cases, it scales better in time (i.e. it is cheaper than random directions). Moreover, note that the highest cost in this procedure is the outer product which can be efficiently computed in GPU (and thus the time cost can be reduced by exploiting it). We will extend Appendix D including this table and other details. Moreover, in order to complete the answer to your question (and the Reviewer's KbSU fourth question) we have to compute and compare the performance in function values using this algorithm instead of random directions. To do that, we repeated the numerical experiments plotting the computational time in the x-axis, and reported the results in the global response (see Figure 2) as requested by Reviewer kbSU (the choice of the parameters is described in Section 4 and Appendix C). As we can observe, orthogonal directions still provide better performance than random directions. Such results confirm the empirical observations of [11].
**References**
1. A. Genz. Methods for generating random orthogonal matrices.
2. F. Mezzadri. How to generate random matrices from the classical compact groups.
3. K. Choromanski, M. Rowland, W. Chen, and A. Weller. Unifying orthogonal monte carlo methods.
4. A. Hedayat and W. D. Wallis. Hadamard matrices and their applications.
5. Å. Björck. Numerics of gram-schmidt orthogonalization.
6. T. W. Anderson, I. Olkin, and L. G. Underhill. Generation of random orthogonal matrices.
7. A. Barvinok. Approximating orthogonal matrices by permutation matrices.
8. C. Rusu and L. Rosasco. Fast approximation of orthogonal matrices and application to pca.
9. C. Boutsidis and A. Gittens. Improved matrix algorithms via the subsampled randomized hadamard transform.
10. P. Virtanen et. al. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.
11. A. Berahas, L. Cao, K. Choromanski and K. Scheinberg. A Theoretical and Empirical Comparison of Gradient Approximations in Derivative-Free Optimization.