Polynomial Chaos Kriging

[jonathf]

Hello,

I am the maintainer of an uncertainty quantification library chaospy.
A user of my software has asked a question about using chaospy to do so kalled Polynomial Chaos Kriging. See jonathf/chaospy#300 for the details.

Following the idea of trying to implement this, would it be possible to provide Universal Kriging in PyKrige with linear trend polynomial? It would be really cool if it were posible to make PCK using PyKrige and chaospy together.

[MuellerSeb]

Hey @jonathf!
This sounds like a good idea for a work-together! in PyKrige, you can provide a list of callable drift functions for Universal kriging.
This is also possible in GSTools and PyKrige will be relying on that in the future.

• Here is an example for Universal kriging in GSTools.
• Here is an example in PyKrige (with the current version v1.5.1)

Some questions from my side:

• You are talking about trend and following the given links, I assume that you rather want Detrended Kriging, or simply ordinary Kriging applied to the residuals. Is the "trend" model known in your case?
• if the trend model is unknown (you want to fit a polynomial trend, right? I am not familiar with PCK), maybe regression kriging could be interesting: Example
• "linear trend polynomial" sounds like a contradiction to me. What do you mean by that?

So far from my side!
Cheers, Sebastian

[jonathf]

Yeah, a work-together on this was what I had in mind.

So it has been a few years since I did Spacial statistics and Kriging, so I appolgies for being a bit rusty on the termonology.
And yes linear trend polynomials is a contradictory in terms. But in my defence as a statistician, polynomial regression model sum(coeffs*polys) is still linear in terms of the coefficients.

I take my definition of Universal Kriging from Wikipedia. It defines it as having the intial mean field (or trend) in terms of a weighted sum of polynomials. The idea of polynomial chaos kriging (PCK) is that these weights and polynomials are provided from a specific polynomial chaos expansion (PCE) scheme. chaospy can create this PCE. Question is if we can provide this initial mean field to the Kriging model.

Doing de-trending was my initial thought of a way forward as well (though I used the term "modelling the discrepency between the PCE and the data" in the issue thread.) But I imagine that UK and detrended OK yields different results. The PCK paper assumes UK, so they did not go for the de-trending route.

I'll take a look at the resources you have provided and see if I can become wiser.

[MuellerSeb]

If you want to fit the coefficients of the polynomials during kriging, Universal-Kriging is the right choice.
If the coefficients are fitted beforehand, you should subtract the trend from the data and apply Ordinary-Kriging to the residuals.

ASFAIU you come up with a set of polynomials with the so-called "LARS" method, right?
Both, GSTools and PyKrige, provide the option to pass callable functions as drift-functions in Universal kriging, see here:

• gstools.krige.Universal -> drift_functions
• pykrige.uk.UniversalKriging -> functional_drift

[jonathf]

So in principle LARS does two things at the same time:

1. selects a set of polynomials from a candidate pool.
2. caclulate coefficients using regularized polynomial regression.

It is a search question where (1) and (2) are coupled. In practice one usually see these two steps defined together, not as separate steps.

From what you are telling me now, and reading the UQLab docs once more (I don't have access to the paper yet), I am realizing that I was wrong to assume that (1) and (2) are provided to kriging. Instead (1) and (2) is both done (as it is required), but only (1) is provided as input.

In other words, UK and drift-functions sounds exactly like what we are looking for.
I'll play around and see if I can get the UQLab example up and running.

[jonathf]

I replicated the 1-D intro example from UQlabs webpage. Here is a demo I plan to publish:
https://chaospy.readthedocs.io/en/docs/tutorials/gstools_kriging.html
A multi-dimensional model should be straight forward from here.

I've also gotten hold of the PCK paper UQLabs made.
It defines two approaches: sequential and optimal PCK. Based on the LARS implementation in scikit-learn, I was able to do the sequential approach. For the optimal approach, I am unsure if scikit-learn's LARS implementation provides everything needed. I'll look into it.

@MuellerSeb, I have a question.

The paper defines the kriging term of the equation as \sigma^2 Z(x) where Z(x) is zero mean, unit variance and stationary Gaussian process, and \sigma^2 is a Kriging parameter to be estimated from data.
In the code I used gstools.Gaussian(dim=1, var=1) as covariance model.
Am I right to assume that the \sigma^2 is estimated internally in gstools and the choice of var=1 referrs to Z(x)'s variance?

[jonathf]

I've updated the link to be address a 2D problem:
https://chaospy.readthedocs.io/en/docs/tutorials/gstools_kriging.html

There are a few caveats worth noting.
A) Lars is not Lars. It is a custom procedure referenced in the article with its own version of Jackknifing/LOO. Im using sklearn.LarsCV, as the latter requires a lot more work than I have time for.
B) Depending on the expansion output I am getting a lot of "Singular Matrix" errors without explaination. Is there a constrain I whould be aware of that is causing it?
C) I am getting better results when removing the constant polynomial. Is that something I sould be doing anyway?
D) results so far is not very impressive so far. Not sure if that is my implementation, the UQLab theory or both that is a t fault.

[MuellerSeb]

@jonathf

The paper defines the kriging term of the equation as \sigma^2 Z(x) where Z(x) is zero mean, unit variance and stationary Gaussian process, and \sigma^2 is a Kriging parameter to be estimated from data.
In the code I used gstools.Gaussian(dim=1, var=1) as covariance model.
Am I right to assume that the \sigma^2 is estimated internally in gstools and the choice of var=1 referrs to Z(x)'s variance?

I am a bit puzzled by this formula. I think it should be \sigma Z(x), since this is the correct formula to have a gaussian distribution with zero mean and a given standard deviation (See link). Assuming this, sigma^2 refers to var in GSTools and is fixed by setting it in the model. On the other hand, you also get the kriging variance, which can be interpreted as a confidence interval, since it is depending on the distance to the conditional points: Example. This variance is estimated.

Regarding your "Singular Matrix" problem: I fixed this by using the pseudo-inverse in the kriging system: GeoStat-Framework/GSTools#97. This is already present in the develop branch of GSTools and will be part of the next release: GeoStat-Framework/GSTools#110. Hopefully, we will release it in december.

You could also move this discussion to an issue in GSTools, since you apparently went with GSTools' kriging routines (which is also more future prove, since PyKrige will be refactored in order to use these under the hood in the future).

[jonathf]

Yeah the \sigma^2 thing is weird. I am wondering if it is a typo.

Good to know about the Singular Matrix issue. I will explore.

At this point, the way forward is to implement the "custom" LOO method. If/When I have some spare time, I might give it a try. I'll see where we get.
If this approach succeeds and it starts producing better fits as the paper claims.

Next time I have something to report, I'll start a new thread in the gstools repo.