Orthogonalize polynomials with respect to Gagliardo seminorm?

Question

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?

Answer 1

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,
   ip`power[i_Integer?EvenQ] =        (* integral of even powers *)
    Integrate[x^(2 k) , {x, -1, 1}] /. k -> i/2;
   ip`moment[i_Integer?EvenQ, j_Integer?EvenQ] =
    ip`power[i] ip`power[j];          (* the integral of x^i y^j is a product *)
   ]];
ip`moment[_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.(ip`moment @@@ 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}  *)

Answer 2

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)}

(same as in Michael's answer)

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.