2.1 Pure algebraic point of view

The main object of study are homomorphisms between vector spaces: linear transformations, also called linear maps. We will later see that tensors are such maps.

Definition 1 (vector space). A vector space (over a field S) is a commutative group $\mathsf{v}$ equipped with a compatible notion of scaling by elements of S.

A vector space must additionally satisfy a number of laws, including that scaling is a linear operation: $s \triangleleft (x+y) = s \triangleleft x+ s \triangleleft y$ .

The exact nature of this field of scalars (S) has little bearing on the algebraic developmentFootnote ³ , but we assume throughout that they are real numbers. Note that S is itself a vector space, with scaling ${(}\mskip 0.0mu\triangleleft \mskip 0.0mu{)}$ then being scalar multiplication.

Definition 2 (linear map). A function $f: V \!\longrightarrow\! W$ is a linear map iff. for all collections of scalars $c_{i}$ and vectors $\vec v_{i}$ we have

\begin{equation*}f \left(\sum _{i} c_{i} \triangleleft \vec v_{i}\right) = \sum _{i} c_{i} \triangleleft (f(\vec v_{i}))\end{equation*}

For a fixed domain and codomain, linear maps themselves form a vector space.

The eager reader should be warned that, for now, indices are used to range of over arbitrary sets of vectors and scalars (and bound by $\sum$ ), in a usual way. Indices take a special meaning only when we get to coordinates and the Einstein notation (from Sections 2.2 and 2.3).

Since covector spaces are special cases of linear maps, they form vector spaces too. In a similar vein, the set of linear maps $f: S \!\longrightarrow\! W$ is isomorphic to W. (Indeed $f(s) = f (s \triangleleft 1) = s \triangleleft f(1)$ , and thus the vector f(1) in W fully determines the linear function f.)

Definition 4 (bilinear map). A function $f: V \times W \longrightarrow U$ is a bilinear map iff. for all $c_{i}, d_j : S$ , $\vec v_{i} : V$ , and $\vec w_j : W$ we have

\begin{equation*}f \left(\sum _{i} c_{i} \triangleleft \vec v_{i}, \sum _j d_j \triangleleft \vec w_j\right) = \sum _{i,j} c_{i}d_j \triangleleft (f(\vec v_{i},\vec w_j))\end{equation*}

Examples: Here is an attempt at providing an intuition for what is, and is not, a bilinear function. Consider the simplest case of the definition of bilinear map where there is just one vector $\vec v$ as the first argument and one vector $\vec w$ as the second argument to f. We then have $f(\vec v,0) = f(1\triangleleft \vec v,0\triangleleft \vec w) = (1 \times 0)\triangleleft f(\vec v,\vec w) = 0$ . This means that vector addition is not bilinear because $\vec v+0=\vec v \neq 0$ . Similarly, f cannot be first or second projection, because they are also linear, not bilinear.

We also have that we can “move constant factors” between $\vec v$ and $\vec w$ : $f(c \triangleleft \vec v, 1 \triangleleft \vec w) = (c \times 1) \triangleleft f(\vec v,\vec w) = (1 \times c) \triangleleft f(\vec v,\vec w) =f(1 \triangleleft \vec v, c \triangleleft \vec w)$ . In connection with the tensor product this means that even though, for any two vectors $\vec v : V$ and $\vec w : W$ , we can construct a tensor $u = \phi(\vec v,\vec w) : V\otimes W$ which looks like we have embedded a pair, we cannot extract $\vec v$ and $\vec w$ again—they are mixed up together (entangled).

What a bilinear function can (and must) do, as we can see from the definition, when given two linear combinations, is to compute a linear combination based on all pairwise products of the coefficients, without depending on the coefficients themselves.

Order of a tensor

Often, tensors are used in a context where there is a single (atomic) underlying vector space $\mathsf{T}$ which is not just the scalars. Then the complexity of a vector space built from $\mathsf{T}$ can be measured by its order. The order of $\mathsf{T}$ is defined to be 1 and the order of the scalar space is 0. The order of a tensor space $\mathsf{V}\mskip 3.0mu{\otimes }\mskip 3.0mu\mathsf{W}$ is the sum of the order of spaces $\mathsf{V}$ and $\mathsf{W}$ , and this way we can build spaces of arbitrarily large order. The order of a linear map can be defined either as the pair of the orders of its input and output spaces, or as their sum (depending on convention). For example, a linear operator on an atomic vector space has order (1,1) or 2 in the respective conventions. Morphisms of order three or more are properly called tensors. Conversely, tensors of any order (including 0, 1 and 2) are linear maps, of the appropriate domain and codomain. When there is more that one underlying vector space, the order is not enough to characterise a tensor space: the full type needs to be specified, as in Section 3. (Yet this level of complexity won’t be exercised in this paper.)

2.2 Coordinate representations

In practice, the algebraic definitions are not easy to manipulate for concrete problems, thus one most commonly works with coordinate representations instead. (Our goal will be to break free of those eventually.) As a reminder, given a basis $\vec e_{i}$ , any vector $\vec x\in V$ can be uniquely expressed as $\vec x = \sum _{i} x^i\triangleleft \vec e_{i}$ where each $x^i $ coordinate is a scalar. In this way, given a basis $\vec e_{i}$ , a vector space is isomorphic to its set of coordinate representations. Note that a superscript is used for the index of such coordinates. The general convention that governs whether one should write indices in low or high positions is explained in Section 2.3; for now, it is enough to know that they are indexing notations.

Like vectors, linear maps are also commonly manipulated as matrices of coefficients. For a linear map f from a vector space with basis $\vec {d_{i}}$ to a space with basis $\vec e_j$ , each column is given by the coefficients of $f(\vec {d_{i}})$ . Indeed, using $F{_i}{^j}$ to denote the coefficients, we have:

\begin{equation*} f(\vec x) = f\left(\sum _{i} x^i \triangleleft \vec {d_{i}}\right) = \sum _{i} x^i \triangleleft f(\vec {d_{i}}) = \sum _{i} x^i \triangleleft \left(\sum _j F{_i}{^j} \triangleleft \vec e_j\right) = \sum _j \left(\sum _i F{_i}{^j} x^i \right) \triangleleft \vec e_j\end{equation*}

In general, the values of the matrix coefficients $F{_i}{^j}$ depend on the choice of bases $\vec {d_{i}}$ and $\vec e_j$ , but to reduce the number of moving parts one usually works with a coherent set of bases.

Coherent bases

Starting from an atomic vector space T, one can build a collection of more complicated tensor spaces using tensor product, dual, and the unit (the scalar field S). For coordinate representations each such space could, in general, have its own basis, but it is standard to work with a collection of coherent bases. Given a basis $\vec e_{i}$ for a finite-dimensional atomic vector space T, the coherent basis for $T^*$ is the set of covectors $\tilde e^j$ such that $\tilde e^j(\vec e_{i}) = \delta_{i}^j$ . (It is usually called the dual basis.) Likewise, given two coherent bases $\vec {d_{i}}$ and $\vec e_{i}$ respectively for V and W, the coherent basis for $V\otimes W$ is $b_{i,j} = \phi(\vec {d_{i}},\vec e_j)$ , where $\phi$ is given by Definition 5. Note that this basis is indexed by a pair. Accordingly, if the dimension of V is m and the dimension of W is n, the dimension of $V\otimes W$ is $m \times n$ . Additionally, re-associating tensor spaces do not change coherent bases ( $\vec e_{(i,j),k}$ is the same as $\vec {d}_{i,(j,k)}$ , up to applying the corresponding associator). Finally, the scalar vector space has dimension one, and thus has a single base vector, which is coherently chosen to be the unit of the scalar field (the number $1 : S$ ).

Coordinate transformations

Assuming one basis $\vec e_j$ and another basis $\vec {d_{i}}$ for the same vector space V such that $\vec {d_{i}} = \sum _j F{_i}{^j} \vec e_j$ , then the coordinates in basis $\vec {d_{i}}$ for $\vec x$ are $\hat x^j = \sum _{i} F{_i}{^j} x^i $ . We say that the matrix F is the transformation matrix for V given the choice of bases made above, and denote it J(V). Then the transformation matrices for vector spaces built from an atomic space T using the coherent set of bases defined above are given by the following structural rules:

