The solution $x$ (and hence any simple response computed from it) is an analytic function of $\sigma$, hence it can be well approximated locally by rational functions. The response has poles at the eigenvalues of $-A$, which is why rational interpolation is needed. The following interpolation technique should work quite well with a number of factorizations approximately equal to the number of peaks in the responses.
Compute factorizations for three values of $\sigma$, namely the endpoints and the midpoint of your interval, and then calculate the responses of interest and their first two (or more) derivatives with respect to $\sigma$. Note that given the factorization, you can easily get the derivatives by multiple backsolves.
Then use piecewise rational Hermite interpolation, first a single rational, then two rational pieces, one in each subinterval. Assuming that the responses of interest decay to zero for $\sigma\to\infty$ (if not one needs to change the degrees in the interpolant), $d$ derivatives at the endpoints of a subinterval allow one to interpolate by a quotient of a degree $d-1$ polynomial divided by a degree $d$ polynomial. For $d=2$, this captures a single Lorentz peak exactly, for $d=4$ it captures a superposition of two. (If this is not enough, one can also consider multipoint Hermite interpolation, making better use of the computed data, at the expense of more programming work.)
Test the accuracy of the current interpolant by an additional exact calculation at the point where the last two interpolants differed most, using a new factorization and a new rational fit in the two intervals created by the new point, to decide when you can quit.
Unless your spectrum is very complicated, you'll need few factorizations only.