Power and Rayleigh iterations are fundamentally different beasts, because the former uses only matrix-vector products with $A$, while the latter uses matrix-vector products with $(A+zI)^{-1}$. This is a very important point -- neither algorithms require explicit access to the matrix itself, but only certain actions of the matrix.
The choice to chose one over the other is heavily weighted by the application-specific complexity of each operation, which does not always scale according to the figures you have quoted.
Consider the case where $A$ is actually dense, but matrix-vector products can be efficiently computed. For example, gravity is a dense interaction, meaning that every body interacts with every other body. If the strengths of the interactions between $n$ bodies are tabulated as the $n\times n$ matrix $A$, then it will contain $n^2$ nonzero elements. For very large $n$, these $n^2$ elements might not even fit in memory. However, the matrix-vector product $Ax$ can be computed in just $O(n)$ complexity using the fast multipole method, without ever having to form the matrix. Accordingly, it is possible to compute $\lambda_\max(A)$ using power iterations, even as $n$ scales into the billions.
On the other hand, when $A$ is a sparse matrix, there often exist a reordering of the columns of the matrix such that $(A+zI)^{-1}x$ can be computed in near-linear or linear complexity. After an initial symbolic factorization step, only $\approx O(n)$ work is required to compute each Rayleigh iteration. In this case, there is very little reason to even consider power iterations.
Finally, I should note that power iterations is a suboptimal Krylov subspace method. Wherever power iterations can be used, the Lanczos iterations can do better. In fact, using exact arithmetic, Lanczos is guaranteed to converge after $n$ iterations.
