# Series expansion in terms of Hermite polynomials

I am trying to expand a polynomial in terms of orthogonal polynomials (in my case, Hermite). Maple has a nice built-in function for this, ChangeBasis.

Is there a similar function in Mathematica? And if not, where should I look for the algorithm?

-

The Hermite polynomials are orthogonal with respect to the inner product $$\langle f,g \rangle = \int_{-\infty}^{\infty} f(x)g(x)e^{-x^2} \, \mathrm dx.$$ Thus, the $n$-th coefficient can be computed using the inner product of your polynomial with the $n$-th normalized Hermite polynomial.

Example:

p[x_] = 1 + x + x^2 + x^3;
coeffs = Table[
Integrate[HermiteH[n, x]*p[x]*Exp[-x^2], {x, -∞, ∞}]/
Integrate[HermiteH[n, x]^2*Exp[-x^2], {x, -∞, ∞}],
{n, 0, 3}]
(* Out: {3/2, 5/4, 1/4, 1/8} *)

coeffs.Table[HermiteH[n, x], {n, 0, 3}] // Expand
(* Out: 1 + x + x^2 + x^3 *)

-

For polynomials, you don't need to do any integrals to find the expansion. Take a polynomial p and a list basis containing the basis functions. Then define a function that takes these two, identifies the variable x, and solves for the coefficients in basis that make the two polynomials equal in terms of their CoefficientLists:

