Gaussian quadrature improves on Newton-Cotes formulas by allowing the abscissas to vary along with the weights in order to integrate higher order polynomials.Can this idea be extended to numerical differentiation? To wit, can I choose a set $\{h_{i}\}$ such that $f'$ is much better approximated by a weighted sum of evaluations at $x+h_{i}$ than it could be at equally spaced data point? For example, maybe the following relation could hold for some $j$:\begin{align*}f'(x) = \sum_{i=1}^{n}\frac{w_i}{h_i}f(x+h_{i}) + \mathcal{O}(\left\| hf^{(n+j)}\right\|^{n+j})\end{align*}
Gaussian Numerical Differentiation
finite difference;quadrature
Yes. As you may know, numerical differentiation and integration is closely related to (polynomial) interpolation: The idea to approximately differentiate or integrate a given function is to approximate it with a function (often an interpolating polynomial) that can be differentiated or integrated exactly. For example, the standard central difference quotient formulas  for $f(x)$ can be derived from differentiating a quadratic  interpolating polynomial through the points $x-h, x, x+h$.The benefit ist that the error in approximating the derivative or integral is only determined by the error in approximating the function -- which is well understood for polynomial interpolation. In particular, it turns out that uniform interpolation points are a poor choice in general, and interpolation points based on roots of orthogonal polynomials are much better. In the context of quadrature, this leads to the different variants of Gaussian quadrature (Legendre, Chebyshev, Jacobi, Laguerre, Hermite...); in the context of differentiation, this is referred to as spectral collocation. Since it is the basis of spectral methods for solving differential equations, rather than a finite-difference quotient this is usually realized by a differentiation matrix that maps the values of $f$ at a selection of collocation points to the corresponding values of $f'$. (In contrast to the standard finite-difference matrices, these are usually dense; spectral methods are therefore global methods.) The deep relations between these concepts are explained in Trefethen, Lloyd N., Approximation theory and approximation practice, Other Titles in Applied Mathematics 128. Philadelphia, PA: Society for Industrial and Applied Mathematics (SIAM) (ISBN 978-1-611972-39-9/pbk). 305~p. (2013). ZBL1264.41001.In particular, chapter 21 is concerned with spectral collocation, and Theorem 21.1 is precisely the kind of error estimate you wrote.The standard choice for Gaussian quadrature are the roots of the Legendre polynomials, and you can in fact use the same points for differentiation as well (aptly called Legendre collocation); here's a Matlab script that sets up the corresponding differentiation matrix. More commonly used are the roots of Chebyshev polynomials leading to Chebyshev collocation, which is the basis of the remarkable Chebfun Matlab package (and its Julia child, ApproxFun.jl).