Initialization for Spectral Mixture Kernels

[jayhop]

Hello! This discussion refers to #1052.

Since Initialization is crucial for SM kernels, I wanted to look into it more.
I have seen that in gpytorch an initialization strategy is used which is based on the empirical spectrum.
I do understand the first part that computes the empirical spectrum using the discrete fourier transform.

N = train_x.size(-2)
emp_spect = np.abs(fft(train_y.cpu().detach().numpy())) ** 2 / N
M = math.floor(N / 2)

freq1 = np.arange(M + 1)
freq2 = np.arange(-M + 1, 0)
freq = np.hstack((freq1, freq2)) / N
freq = freq[: M + 1]
emp_spect = emp_spect[: M + 1]

But I do not understand the second part.
Why is the spectrum integrated cumulatively or what is inv_spec and why is emp_spec not used directly for initialization?

total_area = np.trapz(emp_spect, freq)
spec_cdf = np.hstack((np.zeros(1), cumtrapz(emp_spect, freq)))
spec_cdf = spec_cdf / total_area

a = np.random.rand(1000, self.ard_num_dims)
p, q = np.histogram(a, spec_cdf)
bins = np.digitize(a, q)
slopes = (spec_cdf[bins] - spec_cdf[bins - 1]) / (freq[bins] - freq[bins - 1])
intercepts = spec_cdf[bins - 1] - slopes * freq[bins - 1]
inv_spec = (a - intercepts) / slopes

[gpleiss]

This initialization script comes from @andrewgordonwilson 's original Matlab tutorial. See the link here: https://people.orie.cornell.edu/andrew/pattern/#Kronecker

[jayhop]

Thank you @gpleiss!