expandPoly[p_, basis_, x_] :=
# /. First@Solve[CoefficientList[#.basis, x] == #2, #] &[
Array["a", Length[#]], #] & @

expandPoly[1 + x + 3 x^2 + 7 x^3, HermiteH[Range[4] - 1, x], x]

(* ==> {5/2, 23/4, 3/4, 7/8} *)


Edit

In response to belisarius: if you already know that you're only interested in a basis of HermiteH, you could incorporate that into the function and do away with the specification of the variable basis as follows:

expandPoly[p_, x_] := # /.
First @ Solve[
CoefficientList[#.HermiteH[Range[Length[#]] - 1, x],
x] == #2, #] &[Array["a", Length[#]], #] & @ CoefficientList[p, x]

expandPoly[1 + x + 3 x^2 + 7 x^3, x]

(* ==> {5/2, 23/4, 3/4, 7/8} *)


Edit 2

With the general function given as the first solution above, you can specify any set of polynomials that is known to form a basis for degree n or larger. This means the basis functions don't have to be orthogonal polynomials at all.

-
(+1) Very nice. I was just doing similar thing with another problem, but didn't think of it here. – Michael E2 Apr 14 '13 at 22:02
How to get rid of that Range@4? – Dr. belisarius Apr 14 '13 at 22:08
@belisarius See edit - I hardcoded the HermiteH part into the function in case you're only interested in those basis functions. Then the counting of terms is automatic. In the first version, I wanted to keep the list basis deliberately general so you can use other polynomials there, too. – Jens Apr 14 '13 at 22:13
Ahh... Algebra. I accidentally studied analysis in graduate school. :) – Mark McClure Apr 14 '13 at 22:17
If the expansion is known finite, one should be able to manage the upper bound automagically. – Dr. belisarius Apr 14 '13 at 23:08

The inner product for the Hermite polynomials, $$\langle f, g\rangle \int_{-\infty}^{\infty} f(x)\,g(x)\,e^{-x^2}\;dx\,,$$ has nice formulas for power functions (where $n=a+b$) and for the Hermite polynomials: \begin{align} \langle x^a, x^b \rangle = \langle x^n, 1\rangle &= \frac{1}{2} \left((-1)^n+1\right)\, \Gamma \left(\frac{n+1}{2}\right)\cr \langle H_n(x), H_n(x) \rangle &= \sqrt{\pi}\,2^n n! \cr \end{align}

These can be used to give a quick change of basis function for polynomials.

hermiteIP[f_, g_, x_] := With[{coeff = CoefficientList[f g, x]},
coeff.Table[1/2 (1 + (-1)^(-1 + n)) Gamma[n/2], {n, Length@coeff}]];

hermiteExpand[poly_, var_] /; PolynomialQ[poly, var] :=
Sum[hermiteIP[poly, HermiteH[n, var], var] H[n, var]/(Sqrt[Pi] 2^n n!),
{n, 0, Exponent[poly, var]}]


I used H[n, x] as a place holder for HermiteH[n, x].

hermiteExpand[(1 + x)^5, x]
(* 39/4 H[0, x] + 95/8 H[1, x] + 25/4 H[2, x] + 15/8 H[3, x] +
5/16 H[4, x] + 1/32 H[5, x]  *)

hermiteExpand[(1 + x)^5, x] /. H -> HermiteH
(* 39/4 H[0, x] + 95/8 H[1, x] + 25/4 H[2, x] + 15/8 H[3, x] +
5/16 H[4, x] + 1/32 H[5, x]  *)

% // Factor
(* (1 + x)^5 *)

-
Probably the best answer. It uses the orthogonality and speeds things up using the closed form. +1 – Mark McClure Apr 14 '13 at 22:24

Whenever I want to convert some polynomial expressed with respect to a certain basis in terms of another polynomial basis. my go-to algorithm is Salzer's algorithm. It's rather fast, since it relies only on recurrences. Here's a specialization of that algorithm for the case of monomial-Hermite conversion:

monomialToHermite[cofs_?VectorQ] := Module[{n = Length[cofs] - 1, a},
a[0, 0] = cofs[[n + 1]]; a[0, 1] = cofs[[n]];
a[1, 1] = cofs[[n + 1]]/2;
Do[
a[0, k + 1] = cofs[[n - k]] + a[1, k];
Do[
a[m, k + 1] = (m + 1) a[m + 1, k] + a[m - 1, k]/2,
{m, k - 1}];
a[k, k + 1] = a[k - 1, k]/2; a[k + 1, k + 1] = a[k, k]/2,
{k, n - 1}];
Table[a[m, n], {m, 0, n}]]


The algorithm as I presented it here uses an implicit two-dimensional array, a, to clearly show off the recurrence. The algorithm can be easily rewritten so that it uses only a pair or so of one-dimensional arrays, but I'll leave out that version for now.

Here's a test of Salzer's method:

monomialToHermite[{1, 1, 3, 7}]
{5/2, 23/4, 3/4, 7/8}

{1, 1, 3, 7}.x^Range[0, 3] == {5/2, 23/4, 3/4, 7/8}.HermiteH[Range[0, 3], x] // Expand
True

CoefficientList[(1 + x)^5, x] // monomialToHermite
{39/4, 95/8, 25/4, 15/8, 5/16, 1/32}

%.HermiteH[Range[0, 5], x] == (1 + x)^5 // Expand
True


(Other instances where I used Salzer's algorithm include this and this.)

-
(Of course, the algorithm is easily implemented even in languages that don't have symbolic capabilities.) – J. M. Apr 15 '13 at 0:44
Late answers are usually underestimated, +2. – Artes Apr 15 '13 at 19:29
Can't argue with that. Thanks @Artes! – J. M. Apr 16 '13 at 5:13
@seb, to clarify: are you expecting the terms of such an expansion to then come out as something like $c_{jk}H_j(x)H_k(y)$? – J. M. Jul 18 at 10:14
J.M., sorry for the unfortunate timing. I deleted my comment, because I remembered that this part is not the bottleneck in my calculation. But yes, I want expressions of the form $c_{kl} H_k(x)H_l(y)$. (In fact this is an XY problem and in the end I'm only interested at the values at $x=y=0$.) – sebhofer Jul 18 at 11:26

SolveAlways can find coordinates of a polynomial with respect to any given basis. You set up an equation, setting the given polynomial equal to a linear combination of your basis polynomials. This approach will work generally with any polynomial that is a linear combination of a given set of (linearly independent) polynomials.

poly = 1 + x + 3 x^2 + 7 x^3;       (* Jens' example *)

params = Table[C[i], {i, 0, Exponent[poly, x]}];
basis = HermiteH[Range[0, Exponent[poly, x]], x];

coeff = params /. First@SolveAlways[poly == params.basis, x]
(*  {5/2, 23/4, 3/4, 7/8}  *)


Update: Extension to ChangeBasis

Here is a version of the Maple ChangeBasis that works on polynomials (but not on infinite sums). It is based on the above idea.

ClearAll[changeBasis];
changeBasis[poly_, basisfamilies : {_, _} ..] :=
Module[{vars, families, degrees, params, basis, coeff},
vars = {basisfamilies}[[All, 1]];
families = {basisfamilies}[[All, 2]];
degrees = Exponent[poly, vars];
params = Outer[C, ##] & @@ (Range[0, #] &) /@ degrees;
basis =
With[{bfn = Times @@ MapThread[#1[#3, #2] &,
{families, vars, Array[Slot, Length[vars]]}]},
Outer[HoldForm[bfn] &, ##] & @@ (Range[0, #] &) /@ degrees
];
coeff = params /. First@SolveAlways[
poly == Total[params*ReleaseHold@basis, -1], vars];
With[{res = Total[coeff*basis, -1]},
Hold[res] /. HoldForm[v_] :> v
]
]


Some examples from the Maple docs. GegenbauerC[n, 1, x] is automatically simplified to ChebyshevU[n, x], which affect the form of the output.

changeBasis[1 + 2 x + 3 x^3, {x, GegenbauerC[#1, 1, #2] &}]
ReleaseHold[%] // Expand
GegenbauerC[n, 1, x]        (* Same as ChebyshevU[n, x] *)
(*
Hold[ChebyshevU[0, x] + 7/4 ChebyshevU[1, x] + 3/8 ChebyshevU[3, x]]

1 + 2 x + 3 x^3
ChebyshevU[n, x]
*)


Multivariable function:

changeBasis[y^2 + x^2 - 1, {x, LaguerreL[#1, 1, #2] &}, {y, ChebyshevU}]
ReleaseHold[%] // Expand
(*
Hold[21/4 (ChebyshevU[0, y] LaguerreL[0, 1, x]) +
1/4 ChebyshevU[2, y] LaguerreL[0, 1, x] -
6 (ChebyshevU[0, y] LaguerreL[1, 1, x]) +
2 (ChebyshevU[0, y] LaguerreL[2, 1, x])]

-1 + x^2 + y^2
*)

-
In the definition of changeBasis i think {x, y} should be {basisfamilies}[[All, 1]] or vars – Coolwater Jul 9 at 19:32
@Coolwater D'oh! Thanks! -- I probably missed replacing it by vars. – Michael E2 Jul 9 at 21:07

As another variation, here's another method based on repeated greedy division:

Reap[Fold[Block[{q, r}, {q, r} = PolynomialQuotientRemainder[#1, #2, x]; Sow[q]; r] &,
x^3 + x^2 + x + 1, HermiteH[Range[3, 0, -1], x]]][[-1, 1]] // Reverse
{3/2, 5/4, 1/4, 1/8}


Check:

%.HermiteH[Range[0, 3], x] == x^3 + x^2 + x + 1 // Expand
True

-