Covariance rather than variance in Fixed "Noise" Likelihood

[CarloGraziani]

Hello.

I have a use case in which it would be convenient to represent experimental procedure uncertainty -- not noise, really, but case-to-case variability -- in the likelihood, by means of a full covariance matrix (rather than a variance vector) that applies to a certain subset of data.

In other words, data would come a segments of N numbers, each with its own, fixed NxN covariance matrix, assigned to through the likelihood function. The various segments are then connected through a more general GP model that joins them through the parameter space in which they are all embedded, and regression, prediction, and all that ensues.

There would be value in having these data covariances contain joint adjustable parameters for training, as well.

It's not clear to me that the Multitask likelihood accomplishes what I want, exactly, as it seems to want to fit some kind of joint likelihood covariance term to all the data, if I understand the documentation correctly. But I'd be happy to be corrected.

Does some such "covariance likelihood" exist as an obscure option in the code, or can it be forced into existence by doing something weirdly clever, or would it require a code change?

Thanks,

Carlo

[CarloGraziani]

On further reflection, I see haw to do this using the model rather than the likelihood: index the segments, and add a model component that has a Kronecker dependence on segment indices.

I do still think that there would be conceptual value in adding this functionality to the Likelihood part of the API, because it generalizes the notion of experimental uncertainty, currently only capable of representing uncorrelated measurement noise, to a broader class of experimental uncertainty that certainly is relevant to real experimental situations -- degrees of systematic case-by-case irreproducibility.

[CarloGraziani]

I may have confused the issue by mentioning multitask. This is not, strictly speaking a multitask-related modeling problem.

An experimental situation in which this problem might arise is one in which one must take repeated Raman spectral measurements of chemical mixtures at different concentrations. Even if a spectral measurement at a certain concentration is repeated, it is very difficult (for practical experimental setups, impossible) to repeat the same experiment within noise, because one cannot place the sample exactly at the same spot, because of variations in laser power, because of ambient temperature variations, because of quirks of the detector response. So the spectrum that one obtains in repeated measurements differ not (only) by noise, but by a draw from a random function, which it is not unreasonable to model using a zero-mean GP over wavenumbers, uncorrelated with other spectra that one might measure in other experiments.

This is what one would want to generalize in the likelihood. It isn't really multitask-related, so much as array-index sensitive. There's an "outer" array index that plays the role of the unique experiment identifier, and gets a Kronecker dependence in the additive likelihood term; and an "inner" index that plays the role of wavenumber -- in this particular example -- and that in general lends itself to kernel matrix generation. Repeatable experimental settings (concentrations here) are a function of the "outer" index. There are also parameter dependencies -- the hyperparameters of the individual experimental spectral covariances, in the example -- which can in general be global, or (for example) functions of the concentrations with global adjustable parameters.

Not sure how coherent that last reads. The basic idea is just to generalize measurement point noise -- random iid draws from a univariate normal -- to experimental process uncertainty -- random iid draws from a GP.

[gpleiss]

Gotcha. Is the experimental process uncertainty fixed or inferred?

[CarloGraziani]

We're trying to infer it in this case, or at least to model the uncertainty so as to get a valid overall model and extract the latent chemistry model. There's no real hope in this case of a first-principles model of detector behavior or laser power variations, or distribution of sample placement on the stage, for example. So having the ability to parametrize that stuff empirically, and then attempt to cross-validate that the model is adequate is probably the best that one can do.

[gpleiss]

If you are trying to infer the experimental process uncertainty you'll need to resort to some sort of approximate inference. In some sense you'll have a 2-layer DeepGP: the first GP "layer" models your latent process and your second GP "layer" will model your experimental process uncertainty. You could use this model already with GPyTorch (note that you'll have to add a GaussianLikelihood on top of this model which then also would include independent observational noise, but you could fix the noise parameter of this likelihood to be really small if you don't really want to include this).

[CarloGraziani]

Sorry, I misinterpreted your question. We don't want to infer experimental procedure uncertainty: In fact, the experimental procedure uncertainty is "nuisance" here, in exactly the same way that measurement noise is "nuisance". We want a model term for experimental procedure uncertainty that allows us to get it out of the way so that we can come to terms with the latent model describing the actual chemistry.

