We simulate the spectral method and compare its empirical performance with the asymptotic limit given in [1]}. In our experiment, the measurement model is set to be \(y_i = \tanh \big (\vert a_i^{\hspace{-0.83328pt}\mathsf {T}}\xi \vert \big )\) . We compute the leading eigenvector \(x_1\)  by using the Krylov-Schur algorithm [2]}, which involves the repeated multiplication of \(D\)  with some vectors. With the forms of \(D\)  and \(A\)  given above, this algorithm can again be regarded as a special case of the general dynamics in (REF ). We use the HD algorithm for the simulation and show the results in Figure REF  for two different matrix dimensions: \(n = 10^3\)  and \(n = 10^5\) . Observe that, at \(n = 10^3\) , there is still noticeable fluctuations between the actual performance of the spectral method (shown as green dots in the figure) and the theoretical prediction (the blue line). To get a better match, the standard practice is to do many independent trials (2000 in our experiment) and average the results. This gives us the green curve in the figure. Averaging can indeed reduce statistical fluctuations, but there are still strong finite size effects, especially near the phase transition point. This is a case where the capability of the proposed HD algorithm to handle large matrices becomes particularly attractive: when we increase the dimension to \(n = 10^5\) , the empirical results match the theoretical curve very closely in any (typical) trial, with no need for averaging over repeated simulations. In terms of runtime, it takes the HD algorithm less than 4 seconds on average to obtain an extremal eigenvalue/eigenvector of \(D\)  for \(n = 10^5\) .