\begin{align*}J(V\otimes W)&= J(V)\otimes J(W)&\text{Kronecker product of matrices}\\J(V^*)&= J(V)^{-1}&\text{matrix inverse}\\J(S)&= 1&\text{scalar unit}\end{align*}

Furthermore, the matrix representation G of a linear map $g : V \!\longrightarrow\! W$ is transformed to: $J(W) \cdot G \cdot J(V^*)$ , where $(\cdot)$ is matrix multiplication. These are the “transformation rules” that Rowland & Weisstein (Reference Rowland and Weisstein2023) allude to in the above quote.

2.3 Einstein notation

The previous section showed how to deal with concrete matrices, using a concrete choice of bases. The next step is to manipulate symbolic expressions involving matrices. The language of such expressions (together with a couple of simple conventions) is colloquially referred to as Einstein notation.

In this notation, every index ranges over the dimensions of an atomic vector space.Footnote ⁴ Consequently, the total number of free (non repeated) indices indicates the order of a tensor expression in index notation. An index can be written as a subscript (and called a low index) or as a superscript (and called a high index).

The location (high or low) of the index is dictated by which coordinate transformation applies to it. That is, if a high index ranges over the dimensions of V, then J(V) applies, whereas $J(V^*)$ applies for a low index. Additionally, every reference to a tensor is fully saturated, in the sense that a symbolic tensor is always applied to as many indices as its order. Thus, for instance, $x^i$ denote (the components of) a vector, and $y_j$ denote (the components of) a covector. The expression $t{_i}{^j}$ refers to (components of) a linear transformation of order (1,1). In the absence of contraction (see below), multiplication increases the order of tensors. For instance, $x^i\,y_j$ also has order (1,1). In general, if t and u are expressions denoting tensors of order m and n, respectively, then their product $t\,u$ denotes a tensor of order $m+n$ .

Contraction

In Einstein notation, the convention is that, within a term, a repeated index is implicitly summed over. (In terms familiar to this journal: such indices are implicitly bound by a summation operator.) Because summation is a linear operator, within a term all the well-scoped locations of the summation operator are equivalent—so it makes a lot of sense to omit them. Additionally, when an index is repeated, it must be repeated exactly twice; once as a high index and once as a low index. Mentioning an index twice is called contraction. Viewing tensors as higher dimensional matrices of coefficients, contraction consists in summing coefficients along a diagonal. Therefore, a contraction reduces the order of the tensor by two.Footnote ⁵ The indices which are contracted are sometimes called “dummy” and those that are not contracted are called “live” (In terms familiar to the functional programming community, dummies are bound variables and live indices are free variables.) To be well-scoped, every term in a sum must use the same live indices. For instance, the expression $t{_l}{^j}u{_m}{^k}v{_i}{^l}{^m} + v{_i}{^j}{^k}$ denotes a tensor of order (1,2). Its live indices are ${}_i,{}^j,{}^k$ , and the indices m and l are dummies.

At this stage, one can see the Einstein notation as a convenient way to notate expressions which manipulate coordinates of tensors. The high/low index convention makes it clear which transformations apply. Yet it may be mysterious why indices must be repeated exactly twice, and why (live) indices cannot be omitted from a term. The answer lies in the following observation. Even though the Einstein notation may originate as a convenient way to express coordinates, it really is intended to describe algebraic objects. The physicists Thorne & Blandford (Reference Thorne and Blandford2015) put it this way:

(A “slot” is a component of input or output tensor space.) The key of this “mind flip” is that live indices correspond to inputs (or outputs) of linear functions, and contraction corresponds to connecting inputs to outputs. The main contribution of this paper is to work out this connection in full as a pair of two edsls.

3.1 Matrix instances

An important instance of the compact closed category structure is the category of matrices of coefficients, which we encountered in Section 2.2. In our host functional language, we define them as a function from (both input and output) indices to coefficients (of type $\mathsf{S}$ ):Footnote ⁸

To emphasise the dependency on the basis, we use a subscript when referring to a specific matrix category morphism, as in $M_{\mathsf{b}}$ where $\mathsf{b}$ is a reference to the choice of basis. In the Haskell implementation, this basis is represented by a phantom type parameter.

The identity morphism is the identity matrix, and composition is matrix multiplication:

In this instance, the objects are identified with sets that index the bases of the vector spaces that they stand for. These sets are assumed to be enumerable and bounded (so we have access to their $\mathsf{inhabitants}$ ) and we can compare indices for equality.Footnote ⁹ The instance of the smc structure for matrix representations in coherent bases is then:

Because objects index the bases of the corresponding vector spaces, tensor products are represented as usual pairs. In the above definition, the right-hand sides are Haskell code. This means that the asterisk operator ${(}\mskip 0.0mu{*}\mskip 0.0mu{)}$ denotes multiplication between scalars as components of matrices. In contrast, the operator $(\star)$ defined in Section 4 denotes multiplication between abstract scalar (order-0) tensor expressions (independent of the chosen tensor representation).

With the coherent choice of bases, $\eta$ and $\epsilon$ are simply realised as the identity.

The object ${\mathsf{a}^*}$ has the same dimensionality as $\mathsf{a}$ , so in our Haskell encoding we use $\mathbf{newtype}$ for it, . For concision (and following tradition), we use an asterisk as a shorthand, so ${\mathsf{a}^*}$ stands for the type $\mathsf{DualObject}\mskip 3.0mu\mathsf{a}$ .

Coordinate representation functors. It is worth mentioning that the transformation functions between linear maps $\mathsf{L}$ and their representations $M_{\mathsf{e}}$ in a given basis $\mathsf{e}$ are a pair of functors which are the identity on objects and just change the morphisms:

Furthermore, this pair defines an isomorphism between the respective compact-closed categories. Therefore, even though different representations form different categories, one can always transform one to another. The transformation between systems of coordinates, usually presented using transformation matrices (see Section 2.2) can be understood as the composition of $\mathsf{fromCoordinates}$ in the source basis, and $\mathsf{toCoordinates}$ in the target basis:

5.1 Tensor fields

In this context a field means that a different value is associated with every position on a given manifold.Footnote ¹⁴ We denote the position parameter by $\mathbf X$ . Hence, the scalars (from Definition 1) are no longer just real numbers, but rather real-valued expressions depending on $\mathbf X$ .Footnote ¹⁵ For instance, $\mathbf X$ could be a position on the surface of the earth, and a scalar field could be the temperature at each such point.

A vector field associates a vector with every position; for instance, the wind direction. The perhaps surprising aspect is that each of these vectors may inhabit a different, local, vector space, which can be thought of as tangent to the manifold at the considered point. So in our example, we assumed that the wind is parallel to the earth surface. Hereafter, we assume such a local space for each category z, and call it $T_{z}$ , leaving the dependency on position $\mathbf X$ implicit. Even though in the typical case the local vector space is different at each position, it keeps the same dimensionality. Therefore, as an object, it is independent of $\mathbf X$ . In Haskell terms, $T_{z}$ is an associated type; see Section 5.2.

When we deal with matrix representations and want to perform computations with them, we need a way to identify the position $\mathbf X$ . For a general manifold, this is difficult to do, but we will restrict our scope to the case where a single coordinate system is sufficient. We also need a basis at each position, which gives a meaning to the entries in a tensor matrix representation (the meaning of these coordinates change with position). Furthermore, different choices of coordinate system are possible for the same manifold. The coordinates used to identify the position will be referred to as the global coordinates, while the coordinates of a tensor will be referred to as local coordinates. This terminology is not usual in mathematical praxis. However, we find that making this distinction is useful to lift ambiguities.Footnote ¹⁶

