6
$\begingroup$

Context

As a mean to understand the growth of structure in the universe, I am interested in characterising the curvature of random fields such as this one:

enter image description here

For this purpose I start with a PDF of the eigenvalues of the second derivative of the field (a measure of the local curvature). For a Gaussian Random Field this PDF reads

Pdf1= 2 Sqrt[2/π] (x1-x2) Exp[1/2 (-x1 (3 x1-x2)-x2 (3 x2-x1))]

and looks like this: 

rg = {{x1, -Infinity, Infinity}, {x2, -Infinity, x1}};
rgn = {{x1, -2, 2}, {x2, -2, x1}};
ContourPlot[Pdf1, Sequence @@ rgn // Evaluate,  ImageSize -> Small]

enter image description here

while when the field becomes non Gaussian (like the above), the PDF can be quite different:

enter image description here

My purpose is to use orthogonal polynomials to represent this Non Gaussian PDF.

Attempt

I have defined a scalar product 

Clear[int]; int[a_, b_] := 
 Integrate[ a b  Pdf1, Sequence @@ rg // Evaluate]

and a numerical integration version of it

 Clear[nint]; nint[a_, b_] := 
 NIntegrate[ a b  Pdf1, Sequence @@ rg // Evaluate]

I define my 2D polynomial

p = 2; pol0 = 
 Table[Table[x1^i x2^(p1 - i), {i, 0, p1}] // Flatten // Union,
    {p1, 0, p}] // Flatten // Join

{1, x1, x2, x1^2, x1 x2, x2}

 pol = Orthogonalize[pol0, int[#1, #2] &]

They look like this:

Map[ContourPlot[# Pdf1  , Sequence @@ rgn // Evaluate, 
   PlotPoints -> 15, PlotRange -> All] &, pol]

enter image description here

If I do the same thing numerically

 pol = Orthogonalize[pol0, intn[#1, #2] &,  Method -> "Reorthogonalization"];

I get the same answer.

But

If I try and find higher order polynomials numerically

 p = 6; pol0 = 
 Table[Table[x1^i x2^(p1 - i), {i, 0, p1}] // Flatten // Union,
    {p1, 0, p}] // Flatten // Join;
pol = Orthogonalize[pol0, intn[#1, #2] &,  Method -> "Reorthogonalization"];

I get this

(* {1.,1.81473 x1-0.804133,-1.53835 x1+2.37903 x2+1.73585,2.33148 x1^2-2.23663 x1+0.170409 x2-0.0991482,-2.84065 x1^2+4.36636 x1 x2+4.97902 x1-2.46156 x2-1.87669,1.84722 x1^2-5.47403 x1 x2-4.89841 x1+4.11307 x2^2+6.90647 x2+2.25076,0,3.15107 x1^2 x2+1.83935 x1^2-3.44839 x1 x2-3.23326 x1+0.212753 x2^2+0.21619 x2+0.676981,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0} *)

In other words the loss of accuracy in the orthogonalisation leads to zero higher order polynomials.

I have also tried Gramm Schmitt by hand (from somewhere on SE)

gs[vecs_, ip___] := Module[{ovecs = vecs},
Do[ovecs[[i]] -= Projection[ovecs[[i]], ovecs[[j]], ip], {i, 2, 
Length[vecs]}, {j, 1, i - 1}]; ovecs]; 
pol1 = gs[pol0, Function[
    NIntegrate[ ##  Pdf1, Sequence @@ rg // Evaluate]]];

but it yields the same loss of accuracy.

I have also tried Eigenvectors of the matrix of scalar products

   mat = ParallelTable[int[i, j], {i, pol0}, {j, pol0}];
   eigs = Eigensystem[mat // N];

But the orthogonal polynomials are not of increasing order:

    Map[ContourPlot[#  Pdf1, Sequence @@ rgn // Evaluate, 
   PlotPoints -> 15, PlotRange -> All] &, 
 pol = eigs[[2]].pol0/Sqrt[eigs[[1]]] // Chop]

enter image description here

Note that while orthogonal they are not orthogonal 'in the same direction'.

Question

How can I compute the higher order orthogonal polynomials accurately?

Side Question

One option would be to stick to symbolic evaluation of the scalar product, but it seems to take forever for higher order polynomials.

Is it possible to tell Orthogonalize to not Normalise the polynomials which symbolically seems to take longest?

Improve this question
$\endgroup$
10
  • $\begingroup$ May I ask why you try stick to a triangle? While orthogonal polynomial bases on triangles are known, it is much easier to deal with orthogonal polynomials on a square or a cube. Moreover, polynomials are not square-integrable over the whole $\mathbb{R}^2$ so Mathematica cannot compute the integrals correctly. I suggest to try (tensorproducts of) Hermite polynomials, multiplied by $\mathrm{e}^{-\frac{x^2}{2}}$; they are orthogonal right from the start. $\endgroup$
    – Henrik Schumacher
    Commented Dec 12, 2018 at 8:49
  • $\begingroup$ Because I can always rank my 2 eigenvalues so that x1>x2 (?). These polynomials are square-integrable over ℝ2 here because of the Gaussian part of the PDF? I have used in the past Hermite like polynomials (e.g. arxiv.org/pdf/1107.1863.pdf) but here I want to find the new polynomials for the Kernel of the eigenvalues. Or are you saying I should start with the Hermitte before orthogonalisation ? $\endgroup$
    – chris
    Commented Dec 12, 2018 at 8:51
  • $\begingroup$ Alternatively I could remove the x1>x2 constraint but then my integrant would have |x1-x2| instead of (x1-x2) $\endgroup$
    – chris
    Commented Dec 12, 2018 at 8:57
  • $\begingroup$ Oh, I overlooked Pdf1 in the integrand. Sorry. Towards the eigenvalues: Yes, you can order them, but I tried to say that computations might be easier without doing so. Of course, you would have to replace Pdf1 by the correct distribution on the whole $mathbb{R}^2$ (by multiplying by 1/2 and throwing in some Abs). This would allow you to use Hermite polynomials. $\endgroup$
    – Henrik Schumacher
    Commented Dec 12, 2018 at 8:58
  • 1
    $\begingroup$ @chris: PrecisionGoal and AccuracyGoal are a bit different: If set to a low value they give NIntegrate a bail out option to stop early, but they don't change the number precision NIntegrate uses for intermediate values in the internal calculations. By default NIntegrate will use fast $MachinePrecision numbers to perform those calculations but this (double) precision might not be enough to prevent catastrophic cancellation for high order oscillating polynomials. WorkingPrecision tells Mma to use arbitrary precision numbers internally, which should solve the cancellation issue. $\endgroup$
    – Thies Heidecke
    Commented Dec 12, 2018 at 13:16

1 Answer 1

Reset to default
5
$\begingroup$

Using Eigensystem is actually a great idea. The following should generate an orthonormal basis basis. Notice however that NIntegrate has its problems to compute the integrals accurately.

B = Threshold[Table[nint[a, b], {a, pol0}, {b, pol0}], 100 $MachineEpsilon];
{λ, U} = Eigensystem[B];
U = (U/Sqrt[λ]);
basis = U.pol0;

Edit

Having the Gram matrix B already, Orthogonalize seems to work quite well:

basis = Orthogonalize[IdentityMatrix[Length[B]], #1.B.#2 &].pol0

{1., -0.804133 + 1.81473 x1, 1.73585 - 1.53835 x1 + 2.37903 x2, -0.0991483 - 2.23663 x1 + 2.33148 x1^2 + 0.170409 x2, -1.87669 + 4.97902 x1 - 2.84065 x1^2 - 2.46156 x2 + 4.36636 x1 x2, 2.25076 - 4.89842 x1 + 1.84722 x1^2 + 6.90648 x2 - 5.47404 x1 x2 + 4.11307 x2^2}

and for p=4 enter image description here

Improve this answer
$\endgroup$
14
  • $\begingroup$ Thanks ! I think that's exactly what I need. It seems obvious retrospectively as always :-) $\endgroup$
    – chris
    Commented Dec 12, 2018 at 10:30
  • $\begingroup$ You're welcome! $\endgroup$
    – Henrik Schumacher
    Commented Dec 12, 2018 at 11:02
  • $\begingroup$ It seems to me this could be an option of the Orthogonalize function? Its fairly efficient. $\endgroup$
    – chris
    Commented Dec 12, 2018 at 13:02
  • $\begingroup$ I have a supplementary pb. For large matrices (> 60x60), even if the matrix is known symbolically, the basis starts to show inaccuracies. I have tried N[#1.B.#2,40]& say but it does not seem to help. Any suggestion? Thanks! $\endgroup$
    – chris
    Commented Jan 9, 2019 at 19:02
  • $\begingroup$ @chris Hm. I am actually clueless here. In fact, I haven't had time yet to wait for the Gram matrix B of size $60 \times 60$ to evaluate. It's clear that nint without further options will have very limited accuracy, but you say you computed B symbolically? That must have taken forever. But even with symbolic Gram matrix, I am not surprised that you will have accucary problems. $\endgroup$
    – Henrik Schumacher
    Commented Jan 10, 2019 at 9:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.