Betti numbers

5. Betti numbers#

The two-dimensional de Rham complex of lowest order finite element spaces is

\[\tag{1}\newcommand{\curl}{\mathrm{curl}}\newcommand{\grad}{\mathrm{grad}} 0 \longrightarrow V_{h1} \stackrel{\grad}{\longrightarrow} N_{h0} \stackrel{\curl}{\longrightarrow} W_{h0} \longrightarrow 0 \]

where \(V_{hp}\), with \(p=1\), is the lowest order Lagrange finite element space defined in a prior notebook, \(N_{hp}\) with \(p=0\) is the lowest order Nédélec space from another prior notebook, and \(W_{hp}\) is the lowest order DG space of piecewise constants (also defined in earlier). We know that (1) is a complex, have studied its cohomology, and know that its Betti numbers are topological invariants.

In this notebook we will see practical manifestations of these theoretical results. In fact, the results not only hold for the lowest order spaces in (1), but also for general orders, as will be seen numerically in the experiments below. In order to understand the cohomologies, by definition, we must understand the ranges and null spaces of each operator in the complex. In the two-dimensional (2D) case, this means understanding the ranges and kernels of the gradient and curl operators acting on the finite element spaces. We start with the curl operator.

5.1. Range of the curl#

Considering a general degree \(p\), observe that the application of the 2D curl operator (discussed at length previously) to a function in the Nédélec space \(N_{hp}\) results in a function with no continuity restrictions across element interfaces in general. Since the application of curl also reduces the polynomial degree within each element, the result is a piecewise polynomial of degree not more than \(p\), i.e., a function in the DG space \(W_{hp}\). Thus, we may view \(\curl\) as an operator from the Nedelec space to the DG space,

\[ \curl : N_{hp} \to W_{hp}. \]

How much of the DG space is filled with the range of \(\curl\)?

To answer this, we will need to understand the range of the above map. Since \(\curl\) may act on smooth function spaces or on finite element spaces like \(N_{hp}\), it is a good idea to indicate the specific domain we consider. So we will write \( \text{`` }\curl: N_{hp} \text{ ''} \) to indicate the \(\curl\) when it’s viewed as an operator acting (only) on \(N_{hp}\) functions. We now proceed to compute the dimension of its range, i.e., the number

\[\newcommand{\rank}{\mathop{\mathrm{rank}}} \rank(\curl: N_{h, p}). \]

To do so, we will need a representation of the operator \(\curl: N_{hp}\) as a matrix obtained using the basis functions of \(N_{hp}\) and \(W_{hp}\). Of course, we will also need a domain and a mesh.

Let us punch some holes into a square to make the domain a bit more interesting than just a square and proceed.

import ngsolve as ng
from ngsolve.webgui import Draw
from netgen.occ import OCCGeometry, WorkPlane #, Rectangle, Glue
from ngsolve import GridFunction, curl

wp = WorkPlane()
sqr = wp.Rectangle(2, 2).Face()
rgthole = wp.MoveTo(1.5,1).Rectangle(1/6, 1/6).Face()                 
lfthole = wp.MoveTo(0.5,1).Rectangle(1/6, 1/6).Face()      
dom = sqr-rgthole -lfthole
Draw(dom);
mesh = ng.Mesh(OCCGeometry(dom, dim=2).GenerateMesh(maxh=1/4))

The Nedelec and DG spaces on this mesh for a given degree \(p\) are set as follows.

p = 4  
N = ng.HCurl(mesh, order=p+1-max(1-p,0), type1=True)
D = ng.L2(mesh, order=p)

Now let us think about how to make a matrix representation of the operator \(\newcommand{\curl}{\mathrm{curl}}\newcommand{\grad}{\mathrm{grad}}\curl: N_{h p} \to W_{hp} .\)

One way to obtain this matrix is to take each shape function in \(N_{hp}\), apply curl, and then expand the result in terms of the shape functions of the DG space. However, since the application of curl in NGSolve does not output an object in the DG space, we would need to do an additional step of using the Set method of the DG space: if you provide it a function that you know should be contained within the DG space, then the Set method interpolates it with zero error, creating the interpolant as a GridFunction object. Such objects can be queried for their vector of coefficients of basis expansions, which become entries of the curl matrix. An implementation would look like this:

