Tensors

Tensor Representation and Order

Column-Centric Approach

SQUINT uses a unique approach to represent n-dimensional data. From here on we’ll refer to this approach as column-centric. Column-centric is opposed to the more common row-centric approach used in almost every other language that deals with multidimensional data. In the opinion of the author, even though this may at first feel unnecessarily counter to the established conventions, the column-centric approach is more intuitive when viewed from a purley mathematical perspective and leads to simpler implementations of mathematical and physical concepts in many cases.

Note

It is important to note that the column-centric and row-centric approaches are orthogonal concepts from memory layouts. In fact, there are analogous memory layouts in both approaches to the traditional memory layouts row-major and column-major.

Representation

In SQUINT’s column-centric approach:

1. A 1st order tensor (vector) is represented as a column with shape m

2. A 2nd order tensor (matrix) has shape m × n (a single row being of shape 1 × n)

3. A 3rd order tensor has shape m × n × l

4. Higher-order tensors follow this pattern, adding dimensions to the right

This means that indexing a tensor is done by specifying the row index first, followed by the column index, and so on for higher-order tensors. The shape of a tensor is defined by the number of elements in each dimension, starting from the leftmost dimension.

\[\begin{split}\begin{array}{|c|c|c|} \hline \text{Tensor Order} & \text{Representation} & \text{Shape} \\ \hline \text{1st order (vector)} & \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{m-1} \end{bmatrix} & m \\ \hline \text{2nd order (matrix)} & \begin{bmatrix} a_{00} & a_{01} & \cdots & a_{0,n-1} \\ a_{10} & a_{11} & \cdots & a_{1,n-1} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m-1,0} & a_{m-1,1} & \cdots & a_{m-1,n-1} \end{bmatrix} & m \times n \\ \hline \text{3rd order} & \begin{bmatrix} \begin{bmatrix} a_{000} & a_{010} & \cdots & a_{0,n-1,0} \\ a_{100} & a_{110} & \cdots & a_{1,n-1,0} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m-1,0,0} & a_{m-1,1,0} & \cdots & a_{m-1,n-1,0} \end{bmatrix} \\ \begin{bmatrix} a_{001} & a_{011} & \cdots & a_{0,n-1,1} \\ a_{101} & a_{111} & \cdots & a_{1,n-1,1} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m-1,0,1} & a_{m-1,1,1} & \cdots & a_{m-1,n-1,1} \end{bmatrix} \\ \vdots \\ \begin{bmatrix} a_{0,0,l-1} & a_{0,1,l-1} & \cdots & a_{0,n-1,l-1} \\ a_{1,0,l-1} & a_{1,1,l-1} & \cdots & a_{1,n-1,l-1} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m-1,0,l-1} & a_{m-1,1,l-1} & \cdots & a_{m-1,n-1,l-1} \end{bmatrix} \end{bmatrix} & m \times n \times l \\ \hline \end{array}\end{split}\]

Indexing

• 1st order tensor: t[i] where i is the row index

• 2nd order tensor: t[i, j] where i is the row index and j is the column index

• 3rd order tensor: t[i, j, k] where i is the row index, j is the column index, and k is the depth index

Comparison with Row-Centric Approaches

To visualize the difference between column-centric and row-centric approaches for tensor representation, consider the following diagram:

Column-centric vs Row-centric Tensor Representation

Column-centric vs Row-centric Tensor Representation

This diagram illustrates how 1st and 2nd order tensors are represented in both column-centric and row-centric approaches. Note how the column-centric approach emphasizes columns as the primary structure, while the row-centric approach emphasizes rows.

The column-centric approach has memory layouts that are analogous to the row-centric layouts but with a different ordering of dimensions. Consider a simple example with the following tensors and their representations below:

\[\begin{split}\begin{equation} \begin{bmatrix} a \\ b \end{bmatrix}\tag{A} \end{equation}\end{split}\]

\[\begin{split}\begin{equation} \begin{bmatrix} a & c \\ b & d \end{bmatrix}\tag{B} \end{equation}\end{split}\]

\[\begin{split}\begin{equation} \begin{bmatrix}\begin{bmatrix} a & c \\ b & d \end{bmatrix} \\ \begin{bmatrix} e & g \\ f & h \end{bmatrix}\end{bmatrix}\tag{C} \end{equation}\end{split}\]