While the choice of basis field is arbitrary from an algebraic perspective, some choices of basis will make certain computations easier than others. Given a system of coordinates to identify positions in the manifold, there is a canonical way to define the local basis field. It is to let the base vectors be the partial derivatives of the position $\mathbf X$ with respect to each global coordinate. This yields base vectors which are tangent to coordinate lines in the manifold. In Figure 5, one example follows Cartesian coordinate lines, and the other polar coordinate lines. In the polar case, we have the base vector fields $(\mathbf e_\rho, \mathbf e_\theta)$ Footnote ¹⁷ as basis for $T_{M_{\mathsf{p}}}$ , with $\mathbf e_\rho = \partial \mathbf X / \partial \rho$ and $\mathbf e_\theta = \partial \mathbf X / \partial \theta$ .

Fig. 5: Possible fields of bases for a local space field covering the Euclidean plane. In both examples, basis vectors are tangent to coordinate lines; either Cartesian or polar coordinates. In the second instance, the basis vectors are undefined at the origin.

5.2 Metrics and index juggling

An important additional structure that one can add to vector spaces are associated metrics. The (covariant) metric, noted g, is what defines the inner product of (local) vectors. It is a tensor field that, when given two vectors as input, returns a scalar, their inner product. One defines the contravariant metric g’ as the inverse of the covariant metric.Footnote ¹⁸ This can be specified as $g{_i}{_k}g'{^k}{^j}$ = $\delta_{i}^j$ .

To capture this algebraic structure, we distinguish the vector space for which we define the metric as an associated type, $T_{z}$ , of a new category class $\mathsf{MetricCategory}$ .

The various notations for metrics are shown in Figure 6.

Fig. 6: Tensor field primitives in various notations, and their types.

Then we can embed the metric morphism in Albert as follows:

The metric tensor is symmetric—a fact which can be expressed as its antisymmetric part () being zero. The coefficients of the matrix representation of the metric are given by the inner product of the basis vectors. For the usual Cartesian basis, the metric is represented by the identity matrix. For the basis $(\mathbf e_\rho, \mathbf e_\theta)$ defined above, the matrix representation of the metric is

\begin{equation*}g_\mathsf{Polar} =\left[\begin{array}{cc}\rho^{2}&0\\0&1\end{array}\right]\end{equation*}

With this knowledge, we can define an instance for coordinate representations of the $\mathsf{MetricCategory}$ class:

The scalars associated with this additive category are expressions where the variables $\rho$ and $\theta$ occur, and the $\mathsf{variable}$ function embeds a coordinate name into such expression types:

Thus, the $\mathsf{Polar}$ data type serves triple duty. First, it is the set of variables used in such expressions. Second, it serves to identify the meaning of coordinates as a parameter in the type. Third, it is used to construct a representation of the atomic vector space. For the matrix category, an object is a set indexing the base vectors, so $\mathsf{Atom}$ is a simple wrapper around the coordinate type:

One can contrast the type of metrics with that of $\eta$ and $\epsilon$ . First, only the tangent space $T_{z}$ has a metric, while every object has a dual. Indeed, $\eta$ and $\epsilon$ are both realised as $\delta$ in the tensor instance, and are an identity operation in Einstein notation. They have no geometric significance and thus can be defined generically. In contrast, the metric depends on the geometric properties of the space $T_{z}$ that they operate upon. Second, the metrics do not dualise objects, whereas $\eta$ and $\epsilon$ do. Accordingly, our diagrammatic notation does not make any special mark on the input/output wires of metrics. Yet, metrics satisfy a variant of the snake laws: is equal to .

By combining the compact closed categorical structure with metrics, one can construct the following two morphisms (in Roger):

In diagram form $\mathsf{juggleUp}$ is . The existence of such morphisms have a direct consequence on Einstein notation: any subscript index $\mathsf{P}\mskip 3.0muz\mskip 3.0mu\mathsf{r}\mskip 3.0mu\mathsf{a}$ can be raised into $\mathsf{P}\mskip 3.0muz\mskip 3.0mu\mathsf{r}\mskip 3.0mu{\mathsf{a}^*}$ , and vice versa. In other words, any index can be used as super or subscript as needed. As a first example, given a vector $v^{i}$ , the covector $v_i$ is a shorthand for $v^kg_{ki}$ . It is important to note that the position of the indices changes the value of the expression: the actual numbers (the components) of a vector $v^{i}$ are different from those of the corresponding covector $v_i$ unless the metric representation is the identity (as in global Cartesian coordinates). As another example, assuming a tensor $\mathsf{t}\mskip 3.0mu{:}\mskip 3.0muT_{z}\mskip 3.0mu{\otimes }\mskip 3.0muT_{z}\,\mskip 1.0mu\overset{z}{\leadsto }\mskip 1.0mu\,T_{z}$ , the expression $t_{ijl}$ is a shorthand for $t_{ij}{}^{k} g_{kl}$ .

In the mathematical praxis, raising and lowering indices are also referred to as index juggling. In Albert, index juggling is realised explicitly: our recent example can be written $\mathsf{t}\mskip 3.0mu_{\mathsf{i}}\mskip 3.0mu_{\mathsf{j}}\mskip 3.0mu{(}\mskip 0.0mu\mathsf{raise}\mskip 3.0mu_{\mathsf{l}}\mskip 0.0mu{)}$ . The raising and lowering functions have the following types and are implemented by embedding $\mathsf{juggleUp}$ and $\mathsf{juggleDown}$ :

5.3 Change of global coordinate system

As we have seen, Einstein notation is carefully set up to work without reference to the system of coordinates used to identify positions in the manifold. In practice, this means that one can reason about tensors, their derivatives (Section 5.4), etc. without worrying about the choice of coordinates. It is only at the very last stage that one chooses a system of global coordinates where the data of the problem is easy to express and calculate with. We have seen that one can convert between bases by applying the appropriate transformations on the representation matrices (see Section 2.2). For canonical bases (defined as partial derivatives of the position), the $\mathsf{transform}$ function is obtained by composing with Jacobian. In our running example, we have at our disposal the Jacobian J(T) and its inverse $J(T^*)$ to convert between polar and Cartesian global coordinates with $J(T^*) = $ $\left[\begin{array}{cc}\cos{\theta}&\sin{\theta}\\\rho^{-1}\sin{\theta}&-\rho^{-1}\cos{\theta}\end{array}\right]$ .

As an illustration, let $\tilde v_{\mathsf{Polar}}$ = $\left[\begin{array}{cc}-\rho^{-1}&0\end{array}\right]$ be the local coordinates of a co-vector field in the polar tangent basis $\tilde v_{\mathsf{Polar}}$ ${:}\mskip 3.0muT_{M_{\mathsf{Polar}}}\mskip 1.0mu\overset{M_{\mathsf{Polar}}}{\leadsto }\mskip 1.0mu\mathbf{1}$ . Then the local coordinates of the co-vector field in the Cartesian tangent basis are given by ${v_{\mathsf{Cartesian}}}_i = J(T)_i^j {v_{\mathsf{Polar}}}_j$ , and we have: $\tilde v_{\mathsf{Cartesian}}$ = $\left[\begin{array}{cc}-\rho^{-1}\cos{\theta}&-\rho^{-1}\sin{\theta}\end{array}\right]$ .

5.4 Spatial derivative: Levi-Civita connection

The local vector space can change from position to position, but it is assumed to vary smoothly. One says that a connection is defined between neighbouring spaces. This means that we can take the derivative of tensor fields with respect to position. Two cautionary remarks are in order regarding our use of the phrase spatial derivative. First, the terminology is meant to include hyperbolic geometries, thus space-time (Minkowski space) is covered. Second, there are other notions of spatial derivatives, but here we will only consider the case of the Levi-Civita connection, more precisely, which has additional properties: see footnote 20.

Considering the simplest case, the derivative of a scalar field $\mathsf{s}$ is a covector field: its gradient. That is, given a direction vector $\vec v$ , the gradient will return the slope of $\mathsf{s}$ in the $\vec v$ direction. In general, the spatial derivative of a tensor takes a vector argument and returns the variation in this direction. So the derivative of a tensor field is itself a tensor field, whose (covariant) order is one more than its argument, and for this reason, the spatial derivative of tensors is called the covariant derivative. We capture this in the following class, which signals the presence of a connection (and spatial derivative) in a category.

