A least squares approach to Whitney forms
Ludovico Bruni Bruno Giacomo Elefante


No Train but Gain:
Language Arithmetic for training-free Language Adapters enhancement
    [
    April 29, 2024
======================================================================================




    In this work we describe and test the construction of least squares Whitney forms based on weights. If, on the one hand, the relevance of such a family of differential forms is nowadays clear in numerical analysis, on the other hand the selection of performing sets of supports (hence of weights) for projecting onto high order Whitney forms turns often to be a rough task. As an account of this, it is worth mentioning that Runge-like phenomena have been observed but still not resolved completely. We hence move away from sharp results on unisolvence and consider a least squares approach, obtaining results that are consistent with the nodal literature and making some steps towards the resolution of the aforementioned Runge phenomenon for high order Whitney forms.






§ INTRODUCTION


Interpolation and approximation are crucial aspects of the numerical analysis of both physical- and engineering-oriented problems <cit.>. These techniques make it possible to translate observations and measurements into mathematical objects on which computations may be performed. They are hence a ground step of all numerical methods <cit.>.



Interpolating a scalar function means representing it by a simpler object, 
typically a polynomial, 
asking that they share the same values (measurements) on a set of nodes. This is known as Lagrange interpolation <cit.> and applies to physical quantities that have a natural pointwise meaning, such as the temperature. The reliability of this process is essential and in this direction much work has been done to avoid undesired behaviours such as the Runge phenomenon <cit.>. 

When vector fields (e.g., forces) or tensor fields (e.g., metric tensors) come at play the situation becomes more complicated.  
These quantities are measured via diffused data, such as circulations along lines and fluxes across surfaces, and one hence has to replace nodal evaluations by currents <cit.>. 
This is formalised in the mathematical language of differential forms <cit.>, which one may intuitively think of as objects that can be integrated. The domain of integration must be consistent with the order of the form: for example, the electric field has a well defined line integral, whereas the magnetic field has a meaningful surface integral <cit.>, unveiling the intrinsic difference between such objects. 
To retrieve a concrete extension of Lagrange interpolation to differential forms, 


one may invoke the theory of weights <cit.>. 


Weight-based interpolation 
consists in associating a differential form in an infinite dimensional space (such as that of smooth forms Λ^k, see <cit.>) with a differential form in a finite dimensional space (such as that of high order Whitney forms 𝒫_r^- Λ^k, see <cit.>) by matching their integrals 
<cit.>. The main drawback of this approach, which is inherited from multivariate Lagrange interpolation, lies in the complexity of determining whether a specific collection of supports induces a well defined interpolator. Sharp proofs in this direction are often heavy and confined to very specific cases <cit.>; this opens the gap between redundance <cit.> and minimality <cit.>. 

In place of exact interpolation, one may consider least squares approximation <cit.>. Instead of matching the degrees of freedom, one looks for the object that minimises a distance (typically the 2-norm) from the observations. This is a classical technique in the approximation of functions <cit.>. In the context of  differential forms it found applications to fluids <cit.> and electromagnetism <cit.>. In this work we pursue this approach, giving a general way for constructing a least squares framework for high order Whitney forms. We observe interesting behaviours:

    
  * if the total polynomial degree is fixed and hence the redundant and the minimal sets under consideration are chosen as in <cit.> and <cit.>, performances are very similar between least squares approximation and interpolation. This happens both with an easy to interpolate differential forms (see Fig. <ref>, left) and with Runge-type differential forms, such as those investigated in <cit.>;
    
  * if the total polynomial degree is fixed, increasing the number of observations over a certain quantity (depending on the degree itself) does not improve the least squares approximation error, which stabilises below the interpolation error (see Fig. <ref>, left). This is consistent with nodal observations <cit.>; as a consequence, this is not a possible solution to the Runge phenomenon;
    
  * combining the above strategies, carefully increasing the number of observations along with the polynomial degree leads also to convergence of Runge-like forms even for uniform (in the sense of <cit.>) distribution of supports (see Fig. <ref>, right), thus defeating this undesirable effect.


The paper is organised as follows. In Section <ref>, we recall the main features of Whitney forms and an algorithm for the computation of bases. We introduce the concept of weight and the relative interpolation and least squares approximation problems. In Section <ref>, we translate the problem to the more common language of vector calculus and explain how usual quadrature formulae can be applied in this context. Section <ref> contains numerical experiments and verification. In Section <ref>, we gather conclusions and overview possible further developments.




§ APPROXIMATION BY WHITNEY FORMS
 

Differential forms are nowadays a widely spread language in physics and numerical methods, see e.g. <cit.>. A peculiar finite dimensional space that emerges in applications is that of Whitney forms <cit.>, also known as trimmed polynomial differential forms. They play a relevant role in computational electromagnetism <cit.> (see also the whole collection of Bossavit's Japanese papers, such as <cit.>) and, more generally, in the construction of families of finite elements <cit.>. In particular, they parametrise (in the sense of <cit.>) the first family of Nédélec finite elements <cit.>. What makes Whitney forms special is their deep relationship with the underlying geometry of the element.



 §.§ From low-order to high-order Whitney forms


Let T ⊂^n be an n-simplex spanned by vertices {v_i }_i=0^n and endowed with barycentric coordinates {λ_i }_i=0^n, i.e., the collection of affine polynomials such that λ_i (v_j) = δ_i,j. We denote by F_σ be the k-subsimplex of T associated with the increasing mapping σ: { 0, …, k }→{0, …, n }, that is the k-simplex generated by {v_σ(0), …, v_σ(k)}, see e.g. <cit.>. If we let 𝔖_k^n ≐{σ: { 0, …, k }→{0, …, n }} denote the set of all above increasing mappings, a basis for degree 1 (or lowest order) Whitney forms is given by the “Lagrange basis” dual to integration over F_σ, as σ ranges over 

    ω_σ:   ∫_F_σ'ω_σ = δ_σ, σ' ,

possibly up to a rescaling factor k!, see <cit.>. An explicit representation of Whitney forms is thus

    ω_σ = ∑_i=0^k (-1)^i λ_σ(i)λ_σ(0)∧…∧λ_σ(i)∧…∧λ_σ(k),

the hat denoting that the underlying term has been removed. We set 

    _1^- Λ^k (T) ≐span{ω_σ|σ∈𝔖_k^n } ,

and from (<ref>) we deduce that the dimension of this space is n+1k+1. This implies that _1^- Λ^k (T) is a proper subspace of the space of all degree 1 polynomial differential forms _1 Λ^k (T) unless k = 0.



There are several ways to extend the definition of the above space to r > 1, see e.g. <cit.>, <cit.> and <cit.>. If _r Λ^k (T) denotes the space of k-forms whose coefficients are n-variate degree r polynomials, we may think of the space 𝒫_r^- Λ^k (T) as the image of the product ℙ_r-1 (T) ×_1^- Λ^k (T) →_r Λ^k (T), namely

    𝒫_r^- Λ^k (T) = {η∈Λ^k (T) |η = p ·ω,   p ∈ℙ_r-1 (T),  ω∈_1^- Λ^k (T) } .

Since the dimensional inequality

    𝒫_r^- Λ^k (T) = r+k-1kn+rn-k≤n+rnnk = 𝒫_r Λ^k (T)

holds, as proved in <cit.>, the inclusion _r^- Λ^k (T) ⊆_r Λ^k (T) is strict unless k = 0.

The homogenization formula for a degree r polynomial allows to write a basis of ℙ_r-1 (T) in terms of monomials in barycentric coordinates and multi-indices ℐ (n+1, r-1) of length n +1 and weight r-1, that means λ^α≐∏_i=0^n λ_i^α_i, with | α | ≐∑_i=0^n |α_i | = r-1. Merging this with (<ref>) and (<ref>), it follows that 𝒫_r^- Λ^k (T) is generated by forms λ^αω_σ, with α∈ℐ (n+1,r-1). 
Also, (<ref>) shows that this is necessarily a system of generators and not a basis for 𝒫_r^- Λ^k (T) for k > 0 and r > 1. To retrieve a basis we shall invoke, for instance, <cit.>.


 
    For each σ∈𝔖_k^n, define ℐ_σ≐{α∈ℐ(n+1,r-1) |α_i = 0   ∀ i < σ(0) }. The set
    
    {λ^αω_σ|σ∈𝔖_k^n, α∈ℐ_σ}

    is a basis for _r^- Λ^k (T).


Extensions of Whitney forms to other geometries are also considered in literature <cit.>.



 §.§ Weights
 

We saw in (<ref>) that integrals over k-subsimplices of T give natural degrees of freedom for _1^- Λ^k (T). Since they are too few when r > 1, i.e. in the case of _r^- Λ^k (T), enriched collection of supports 𝒮≐{ s_1, …, s_N }⊂ T are needed. These objects are called small simplices[Historically, only the set s_(α, σ), here defined in the sequel, was named so <cit.>.] 
and integration of a differential k-form ω over elements of the set 𝒮 is called a weight <cit.>:

    ∫_s_iω .

The collection 𝒮 of small simplices is unisolvent for 𝒫_r^- Λ^k (T) if the Vandermonde matrix

    V_i,j≐∫_s_iω_j,

written with respect to any basis {ω_j }_j=1^N for 𝒫_r^- Λ^k (T), is full rank. Having at hand an invertible Vandermonde matrix (<ref>) allows us to compute the Lagrange basis dual to weights associated with 𝒮, see <cit.>, and hence to expand the interpolated as Πω = ∑_s_i ∈𝒮( ∫_s_iω) ω_s_i,
where ω_s_i is such that ∫_s_jω_s_i = δ_i,j. The coefficients of such an expansion are now measurements of ω on 𝒮. We address to <cit.> for a computational-oriented discussion about the duality between bases and degrees of freedom.

Generically, any collection of simplicial supports is unisolvent. Since any s_i ∈𝒮 is parametrised by its k+1 vertices, any collection 𝒮 is an element of the Euclidean space (ℝ^k+1)^N. For any basis {ω_j }_j=1^N of 𝒫_r^- Λ^k (T), via the Lasserre theorem <cit.>, the determinant of the Vandermonde matrix (<ref>) defines a polynomial function of the vertices of the supports {s_i}. Hence, the complement of the non-unisolvence condition V = 0 defines a dense open set of (ℝ^k+1)^N. As a consequence, we deduce that even a non-unisolvent family may be slightly perturbed to obtain a unisolvent one. 











Two relevant unisolvent sets are studied in literature. Following the original definition (see also <cit.>), a small k-simplex s_(α, σ) is the datum of a pair (α, F_σ) and a contracting factor r, so that

    s_(α, σ) = 1/r( F_σ + [v_0 | … | v_n] α^T ) .

This construction is detailed in <cit.>. 
It is immediate to observe that small k-simplices are completely supported in T. 




The following result is proved in <cit.>.

 
    Let r ∈ℕ. For any k ≥ 0, the set
    
    X_r^k (T) ≐{ s_(α, σ)|σ∈𝔖_k^n,  α∈ℐ (n+1,r-1) }

    induces weights that determine _r^- Λ^k (T).


There is a natural bijection between the set X_r^k (T) and the system of generators λ^αω_σ, with α∈ℐ (n+1,r-1) and σ∈𝔖_k^n. As a consequence, # X_r^k (T) ≥_r^- Λ^k (T), and the set of Proposition <ref> is minimal (hence unisolvent in a strict sense) if and only if k = 0, n or r = 1. To restore minimality one may perform the same selection as in Theorem <ref>. This is proved in <cit.>.

 
    Let r ∈ℕ. For any k ≥ 0, the set
    
    X_r,min^k (T) ≐{ s_(α, σ)|σ∈𝔖_k^n,  α∈ℐ_σ (n+1,r-1) }

    induces unisolvent weights for _r^- Λ^k (T).


The sets of Proposition <ref> and Theorem <ref> coincide for k = 0 and k = n. In Figure <ref>, the difference between X_2^k (T) and X_2,min^k (T) consists only of the small edge ( E_0, (1,0,0) ), since E_0 = [v_1, v_2] is associated with the permutation σ (0,1) = (1,2) and α_0 = 1, being α = (1,0,0). In the three dimensional framework their difference is studied in <cit.> for the case of edges k = 1 and in <cit.> for the case of faces k = 2. Despite these differences, both the above sets are relevant: the first one makes it possible to keep a rich geometrical and algebraic structure on small simplices <cit.>, whereas the second gives a square and invertible Vandermonde matrix <cit.>.


The possibility of geometrically choosing small simplices makes weights 
a natural extension of Lagrange interpolation to differential forms. Consistently, the choice of these supports severely affects the stability of the corresponding interpolation operator, leading to Runge-like phenomena <cit.>. This is explained in terms of the generalised Lebesgue constant <cit.>, although the analysis of the theoretical and numerical performances of minimal sets different from X_r,min^k (T) are mostly conjectural.



 §.§ Least squares approximation
 

Given the complexity of explicitly designing minimal sets of weights, in place of looking for the above described interpolated Πω, one may look for a Whitney form Πω∈_r^- Λ^k (T) that minimises the quantity

    ‖ r ‖^2 = ∑_i=1^M ( ∫_s_i( Πω - ω) )^2,    s_i ∈𝒮.

The resulting operator Π : Λ^k (T) →𝒫_r^- Λ^k (T) is a projector <cit.>, i.e. Π(Πω)= Πω. Assume 𝒮 is rich and general enough, so that #𝒮 = M ≥ N ≐_r^- Λ^k (T). Since Πω∈_r^- Λ^k (T), it may be written as Πω = ∑_j=1^N a_j ω_j, where in this case coefficients a_j are not anymore measurements over 𝒮 but are chosen to minimise (<ref>). 


The residual ‖ r ‖^2 may be therefore thought as a function of variables { a_i }_j=1^N which we look to minimize. Let W_i,j≐∫_s_iω_j. We hence impose, for each j = 1, …, N,

    0 = ∂‖ r ‖^2/∂ a_j = 2 ∑_i=1^N W_i,j( ∑_ℓ=1^N W_i,ℓ a_ℓ - f_ℓ),




Up to the scaling constant, the right hand side of (<ref>) is the j-th component of ( W a - f)^T W. As consequence, collecting the 
N equations of (<ref>) yields the rectangular linear system of normal equations

    W^T W a = W^T f .


Theorem <ref> comes now handy. If 𝒮 contains a unisolvent set (for instance considering the set described in Proposition <ref>), then rank (W) = rank (W^T W) = N and (<ref>) has a unique solution

    a = W^†f ,

where W^† = (W^T W)^-1 W^T is the Moore-Penrose pseudo-inverse of W. 



 §.§ Weighted least squares


When all the data to be fitted may be not equally relevant or reliable, least squares approximation can be appropriately modified by introducing correction coefficients. This leads to weighted least squares <cit.> (this is an unfortunate coincidence of nomenclature: this concept of weight has nothing to do with the degrees of freedom for Whitney forms). This transforms (<ref>) into

    ‖ r ‖^2 = ∑_i=1^M b_i ( ∫_s_i( Πω - ω) )^2,    s_i ∈𝒮 .

Defining B as the matrix with constant diagonal entries B_i,i = b_i, so B_i,i^1/2 = √(b_i), one looks for

    a∈^narg  min‖ B^1/2( f - W a) ‖ .

The analysis may be carried exactly as that of Section <ref>, and normal equations (<ref>) read W^T B W a = W^T B f.

If we select b_i = |s_i|^-2, the square of the reciprocal of the measure |s_i| of the support of integration s_i, we normalise the (physical) weights introduced in Section <ref>. This reliefs measurements from being influenced by the size of the underlying domain.


    When dealing with interpolation of Whitney forms, the norm
    
    ‖ω‖_0 ≐sup_c ∈𝒞^k (T)1/|c||∫_c ω|

    is taken into account, 
    for its role in the generalised Lebesgue constant, see <cit.>. This norm is designed to provide local information, as the supremum is sought over the collection of k-chains supported in T, denoted by 𝒞^k (T). 
    The space 𝒞^k (T) is sometimes replaced by the set 𝒮^k (T) of simplices supported in T, as in <cit.>. 
    Considering weighted least squares and observing that
    
    √(∑_s_i ∈𝒮1/|s_i|^2( ∫_s_iω - Πω)^2)≥sup_s_i ∈𝒮1/|s_i||∫_s_iω - Πω| ,

    we immediately deduce that a small weighted least squares approximation error over a rich collection of supports suggests a low error in the norm ‖·‖_0.




 §.§ Solving the system


In practical situations, the vector f is a collection of measurements and hence carries some error <cit.>. The propagation of this error depends on the numerical methods used to solve the linear system, see <cit.>. Via the SVD decomposition, one sees that the condition number κ_2 (W) ≐σ_max (W)/σ_min (W) of W is a relevant quantity for the stability of the problem.

Fig. <ref> shows an interesting behaviour of κ_2 (W) as a function of the supports X_r^k (T). In such a figure, the approximating polynomial degree is fixed to 4 and the set of supports 𝒮 = X_r^k (T) ranges in r = 2, …, 20. The red line represents the conditioning of the interpolatory Vandermonde matrix V. These data are consistent with the literature <cit.>. The basis for 𝒫_4^- Λ^k (T) is that described in Theorem <ref>. 
The peak of the conditioning is obtained when the least squares approximation problem is close to the interpolation problem. As 𝒮 is enriched, the conditioning drops and the solution remains stable. An interesting discussion about stability and accuracy of least squares problem is given in <cit.>. 




The low conditioning 
shown in Fig. <ref> and the relatively small dimension of the system prescribed by Eq. (<ref>) suggest that we shall not observe a sensible difference in dependence of the method adopted for the resolution of the linear system. For ease, we hence opt for a QR factorisation. A comparison of methods adopted for solving least squares problems is given in <cit.>.






§ COMPUTATION OF HIGH ORDER WEIGHTS
 

The construction of the matrix W requires the computation of weights for Whitney forms. There is a nice amount of literature about this problem. For the lowest degree Whitney forms, i.e., for elements of the space _1^- Λ^k (T), no quadrature rules are involved <cit.>. For their high order counterparts, closed formulae may be obtained by combining localisation techniques <cit.> and a wise use of “local” and “global” barycentric coordinates <cit.>. Alternatively, quadrature formulae may be exploited. We expand this latter idea.






 §.§ Edge least squares
 

In the edge framework, the matrix W and the vector f contain line integrals of 1-forms. 
While their meaning is clear, their computation requires a coordinate representation. 
Coordinates should be taken at convenience in order to easily handle the corresponding quadrature formulae.

Let T be an n-simplex and let ω∈Λ^1 (T). We have a coordinate expansion

    ω = ∑_i=1^n f_i (x_1, …, x_n)  x_i .

Exploiting the duality x_i ( ∂/∂ x_j) = δ_i,j between 1-forms and vector fields, we may associate with ω the vector field F_ω = (f_1, …, f_n ). In coordinates we then have

    ∫_e_iω = ∫_e_iω|_e_i = ∫_e_i j^* ω =∫_e_i F_ω·t_e_i e_i ,

where t_e_i is the unit vector tangent to e_i = [p_i, q_i], i.e., t_e_i = q_i- p_i/‖q_i- p_i ‖, and j^* : Λ^1 (T) →Λ^1 (e_i) is the pullback induced by the inclusion map j : e_i ↪ T. Since j^*ω 
is a univariate function (in an appropriate chart), we may approximate the above integral with the aim of a univariate quadrature rule:

    ∫_e_i F_ω·t_e_i e_i ≈∑_j=1^m c_j (F_ω·t_e_i) (x_j) ,

being c = (c_1, …, c_n) the vector of weights of the quadrature rule. The evaluation nodes x_j, for j = 1, …, m, belong to e_i, i.e. to the segment p_i + γ (q_i - p_i), with γ∈ [0,1]. Writing this relation componentwise, one explicitly retrieves the evaluation nodes from the one dimensional case, independently of the dimension of the ambient space. If m is chosen accordingly to the degree r, the quadrature formula (<ref>) is exact for degree r Whitney 1-forms, hence entries of the whole matrix W carry no approximation error.















Equation (<ref>) provides a way to approximate circulations. Edge least squares for Whitney forms, for polynomial degree r = 2, are invoked in literature in the context of incompressible flows <cit.>.



 §.§ Face least squares
 

We treat this section in the evocative framework of 2-simplices supported in a tetrahedron T. Remark <ref> extends the idea to any simplex of codimension 1 in an n-dimensional framework. In a tetrahedron, the matrix W and the vector f contain surface integrals of 2-forms ω. In coordinates, we have

    ω = f_z (x,y,z)  x ∧ y + f_y (x,y,z)  x ∧ z + f_x (x,y,z)  y ∧ z ,

and we want to compute

    ∫_s_iω = ∫_s_iω|_s_i = ∫_s_i j^* ω ,

being j : s_i ↪ T the inclusion map and j^* : Λ^2 (T) →Λ^2 (s_i) its pullback. 
We may identify the 2-form ω with a vector field by means of the Hodge star operator ⋆. Letting F_ω = (f_x, f_y, f_z), we have

    j^* ω = F_ω·n_s_i s_i ,

where n_s_i is the outward unit normal to s_i, which is well defined since s_i has codimension 1. If {v_0, v_1, v_2 } are the vertices of s_i, such a vector is readily computed as n_s_i = (v_1-v_0) × (v_2-v_0)/‖ (v_1-v_0) × (v_2-v_0) ‖. Hence, we may exploit a bivariate quadrature rule of weights c = (c_1, …, c_n):

    ∫_s_iω = ∫_s_i F_ω·n_s_i s_i ≈∑_j=1^m c_j (F_ω·n_s_i) (x_j) ,

since j^* ω is a bivariate function in an appropriate chart. 

The nodes x_j, for j = 1, …, m, have the same barycentric coordinates on s_i that the quadrature nodes have on the reference simplex. 
Again, (<ref>) is an equality when ω is a Whitney form and m is chosen appropriately with respect to their degree r, so W carries no approximation error. 

 
    The above reasoning extends to the multi-dimensional framework as follows. First, note that an (n-1)-form reads in coordinates as
    
    ω = ∑_i=1^n f_i (x_1, …, x_n)  x_1 ∧…∧ x_i∧…∧ x_n .

    
    Defining the field F_ω = (f_1, …, f_n) and computing the pullback j^*: Λ^n-1(T) →Λ^n-1 (s_i), one finds
    
    ∫_s_iω = ∫_s_i F_ω·n_s_i s_i ,

    where s_i is the (n-1)-simplex spanned by vertices {v_0, …, v_n-1}. To compute n_s_i we may define u_j≐v_j - v_0, and ask that n_s_i⊥u_j for j = 1, …, n-1, i.e., that u_j^T ·n_s_i = 0. Defining the rectangular matrix U = [ u_1^T | … | u_n-1^T ], all normal vectors are solutions of the linear system U n_s_i = 0. The desired one is obtained by normalising and imposing that [U^T  |  n_s_i] > 0. 














§ NUMERICAL EXPERIMENTS
 


For computational purposes, we express (<ref>) and (<ref>) 
in terms of the canonical basis {e_1, …, e_n }.

 For k = 1, i.e., the edge case, the identification of 1-forms and vector fields given in Section <ref> is enough to obtain a Cartesian representation of Whitney forms. 
We thus get

    λ_i = e_i  if  i > 0,   λ_0 = - ∑_i=1^n e_i .


For k = n-1, to which the case of faces in a tetrahedron belongs, following the identification between (n-1)-forms and vector fields of Section <ref> one gets

    λ_1 ∧…∧λ_i∧…∧λ_n = e_i .

If λ_0 appears in the expansion, such an expression shall be replaced by λ_0 = - ∑_i=1^n λ_i and 
rule (<ref>) now applies, paying attention to signs that change due to alternation.

Table <ref> summarises cases k = 1 in a triangle T ⊂^2 and k = 2 in a tetrahedron T ⊂^3. Matlab demos and codes are available at
<https://github.com/gelefant/LSwhitney>.


 


    To choose the correct degree for the quadrature formula, we observe that when the components of the field F_ω are polynomials of degree at most r, the (multivariate) function F_ω·n_s_i is a (multivariate) polynomial of degree at most r.




 §.§ Choice of supports of approximation
 




When the form to reconstruct is “easy” to capture (in the sense of <cit.>), interpolation in 𝒫_r^- Λ^k (T) on the set X_r,min^k (T) and least squares approximation in 𝒫_r^- Λ^k (T) on the set X_r^k (T) offer close results, as shown in Fig. <ref>, left. 
To improve approximation features of the least square approach we shall enrich the set of supports 
𝒮 and the following question arises somewhat naturally: is the set X_r^k (T) of Proposition <ref> feasible, if we replace r by a large R? 
To obtain an affirmative answer,  the choice of such parameter R requires some attention. We present two possibilities:

    
  * R is significantly larger than the largest degree considered r_max, for instance R = j · r_max with j = 3: this allows to generate the set 𝒮 once and for all. In presence of a loop on the degree r, for the lowest degrees it may give an unnecessarily large computational cost.
    
  * R is a function of r, for instance R = r + j, with j = 10: since the set X_r^k (T) depends on r, at each iteration of the loop the collection 𝒮 must be recomputed.





Fig. <ref> shows the sensitivity of the approximation process to this choice. If the form to be reconstructed is “good”  such as that defined in (<ref>), even for R = r approximation and interpolation both produce satisfactory and comparable results, see Fig. <ref>, left. In contrast, if we consider the Runge 1-form (<ref>), the value of R must be carefully selected. In Fig. <ref>, right, we fix R = r + 10. When r is low enough (so that when we project onto 𝒫_r^- Λ^k (T) we have R sensibly grater than r), the approximated seems to converge towards the Runge form. However, as r increases, this convergence deteriorates. Errors are computed with respect to ‖·‖_0. This example is detailed in Section <ref>.













 §.§ On the real line


The case k = 0 corresponds to the classical nodal least square interpolation, see e.g. <cit.> and <cit.>. We consider and extend a peculiar aspect of least squares, which is treated in literature. We consider the Runge <cit.> function f: [-1,1] → given by

    f(x) = 1/1+25x^2 .

Such a function is well known to be not suited for interpolation at equispaced nodes. In fact, if r ≥ 1 is the polynomial degree of the interpolant p_r (x) at equispaced nodes x_i = -1 + 2i/r, the quantity ‖ p_r (x) - f(x) ‖_∞ is well known to diverge as r increases (see, e.g., <cit.>).

In <cit.>, authors try to find a least square solution to the Runge phenomenon induced by interpolation on equispaced nodes. They fix a low polynomial degree r, refine the grid x_i = -1 + 2i/M with M sufficiently large and see if the polynomial p_r : [-1,1] → that solves the least square problem

    ∑_i=1^M ‖ p_r(x_i) - 1/1+25(x_i)^2‖_2

is in fact close enough to the Runge function (<ref>). The answer is negative.

The case k = 1 is that of segmental interpolation (or histopolation <cit.>). 
We consider again the function (<ref>) and, in place of considering nodal evaluations on equispaced points x_i = -1 + 2 i /r, we consider as degrees of freedom its uniform weights,

namely integrals of (<ref>) on s_i = [ -1+2i/r, -1+2(i+1)/r], for i = 0, …, r-1. That the corresponding interpolating polynomial do not converge towards f is known <cit.>. We thus set up (<ref>) to this framework, considering a large M, uniform segments s_i = [ -1+2i/M, -1+2(i+1)/M], and we look for a polynomial p_r : [-1, 1] → of degree r such that

    ∑_i=1^M ‖1/|s_i|∫_s_i p_r(x)  x - 1/|s_i|∫_s_i1/1+25(x_i)^2 x ‖_2

is the least possible. Since any family of r + 1 segments in the set { s_i }_i=1^M is unisolvent for ℙ_r (see <cit.> or <cit.>), the solution to the least squares problem (<ref>) is unique. Consistently with the nodal case, the left hand panel of Figure <ref> shows that the increment in the number of segments, when the degree is fixed (in the picture, r = 5), does not push toward convergence. To defeat the Runge phenomenon and hence obtain convergence, it is necessary to combine an increment of the degree with a large number of supports, as shown in the right hand panel of Figure <ref>.


    
An explicit implementation strategy for this framework can be derived from <cit.>.



 §.§ Segments in the triangle
 

We perform two numerical tests to test the effectiveness of the least square approximation of Whitney 1-forms. First, we consider the form

    ω_1 ≐ e^x+y x + sin(xy)  y .

We compute its interpolated with respect to weights of the set X_r,min^1 (T) of Theorem <ref> and its approximated with respect to least squares over the set X_r^1 (T) of Proposition <ref>. We then take a rich collection of 7056 segments supported in T and here estimate the error with respect to (<ref>). Results of the two techniques are comparable, see Fig. <ref>, left.

We then perform the same test on a Runge-like 1-form. The construction of these forms is based on generalisations of the univariate Runge function <cit.> to the multivariate framework <cit.>. 
The differential of any function suffering from the Runge phenomenon is in fact a 1-form that presents the same behaviour with respect to the norm (<ref>), see <cit.>. A neat strategy in the choice of supports to counter this effect has not been found yet, and only mild solutions have been obtained with interpolation strategies (see <cit.> and <cit.>).


We follow the above strategy and consider the function

    f(x,y) = 1/1 + α( x - 1/3)^21/1 + α( y - 1/3)^2 ,


where α is sufficiently large to observe the phenomenon (following <cit.> and <cit.>, in the experiments we take α = 100) and the shift of factor 1/3 centers the function to the standard simplex. Then


    f = 2 α(x-1/3)/( 1+ α(x-1/3)^2 )^2 ( 1 + α( y - 1/3) ) x + 2 α(y-1/3)/( 1+ α(x-1/3) ) ( 1 + α( y - 1/3)^2 )^2 y .

The complexity of finding a collection of weights that gives converge of the interpolated is discussed in <cit.> and <cit.>. To be more definite, in the numerical tests available in literature, what is achieved is a non-increasing trend of the error as the degree increases. Fig. <ref>, right, consistently with the behaviour observed on the real line, shows that a combination of least squares and redundant simplices of Proposition <ref> provide a possible solution to mitigate the Runge effect. In our experiments, for any degree r =1,…,r_max, as supports for the least squares approximation we take X_3r_max^1 (T).





 §.§ Faces in the tetrahedron



To implement the quadrature formula (<ref>) we start from a two-dimensional quadrature formula on a two-dimensional simplex. 
We embed this triangle in ^3 by adding a third coordinate and setting it to zero. Call this T. We shall map T to  any 2-face T_i = {y_0, y_1, y_2 }⊂ T ⊂^3. Put

    v_j = x_j - x_0  for  j = 1,2 ,   n = v_1 ×v_2

and 

    w_j = y_j - y_0  for  j = 1,2 ,   n_i = w_1 ×w_2 .

Since the above are both bases for ^3, the affinity φ_i(x) = A_i x + b_i that maps T to T_i (hence quadrature nodes) is completely determined by the set of equations

    A_i v_j = w_j    j = 1, 2 
    
        A_i n = n_i 
    b_i = y_0 - A_i x_0 .


The matrix A_i and the vector b_i must be determined for each triangle T_i. To avoid the resolution of this linear system, we start with a quadrature formula on the standard triangle in ^2. 
Under this hypothesis, the reference triangle reads T = {0, e_1, e_2 }, whence v_1 = e_1, v_2 = e_2 and n = e_3. Thus

    A_i = [y_1 - y_0  |  y_2 - y_0  |  n_i]  and b_i = y_0 .


We again perform two numerical experiments. In the first we consider the 2-form


    ω_2 ≐ e^z+xy x ∧ y + sin (xy) √(z+1) x ∧ z + cos(x+y)  y ∧ z

and compare the interpolated with respect to the set X_r,min^2 (T) with the approximated obtained from the set X_r^2 (T), in terms of the 0-norm (<ref>). Results are close, as shown in Fig. <ref>.

For the second case, we construct a Runge-like 2-form. Although they are not explicit in literature, we may exploit the identification given in section <ref> and take the Runge 2-form ω_2 whose coefficients (up to reordering) coincide with the Runge 1-form obtained by differentiating

    f(x,y,z) = 1/1 + α( x - 1/4)^21/1 + α( y - 1/4)^21/1 + α( z - 1/4)^2 .


The shift centres the function with the barycenter of the standard tetrahedron T and α = 10, as in Section <ref>. Consistently with the segmental framework, the interpolated over the set X_r,min^2 (T) diverges, and so does the least squares approximation obtained on X_r^2 (T). In contrast, a combination of the increment of the polynomial degree r, with r=1,…,r_max, and a richer least square approximation, obtained over the set X_r_max+10^2 (T) shows convergence towards the Runge 2-form.

In all the above numerical test the collection of supports used to estimate the norm of the error contains 13680 triangles. 





§ CONCLUSIONS AND OPEN PROBLEMS
 

In this paper, we provided a framework for the construction of Whitney forms based on least squares weights. This bypasses the hard task of studying unisolvence of sets of k-simplices. When the sets of redundant supports are entangled with the corresponding minimal sets, the behaviour of the least squares based approximant is close to the interpolant. In contrast, when these former sets are sensibly enriched, the approximated and the interpolated form diverge one from the other. If this is not particularly relevant on easy to capture objects, it becomes relevant when dealing with Runge-like forms. As a consequence, a wise combination of the total polynomial degree and the enrichment of the set supports may be used to tame the Runge phenomenon.

Some deep questions, which we aim to treat in the future, remain open. First of all, we presented interpolation and approximation as two antithetic techniques. This is not completely fair, as blended methods may be applied: for instance, we could ask a polynomial to fit some data and minimise the distance from others, as in <cit.>. This is for sure a relevant development, which requires also a non trivial understanding of the geometry of the problem.

We selected a fixed set of supports of integration in the construction of the matrix W (see discussion in Section <ref>). This may be expensive. To contain the computational cost, it is of interest to retrieve higher-degree matrices by exploiting lower-degree ones. This depends on the structure of the space 𝒫_r^- Λ^k (T) and its decomposition. This topic is of self-standing interest, see e.g. <cit.>.

Further, in an academic problem, the error of the right hand side can be handled by considering a sufficiently strong quadrature rule and all sizes at play are controlled. In a real physical or engineering problem the number of observations might sharply grow and the vector f may contain a non trivial noise. The reconstruction of the approximant requires a thorough study of the numerical linear algebra of the problem: stable methods such as QR or SVD factorisation bear a large computational cost; 
in contrast, faster methods might overlook the propagation of the error. For large systems, the analysis of stabilization methods and involved techniques is thus of interest.
















§ ACKNOWLEDGEMENTS



The authors thank Mariarosa Mazza for her precious remarks on the topic and Alvise Sommariva for his suggestions and repository of quadrature formulae on triangles available at <cit.>. This research has been accomplished within Rete ITaliana di Approssimazione (RITA), the thematic group on Teoria dell'Ap­pros­si­ma­zio­ne e Applicazioni(TAA) of the Italian Mathematical Union. The authors are members of the Gruppo Nazionale Calcolo Scientifico-Istituto Nazionale di Alta Matematica (GNCS-INδAM), which partially supported this research. The first author is funded by INδAM and supported by Università di Padova. The second author has been supported by Fondazione CRT, project 2022 Modelli matematici e algoritmi predittivi di intelligenza artificiale per la mobilità sostenibile.