It's OK if the term is approximate/empirical. However, we do want to validate any model we try, say by leaving out spectra from training and checking whether their Mahalanobis chi-squared in their predictive distributions is reasonable. But that is admittedly different from inferring the procedure uncertainty.

[gpleiss]

After our discussion above, I think what you're looking for is FixedNoiseGaussianLikelihood, which is a heteroskedastic likelihood for training inputs.

[CarloGraziani]

I don't see why you think so. That's for latent point data, in which each datapoint is compounded with univariate zero-mean normal error representing measurement uncertainty. What I need is latent sets of points (e.g. spectra) compounded by zero-mean internally-correlated normal errors (offering some modeling possibilities) representing experimental procedure uncertainty. The FixedNoiseGaussianLikelihood doesn't seem to offer the possibility of non-diagonal noise blocks that apply to segments of the data. That is the feature request here.

It is actually something that can also be done at the model level, so I'm not beating the drums to have this implemented. I just think it's a nice generalization of ordinary uncorrelated noise, and as I said it has a sensible experimental interpretation. In a couple of months I expect to have a paper out on it (OK, such statements always offer hostages to fortune, I know...), and I can share it then.

[Balandat]

I implemented support for this locally at some point (years ago). I can try to dig up what I did then.

[Balandat]

Was literally just talking about this with @SebastianAment, seems like the time is ripe to unearth this.

[CarloGraziani]

I've run into a problem implementing what I want using the model interface.

What I would like, to represent procedure uncertainty, is very like batch independent multi-output GP, with each task representing a single experiment. However, I don't want the parameters controlling the kernel/mean in each task (i.e. producing the task's diagonal covariance block) to be independent of each other. Instead, I would like, say, outputscale[k] = f(k,p), where k is the task number, and f(k,p) is a user-supplied parametrized function of a parameter vector p. f() represents a dependency on the experimental conditions corresponding to task k. I would like to be able to fit p. And to do this for the other block kernel parameters besides outputscale as well.

In other words, the model has no correlation between blocks, but instead of having N_task * N_k parameters to fit, it only has to optimize the sum of the log-likelihoods with respect to the N_p parameters in the vector p.

I'm unsure how to accomplish this. Perhaps a new constraint class?

[gpleiss]

@CarloGraziani outputscale[k] = f(k,p) is already implicitly accomplished by the MultiTask kernel, or by the LMC model (essentially the same thing as MultiTask kernel). These models assume each task is a linear combination of multiple latent GPs. The coefficients of the linear combination essentially act like output scales as you describe.

See https://arxiv.org/abs/1106.6251 for a more thorough description.

[CarloGraziani]

Hi Geoff. Thanks for the response. I am familiar with that paper. I've also been reading the MT/LMC code rather carefully lately. But I'm afraid that the functionality that I want is not what you describe. For example, if I have 40 experiments described by 40 tasks, I want the 40 outputscales to be functions of (say) 3-4 optimizable parameters. I don't want 40 independent diagonal matrix parameters + 40*rank independent matrix factor parameters.

I actually now have subclassed versions of index_kernel.py, multitask_kernel.py, and lcm_kernel.py that should do what I want, and I'm in the process of testing them. Once they're working, I can make them available to the project by a PR, together with a test notebook, if there's any interest.

And, on a related note: I noticed that the missing values code only relates to Pyro functionality. This is quite sensible in the case of scalar GPs, since a missing value is simply a latent value with no effect. But for MT GPs, it can happen that the training data for an input setting is missing one or more task values, whereas others are available. One does not want to throw all the data away for that input just because a few output values are missing. There is a natural way to handle this in GP modeling, which is to mask those values, and drop the corresponding rows and columns of the covariance and mean. Since I need this too (some experiments can have incomplete results), I've subclassed multitask_multivariate_normal.py and multitask_gaussian_likelihood.py to support the required masking. That's also in testing.

[gpleiss]

@CarloGraziani sounds good. #2288 should hopefully add better support for missing values in the multitask setting.

[CarloGraziani]

Excellent, @gpleiss , @Turakar . Great minds think alike, but ours do too 😄. The PR produced by @Turakar looks much more complete than anything that I could have supplied for missing values.