Accordingly, in Einstein notation, the covariant derivative uses an additional (lower) index.

The actual implementation is too involved to present just yet, and deferred to Section 7. In the diagrammatic notation, we represent the covariant derivative as a thick box around the tensor whose derivative is taken. This box adds an input wire for the additional input vector, and propagates the wires of the tensor which it encloses, reflecting the propagation of the types $\mathsf{a}$ and $\mathsf{b}$ in the type of $\nabla$ . In Figure 6, we illustrate with a derivative applied to a tensor of order (1,1), but it can be applied to any morphism, with arbitrary domain and codomain—and hence to tensors of arbitrary order. For instance, the covariant derivative of a tensor $\mathsf{t}$ of order (2,2) is written $\nabla{_i}t{_j}{_k}{^l}{^m}$ or , but (still) ${\nabla}\mskip 3.0mu\mathsf{t}$ in Roger.

Example: Laplacian. As a simple example, we express the Laplacian (the divergence of the gradient) of a given scalar field P. Its physical meaning has no bearing on our development; but if P is an energy potential, its gradient is the corresponding force field, and its divergence the density of an associated charge or mass. We write P for the potential field, or as a diagram. The scalar character of this field is represented by the lack of indices, or the use of the unit object for its domain and codomain. Its gradient is denoted $\nabla{_i}P$ , or . Again, indices or objects indicate that we have an order (1,0) (covector) field. At each point, the local covector is a linear function that takes a direction (vector) into the slope in that direction (a scalar).

The Laplacian is a linear combination of second order derivatives. To compute the linear combination, we need to apply contraction, which always expects an upper and a lower index, but we have two covariant (lower) indices, so we must raise one index by multiplying with the contravariant metric. The Einstein notation is $g'{^i}{^j}\nabla{_i}\nabla{_j}P$ , the corresponding diagram is , and in Roger it is ${(}\mskip 0.0mu{\nabla}\mskip 3.0mu{(}\mskip 0.0mu{(}\mskip 0.0mu{\nabla}\mskip 3.0mu\mathsf{P}\mskip 0.0mu{)}\mskip 3.0mu{\circ }\mskip 3.0mu\rho\mskip 0.0mu{)}\mskip 0.0mu{)}\mskip 3.0mu{\circ }\mskip 3.0mu\mathsf{g'}$ . The Einstein notation, while already economical, can be made even more concise by using the index juggling convention: $\nabla_{i}\nabla^i P$ . In Albert, the same thing would be written

Here the high index from the contraction must be lowered because $\mathsf{deriv}$ can only take a low index as its first argument.

5.4.1 Laws of covariant derivatives

As one might expect, the covariant derivative satisfies the product law for derivatives:

(5.1) \begin{equation}\nabla_{i} (t u) = (\nabla_{i} t) u + t (\nabla_{i} u)\end{equation}

The above formulation is concise, but using it as basis for implementation can be tedious, because one needs to track free and bound variables. A pitfall is that t and u stand for arbitrary expressions, and index variables may occur free in them. Therefore, a specific implementation difficulty is that one needs to preserve the linearity of index variables at the level of the host language.Footnote ¹⁹ Thus, we find that the morphism edsl Roger is a better implementation vehicle in this case. The laws become noticeably more verbose, but not horribly so, and dispense of tracking free variables. In this notation, the product law is expressed as two cases, one for each of the composition ${(}\mskip 0.0mu{\circ }\mskip 0.0mu{)}$ and tensor ${(}\mskip 0.0mu{\otimes }\mskip 0.0mu{)}$ operators:

\begin{align*}{\nabla}\mskip 3.0mu{(}\mskip 0.0mu\mathsf{t}\mskip 3.0mu{\circ }\mskip 3.0mu\mathsf{u}\mskip 0.0mu{)}&=\mathsf{t}\mskip 3.0mu{\circ }\mskip 3.0mu{(}\mskip 0.0mu{\nabla}\mskip 3.0mu\mathsf{u}\mskip 0.0mu{)}\mskip 3.0mu{+}\mskip 3.0mu{(}\mskip 0.0mu{\nabla}\mskip 3.0mu\mathsf{t}\mskip 0.0mu{)}\mskip 3.0mu{\circ }\mskip 3.0mu{(}\mskip 0.0mu\mathsf{id}\mskip 3.0mu{\otimes }\mskip 3.0mu\mathsf{u}\mskip 0.0mu{)}\\{\nabla}\mskip 3.0mu{(}\mskip 0.0mu\mathsf{t}\mskip 3.0mu{\otimes }\mskip 3.0mu\mathsf{u}\mskip 0.0mu{)}&={(}\mskip 0.0mu{\nabla}\mskip 3.0mu\mathsf{t}\mskip 3.0mu{\otimes }\mskip 3.0mu\mathsf{u}\mskip 0.0mu{)}\mskip 3.0mu{\circ }\mskip 3.0mu\bar{\alpha}\mskip 3.0mu{+}\mskip 3.0mu{(}\mskip 0.0mu\mathsf{t}\mskip 3.0mu{\otimes }\mskip 3.0mu{\nabla}\mskip 3.0mu\mathsf{u}\mskip 0.0mu{)}\mskip 3.0mu{\circ }\mskip 3.0mu\alpha\mskip 3.0mu{\circ }\mskip 3.0mu{(}\mskip 0.0mu\sigma \mskip 3.0mu{\otimes }\mskip 3.0mu\mathsf{id}\mskip 0.0mu{)}\mskip 3.0mu{\circ }\mskip 3.0mu\bar{\alpha}\end{align*}

As before, we find the corresponding diagrams more readable:

(5.2)

(5.3)

The derivative of all constant morphisms is zero: $\nabla\mskip 3.0mu\mathsf{id}\mskip 3.0mu{=}\mskip 3.0mu\mathrm{0}$ , $\nabla\mskip 3.0mu\rho\mskip 3.0mu{=}\mskip 3.0mu\mathrm{0}$ , $\nabla\mskip 3.0mu\alpha\mskip 3.0mu{=}\mskip 3.0mu\mathrm{0}$ , etc. This property also holds for the (co)metric tensors:Footnote ²⁰ $\nabla\mskip 3.0mug\mskip 3.0mu{=}\mskip 3.0mu\mathrm{0}\mskip 3.0mu{=}\mskip 3.0mu\nabla\mskip 3.0mug'$ . Together, the above laws fully specify the structural behaviour of the derivative and the implementation of the covariant derivative falls out from them with no additional difficulty. However, the important case of the derivative of a tensor in a coordinate category remains to be addressed. This question leads us to the concept of affinity.

5.4.2 Partial derivatives, christoffel symbols and affinities

One may be tempted to think that the coefficient representation of the covariant derivative is the index-wise derivative of the coefficient representations. While this is true if the metric is the identity everywhere on the manifold, it is not the case in general.

Conventionally, one speaks of the “partial derivative” for the index-wise derivative of coefficients and retains the term “covariant derivative” for the spatial derivative. (The notations are shown in Figure 6, but note that $\partial_{i}t$ stands for the partial derivative of the expression $\mathsf{t}$ with respect to the $i\text{th}$ coordinate, sometimes also written $\partial t/\partial x^i $ .) As usual, we express the availability of this new operation by means of a class:

The Levi-Civita connection is unique for a given metric, but we still require the user to provide an implementation. For categories with coordinates, discussed in the next subsection, a canonical implementation of $\nabla$ is provided.

We make the concept available in Albert like this:

