Smooth low-rank approximation of random field with karhunen loeve

[baxmittens]

Hello all,

first of all, I think you provided an excellent package here. Thank you for your work. I have a problem using it. I'll try to describe.

The following code

using GaussianRandomFields
using Plots
import PlotlyJS
cov = CovarianceFunction(1, Exponential(.5))
pts_10 = range(0, stop=2, length=10)
pts_100 = range(0, stop=2, length=100)
grf_KL_10_p_10 = GaussianRandomField(cov, KarhunenLoeve(10), pts_10) # KL with 10 modes and 10 points
grf_KL_10_p_100 = GaussianRandomField(cov, KarhunenLoeve(10), pts_100) # KL with 10 modes and 100 points
grf_KL_100_p_100 = GaussianRandomField(cov, KarhunenLoeve(100), pts_100)# KL with 100 modes and 100 points
PlotlyJS.plot(PlotlyJS.GenericTrace[
        PlotlyJS.scatter(y=grf_KL_10_p_100.data.eigenfunc[:,1], x=pts_100,name="KL_10_p_100"), # this looks strange
        PlotlyJS.scatter(y=grf_KL_10_p_10.data.eigenfunc[:,1], x=pts_10,name="KL_10_p_10"),
        PlotlyJS.scatter(y=grf_KL_100_p_100.data.eigenfunc[:,1], x=pts_100,name="KL_100_p_100"),
        ]
        ) 

produces this image (KL_10_p_10 and KL_100_p_100 are looking fine) which looks strange.

The following code