def curlmatrix(N, D):
    """Return the matrix representation of curl: N -> D"""
    curlmat = []
    shapeN = GridFunction(N)
    curlD = GridFunction(D)
    for i in range(N.ndof):
        shapeN.vec[:] = 0;  shapeN.vec[i] = 1   # i-th shape function 
        curlD.Set(curl(shapeN))                 # interpolation error is 0
        curlmat.append(list(curlD.vec)) 
    return np.array(curlmat).T

However, the code above is inefficient. As you know, python loops are slow. Moreover, inside the above loop, the Set does a global computation when what was really needed is just a local computation on the support of each shape function. Thankfully, NGSolve provides a ConvertOperator which accomplishes the same thing faster. (Go ahead and check that the above curlmatrix function gives the same matrix as the few lines below; and check the speed difference too!)

from ngsolve import curl 
from ngsolve.comp import ConvertOperator
from scipy.sparse import coo_matrix

# Compute the matrix representation of curl: N -> D

u = N.TrialFunction()
Crl = ConvertOperator(N, D, trial_proxy=curl(u))
i, j, val = Crl.COO()
curlmat = coo_matrix((val, (i, j))).toarray()

Here we have used

  • ConvertOperator(X, Y, trial_proxy=M) to make the sparse matrix of the map \(M:X \to Y\), with \(M, X, Y\) set to \(\curl,\) \(N_{hp}\) and \(W_{hp}\), respectively,

  • converted the resulting ngsolve sparse matrix to the standard coordinate sparse format COO,

  • and used scipy’s coo_matrix to pass the matrix to scipy.

Since our example is quite small, we have also converted the sparse matrix to the dense matrix curlmat for convenience (but we should not do this for large problems). Now that we have the matrix representation curlmat of the linear operator \(\newcommand{\curl}{\mathrm{curl}} \curl: N_{hp}\), computing its rank can be done by any number of linear algebra techniques. Here is one way, which proceeds by counting the number of nonzero singular values.

from scipy.linalg import svdvals

s = svdvals(curlmat)
print('Rank of curl =', sum(s>1e-15))

Rank of curl = 2400

Of course, the rank should be at most the dimension of the codomain, which is immediately retrieved from the ndof attribute of the finite element space:

D.ndof

2400

This equals the same number output in the previous cell. Thus we are computationally led to the interesting discovery,

\[\newcommand{\rank}{\mathop{\mathrm{rank}}} \newcommand{\curl}{\mathrm{curl}} \rank(\curl: N_{hp}) = \dim W_{hp}. \]

i.e., the 2D operator \(\curl : N_{hp} \to W_{hp}\) is surjective!

Do feel free to hit the above code with other meshes of your choice. You will end up with the same observation.

5.2. Null space of the curl#

Having discussed the range of curl, let us now turn to the kernel (or the null space) of the same linear operator \(\newcommand{\curl}{\mathrm{curl}} \curl: N_{hp} \to W_{hp}\). Recalling the rank-nullity theorem from linear algebra, we can immediately compute the dimension of the kernel using the above-computed dimension of the range:

\[\begin{split} \begin{aligned} \dim(\ker(\curl:N_{hp})) & = \dim(N_{hp}) - \rank(\curl:N_{hp}) \\ & = \dim(N_{hp}) - \dim W_{hp} \end{aligned} \end{split}\]
dimkercurl = N.ndof - D.ndof
dimkercurl

2100

5.3. Range and null space of the gradient#

Let us now turn to the dimension of the subspace of gradient fields in \(N_{hp}\), i.e., the rank of

\[ \newcommand{\curl}{\mathrm{curl}}\newcommand{\grad}{\mathrm{grad}} \grad : V_{h, p+1} \to N_{hp}. \]

Of course, since \(\curl \circ \grad = 0\), the range of the gradient map is contained in the kernel of curl. However, we would like to compute the exact dimension of its range.

The dimension of \(\grad V_{h, p+1}\), is easy to compute because we know that on a connected domain (such as the domain we are currently working with), the kernel of \(\grad: V_{h, p+1}\to N_{h, p}\) consists of the one-dimensional space of constant functions. Hence,

