PolynomialBasis is not orthogonal?

[oarcelus]

I am running the following test

distribution = JointIndependent(
    [
        Uniform(1e-16, (1e-13 - 1e-16)),
        Uniform(1e-15, (1e-13 - 1e-15)),
        Uniform(1e-12, (1e-7 - 1e-12)),
        Uniform(1e-12, (1e-7 - 1e-12)),
        Uniform(0.0001, (0.001 - 0.0001)),
        Uniform(0.0001, (0.001 - 0.0001)),
        Uniform(0.001, (0.01 - 0.001)),
        Uniform(0.001, (0.01 - 0.001)),
        Uniform(0.01, (30.0 - 0.01)),
        Uniform(0.1, (400.0 - 0.1)),
        Uniform(1.0, (6.0 - 1.0)),
        Uniform(1.0, (6.0 - 1.0)),
        Uniform(0.02, (10.0 - 0.02)),
        Uniform(0.02, (10.0 - 0.02)),
        Uniform(800, (1600 - 800)),
    ]
)

poly = TotalDegreeBasis(distribution, 2)
distribution_r = cp.J(*[cp.Uniform(-1, 1) for _ in range(15)])
X, W = cp.generate_quadrature(2, distribution_r)

for i, p1 in enumerate(poly.polynomials):
    for j, p2 in enumerate(poly.polynomials):
        if i <= j:
            p1_eval = p1.evaluate(X.T)
            p2_eval = p2.evaluate(X.T)

            res = np.sum(p1_eval * p2_eval * W)
            print(f"Inner product of P_{i} and P_{j}: {res:.5e}")

The result I get is

Inner product of P_0 and P_0: 1.00000e+00
Inner product of P_0 and P_1: -1.71779e+00
Inner product of P_0 and P_2: -1.75087e+00

Am I doing something wrong? I would expect inner products between different Legendre polynomials to be close to 0 when doing a quadrature between -1 and 1.

[oarcelus]

Rushing is a danger. I was generating quadratures on the wrong distribution. I changed to random sampling and, appart from some numerical error I indeed see the diagonal elements very close to 1 and the rest 0.

distribution = JointIndependent(
    [
        Uniform(1e-16, (1e-13 - 1e-16)),
        Uniform(1e-15, (1e-13 - 1e-15)),
        Uniform(1e-12, (1e-7 - 1e-12)),
        Uniform(1e-12, (1e-7 - 1e-12)),
        Uniform(0.0001, (0.001 - 0.0001)),
        Uniform(0.0001, (0.001 - 0.0001)),
        Uniform(0.001, (0.01 - 0.001)),
        Uniform(0.001, (0.01 - 0.001)),
        Uniform(0.01, (30.0 - 0.01)),
        Uniform(0.1, (400.0 - 0.1)),
        Uniform(1.0, (6.0 - 1.0)),
        Uniform(1.0, (6.0 - 1.0)),
        Uniform(0.02, (10.0 - 0.02)),
        Uniform(0.02, (10.0 - 0.02)),
        Uniform(800, (1600 - 800)),
    ]
)
poly = TotalDegreeBasis(distribution, 1)

X = distribution.rvs(1000000)
M = poly.evaluate_basis(X)
for i in range(M.shape[1]):
    for j in range(M.shape[1]):
        if i <= j:
            res = np.sum(M[:, i] * M[:, j] / 1000000)

            print(i, j, res)