xi = randn(randdim(grf_KL_10_p_100))
PlotlyJS.plot(PlotlyJS.GenericTrace[
        PlotlyJS.scatter(y=sample(grf_KL_10_p_100,xi=xi),x=pts_100,name="KL_10_p_100"),
        PlotlyJS.scatter(
            y=(repeat(xi',100,1).*grf_KL_100_p_100.data.eigenfunc[:,1:10])*sqrt.(grf_KL_100_p_100.data.eigenval[1:10]), 
            x=pts_100, 
            name="Random field generated by 10 modes from KL_100_p_100"
        )
    ]
    )

produces this image

of how the random field generated from KL_10_p_100 looks like and versus how I think the random field should look like.

Am I using the package wrong or is there some strange?

Additionally: I am not sure if things like

GaussianRandomField(cov, KarhunenLoeve(500), pts_100)

should work in general. Why can I generate a KL with 500 terms while I only provided a discretization with 100 DOFs?

Greetz Max

[PieterjanRobbe]

Hi Max,

What you're observing are numerical issues with the computation of the eigenvalues and eigenfunctions/vectors. We use the Nyström method to solve the generalized eigenvalue problem. A poor choice of quadrature points can lead to very strange results, as you observe. There's a note about this in the documentation for KarhunenLoeve:

julia> ?KarhunenLoeve
[...]
  │ Note
  │
  │  Techniqually, the KL expansion is computed using the Nystrom method. For nonstructured grids, we use a bounding box approach. Try
  │  using the Spectral method if this is not what you want.

  │ Warning
  │
  │  To avoid an end effect in the eigenvalue decay, choose nq^d ≥ 5n.
[...]

The "end effect" causes the unexpected shape of the eigenfunction (and sample of the random field) computed with only 10 terms (= 10 quadrature points in 1D), and evaluated in 100 points. To avoid this, you can specify the keyword nq = 50, as recommended in the documentation. This should solve both of the issues you're observing.

You can also observe convergence with respect to the eigenvalues with increasing nq as follows:

julia> p = plot()

julia> for f in 1:5
           plot_eigenvalues!(GaussianRandomField(cov, KarhunenLoeve(10), pts_100, nq=10f), label="nq=$(10f)")
       end

julia> display(p)

This gives the following result:

Finally, yes, you can generate a random field with more terms than discretization points, thanks to the Nyström method. In your case, using

grf_500_terms = GaussianRandomField(cov, KarhunenLoeve(500), pts_100)

the eigenfunctions will be generated in 500 equidistant points. The first 100 eigenvalues of that random field should agree with the eigenvalues of a field with 100 terms (provided that nq is large enough)

grf_100_terms = GaussianRandomField(cov, KarhunenLoeve(100), pts_100)
grf_100_terms_nq_500 = GaussianRandomField(cov, KarhunenLoeve(100), pts_100, nq=500)

julia> sum((grf_500_terms.data.eigenval[1:100] - grf_100_terms.data.eigenval).^2)
0.0012606162038396264

julia> sum((grf_500_terms.data.eigenval[1:100] - grf_100_terms_nq_500.data.eigenval).^2)
3.084966813077784e-11

[baxmittens]

Hi Pieterjan,

thank you for your quick response.

The nq keywords fixed the issues for me. I missed that while reading the docs.

I still cannot wrap my head around how things like

grf500 = GaussianRandomField(cov, KarhunenLoeve(500), pts_100, nq=500)

should work since the resulting matrix grf500.data.eigenfunc is of size 100x500, and with 100 points in space, it is not possible to represent the 500th mode properly. But for my use-case, I will always use way fewer modes than discretization points.

So thank you again. I close this one.

Greetz max

[PieterjanRobbe]

Great! You may be interested in reading section 3.2 in the reference below for more details on how this works.

Betz, W., Papaioannou, I., & Straub, D. (2014). Numerical methods for the discretization of random fields by means of the Karhunen–Loève expansion. Computer Methods in Applied Mechanics and Engineering, 271, 109-129.

[baxmittens]

@PieterjanRobbe I read this and I know the theory in general. But I don't see why this should work.

I played around with it a little bit

nq = 4000
grf_KL_500_p_100 = GaussianRandomField(cov, KarhunenLoeve(200), pts_100, nq=nq) # KL with 500 modes and 100 points
grf_KL_500_p_500 = GaussianRandomField(cov, KarhunenLoeve(200), pts_500, nq=nq) # KL with 500 modes and 500 points
grf_KL_500_p_1000 = GaussianRandomField(cov, KarhunenLoeve(200), pts_1000, nq=nq) # KL with 500 modes and 1000 points2
grf_KL_500_p_2000 = GaussianRandomField(cov, KarhunenLoeve(200), pts_2000, nq=nq) # KL with 500 modes and 2000 points

I set nq to different values ranging from 200 to 10_000 The outcome was like I did expect, independent of nq (provided it was sufficiently high, in general)

Here is the 60th eigenfunction for the case generated with a discretization with 100 points

while it should look like something like this

the 200th for 100 points looks like this

while for 2000 points it looks like this

I could imagine that both solutions are approximately correct in a sense of integrals but only the last one is a proper solution.

I would guess that in 1d, you need at least twice as many discretization points in space as KL-modes you want to display (something connected to the sampling theorem) in order to represent your eigenfunction in a sensible way.

But that is only a guess. What do you think?

[PieterjanRobbe]

That's just an aliasing effect. You're plotting discrete approximations to the continuous eigenfunction, represented with only a limited number of grid points. Here's an example representing the 6th eigenfunction with 5 and 100 grid points:

grf_KL_6_p_5 = GaussianRandomField(cov, KarhunenLoeve(6), range(start=0, stop=1, length=5), nq=100)
grf_KL_6_p_100 = GaussianRandomField(cov, KarhunenLoeve(6), range(start=0, stop=1, length=100), nq=100)
plot_eigenfunction(grf_KL_6_p_100, 6, label="100 points", legend = :outertopright)
plot_eigenfunction!(grf_KL_6_p_5, 6, label="5 points", markershape=:circle, markerstrokewidth=0)

gives

(Note: you may need to generate the random field several times, the orientation of the eigenfunctions is not guaranteed to be the same).

[baxmittens]

That is exactly what I was describing. But this is not just an aliasing effect. As you can clearly see, you are "missing frequencies". Your discretization is just not able to represent those high-frequency modes. Your resulting random field will not necessarily satisfy the desired properties. Surely, this is not an error or something like that on your side. That is just a thing the user has to be aware of.

[baxmittens]

The i-th (i=0,...,n) eigenfunction for a KL-expansion with an exponential kernel in 1d has the property that it has exactly i zero crossings. A property the coarse-scale solution is missing here, and thus, it is no proper approximate solution. I would like at least a warning in that case...or a hint in the documentation if you do not want that under any circumstances.

[PieterjanRobbe]

You can consider creating a PR with a proposal for a documentation change.

[baxmittens]

I'll do it on one of the following days. It is really just a minor thing but I think it could be a useful addition for some users. have a great day and thank you for responding once again.