Tensor Memory Layout Comparison

Tensor	Row-Centric Row-Major	Row-Centric Column-Major	Column-Centric Row-Major	Column-Centric Column-Major
A	Memory: [a, b] Strides: [1, 1] Shape: [m, 1] Indexing: t[row, column]	Memory: [a, b] Strides: [1, 1] Shape: [m, 1] Indexing: t[row, column]	Memory: [a, b] Strides: [1] Shape: [m] Indexing: t[row]	Memory: [a, b] Strides: [1] Shape: [m] Indexing: t[row]
B	Memory: [a, c, b, d] Strides: [n, 1] Shape: [m, n] Indexing: t[row, column]	Memory: [a, b, c, d] Strides: [1, m] Shape: [m, n] Indexing: t[row, column]	Memory: [a, c, b, d] Strides: [n, 1] Shape: [m, n] Indexing: t[row, column]	Memory: [a, b, c, d] Strides: [1, m] Shape: [m, n] Indexing: t[row, column]
C	Memory: [a, c, b, d, e, f, g, h] Strides: [m*n, n, 1] Shape: [l, m, n] Indexing: t[depth, row, column]	Memory: [a, e, b, f, c, g, d, h] Strides: [1, m, m*n] Shape: [l, m, n] Indexing: t[depth, row, column]	Memory: [a, e, c, h, b, f, d, g] Strides: [m*n, n, 1] Shape: [m, n, l] Indexing: t[row, column, depth]	Memory: [a, b, c, d, e, f, g, h] Strides: [1, m, m*n] Shape: [m, n, l] Indexing: t[row, column, depth]

Note

In the 1D case (A), all representations are equivalent with the only difference being an additional dimesion of 1 for the columns needs to be added to the shape of the row-centric views in order to represent a column. In the 2D case (B), the column-major and row-major layout are equivalent to their corresponding layout in the other approach. In the 3D case (C), the column-centric approach adds the new dimension to the right, while the row-centric approach adds the new dimension to the left. This difference is reflected in the memory layout and indexing.

SQUINT uses the column-centric approach with the column-major layout by default since this is the most straightforward view when we maintain the idea of columns as the fundamental building blocks of tensors. You can specify any sequence to represent the strides and shape of the tensor, which allows you to use any approach with any memory layout you prefer.

Tensor Construction

SQUINT provides several ways to construct tensors:

1. Using initializer lists:

// Alias types use the default column-centric approach and column-major layout
mat2x3 A{1, 4, 2, 5, 3, 6};
// A represents:
// [1 2 3]
// [4 5 6]

2. Factory methods:

auto zero_matrix = mat3::zeros();
auto ones_matrix = mat4::ones();
auto identity_matrix = mat3::eye();
auto random_matrix = mat3::random(0.0, 1.0);
// and more ...

3. Element-wise initialization:

mat3 custom_matrix;
for (size_t i = 0; i < 3; ++i) {
    for (size_t j = 0; j < 3; ++j) {
        custom_matrix(i, j) = i * 3 + j;  // Note the use of () for element access
    }
}

4. Construction from other tensors or views:

mat3 original{{1, 2, 3, 4, 5, 6, 7, 8, 9}};
mat3 copy(original);

mat4 big_matrix = mat4::random(0.0, 1.0);
mat3 sub_matrix(big_matrix.subview<3, 3>(0, 0));

5. Dynamic tensor construction:

std::vector<size_t> shape = {3, 4, 5};
dynamic_tensor<float> dynamic_tensor(shape);
dynamic_tensor<float> filled_tensor(shape, 1.0f);

6. Tensor construction with quantities:

vec3_t<length_t<double>> position{
    units::meters(1.0),
    units::meters(2.0),
    units::meters(3.0)
};

Basic Operations

SQUINT supports a wide range of operations for tensors:

auto C = A + B;  // Element-wise addition
auto D = A * B;  // Matrix multiplication
auto E = A * 2.0;  // Scalar multiplication
auto F = A / B;  // General least squares or least norm solution

// Element access (note the use of () for multi-dimensional access)
auto element = A(1, 2);  // Access element at row 1, column 2

// Iteration (column-major order by default)
for (const auto& element : A) {
    // Process each element
}

