Indices and Equations

One of the main features of xerus is the ability to write arbitrary Tensor contraction in a Einstein-notation like form. This allows to write very readable sourcecode that closely resembles the mathematical formulas behind the algorithms.

Indices

To write indexed equations, we first have to declare the variables that we will use as indices.

xerus::Index i,j,k,l;

i,j,k,l = xerus.indices(4)

Simple Equations

The most basic equations utilizing indexed expressions are reshufflings and contractions. Assuming A is a tensor of degree 2 and bof degree 1 we can write:

// transposing a matrix
A(i,j) = A(j,i);

// a simple matrix-vector product
xerus::Tensor c;
c(i) = A(i,j) * b(j);

# transposing a matrix
A(i,j) << A(j,i)

# a simple matrix-vector product
c = xerus.Tensor()
c(i) << A(i,j) * b(j)

In analogy to the Einstein notation, such an expression will contract those modes on the right hand side that are indexed by the same index and assign the result to the left hand side. If necessary it is reshuffled first to obtain the index order denoted on the left hand side.

The left hand side (c in the above example) is not required to have the right degree or dimensions. The type of c on the other hand does change the meaning of the equation and so has to be set correctly. E.g. if c is a xerus::TensorNetwork, no contraction is performed as A(i,j)*b(j) is in itself a valid tensor network. See also the tutorials on TT-Tensors and Tensor Networks for details on the respective assignments.

try {
    c(j) = A(i,j) * b(j); // runtime error!
} catch(xerus::misc::generic_error &e) {
    std::cout << "something went wrong: " << e.what() << std::endl;
}

try:
    c(j) << A(i,j) * b(j) # runtime error!
except xerus.generic_error as err:
    print("something went wrong:", err)

// do NOT do this!
auto evil = A(i,j) * b(j);

# do NOT do this!
evil = A(i,j) * b(j)

Summation, subtraction, multiplication by scalars and calculating a norm all work as one would expect in such equations.

// perform a single gradient step of stepsize alpha
// x' = x + \alpha * A^T * (b - A*x)
x(i) = x(i) + alpha * A(j,i) * (b(j) - A(j,k)*x(k));

# perform a single gradient step of stepsize alpha
# x' = x + \alpha * A^T * (b - A*x)
x(i) << x(i) + alpha * A(j,i) * (b(j) - A(j,k)*x(k))

High-Dimensional Equations

Writing high-dimensional equations in Einstein notation with individual indices can be cumbersome. What is more: in general functions in our code we might not even know the degree of the tensors arguments before runtime. To still be able to use indexed equations in these settings xerus uses multi-indices that can span anything between 0 and all modes of a tensor.

To denote such multi-indices, the regular indices defined as above have to be modified by operators inside the indexed expressions.

i^d // an index that spans d modes

i&d // an index that spans all but d modes
i/n // an index that spans degree()/n modes

i^d  # an index that spans d modes
i**d # an index that spans d modes
i&d  # an index that spans all but d modes
i/n  # an index that spans degree()/n modes

As it might be difficult to tell how many modes the indices on the left hand side of the equation will span, it is not necessary to declare them as multi-indices.

// contract the last mode of A with the first of B
C(i,k) = A(i&1, j) * B(j, k&1);

// C is now of degree A.degree()+B.degree()-2

# contract the last mode of A with the first of B
C(i,k) << A(i&1, j) * B(j, k&1)

# C is now of degree A.degree()+B.degree()-2

The division i/n is useful for example to write equations with high dimensional operators such as TT-Operators for which the indices are ordered per default such, that the application of the operator can be written in analogy to matrix-vector products as:

// assumes e.g.: TTTensors u,v; TTOperator A; of any compatible degree
u(i&0) = A(i/2, j/2) * v(j&0);

# assumes e.g.: TTTensors u,v; TTOperator A; of any compatible degree
u(i&0) << A(i/2, j/2) * v(j&0)

Notice how the declaration of multi-indices with & or / are always relative to the current tensor. As such, an index (such as j in the above example) can appear with different modifiers. The important thing is, that each occurence of every index results in an equal number of modes spanned and that these modes are of equal dimensions (in the example the placement of j requires that the degree of v is equal to half the degree of A and that its dimensions are equal to the second half of dimensions of A).

Blockwise Construction of Tensors

A common use for indexed expressions is to construct tensors in a blockwise fashion. In the following example we were able to calculate the tensor comp whenever the first index was fixed, either by numerical construction (A and B) or by showing mathematically, that it is then equal to a well known tensor (here the identity matrix). The full tensor can thus be constructed with the help of the named constructors of xerus::Tensor (see the Tensor tutorial) as follows.

// construct comp s.th.:
// comp(0, :,:) = A+identity
// comp(1, :,:) = B+identity
// comp(2, :,:) =   identity
comp(i, j^2) = 
      xerus::Tensor::dirac({3}, 0)(i) * A(j^2) 
    + xerus::Tensor::dirac({3}, 1)(i) * B(j^2)
    + xerus::Tensor::ones({3})(i) * xerus::Tensor::identity({64,64})(j^2);

# construct comp s.th.:
# comp(0, :,:) = A+identity
# comp(1, :,:) = B+identity
# comp(2, :,:) =   identity
comp(i, j**2) << \ 
      xerus.Tensor.dirac([3], 0)(i) * A(j**2) \
    + xerus.Tensor.dirac([3], 1)(i) * B(j**2) \
    + xerus.Tensor.ones([3])(i) * xerus.Tensor.identity([64,64])(j**2)