How to create tensors at runtime?

Question

I have a small personal CAS project around linear algebra which involves, among other things, manipulating tensors. When I talk about tensors, I mean them according to the physical / mathematical definition, so this is of little importance for the rest of the question, but I would like to make it clear that what I want to manipulate is not simply arrays of any dimensions (a tensor is not a multidimensional matrix, even if I will later contradict myself, but for good reasons).

Since I am trying to create a CAS, I have a DSL in which the user can enter expressions and define objects, such as tensors. The user can define a tensor quite symbolically, without specifying values for components, in which case the problem I am encountering does not arise.

In a very general way, a user who wishes to specify the components of a tensor will do so component by component, so the component of a tensor of order $N$ will be of the form $C_{i, j, ..., N}$ (if for the sake of simplicity we consider here only covariant components). For example, if I define a 3 times covariant tensor $T_{ijk}$ and I want to define its component $T_{3,2,1}$ with any expression $x$, I would write $T_{3, 2, 1} = x$. In practice, one will often use matrices to do this (or a new tensor will arise from another operation, as a tensor product), but I want to be as general and abstract as possible, if possible.

So that's where my problem lies. If a user, for whatever reason, decides to create (or generate by a particular operation) a tensor of order N, the DSL must be able to represent internally the different components that constitute it in an adequate way. Even if it is not mathematically correct to define a tensor as a multi-dimensional array, at a given time, it must be manipulated numerically and at runtime, so that it is possible for my implementation to access each component from a given position $i, j, ..., N$ where N is not known at compile-time.

I am therefore looking for a way to define a kind of multi-dimensional array at runtime, i.e., for example, do something like (C#) int[, ,] arr3d; int[, , ,] arr4d; but dynamically. I've been programming long enough to know that this is not possible (at least not in a clean way).

My idea is then to use a geometric approach and to consider only one list of elements representing all the components. The position, in the geometrical approach, of each of the terms of this list individually would then be defined by another list [i, j, ..., N] whose general form would be fixed (at runtime) at the creation of the tensor. To clarify my thinking, here is a little pseudo-code that defines the useful types:

/// The position of a component represented by a list i, j, ..., N
type Position = int list

/// The type for a single component
type 'K Component = { position: Position; value: 'K }

/// The numerical components of a tensor
type 'K Components = 'K Component list

For example, a simple scalar would then have the components $(i)$ where $i$ can only be 0 (otherwise it is a vector), so I would make

let myScalar (x: int) =
    let pos : Position = [0]
    let components : Component<int> = { position=pos; value=x }
    let result : Components<int> = [ components ]
    
    result

Another example, a vector:

let myVector (len: int) (values: int list) =
    assert (values.Length = len)

    // "a vector of size len", so "len" is a "border"
    let pos = [len]

    let result : Components<int> = []

    for i=0 to len do
        let c = { position=[i]; value=value[i] }
        result += c

    result

Another example, a matrix:

let myMatrix (n: int) (m: int) (values: int list) =
    assert (n * m = values.Length)

    // like a square, the borders are "n" and "m"
    let pos = [n; m]

    let result : Components<int> = []

    for i=0 to n do:
        for j=0 to m do:
            let c = { position=[i; j]; value=... }

    .....

This is obviously pseudo code (F#), but I really hope you get the idea... The problem with this approach is

1. it becomes difficult to determine "edges" at N-dimensions.
2. it is not possible to define new for loops at runtime.

Since then, I have been thinking about this problem, but I have not found any way to meet my expectations, so I summarise them:

• I want to be able to represent the components of a tensor for non-fixed positions at compile-time.
• I want to be able to access the $(i, j, ..., N)$-th component by mentioning it by means of a list (if it had been fixed at compile-time, I would have remapped the list by a tuple of course).

My question is then "simply" how to do what I have just explained, by completing my geometrical approach for example, or by pointing me to a completely different one.

PS: in fact I talk about tensors in my question to give a real context of my problem, if I had simply asked for an array with indeterminate dimensions at compile-time, I would have been advised to review my algorithms. By the way, if you have a better idea to represent such an object, I'm interested! I did a lot of research on the Internet, but I didn't find anything...

PSS: I wasn't very formal in the explanations I gave about tensors, but who cares, I don't do math here, I only implement them like I can on a computer that even no know about the existance of the number 3.

Answer 1

If you have an $N_1 \times \cdots \times N_d$ tensor, you can index the (zero-based) $(i_1,\ldots,i_d)$'th position using the formula $$ \sum_{j=1}^d i_j \prod_{k=j+1}^d N_k, $$ which you can compute in $O(d)$ if you're careful.

If you want to go over all positions, you can go over $I \in \{0, \ldots, \prod_{i=1}^d N_i - 1\}$, and extract the individual indices as follows: $$ i_j = \left\lfloor I\middle/\prod_{k=j+1}^d N_k \right \rfloor \bmod N_j. $$ Once again, you can compute this in $O(d)$ if you're careful.

The two operations are inverses of one another, and they are known as ranking and unranking.

There are other ways to over all positions. For example, you could start at the position $(0,\ldots,0)$, which you then increment as follows: start at $j = d$; increment $i_j$ modulo $N_j$; if there was an overflow, decrease $j$. If $j$ ever gets to zero, you have finished going over all entries.

---

In your case, you might also consider to write specialized code for tensors of small dimension, since these are the ones most likely to be encountered in practice. If you know the dimension, there is no problem with writing the nested loops as usual.