// Iteration of rows
for (const auto& row : A.rows()) {
    // Process each row
}

// Iteration of views
for (const auto& view : A.subviews<2,3>()) {
    // Process each view
}

For matrix multiplication, the operation performed is:

\((AB)_{ij} = \sum_{k=1}^n A_{ik}B_{kj}\)

Views and Reshaping

SQUINT provides powerful view and reshaping capabilities:

auto view = A.view();  // Create a view of the entire tensor
auto subview = A.subview<2, 2>(0, 1);  // Create a 2x2 subview starting at (0, 1)
auto reshaped = A.reshape<6>();  // Reshape to a 1
D tensor
auto transposed = A.transpose();  // Transpose the tensor
auto permuted = A.permute<1,0>(); // Permutation of the tensor

// For dynamic tensors
auto dynamic_reshaped = dynamic_tensor.reshape({6, 4});
auto dynamic_transposed = dynamic_tensor.transpose();

Linear Algebra Operations

SQUINT provides comprehensive linear algebra operations:

• Solving Linear Systems:

auto result = solve(A, b);  // Solves Ax = b for square systems

This solves the system of linear equations:

\(Ax = b\)

A will be overwritten with the LU decomposition of A and b will be overwritten with the solution x.

• Least Squares / Least Norm Solution:

auto result = solve_general(A, b);  // Solves Ax = b for non-square systems

\(Ax = b\)

The system is solved in the least squares sense, where A is an m x n matrix with m >= n and in the least norm sense when m < n.

A will be overwritten with the QR decomposition of A and b will be overwritten with the solution x.

Note

b must have enough rows to store the solution.

• Matrix Inversion:

auto inverse = inv(A);  // Computes the inverse of a square matrix

The inverse \(A^{-1}\) satisfies:

\(AA^{-1} = A^{-1}A = I\)

• Pseudoinverse:

auto pseudo_inverse = pinv(A);  // Computes the Moore-Penrose pseudoinverse

For a matrix \(A\), the Moore-Penrose pseudoinverse \(A^+\) satisfies:

\(AA^+A = A\) \(A^+AA^+ = A^+\) \((AA^+)^* = AA^+\) \((A^+A)^* = A^+A\)

• Matrix Determinant:

auto determinant = det(A);

The determinant of a square matrix \(A\) is denoted as \(\text{det}(A)\).

Vector Operations

• Cross Product (for 3D vectors):

auto cross_product = cross(a, b);

For vectors \(a = (a_x, a_y, a_z)\) and \(b = (b_x, b_y, b_z)\):

\(a \times b = (a_y b_z - a_z b_y, a_z b_x - a_x b_z, a_x b_y - a_y b_x)\)

• Dot Product:

auto dot_product = dot(a, b);

For vectors \(a\) and \(b\):

\(a \cdot b = \sum_{i=1}^n a_i b_i\)

• Vector Norm:

auto vector_norm = norm(a);

The Euclidean norm of a vector \(a\) is:

\(\|a\| = \sqrt{\sum_{i=1}^n |a_i|^2}\)

Matrix Operations

• Trace:

auto matrix_trace = trace(A);

The trace of a square matrix \(A\) is:

\(\text{tr}(A) = \sum_{i=1}^n A_{ii}\)

Statistical Functions

• Mean:

auto tensor_mean = mean(A);

For a tensor \(A\) with \(n\) elements:

\(\text{mean}(A) = \frac{1}{n} \sum_{i=1}^n A_i\)

Tensor Contraction

• Tensor Contraction:

// for dynamic shape tensors
auto contraction_pairs = std::vector<std::pair<size_t, size_t>>{{I, J}};
auto contracted = contract(A, B, contraction_pairs);
// for fixed shape tensors
auto contracted = contract(A, B, std::index_sequence<I>{}, std::index_sequence<J>{});

For tensors \(A\) and \(B\), the contraction over indices \(i\) and \(j\) is:

\((A \cdot B)_{k_1...k_n l_1...l_m} = \sum_{i,j} A_{k_1...k_n i} B_{j l_1...l_m}\)

Einstein Summation (einsum)

SQUINT provides an implementation of Einstein summation convention through the einsum function. This powerful method allows for concise expression of many tensor operations including, but not limited to, transpose, matrix multiplication, and trace.

Syntax

For dynamic shape tensors:

auto result = einsum("subscript_string", tensor1, tensor2);

For fixed shape tensors:

auto result = einsum<InputSubscripts1, InputSubscripts2, OutputSubscripts>(tensor1, tensor2);

Where:

• subscript_string is a string specifying the operation in Einstein notation

• InputSubscripts1, InputSubscripts2, and OutputSubscripts are std::index_sequence representing the input and output subscripts

• tensor1 and tensor2 are the input tensors

For operations on a single tensor, use:

For dynamic shape tensors:

auto result = einsum("subscript_string", tensor);

For fixed shape tensors:

auto result = einsum<InputSubscripts, OutputSubscripts>(tensor);

Examples

1. Matrix Multiplication:

Dynamic shape:

auto result = einsum("ij,jk->ik", A, B);

Fixed shape:

auto result = einsum<seq<I, J>, seq<J, K>, seq<I, K>>(A, B);

2. Dot Product:

Dynamic shape:

auto result = einsum("i,i->", A, B);

Fixed shape:

auto result = einsum<seq<I>, seq<I>, seq<>>(A, B);

3. Outer Product:

Dynamic shape:

auto result = einsum("i,j->ij", A, B);

Fixed shape:

auto result = einsum<seq<I>, seq<J>, seq<I, J>>(A, B);

4. Trace:

Dynamic shape:

auto result = einsum("ii->", A);

Fixed shape:

auto result = einsum<seq<I, I>, seq<>>(A);

5. Diagonal:

Dynamic shape:

auto result = einsum("ii->i", A);

Fixed shape:

auto result = einsum<seq<I, I>, seq<I>>(A);

6. Permutation:

Dynamic shape:

auto result = einsum("ijk->kji", A);

Fixed shape:

auto result = einsum<seq<I, J, K>, seq<K, J, I>>(A);

The einsum function provides a unified interface for many tensor operations, allowing for concise and readable code. It automatically handles the necessary contractions and permutations based on the specified subscripts.

Tensor Error Checking

SQUINT provides optional error checking for tensors, which is separate from and orthogonal to error checking for quantities. When enabled, tensor error checking primarily focuses on bounds checking and additional shape checks at runtime, especially for dynamic tensors.

Enabling Error Checking

Error checking for tensors can be enabled by specifying the error_checking::enabled policy when declaring a tensor:

using ErrorTensor = squint::tensor<float, dynamic, dynamic, error_checking::enabled>
ErrorTensor t({2,3}, std::vector<float>{1, 4, 2, 5, 3, 6});

Types of Checks

When error checking is enabled for tensors, SQUINT performs the following types of checks:

1. Bounds Checking: Ensures that element access is within the tensor’s dimensions.

   // This will throw std::out_of_range
   t(2, 0);
   t(0, 3);
   

2. Shape Compatibility: Verifies that tensor operations are performed on compatible shapes.

   ErrorTensor a({2, 3});
   ErrorTensor b({3, 4});
   ErrorTensor c({2, 4});
   
   // This will compile and run correctly
   auto result1 = a * b;
   
   // This will throw a runtime error due to incompatible shapes
   auto result2 = a * c;
   

3. View Bounds: Ensures that tensor views and reshaping operations are within bounds.

   // This will throw if the subview exceeds the tensor's bounds
   auto subview = t.subview({2,2}, {1, 2});
   

Performance Considerations

While error checking provides additional safety, it does come with a performance cost. In performance-critical code, you may want to disable error checking:

using FastTensor = squint::tensor<float, squint::shape<2, 3>, squint::strides::column_major<squint::shape<2, 3>>, squint::error_checking::disabled>;
FastTensor ft{1, 4, 2, 5, 3, 6};

// No bounds checking performed, may lead to undefined behavior if accessed out of bounds
auto element = ft(1, 1);

Error Checking and Quantities

It’s important to note that tensor error checking is independent of quantity error checking. You can have tensors of quantities with different error checking policies:

// Tensor with error checking, containing quantities without error checking
tensor<length_t<double>, shape<3>, strides::column_major<shape<3>>, error_checking::enabled> t1;

// Tensor without error checking, containing quantities with error checking
tensor<quantity<double, dimensions::L, error_checking::enabled>, shape<3>, strides::column_major<shape<3>>, error_checking::disabled> t2;