Compute the L2 norm of a Field

Python usage
1

Hi!

We have a Field with index parameter in 2D and scalar value and we want to compute its L^2 norm:

\|X\|_{L^2}^2 = \int_\mathcal{D} X(\omega,{\bf t})^2 d{\bf t}

Is there a method in OT to compute this? Otherwise, it is necessary to take into account the function values and the volume of the vertices with Mesh.computeWeights(), which is doable more requires some computations.

Best regards,

Michaël

2

It is quite tricky. I can see several options:

  • Approximate the true squared L^2 norm using a quadrature rule in the spirit of the midpoint rule (it is exactly this rule for 1D meshes). For that, use mesh.computeWeights()
  • Compute the the exact squared L^2 norm of a piecewise interpolation of the field given its values at the vertices. For that, use mesh.computeP1Gram()

Here is a script to compare these options, as well as reference values obtained using adaptive quadratures of the model or the P_1 interpolation:

import openturns as ot

integrator = ot.GaussKronrod(65536, 1e-15, ot.GaussKronrodRule(ot.GaussKronrodRule.G7K15))

# The reference function
model = ot.SymbolicFunction("x", "sin(x+x^2/4)")
# An helper function
square = ot.SymbolicFunction("x", "x^2")

N = 10
domain = ot.Interval(-2.0, 2.0)
mesh = ot.IntervalMesher([N]).build(domain)

# Here we deform the mesh in order to have a nonuniform sampling of the domain
deform = ot.SymbolicFunction("x", "x+0.1*sin(pi_*x)")
mesh.setVertices(deform(mesh.getVertices()))

# g = mesh.draw()
# g.add(model.draw(domain.getLowerBound(), domain.getUpperBound()))
# Show(g)

###################################################
# Exact squared L2 norm using adaptive quadrature #
###################################################
refL2NormSquared = integrator.integrate(ot.ComposedFunction(square, model), domain)[0]
print("refL2NormSquared   =%.12g" % refL2NormSquared, "rel. err=%.4e" % abs(1-refL2NormSquared/refL2NormSquared))
##############################################################################
# Mid-point approximation of the true squared L2 norm using the mesh weights #
##############################################################################
values = model(mesh.getVertices())
field = ot.Field(mesh, values)
fieldSquared = ot.ValueFunction(square, mesh)(field)
weightL2NormSquared = fieldSquared.asPoint().dot(mesh.computeWeights())
print("weightL2NormSquared=%.12g" % weightL2NormSquared, "rel. err=%.4e" % abs(1-weightL2NormSquared/refL2NormSquared))
########################################################################
# Exact squared L2 norm of a P1 interpolation using the P1 Gram matrix #
########################################################################
gram = mesh.computeP1Gram()
gramL2NormSquared = values.asPoint().dot(gram*values.asPoint())
print("gramL2NormSquared  =%.12g" % gramL2NormSquared, "rel. err=%.4e" % abs(1-gramL2NormSquared/refL2NormSquared))
#########################################################################
# Exact squared L2 norm of a P1 interpolation using adaptive quadrature #
#########################################################################
myFieldSquaredFunction = ot.ComposedFunction(square, ot.Function(ot.P1LagrangeEvaluation(field)))
quadL2NormSquared = integrator.integrate(myFieldSquaredFunction, domain)[0]
print("quadL2NormSquared  =%.12g" % quadL2NormSquared, "rel. err=%.4e" % abs(1-quadL2NormSquared/refL2NormSquared))

I get the following results:

refL2NormSquared   =1.82317483123 rel. err=0.0000e+00
weightL2NormSquared=1.84393295915 rel. err=1.1386e-02
gramL2NormSquared  =1.74962050461 rel. err=4.0344e-02
quadL2NormSquared  =1.74962050461 rel. err=4.0344e-02

As you can see, there are several ways to do the job. The cheapest option (and here, the most accurate one) is based on the mesh weights. If your mesh has n vertices in dimension d and m simplexes, it needs a storage of size n for a computational cost \mathcal{O}(md). The option based on the Gram matrix needs a storage of size n^2 (no sparse representation of the Gram matrix) and a computational cost \mathcal{O}(md^2). Note that this last option is essentially a one-liner.

1 Like
3

thanks, I’m coding the norm for the mesh weights then