In the general case, to compute the covariant derivative, one must account for the variation of the basis. Therefore, the partial derivative must be corrected by a so called affinity term. The partial derivative accounts for the variation of (the representation of) the tensor field itself as the position varies, while the affinity term accounts for the variation of the basis. The variation of the basis is measured by the Christoffel symbol, denoted $\Gamma$ and of type ${(}\mskip 0.0muT_{z}\mskip 3.0mu{\otimes }\mskip 3.0muT_{z}\mskip 0.0mu{)}\,\mskip 1.0mu\overset{z}{\leadsto }\mskip 1.0mu\,T_{z}$ (Figure 6). Different choices of local basis field for the same manifold will yield different values for it. (Therefore, even though $\Gamma$ is a morphism in a matrix category, and even though it can be transformed to another basis by multiplication with Jacobians, this transformed version will not be the Christoffel symbol for the new basis. This fact is sometimes expressed in textbooks as “ $\Gamma$ is not a tensor.”)

The Christoffel symbol is often treated abstractly, but it is determined by the (coefficient representation of) the metric:

(5.4) \begin{equation}\Gamma{_i}{_j}{^k} = \frac{1}{2}g'{^l}{^k}\partial{_j}g{_i}{_l} + \frac{1}{2}g'{^m}{^k}\partial{_i}g{_j}{_m} - \frac{1}{2}g'{^n}{^k}\partial{_n}g{_j}{_i}\end{equation}

The Christoffel symbol is always symmetric in its first two indices, which is equivalent to asserting that the following diagram is zero: . In the implementation, we make $\Gamma$ a method of the class $\mathsf{CoordinateCategory}$ , but with a default definition in terms of Equation (5.4). We make it available in Albert by embedding the morphism $\Gamma$ as follows:

For 2-dimensional atomic vector space, the Christoffel symbol has $2^3=8$ components. For our running example of the polar coordinate system and associated canonical tangent space, we can write $\Gamma$ as two $2 \times 2$ symmetric matrices, as follows:

\begin{align*}\Gamma^{\rho} =\left[\begin{array}{cc}0&0\\0&-\rho\end{array}\right] \quad\Gamma^{\theta} =\left[\begin{array}{cc}0&1/\rho\\1/\rho&0\end{array}\right]\end{align*}

In textbooks on tensors, for instance that of Lovelock & Rund (Reference Lovelock and Rund1989), one often sees the relation between covariant and partial derivatives expressed as a family of equations, depending on the order of the tensor whose derivative is taken:

\begin{align*}\nabla_{i} T& = \partial_{i} T\\\nabla_{i} T_j& = \partial_{i} T_j - \Gamma_{ij}{^k} T_k \\\nabla_{i} T^j& = \partial_{i} T^j + \Gamma_{ik}{^j} T^k\\\nabla_{i} T_j{^k}& = \partial_{i} T_j{^k} - \Gamma_{il}{^k} T_j{^l} + \Gamma_{ij}{^l} T_l {^k}\end{align*}

etc. Using this sort of definition is particularly error prone (can you spot quickly whether there is a mistake in the last line?). In contrast, Roger captures all cases in one go:

(5.5)

The complexity is pushed down into $\mathsf{affinity}$ , which is invoked once for the domain and once for the codomain of $\mathsf{t}$ . The affinity is a family of morphisms $\mathsf{affinity}\mskip 3.0mu{::}\mskip 3.0muT_{M_{\mathsf{b}}}\mskip 3.0mu{\otimes }\mskip 3.0mu\mathsf{a}\,\mskip 1.0mu\overset{M_{\mathsf{b}}}{\leadsto }\mskip 1.0mu\,\mathsf{a}$ for any object $\mathsf{a}$ constructed from the local vector space $T_{M_{\mathsf{b}}}$ and refers to a specific basis for it, using some coherent canonical choice of bases (with $\mathsf{b}$ as a basis for $T_{M_{\mathsf{b}}}$ ). The affinity for arbitrary vector spaces is defined by induction on the structure of the corresponding object. The affinity for a product object is the sum of affinities for each of the components, leaving the other component untouched. The affinity for a dual object is the negative affinity of the underlying object, with input and output suitably swapped. This can be coded by a type-dependent set of equations:Footnote ²¹

Here, Albert is much more concise than Roger (we omit the corresponding expressions entirely). The corresponding graphical notation is shown in Figure 7 for each case (but using atomic types in place of a proper induction).

Fig. 7: Affinities $\mathsf{aff}\mskip 3.0mu\mathsf{a}$ at various types $\mathsf{a}$ (shown in the third row).

5.5 Example: Computing the Laplacian in an arbitrary coordinate system

Returning to the example of the potential field, Equation (5.5) and $\mathsf{aff}$ tell us that its covariant derivative is equal to its partial derivative, regardless of the value of $\Gamma$ : $\nabla{_i}P$ = $\partial{_i}P$ . Indeed, its domain and codomain are both $\mathbf{1}$ , and therefore the affinities are both zero. (To compute the value of $\nabla{_i}P$ in a given coordinate system, one would still need to multiply by the Jacobian as indicated earlier.)

However, when computing the second derivative, a non-zero affinity arises (because the 1st derivative has a non-unit domain.) One has therefore:

(5.6) \begin{equation}\nabla_{i}\nabla^i P =g'{^i}{^j}\nabla{_i}\nabla{_j}P=g'{^i}{^j}\partial{_i}\partial{_j}P - g'{^k}{^l}\Gamma{_k}{_l}{^m}\partial{_m}P\end{equation}

In diagram notation:

As an illustration, let us compute the Laplacian of the scalar field defined as growing with the negated logarithm of the distance to the origin: $P = -\log(\rho)$ .

Its covariant derivative (or gradient) is given by the partial derivatives in polar coordinates (because $\mathsf{aff}\mskip 3.0mu\mathbf{1}\mskip 3.0mu{=}\mskip 3.0mu\mathrm{0}$ , as already mentioned). The components in the polar tangent basis are

\begin{equation*}\nabla P = \partial P - 0 + 0 = \left[\begin{array}{cc}\partial P/\partial\rho&\partial P/\partial\theta\end{array}\right]=\left[\begin{array}{cc}-\rho^{-1}&0\end{array}\right]\end{equation*}

We have seen above that, by composing with the appropriate Jacobian, the components of the gradient in the Cartesian tangent basis are $\left[\begin{array}{cc}-\rho^{-1}\cos{\theta}&-\rho^{-1}\sin{\theta}\end{array}\right]$ .

The second derivative ( $\nabla\nabla P$ ) is a second-order tensor, and computing it manually in the Cartesian basis is error-prone. In contrast, because P has radial symmetry, in the polar tangent basis its partial derivative $\partial P$ has a simple expression (the co-vector $\left[\begin{array}{cc}-\rho^{-1}&0\end{array}\right]$ ). Thus, only the $\rho-\rho$ -component of the second partial derivative is non-zero: $\partial^2P/\partial\rho^2 = \rho^{-2}$ . To compute the covariant derivative, we also need the affinity term, which is obtained by multiplying the gradient $\nabla P$ by $\Gamma$ , then multiply by the contravariant metric (recall Equation 5.6). Because only the $\mathbf e_\rho$ coefficient of $\nabla P$ is non-zero, is suffices to multiply this component ( $-\rho^{-1}$ ) by $\Gamma^\rho$ to obtain the affinity term $\left[\begin{array}{cc}0&0\\0&1\end{array}\right]$ . Finally, we get $\nabla\nabla P$ (in polar tangent coordinates) from Equation (5.5): the partial derivative minus the affinity we computed above plus the second affinity which is zero (because the output space is $\mathbf{1}$ ) which gives us $\left[\begin{array}{cc}\rho^{-2}&0\\0&-1\end{array}\right]$ .

The Laplacian is obtained by 2-way contraction with the contravariant metric. One way to do this is to multiply the above by the contravariant metric, to get $\left[\begin{array}{cc}\rho^{-2}&0\\0&-\rho^{-2}\end{array}\right]$ , and take the trace, which is zero for $\rho>0$ . (Note that even if it were non-zero there was no need to multiply by any Jacobian because the divergence field is scalar.) The fact that the Laplacian of this scalar field is zero means that the charge density is zero away from the origin: in a two-dimensional space the potential of a point charge at the origin is proportional to this scalar field P.

