It has been recently established that for second-order consensus dynamics with additive noise, the performance measures, including the vertex coherence and network coherence defined, respectively, as the steady-state variance of the deviation of each vertex state from the average and the average steady-state variance of the system, are closely related to the biharmonic distances. However, direct computation of biharmonic distances is computationally infeasible for huge networks with millions of vertices. In this article, leveraging the implicit fact that both vertex and network coherence can be expressed in terms of the diagonal entries of pseudoinverse $\boldsymbol {{L}}^{2\dagger }$ of the square of graph Laplacian, we develop a nearly linear-time algorithm to approximate all diagonal entries of $\boldsymbol {{L}}^{2\dagger }$ , which has a theoretically guaranteed error for each diagonal entry. The key ingredient of our approximation algorithm is an integration of the Johnson–Lindenstrauss lemma and Laplacian solvers. Extensive numerical experiments on real-life and model networks are presented, which indicate that our approximation algorithm is both efficient and accurate and is scalable to large-scale networks with millions of vertices.