\[ \ker( \grad : V_{h, p+1} ) = 1 \]

so by the rank-nullity theorem,

\[ \dim \grad (V_{h, p+1}) = \dim V_{h, p+1} - 1. \]
V = ng.H1(mesh, order=p+1)
rankgrad = V.ndof -1
rankgrad

2098

5.4. Betti numbers#

For any cochain complex, recall that the Betti numbers are defined as the dimensions of the cohomology spaces, i.e., the quotient spaces of kernels modulo ranges. In our case, for the complex

\[\tag{2}\newcommand{\curl}{\mathrm{curl}}\newcommand{\grad}{\mathrm{grad}} 0 \longrightarrow V_{h, p+1} \stackrel{\grad}{\longrightarrow} N_{hp} \stackrel{\curl}{\longrightarrow} W_{hp} \longrightarrow 0 \]

viewing \(d^0 = \grad\) and \(d^1 = \curl\), observe that

  • the space of 0-cocycles \(\mathcal{Z}^0 = \ker (\grad: V_{h, p+1})\),

  • the space of 0-coboundaries \(\mathcal B^0 = \{ 0\}\),

  • the space of 1-cocycles \(\mathcal Z^1 = \ker (\curl)\),

  • the space of 1-coboundaries \(\mathcal B^1 = \grad (V_{h, p+1})\),

We have already seen how to compute these ranges and null spaces. Recall that the \(k\)th Betti number \(b^k\) of the complex is the dimension of the \(k\)th cohomology space \(\mathcal H^k = \mathcal Z^k / \mathcal B^k\). Thus, in our case,

  • 0th Betti number \(b^0\) equals \(\dim \ker (\grad: V_{h, p+1})\)

  • 1st Betti number \(b^1\) equals \(\dim (\ker (\curl)) - \dim (\grad (V_{h, p+1}))\).

b0 = V.ndof - rankgrad
b0

1

b1 = dimkercurl - rankgrad
b1

2

In calculus, we learned that every curl-free vector field can be written as a gradient of a scalar function, provided the domain has the topological property of being simply connected. However, we punched two holes in our domain, so our domain is not simply connected. Hence we cannot assert that all curl-free functions on our domain are gradients. Those that are not span a space of dimension \(b^1\), the number we just computed above.

5.5. The emergence of topology#

From the study of simplicial complexes, we have seen that Betti numbers are topological invariants. In particular,

  • the 0th Betti number \(b^0\) equals the number of connected components of the domain, and

  • the 1st Betti number \(b^1\) equals the number of holes in the domain.

Clearly, these match our computational observations above.

We have thus seen a practical manifestation of how fundamental topological invariants are naturally encoded in the finite element spaces of the de Rham complex.

Questions for discussion:

  • How would the above computations change if we punched more holes in the domain?

  • What if we changed the mesh size?

  • What if we changed the polynomial degree \(p\)?

  • Above, we have computed

    \[\newcommand{\grad}{\mathop{\mathrm{grad}}} b^1 = \dim (\ker(\curl: N_{hp})) - \dim(\grad V_{h, p+1}). \]

    and seen its topological significance. What about the number obtained if you replace the Nédélec space \(N_{hp}\) above to the space of vector fields in the Cartesian product of Lagrange spaces \(V_{h, p+1} \times V_{h, p+1}\)? (Hint/suggestion: Use ng.VectorH1 and ng.Grad.)

  • The above experiments considered a connected domain. What happens if do not restrict to connected domains? Understanding this, write a python function that takes as input a mesh and outputs the Betti numbers \(b^0\) and \(b^1\). Display your results from meshes of connected and disconnected domains, with and without holes.

5.6. Summary#

We have seen

  • matrix representations of the operators in the de Rham complex,

  • use of ConvertOperator for efficient change of basis,

  • how to make a scipy sparse matrix from ngsolve’s sparse matrix,

  • computation of ranks and null spaces of operators,

  • how to compute Betti numbers of cohomology spaces

  • the topological invariants of the underlying domain reflected in the Betti numbers, and

  • how topology is encoded in a specific sequence of finite elements.

Contents