Beyond its illustrative benefits, a takeaway from this example is that all computations were free of trigonometry, and involved many zeros, thanks to the choice of a global coordinate system which had a zero gradient in one of the axes. Even though each step in the computations was presented for illustrative purposes, the $\mathsf{MetricCategory}\mskip 3.0muM_{\mathsf{Polar}}$ instance meant that we could run all the above computations as Haskell programs.

5.6 Tensor calculus summary

Athis point we have presented all the classes whose morphisms constitute the combinators of Roger. Their relationships are summarised in Figure 8.

Fig. 8: Inheritance relationships between tensor structures. In the praxis, every instance of a category with a metric $\mathsf{Additive}$ is also additive. However, diagrams don’t support a good representation for addition, so the dotted line is implemented as a subclass relationship in our library.

When it comes to Albert, we have two separate sub-languages. First, we have a number of combinators which only manipulate indices, shown in Figure 9. Indices can be split, merged, raised and lowered. Indices for unit vector-spaces are unimportant and can be created or discarded at will. Second, we have a number of combinators which nominally manipulate scalar-valued expressions (addition, multiplication, embedding of constant tensors, Kronecker delta, and contraction, etc; see Figure 10). Various combinators require various amounts of structure in the underlying category z.

The semantics in terms of morphisms is provided by the $\mathsf{tensorEval}$ function, and $\mathsf{tensorEval}_{1}$ for the special case of closed tensor expressions. For tensor fields, primitives for metrics and derivatives are also available.

6.1 Code for Equation (6.1)

A contraction of the Riemann curvature occurs twice in Equation (6.1), and as such deserves to be extracted as an intermediate definition. It is called the Ricci tensor.

This tensor can be contracted one more time to obtain what is called the scalar (or Gaussian) curvature:

From there, the left-hand-side of Equation (6.1) is defined as follows:

We can then convert it to a morphism:

In the above, the operator & is the linear post-fix application:

6.2 Point-mass example (the Schwarzschild metric)

To complete our test case for Albert, in this section, we define the metric which describes the gravitational effects of a point-sized mass: the Schwarzschild metric. We then verify that this metric satisfies of Equation (6.1) by evaluating both sides and checking that they match.

The first step is to define the coordinate system. Roughly speaking, Schwarzschild coordinates are spherical coordinates with an extra component for time.

The point mass will be located at the origin ( $\rho=0$ ) at every point in time. In other words, the mass is at rest in this coordinate system. Hence, because we have a point-mass, the $T{_i}{_j}$ tensor is zero everywhere except at the origin, where it is infinite.

As for all coordinate representations, the compact-closed category structure for $\mathsf{Matrix}\mskip 3.0mu\mathsf{Spherical}$ falls out directly from Section 3.1. Likewise for the partial derivative. The only missing piece of the structure is the metric. The Schwarzschild metric is defined in terms of the considered mass M or alternatively by the Schwarzschild radius $r_s$ , the two being connected by the equation $r_s = \kappa Mc^2/4\pi$ . We use the parameter $r_s$ in the rest of the section. The metric is then given by tabulating the following function as a matrix.

We refer the reader to a course in general relativity for the physical meaning of it. We can then define the $\mathsf{Matrix}\mskip 3.0mu\mathsf{Spherical}$ instance as follows:

At this point, we can directly evaluate $\mathsf{grLhsM}$ with $z\mskip 3.0mu{=}\mskip 3.0mu\mathsf{Matrix}\mskip 3.0mu\mathsf{Spherical}$ , and obtain a 4 by 4 matrix of symbolic expressions depending on $\mathsf{Spherical}$ coordinate variables. We find that it simplifies to zero everywhere it is defined.Footnote ²³ Thus, we can verify that the Schwarzschild metric satisfies the general relativity equation.

7.1 Complication: Supporting derivatives

The main issue with the above implementation sketch arises when trying to implement the covariant derivative of a tensor expression $\mathsf{t}$ as a sum of its partial derivative and affinity terms. Indeed, the affinity terms of a tensor expression $\mathsf{t}$ depends on the number (and type) of free index variables in it. But this information is not made available by the library of Bernardy & Spiwack (Reference Bernardy and Spiwack2021) which we base our work upon. Simply put, the type $\mathsf{P}\mskip 3.0muz\mskip 3.0mu\mathsf{r}\mskip 3.0mu\mathbf{1}$ is opaque. To explain the solution that we employ, we need to first peek inside the implementation of this library.

Its principle is that a port (in our application, an index) is represented as a morphism in the (free) Cartesian extension of the underlying category z:

The $\mathsf{encode}$ function simply embeds the morphism in the Cartesian extension. The $\mathsf{decode}$ function takes a linear function $\mathsf{f}\mskip 3.0mu{:}\mskip 3.0mu\mathsf{P}\mskip 3.0muz\mskip 3.0mu\mathsf{r}\mskip 3.0mu\mathsf{a}\mskip 3.0mu{\multimap }\mskip 3.0mu\mathsf{P}\mskip 3.0muz\mskip 3.0mu\mathsf{r}\mskip 3.0mu\mathsf{b}$ and applies it to the identity morphism ( $\mathsf{id}\mskip 3.0mu{:}\mskip 3.0mu\mathsf{P}\mskip 3.0muz\mskip 3.0mu\mathsf{a}\mskip 3.0mu\mathsf{a}$ ), obtaining $\mathsf{f}\mskip 3.0mu{(}\mskip 0.0mu\mathsf{P}\mskip 3.0mu\mathsf{id}\mskip 0.0mu{)}\mskip 3.0mu{:}\mskip 3.0mu\mathsf{CartesianExt}\mskip 3.0muz\mskip 3.0mu\mathsf{a}\mskip 3.0mu\mathsf{b}$ . The non-obvious property proven by Bernardy & Spiwack (Reference Bernardy and Spiwack2021) is that, because $\mathsf{f}$ is linear, $\mathsf{f}\mskip 3.0mu{(}\mskip 0.0mu\mathsf{P}\mskip 3.0mu\mathsf{id}\mskip 0.0mu{)}$ is always equivalent to a morphism in the underlying category z; it never needs to refer to projections or duplications. The function which recovers this morphism has the following signature, but remember that it will crash if the input morphism makes essential use of the Cartesian structure:

To be able to support tensor derivatives, we cannot rely on $\mathsf{decode}$ , and we will have to access the internal function $\mathsf{cartesianToMonoidal}$ directly.

With this in mind, we can solve our problem, namely finding an implementation for $\mathsf{R}\mskip 3.0muz\mskip 3.0mu\mathsf{r}$ which lets us track free variables. We do this by embedding a pair of morphisms $\mathsf{t}\mskip 3.0mu{:}\mskip 3.0mu\mathsf{x}\,\mskip 1.0mu\overset{z}{\leadsto }\mskip 1.0mu\,\mathbf{1}$ and $\mathsf{p}\mskip 3.0mu{:}\mskip 3.0mu\mathsf{CartesianExt}\mskip 3.0muz\mskip 3.0mu\mathsf{r}\mskip 3.0mu\mathsf{x}$ (instead of just $\mathsf{p}$ ). Note that the type $\mathsf{x}$ is existentially bound, not a parameter to $\mathsf{R}$ .

Here, $\mathsf{x}$ is the context used by the tensor $\mathsf{t}$ , without any spurious component. The tensor $\mathsf{t}$ carries the payload of the expression. Because it inhabits the non-Cartesian category z, we know that it uses the whole context $\mathsf{x}$ . Even though the type declares that the morphism $\mathsf{p}$ is built from the Cartesian extension of z, the implementation will arrange that it will only use the sub-category of z whose derivative is zero. This is enforced by hiding the $\mathsf{encode}$ function on $\mathsf{P}\mskip 3.0muz\mskip 3.0mu\mathsf{r}\mskip 3.0mu\mathsf{a}$ from the user: we only provide the $\mathsf{raise}$ and $\mathsf{lower}$ functions to manipulate this type, and they are safe because $\nabla g=0$ :

