Create orthonormal expansion from glexindex output

[oarcelus]

I would like to know if there is a way to generate an orthonormal expansion of polynomials from the output of glexindex directly.

[jonathf]

Sure.

Something like the following?

expansion = cp.prod([
  cp.generate_expansion(order, dist[idx])[glexindex[idx]]
  idx in range(len(dist))
], axis=1)

This assumes dist is of type cp.J and glexindex is of type np.ndarray with len(glexindex) == len(dist).

[oarcelus]

Hi, thanks for the input. I have tried to compare this with generate_expansion function directly. And I see there are some differences.

alpha = cp.glexindex(
                start=0,
                stop=p + 1,
                dimensions=dimension,
                cross_truncation=q,
                graded=True,
                reverse=True,
            ).T

polyno = cp.prod(
    [
        cp.generate_expansion(
            config.order, config.distribution[idx], normed=True
        )[alpha[idx]]
        for idx in range(dimension)
    ],
    axis=1,
)

print(polyno)

polyno = cp.generate_expansion(
    p,
    config.distribution,
    normed=True,
    graded=True,
    reverse=True,
    cross_truncation=q,
)

print(polyno)

There is a difference between the two methods.

[1.7320508075688774q0 1.7320508075688774q0 1.7320508075688774q0 1.7320508075688774q0]

[1.0 1.7320508075688774q3 1.7320508075688774q2 1.7320508075688774q1 1.7320508075688774q0]

The coefficients look the same but the polynomial functions are not. Here p = 1 and dimension=4 and q=0.5

[jonathf]

Odd. Ill take a look.

[jonathf]

I was a bit too quick there. Try: axis=0.

[oarcelus]

[1.0 1.7320508075688774q0 1.7320508075688774q0 1.7320508075688774q0 1.7320508075688774q0]

[1.0 1.7320508075688774q3 1.7320508075688774q2 1.7320508075688774q1 1.7320508075688774q0]

I now recover the 1.0 component but still only q0's present. I think since the distribution is indexed distribution[idx] it thinks the problem only has a single random variable?

[jonathf]

Right. They are all in dim 1.

Here is a working example where each expanstion get assigned their own dim:

import chaospy as cp

qs = cp.variable(4)

dist = cp.J(cp.Normal(), cp.Normal(), cp.Normal(), cp.Normal())
alpha = cp.glexindex(
    start=0,
    stop=2,
    dimensions=len(dist),
    cross_truncation=0.5,
    graded=True,
    reverse=True,
).T
print(alpha)

polyno = cp.prod(
    [
        cp.generate_expansion(
            2, dist[idx], normed=True
        )[alpha[idx]](**{f"q0": qs[idx]})
        for idx in range(len(dist))
    ],
    axis=0,
)

print(polyno.round(4))

polyno = cp.generate_expansion(
    1,
    dist,
    normed=True,
    graded=True,
    reverse=True,
    cross_truncation=0.5,
)

print(polyno.round(4))

[oarcelus]

Hi! Now it works well. Thanks!

[oarcelus]

One last question about this. How can we verify that the polynomials are orthonormal?

My distribution is Uniform(-1, 1) as you suggested in a previous github post. The associated orthogonal polynomials are of the Legendre type. However, I see in the tables of polynomial families that the concrete distribution is Uniform(-1,1)/2, so maybe the norms of the polynomials are not 1?

[jonathf]

The quick way is to do this:

print(cp.E(cp.outer(polyno, polyno), dist))

Orthogonal polynomials should produce diagonal matrices here, orthonormal ones should produce identity/eye.

For papers to be more rigorous you would need to generate higher order Gaussian quadratures (X, W) and estimate each pair (i, j):

poly = polyno[i] * polyno[j]
coef = np.sum(poly(X) * W)

This is the same in theory, but has been shown to be a lot more numerically stable as the order increases.

[oarcelus]

Thank you!