Normalization of Jacobi polynomials in UQLab: uq_eval_jacobi versus Matlab function jacobiP

Dear UQWorld community,

I use UQLab to perform polynomial chaos expansions (PCEs). As I aim at further working with the resulting PCE coefficients in Matlab, I try to normalize the Jacobi polynomials in Matlab (function see jacobiP) such that they are consistent with the Jacobi polynomial definition in UQLab. How can I extract the proper normalizing constant?

I found a statement on the normalization of the Jacobi polynomials in the user manual “UQLab user manual Polynomial Chaos Expansions”, page 3, Table 1. What is the definition of the function B(.) there? To my knowledge, the Beta function needs two arguments, so should it be B(a,b) instead of B(a)*B(b)?

I figured out that the UQLab function uq_eval_jacobi gets the parameters (a,b) of the beta distribution as input. Do I have to use the parameters (a-1,b-1) as input to the Jacobi polynomial as it is defined in the Matlab function jacobiP then?

Exemplary, for order n=2, I compared the Matlab output (normalized to unit norm) and the UQLab output of the Jacobi polynomials for x ranging from -1 to 1.

j2_mat(x) = jacobiP(2,a-1,b-1,x)/normConst(2),

where normConst(2)= sqrt((2^(a+b-1)/(2\cdot2+a+b-1)*(Gamma(2+a)*Gamma(2+b)/(Gamma(2+a+b-1)*Gamma(2+1))). This is the definition of J_{a,b,2} in “UQLab user manual Polynomial Chaos Expansions”, page 3, Table 1, for the values a-1 and b-1.

j2_UQLab(x) = uq_eval_jacobi(2, x, [a,b], true).

Why do j2_mat and j2_UQLab lead to different results?

The same scaling procedure worked perfectly fine for the Legendre polynomials by setting normConst(2) = sqrt(1/(2*2+1)) as also specified in “UQLab user manual Polynomial Chaos Expansions”, page 3, Table 1.

What is the normalization constant of the Jacobi polynomials inherent in UQLab?

I would appreciate any help on that issue.

Thanks for taking time to reply.

Hi @klux,

Welcome to UQWorld. Thanks for creating this post here.
Regarding the manual, you are right, the function should be B(a,b) instead of B(a)B(b). In order to extract the constant, I would like to follow the definition of the Jacobi polynomials in Matlab:

\int_{-1}^{1} P^2_n(t;a,b) \, (1-t)^a(1+t)^b \,{\rm d}t = \frac{2^{a+b+1}}{2n+a+b+1}\frac{\Gamma(n+a+1)\Gamma(n+b+1)}{\Gamma(n+a+b+1)n!}.

Because the probability density function of the beta distribution is defined as

f(x;\alpha,\beta) = \frac{x^{\alpha-1} (1-x)^{\beta-1}}{B(\alpha,\beta)}, \quad \text{ for } x\in [0,1],

based on which we will modify the previous equation. First, we choose a = \beta-1 and b=\alpha -1 carry out a change of variable x = \frac{t+1}{2} and thus t = 2x-1. We obtain

\begin{aligned} &\int_{0}^{1} P^2_n(2x-1;\beta-1,\alpha -1) \, (1-(2x-1))^{\beta-1}(1+2x-1)^{\alpha -1} \,{\rm d}(2x-1) \\ =& 2^{\alpha+\beta -1}\int_{0}^{1}P^2_n(2x-1;\beta-1,\alpha -1) \, (1-x)^{\beta-1}x^{\alpha -1} \,{\rm d}x \\ =&2^{\alpha+\beta -1}\int_{0}^{1}P^2_n(2x-1;\beta-1,\alpha -1) \, \frac{x^{\alpha -1}(1-x)^{\beta-1}}{B(\alpha,\beta)}\,B(\alpha,\beta) \,{\rm d}x \\ =&2^{\alpha+\beta -1}\,B(\alpha,\beta)\,\int_{0}^{1} P^2_n(2x-1;\beta-1,\alpha -1) \, \frac{x^{\alpha -1}(1-x)^{\beta-1}}{B(\alpha,\beta)}\,{\rm d}x \\ =&2^{\alpha+\beta -1}\,B(\alpha,\beta)\,\int_{0}^{1} P^2_n(2x-1;\beta-1,\alpha -1) \, f(x;\alpha,\beta)\,{\rm d}x \end{aligned}

Combining this equation with the property of P_n on top of this page, we have

\begin{aligned} 2^{\alpha+\beta -1}\,B(\alpha,\beta)\,\int_{0}^{1} P^2_n(2x-1;\beta-1,\alpha -1) \, f(x;\alpha,\beta)\,{\rm d}x &= \frac{2^{\alpha+\beta-1}}{2n+\alpha+\beta-1}\frac{\Gamma(n+\alpha)\Gamma(n+\beta)}{\Gamma(n+\alpha+\beta-1)n!} \\ \int_{0}^{1} P^2_n(2x-1;\beta-1,\alpha -1) \, f(x;\alpha,\beta)\,{\rm d}x &=\frac{\Gamma(n+\alpha)\Gamma(n+\beta)}{(2n+\alpha+\beta-1)\Gamma(n+\alpha+\beta-1)n! B(\alpha,\beta)} \end{aligned}

The equation above means that L^2 norm of P_n(2x-1;\beta-1,\alpha-1) with respect to the weight function f(x;\alpha,\beta) is

\lVert P_n(2X-1;\beta-1,\alpha-1) \rVert_{2} = \sqrt{\frac{\Gamma(n+\alpha)\Gamma(n+\beta)}{(2n+\alpha+\beta-1)\Gamma(n+\alpha+\beta-1)n! B(\alpha,\beta)}},

which is the normalization constant that you want. Importantly, because we use a = \beta-1 and b=\alpha -1, the Jacobi polynomial is defined by P_n(2X-1;\beta-1,\alpha-1) instead of P_n(2X-1;\alpha-1,\beta-1).

Finally, based on the formula above, we can use the Matlab function jacobiP to obtain the Hermitian basis with respect to the beta distribution f(x; a,b). For example for degree n=2, we have

n = 2;
phitmp = jacobiP(n,b-1,a-1,2*x-1);
normConst = gamma(n+a)*gamma(n+b)/(2*n+a+b-1)/gamma(n+a+b-1)/factorial(n)/beta(a,b);
phiuqlab= phitmp/sqrt(normConst);

Please let me know if this solves the problem.

Xujia

Dear Xujia,

thanks a lot for your detailed answer and taking time to respond! That helps a lot and clarifies the confusion about the normalization constant of the Jacobi polynomials inherent in UQLab compared to the Matlab function jacobiP.
Thanks again for your support.

Best regards!