With this setup the derivative can be computed correctly on a pair $\mathsf{Compose}\mskip 3.0mu\mathsf{t}\mskip 3.0mu\mathsf{p}$ . Indeed, the morphism-level derivative () need only be applied to the $\mathsf{t}$ component of the pair. As for the $\mathsf{p}$ component of the pair, the use of the Cartesian fork operator ${(}\mskip 0.0mu{ \triangle }\mskip 0.0mu{)}$ ensures that the new index is passed to the right component of $\nabla\mskip 3.0mu\mathsf{t}$ .

The transcoding functions between Albert and its categorical semantics then become:

Multiplication remains straightforward to implement:

7.2 Complication: Supporting addition of tensors

The type of the addition operator ( $\mathsf{plus}\mskip 3.0mu{::}\mskip 3.0mu{(}\mskip 0.0mu\mathsf{Bool}\mskip 3.0mu{\rightarrow }\mskip 3.0mu\mathsf{R}\mskip 3.0muz\mskip 3.0mu\mathsf{r}\mskip 0.0mu{)}\mskip 3.0mu{\multimap }\mskip 3.0mu\mathsf{R}\mskip 3.0muz\mskip 3.0mu\mathsf{r}$ ) poses a problem for the above implementation. The issue is that the semantics of $\mathsf{plus}\mskip 3.0mu\mathsf{f}$ depends on both $\mathsf{f}\mskip 3.0mu\mathsf{False}$ and $\mathsf{f}\mskip 3.0mu\mathsf{True}$ . However, $\mathsf{plus}$ is linear in $\mathsf{f}$ , and thus we can call $\mathsf{f}$ only once in the implementation of $\mathsf{plus}$ . The way to work around this problem is to embed a non-linear implementation inside a linear datatype. That is, it suffices to add a constructor $\mathsf{Plus}\mskip 3.0mu{::}\mskip 3.0mu{(}\mskip 0.0mu\mathsf{Bool}\mskip 3.0mu{\rightarrow }\mskip 3.0mu\mathsf{R}\mskip 3.0muz\mskip 3.0mu\mathsf{r}\mskip 0.0mu{)}\mskip 3.0mu{\multimap }\mskip 3.0mu\mathsf{R}\mskip 3.0muz\mskip 3.0mu\mathsf{r}$ to the $\mathsf{R}$ datatype. Operationally, if $\mathsf{f}$ is the argument to $\mathsf{plus}$ , it is stored once in the $\mathsf{R}$ data structure, and it is only at the moment of evaluating $\mathsf{plus}\mskip 3.0mu\mathsf{f}$ to a morphism that we invoke $\mathsf{f}$ . But, at the evaluation point, there is no linearity restriction on the tensor expression of type $\mathsf{R}\mskip 3.0muz\mskip 3.0mu\mathsf{r}$ : it can be used an arbitrary number of times, and, according to the typing rules of Linear Haskell, it means that its payload ( $\mathsf{f}$ ) has itself no usage restrictions. The fact that the evaluation function has no linearity restriction on its tensor expression argument ( $\mathsf{t}$ ) means that $\mathsf{t}$ cannot have any free index variables in it. At first sight, this is a drawback. However, there is a good reason for restricting evaluation to closed tensor expressions: it ensures that index variables do not escape the scope of the tensor expression where they belong.

The rest of the implementation of Albert needs to be modified to support this new constructor, by defining a case for sums. Fortunately, such a case is never difficult to handle: we simply distribute every operation over the operands of the sum.

7.3 Complication: Contraction

Unfortunately, our definition for $\mathsf{R}$ so far does not in fact support contraction. This is because all indices are inputs of the expression, and to connect inputs together we need to use the $\eta$ combinator. To do so, we need an additional $\mathbf{1}$ input—but no such input is explicitly available in the $\mathsf{R}$ type. The workaround for the problem is to add such a unit input explicitly. Hence, the final implementation of the $\mathsf{R}$ type is

The implementation of $\mathsf{contract}$ is then:

The other combinators can simply ignore or thread this unit input. The core of the final implementation is shown in Figure 11.

Fig. 11: Final implementation

8.1 Tensor presentations in introductory texts

We expect many of our readers not to be already familiar with tensors, and therefore they will need to read pedagogical introductions to the topic, as the authors did when preparing this paper. We believe that a warning is in order, because, typically, introductory texts lean heavily on the representations of tensors as (generalised) arrays of coefficients. Consequently, undue importance is enthused to what happens under change of coordinates in the manifold—even though from an algebraic perspective, coordinate systems do not even enter the picture. All the introductory texts that we could find take this approach (Dullemond & Peeters, Reference Dullemond and Peeters2010; Fleisch, Reference Fleisch2011; Grinfeld, Reference Grinfeld2013; Porat, Reference Porat2014; Rowland & Weisstein, Reference Rowland and Weisstein2023). The quotes below are by no means intended as singling out particular authors: this approach is pervasive. These are the kinds of definition that we find:

At the end of the document one finds out that the transformation rules are the multiplication by Jacobians. Other sources take a more pedagogical approach and start with vectors and covectors:

Here y is considered to be some array of numbers, enthused with a property (covariance), which refers to the fact that the inverse of the Jacobian should multiply the coefficients when changing coordinate systems. (Furthermore, $\alpha$ is implicitly lambda bound here—a kind of liberty that most textbooks take in all areas of mathematics.)

For the student already familiar with linear algebra, this can be particularly confusing, because every vector or matrix transforms this way. So why bother highlighting this property? This presentation continues with

In algebraic terms, it is not hard to state that the metric is a linear function of two arguments, but the representational approaches incites the authors to beat around the bush. Besides, the representational approaches defines the metric as the Gram matrix of the base vectors, and it being a tensor requires some proof.

In the algebraic definition this kind of pitfall is avoided. Because of the awkward definition of tensors in the representational view, there is much discussion about what is and what isn’t a tensor. One can find statements like the following:

Algebraically, this is attempting to add objects of different type together. In Albert, the above expression is not well-typed either. (And yet, to confuse matters even more, this expression can be made sense of with the pervasive index juggling conventions.)

But there is a more subtle way in which something might be “not a tensor”:

We believe that this kind of statement is puzzling to the novice. Indeed, the first sentence states that because $\Gamma$ (the Christoffel symbol) is zero in some bases, then it cannot be a tensor (implicit to the above is that no Jacobian can transform a zero tensor to a nonzero tensor). The student might still wonder if this statement applies for a manifold that cannot be straightened (i.e. one that is not flat), and why it can be handled like a tensor anyway.

In the algebraic view that we employ, $\Gamma$ is a morphism like any other. So is it a tensor after all? In fact, the value of $\Gamma$ for any given basis is a tensor. But $\Gamma$ is defined as an expression which explicitly references the basis, and changes with it. The confusing aspect is that the coefficient representation of any tensor changes when the basis changes. But $\Gamma$ changes as a geometric object as the coordinate system changes, which means that it cannot be converted to another basis by the usual jacobian-based transformation while retaining the properties of the Christoffel symbol.

What we find puzzling is that the choice of the representational approach appears to be a conscious one:

It is not to say that the algebraic approach is absent from the literature. While it appears to be chiefly geared towards specialists, it also can be found in textbooks, but we could not find one that explicitly makes the link between all notations. For instance, Bowen & Wang (Reference Bowen and Wang1976) take an algebraic approach, and as such only manipulate the category-oriented (point-free) notations. Using only this notation is sufficient for them because they do not manipulate any complicated tensor expressions, but it also means that readers will have a hard time connecting other notations to the point-free language. Jeevanjee (Reference Jeevanjee2011) does better by describing the tensor algebras, and discussing various representations. However, the equivalence between the various languages is not mentioned.

As mentioned earlier Thorne & Blandford (Reference Thorne and Blandford2015) do mention this equivalence informally. We refer the more advanced reader to Bleecker (Reference Bleecker2005) for a full development in the language of linear algebra and differential forms, without any connection to the underlying categorical structures.

8.2 Einstein notation

The Einstein notation appears to arise as a generalisation of the notation to denote elements of matrices. This indexing notation is so well established that we could not trace where it originates.

