In the general case, if I have a symmetric, invertible (real-valued) matrix where each entry in the matrix has a separate error associated with it (or, in the more general case, where each value has a separate positive and negative error associated with it), how can I rigorously propagate these error values through the eigendecomposition to get error values for each of the eigenvalues and each entry in the eigenvector matrix? For my current application, it is as critical that I avoid overestimating the error as it is that I avoid underestimating it.
The propagation will depend on the diagonalization algorithm--which might be a black box--as well as the multivariate distribution of the errors. Pursuing an analytical solution therefore looks unpromising. Why not just compute an empirical distribution? That is, draw a large number of variants of the original matrix from the hypothesized error distribution and diagonalize them. Study the output distribution of the eigenvectors and eigenvalues.
There are some subtleties, because there will not be a definite matching among the lists of eigenvalues. For instance, in one iteration the sorted eigenvalues might be $(1.0, 0.99, 0.17)$ and in the next they might be $(1.01, 0.98, 0.17)$. Is the $1.01$ in the latter a slight variation of the $1.0$ in the former, or perhaps has the $0.99$ been perturbed into $1.01$ and the $1.0$ into $0.98$? It is impossible to know. Thus, you need to characterize the multivariate distribution of multisets of eigenvalues rather than $n$-tuples of eigenvalues.
The same problem attaches to the eigenvectors, but it gets worse, because there is no unique normalization of the eigenvectors. (They are determined only up to sign.) However, these problems are no different in nature than the ambiguities present in other geometric problems such as characterizing the directions of linear features in a plane (which can be given only up to a multiple of 180 degrees) and so should not present any additional conceptual challenge; they are just going to be a nuisance.
Here is an example of the empirical distributions of the sorted eigenvalues of a 4 by 4 matrix, using 2500 draws from the error distribution.
This scatterplot matrix also shows the lines y=x on each plot to emphasize the constraints imposed by sorting the four eigenvalues.