How to compute autocorrelation with exponential IIR filter efficiently

I have a discrete dataset (called V). I want to compute the autocorrelation of V at multiple time horizons. Normally I know that I can use the IFFT of the PSD. But in my case I need it with exponential kernels.

Let X_s = V * e^(-st) for t >= 0. (Convolute V with exponentially decaying causal filter).

I want to compute <V, X_s> for many s. Each s is a different mix of every frequency so I don’t see how to do this faster than len(V) * num_s time, even with the Fourier transform.

How does one do this more efficiently the way you would with a non-exponential autocorrelation function? In the real world my time series is multivariate and with dimension 3000 and shape(V) = 1,000,000 x 3000. So my output should have num_s * 3000 * 3000 size.

Thanks!
James