However, not every expression which involves accessing indices can be mapped to a tensor. In a nutshell, indices must be treated abstractly and used linearly. Thus, under these conditions, which are exactly those of Albert, there is an equivalence between Einstein notation and the algebraic specification of tensors as morphisms.

Despite our best efforts, some distance remains between Albert and standard Einstein notation. First the Einstein notation leaves the binders of indices implicit, but because Albert is embedded in a functional language, every index variable must be lambda bound. In particular, contraction is written as $\mathsf{contract}\mskip 3.0mu{(}\mskip 0.0mu\lambda\mskip 3.0mu\mathsf{^k}\mskip 3.0mu\mathsf{_k }\mskip 3.0mu{\rightarrow }\mskip 3.0mu{\ldots}\mskip 0.0mu{)}$ . One could imagine a pre-processing step to shorten the notation and more closely approach the Einstein notation, but we find that sticking to the lambda calculus convention avoids a source of confusion.

Second, it should be noted that the use of standalone subscript and superscript variables is a small liberty that we took in typesetting the paper: while the Haskell compiler accepts super and subscript characters in variable names, it disallows them as the first character. Third, the operands of an addition are $\mathbf{case}$ branches. We believe that an alternative syntax should be provided by ghc, because it is generally useful in presence of linear types. Fourth, a pattern syntax for $\mathsf{merge}$ and $\mathsf{split}$ would shorten index manipulation. Pattern synonyms are already available in ghc, but are not currently compatible with linear types. This is a technical shortcoming that could be addressed in Haskell implementations.

While this paper covers all the principles of the tensor notations, we have not aimed for an exhaustive coverage. In particular, we did not discuss the Levi-Civita tensor, which has many uses and is related to antisymmetrisation.

8.3 ${\mathbf{D}\scriptstyle \mathbf{SL}}$ s for array programming and scientific programming

There is overlap between the present work and languages oriented to scientific programming: they can both be used to describe array-oriented computations. Array-oriented programming languages have a long history, perhaps starting with apl (Iverson, Reference Iverson1962). Notable standalone languages include Single-Assignment C (Scholz, Reference Scholz1994) or sisal (Feo et al., Reference Feo, Cann and Oldehoeft1990) and even Matlab (Gilat, Reference Gilat2004). When it comes to array edsls, there is an abundance of libraries available. Limiting ourselves to the Haskell ecosystem, notable examples include RePa (Lippmeier et al., Reference Lippmeier, Keller, Chakravarty, Leshchinskiy and Peyton Jones2010) and accelerate (Chakravarty et al., Reference Chakravarty, Keller, Lee, McDonell and Grover2011).

An important difference between this work and most array-oriented languages is that they focus on representations first, and semantics come only as a means to support optimisations. In contrast, our approach puts categorical semantics at the core. This means that constructing a dataflow representation (typically as a free smc) is natural. This representation can then be optionally optimised and interpreted as operations on matrices and arrays, perhaps relying on the aforementioned libraries or a dedicated compiler (Kjolstad et al., Reference Kjolstad, Kamil, Chou, Lugato and Amarasinghe2017) as backend.

The implementation model that we suggest is to construct a dataflow graph (typically as a free categorical representation), optimise it and then interpret it as operations on matrices. This kind of model was already at the heart of sisal, but has been popularised recently in machine-learning applications by the TensorFlow library (Abadi et al., Reference Abadi, Barham, Chen, Chen, Davis, Dean, Devin, Ghemawat, Irving, Isard, Kudlur, Levenberg, Monga, Moore, Murray, Steiner, Tucker, Vasudevan, Warden, Wicke, Yu and Zheng2016). We note, however, that TensorFlow (and similar machine learning packages) do not offer an index-based notation.Footnote ²⁴ Instead, the programmer must keep track of dimensions using indices. The situation is similar to de Bruijn index representations for lambda terms. Our work, on the other hand, rejects non-tensor (non-linear) primitives, and as such would preclude many useful array operations. These would need to be added as an additional layer.

Another salient feature of the present work is its specific ability to support calculus. Most scientific programming languages provide arrays, but let the user figure out how to model tensor calculus operations. However, there are exceptions. Sussman and Wisdom (Reference Sussman and Wisdom2013) develop a lisp-based dsl for differential geometry. It is a point-free language close in spirit to Roger, even though the categorical structures remain implicit. Diderot (Chiw et al., Reference Chiw, Kindlmann, Reppy, Samuels and Seltzer2012) is specifically oriented towards tensor calculus, with explicit support for indices, as in Einstein notation. Cadabra (Peeters, Reference Peeters2006) is another computer-algebra system supporting tensor calculus expressions.

8.4 Categorical semantics

It is generally more convenient to provide instances of categorical structures rather than handling lambda terms directly. In this paper, we have used matrix and diagram instances, but many other applications exist (Elliott, Reference Elliott2017). A major selling point of categorical semantics is that they avoid the need of manipulating variables explicitly (perhaps as de Bruijn indices). This advantage is already identifiable in the work of Cousineau et al. (Reference Cousineau, Curien and Mauny1985), but Elliott (Reference Elliott2017) has shown how to leverage them in edsls. He does so by providing a compiler plugin that translates lambda terms of the edsl to morphisms in a closed Cartesian category. This last characteristic means that no support for linear types or specific support of smcs is available.

These shortcomings have been addressed by Bernardy & Spiwack (Reference Bernardy and Spiwack2021), who show how to evaluate linear functions to morphisms in an smc, with no compiler modification. While they lay down the foundation for the present work, their library is unfortunately not sufficient for our purposes. The technical additions provided by this paper include the ability to represent tensor derivatives (Section 7.1), tensor addition (Section 7.2) and contraction (Section 7.3).

Categorical semantics are particularly well suited to perform automatic differentiation (AD), as Elliott (Reference Elliott2018) has shown. We have shown that derivatives can be represented, and we expect their symbolic computation can be done with standard techniques. We have not shown how AD can be performed, but because we use categorical semantics, Elliott-style AD is the natural approach to implement our interface for derivatives. In fact, it is the method that we have used to implement the Schwarzschild metric example (Section 6.2).

8.5 Penrose diagram notation

The diagram notation that we have used is a (graphical) dsl in its own right. It is particularly suited for morphisms in smcs (and their extensions; such as compact closed categories). Selinger (Reference Selinger2011) provides a survey for various kinds of categories. The correspondence between definitional and topological equivalence is their main advantage: it means that topological intuitions can be leveraged for formal proofs (Blinn, Reference Blinn2002; Hinze, Reference Hinze2012). Kissinger (Reference Kissinger2012) even built a tool on this premise.

To our knowledge, the diagram notation was in fact first developed for tensor calculus by Penrose (Reference Penrose1971). Later it was generalised to represent morphisms in monoidal categories (Joyal & Street, Reference Joyal and Street1991). For tensors, there is a considerable amount of variation between the diagrammatic notation used by various authors. We have adhered to a standardised subset of this notation. Furthermore, we have ensured that each diagram is built strictly using well-defined building blocks, combined strictly using sequential ${(}\mskip 0.0mu{\circ }\mskip 0.0mu{)}$ and parallel ${(}\mskip 0.0mu{\otimes }\mskip 0.0mu{)}$ composition.

In contrast, the literature on tensor applications does not prescribe a direction for reading diagrams (whereas ours are always read left to right). As long as the various input and outputs of atomic morphisms are clearly identified, there is no problem in not specifying a direction for diagrams of compact closed categories because all directions are equivalent thanks to the snake laws.

In his seminal work, Penrose additionally does not make a graphical difference between the dual ${\mathsf{a}^*}$ and $\mathsf{a}$ . This means, for instance, that metrics are (graphically) indistinguishable from $\eta$ and $\epsilon$ . Most of the time the difference is inconsequential; however, it becomes important when expressing covariant derivatives in terms of partial derivatives and affinities. Indeed, the affinities for $\mathsf{a}$ and ${\mathsf{a}^*}$ are not the same.