# Orthogonalize polynomials with respect to Gagliardo seminorm?

For a function $f\colon [-1,1]\to\mathbb{R}$, the Gagliardo seminorm of $f$ is defined to be

$$|f| = \int_{-1}^1\int_{-1}^1 \frac{(f(x)-f(y))^2}{(x-y)^2}\, \mathrm{d}x\, \mathrm{d} y.$$

Given $(x,x^2,\ldots,x^n)$ for some $n$, I am interested in finding an orthonormal basis of these functions with respect to this seminorm. (We have to start at $x^1$ and not $x^0$ since the seminorm of the latter is zero).

The function Orthogonalize is normally how I do something like this. Indeed, in the help entry for Orthogonalize under Generalizations and Extensions, there is something similar to what I want to do. I can't figure out the notation/syntax for this norm.

What do I need to do?

• To orthogonalize a set of vectors/functions you need an inner product, not just a norm. – Marius Ladegård Meyer Apr 12 '16 at 13:05
• @MariusLadegårdMeyer Yeah, I ran into that too. We could use the polarization identity to construct the inner product from the norm, (hopefully). Or I was thinking that we could use $\langle f,g \rangle = \int_{-1}^1\int_{-1}^1 (f(x)-f(y))(g(x)-g(y)) / (x-y)^2\, \mathrm{d}x\, \mathrm{d} y$ for the inner product. – fred Apr 12 '16 at 13:10

Perhaps this?

ip[f_, g_, v_] := Module[{x, y, int},
int = (((f /. v -> x) - (f /. v -> y)) ((g /. v -> x) - (g /. v -> y)))/(x - y)^2;
Integrate[int, {x, -1, 1}, {y, -1, 1}]
];

Orthogonalize[x^Range[4], ip[##, x] &]
(*
{x/2, 1/2 Sqrt[3/2] x^2, 3/2 Sqrt[5/13] (-((2 x)/3) + x^3),
15/8 Sqrt[21/31] (-((14 x^2)/15) + x^4)}
*)


The inner product ip assumes that f and g are real-valued functions of the parameter v.

Update: For the sake of speed on polynomials

If you're dealing only with polynomials, one can get a speed-up by precomputing the integrals of monomials.

Block[{i, j, k},   (* protect i,j,k since they will be evaluated during the definition *)
Assuming[k ∈ Integers && k >= 0,
ippower[i_Integer?EvenQ] =        (* integral of even powers *)
Integrate[x^(2 k) , {x, -1, 1}] /. k -> i/2;
ipmoment[i_Integer?EvenQ, j_Integer?EvenQ] =
ippower[i] ippower[j];          (* the integral of x^i y^j is a product *)
]];
ipmoment[_Integer, _Integer] := 0;   (* odd powers integrate to zero *)

ipfast[f_, g_, v_] := Module[{x, y, int, coeff, mom},
int = (((f /. v -> x) - (f /. v -> y)) ((g /. v -> x) - (g /. v -> y)))/(x - y)^2;
mom = SparseArray@CoefficientList[Factor@int, {x, y}];  (* get coefficient list *)
coeff = mom["NonzeroValues"];                           (* get coeffs of monomials *)
mom = mom["NonzeroPositions"] - 1;                      (* get powers of monomials *)
coeff.(ipmoment @@@ mom)                               (* integrate and add results *)
]

b = x^Range[5];
res = Orthogonalize[b, ip[##, x] &]; // AbsoluteTiming // First
res2 = Orthogonalize[b, ipfast[##, x] &]; // AbsoluteTiming // First
(*
7.20899
0.021552
*)

res - res2 // Simplify
(*  {0, 0, 0, 0, 0}  *)


I seem to be getting different results depending on whether I use OP's suggested inner product in the comments, or if I use the polarization identity.

Before everything else, however, here is a routine by Velvel Kahan for evaluating the divided difference of a polynomial, based on Horner's method:

polynomialDividedDifference[poly_, {x_, a_, b_}] /; PolynomialQ[poly, x] :=
Module[{d = 0, y = 0},
Do[y = b y + Coefficient[poly, x, k]; d = a d + y,
{k, Max[0, Exponent[poly, x]], 1, -1}];
Expand[d]]


With this, here is the result of using the inner product

$$\langle f,g\rangle=\int_{-1}^1\int_{-1}^1\frac{(f(u)-f(v))(g(u)-g(v))}{(u-v)^2}\,\mathrm du\mathrm dv$$

Orthogonalize[x^Range[4],
Integrate[polynomialDividedDifference[#1, {x, u, v}]
polynomialDividedDifference[#2, {x, u, v}],
{u, -1, 1}, {v, -1, 1}] &]
{x/2, 1/2 Sqrt[3/2] x^2, 3/2 Sqrt[5/13] (-((2 x)/3) + x^3),
15/8 Sqrt[21/31] (-((14 x^2)/15) + x^4)}


Contrast this with the result of using the inner product

$$\langle f,g\rangle=\frac14\left(\left(\int_{-1}^1\int_{-1}^1\frac{(f(u)-f(v)+g(u)-g(v))^2}{(u-v)^2}\,\mathrm du\mathrm dv\right)^2-\left(\int_{-1}^1\int_{-1}^1\frac{(f(u)-f(v)-(g(u)-g(v)))^2}{(u-v)^2}\,\mathrm du\mathrm dv\right)^2\right)$$

based on the polarization identity:

Orthogonalize[x^Range[4],
(Integrate[polynomialDividedDifference[#1 + #2, {x, u, v}]^2,
{u, -1, 1}, {v, -1, 1}]^2 -
Integrate[polynomialDividedDifference[#1 - #2, {x, u, v}]^2,
{u, -1, 1}, {v, -1, 1}]^2)/4 &]
{x/8, (3 x^2)/16, (1036800 (-((719 x)/2880) + x^3))/3838561,
(17920000 (-((3177 x^2)/6400) + x^4))/33294461}


The OP will have to choose which basis is most appropriate for his situation.

• Hmmm, perhaps things are different because we only have a semi-norm instead of a norm. I will have to think about this. – fred Apr 12 '16 at 17:41
• Right, that's why I said that you need to do further investigations on which inner product is the right one for your application. – J. M.'s discontentment Apr 12 '16 at 17:43