Tensors and the tape

This chapter describes the values Mica computes with as mathematics: vectors, matrices and tensors of any rank, whose shapes are part of their types or arrive with their values, the operations between them, the gradient tape that records a computation and walks it backwards to its gradients, and the placement of a tensor on the GPU. Each construct names its rule in the grammar, says what it means, and quotes what the compiler answers when a source breaks the rule. Vectors and matrices, dynamic tensors, tracked values and the first GPU program teach the same ground from the start.

Every example is a complete program the website’s gate compiles with the compiler as it is written here. The output shown beneath an example is the output the program prints; the GPU example needs a card to run, and the chapter says what it computes.

The tensor types

A tensor’s type states its rank and its element type, and a static tensor’s type states its extents as well (tensor-type, dynamic-tensor-type):

A tensor is indexed from zero on every axis, whatever its type, and a constant index outside an extent the type states is refused where it is written, while a variable index outside it stops the program with index_out_of_range. A dynamic tensor starts empty: Resize(address grid, 2, 3) gives a matrix two rows and three columns, every element zero, and Rows, Columns and Length read its extents. A tensor is a value, as a record and an array are: assignment copies the whole of it, and a change to the copy leaves the original as it was. A literal writes a vector whole, v := (1, 2, 3).

Operations

The operators between tensors are the ones of linear algebra:

An operation computes in one element type. An integer literal takes the element type of the tensor it meets, and every other scalar has the element type already: a constant declared without a type is an int32, so a constant that scales a vector of int64 is declared Factor : int64 = 3. Between two vectors, * and / mean nothing: the product of two vectors is Dot, and no operator divides by a matrix. Sum, Transpose, Dot and the rest of the family arrive by imp from the math unit, and the math unit lists them.

Where a static tensor’s extents disagree, the compiler refuses the operation. Where a dynamic tensor’s do, nothing is known until the program runs, so every operation on a dynamic tensor compares the extents first, on every tier, and a disagreement stops the program at the operation with tensor_shape_mismatch: the product of two matrices of two rows and three columns, c := a * b, stops there, at its line and column.

{
    Shapes - vectors and matrices whose shapes are part of their types, and one whose extents arrive at run time.

    'v' and 'u' are vectors of three integers and 'a' a matrix of two rows and three columns: '+' and '⊙' combine two
    vectors cell by cell, a scalar scales one, '*' is the product of a matrix and a vector or of two matrices, and 'Dot'
    is the product of two vectors. 'grid' and 'other' are dynamic matrices, whose product's shape is read from them.
}
program Shapes;

imp
    WriteLn, Resize, Rows, Columns : std;
    Dot, Transpose : math;

const
    RowCount = 2;
    ColumnCount = 3;
    Factor : int64 = 3;
    Half = 0.5;
    VectorLine = "%ls %lld %lld %lld";
    PairLine = "%ls %lld %lld";
    DotLine = "Dot(v, v)      %lld";
    GramLine = "a * Transpose(a)  %lld %lld / %lld %lld";
    ShapeLine = "grid * other is %lld x %lld, its first cell %lf";
    VectorTitle = "v             ";
    SumTitle = "v + v         ";
    ScaledTitle = "3 * (v + v)   ";
    CellTitle = "u ⊙ v         ";
    ProductTitle = "a * v         ";

type
    { A vector of three integers, a vector of two, and the matrices the products answer. }
    Triple = vector[3] of int64;
    Pair = vector[2] of int64;
    Wide = matrix[2, 3] of int64;
    Square = matrix[2, 2] of int64;

var
    v, u : Triple;
    w : Pair;
    a : Wide;
    gram : Square;
    grid, other, product : matrix of float64;
    row, column : int64;

begin
    v := (1, 2, 3);
    for row := 0 to RowCount - 1 do
        for column := 0 to ColumnCount - 1 do
            a[row, column] := row * ColumnCount + column + 1;
    WriteLn(VectorLine, VectorTitle, v[0], v[1], v[2]);
    u := v + v;
    WriteLn(VectorLine, SumTitle, u[0], u[1], u[2]);
    u := Factor * u;
    WriteLn(VectorLine, ScaledTitle, u[0], u[1], u[2]);
    u := u ⊙ v;
    WriteLn(VectorLine, CellTitle, u[0], u[1], u[2]);
    w := a * v;
    WriteLn(PairLine, ProductTitle, w[0], w[1]);
    WriteLn(DotLine, Dot(v, v));
    gram := a * Transpose(a);
    WriteLn(GramLine, gram[0, 0], gram[0, 1], gram[1, 0], gram[1, 1]);
    Resize(address grid, RowCount, ColumnCount);
    Resize(address other, ColumnCount, RowCount);
    for row := 0 to RowCount - 1 do
        for column := 0 to ColumnCount - 1 do
        begin
            grid[row, column] := Half;
            other[column, row] := Half;
        end;
    product := grid * other;
    WriteLn(ShapeLine, Rows(product), Columns(product), product[0, 0]);
end.

It prints:

v              1 2 3
v + v          2 4 6
3 * (v + v)    6 12 18
u ⊙ v          6 24 54
a * v          14 32
Dot(v, v)      14
a * Transpose(a)  14 32 / 32 77
grid * other is 2 x 2, its first cell 0.750000

The tape

A tracked value is a value paired with its place on the program’s gradient tape: tracked before a type makes a variable, a parameter, a field or a result of the pair (value-type). The type behind it is a float, float16, bfloat16, float32 or float64, or a tensor of any rank over one. Ordinary code reads and writes the value half alone: an assignment of a plain value writes it, and an assignment of the tracked variable to a plain one reads it.

A tape window, tape and its statements and end, records every operation on tracked values that runs inside it (tape-window): the operators +, -, * and /, negation, and the recorded verbs of the math unit, which tracked values lists. A window starts a new record, so a value recorded in an earlier window has no gradient in this one. A routine whose parameters or result are tracked records its body when a window calls it, and computes plainly anywhere else.

After the window, Backward walks the record backwards from one scalar, the loss, once, and Gradient then reads the derivative of the loss with respect to a tracked value: a scalar for a scalar and a tensor of the same shape for a tensor. A tracked value the window never recorded has a gradient of zero. Both come from the language’s own unit, imp Backward, Gradient : std;, and the order is checked while the program runs:

The tape records one window at a time, in one task, so a window stands in the body of a routine or of the program, outside every other window and every concurrent block, and a task takes no tracked parameter.

{
    Descent - a weight walked towards its target by the slope the tape computes.

    Each round records the squared error of 'weight' in a tape window, walks the window backwards from the loss, and
    reads the slope of the loss with respect to the weight, which moves the weight a quarter of the slope downhill. The
    last lines record the length of a vector squared, whose gradient is twice the vector. Every value is a sum of
    halves and quarters, so every digit is exact.
}
program Descent;

imp
    WriteLn, Backward, Gradient : std;
    Dot : math;

const
    Rounds = 4;
    Rate = 0.25;
    Target = 3.0;
    Start = 1.0;
    RoundLine = "round %lld: weight %lf, loss %lf, slope %lf";
    VectorLine = "the gradient of Dot(p, p) at (%lf, %lf) is (%lf, %lf)";

var
    weight, loss, length : tracked float64;
    point : tracked vector[2] of float64 := (2.0, 3.0);
    plainWeight, plainLoss, slope : float64;
    plainPoint, slopes : vector[2] of float64;
    round : int64;

{ The squared distance of an estimate from the target, recorded on the tape when a window calls it. }
function SquaredError(estimate : tracked float64) : tracked float64;
begin
    SquaredError := (estimate - Target) * (estimate - Target);
end;

begin
    weight := Start;
    for round := 1 to Rounds do
    begin
        tape
            loss := SquaredError(weight);
        end;
        Backward(loss);
        slope := Gradient(weight);
        plainWeight := weight;
        plainLoss := loss;
        WriteLn(RoundLine, round, plainWeight, plainLoss, slope);
        weight := plainWeight - Rate * slope;
    end;
    tape
        length := Dot(point, point);
    end;
    Backward(length);
    slopes := Gradient(point);
    plainPoint := point;
    WriteLn(VectorLine, plainPoint[0], plainPoint[1], slopes[0], slopes[1]);
end.

It prints:

round 1: weight 1.000000, loss 4.000000, slope -4.000000
round 2: weight 2.000000, loss 1.000000, slope -2.000000
round 3: weight 2.500000, loss 0.250000, slope -1.000000
round 4: weight 2.750000, loss 0.062500, slope -0.500000
the gradient of Dot(p, p) at (2.000000, 3.000000) is (4.000000, 6.000000)

Tensors on the GPU

on gpu behind a tensor’s type places its storage in the memory of an NVIDIA GPU (tensor-type), and behind a float’s name places a single value there (device-scalar-type). The placement moves the storage and nothing else: matrix[2, 2] of float32 on gpu has the shape, the operators and the refusals of the matrix without it, and an operation whose operands are on the device runs on the device. ToDevice copies a value from the processor’s memory to the card’s and ToHost copies it back, and they are the only roads between the two: an operation or an assignment that would mix the two memories is refused. A device tensor holds float32, float64, float16 or bfloat16, and the tape records on the device as it does on the processor.

A program that places a tensor on the GPU builds on a machine without a card, and running it there names what it looked for. Running it needs an NVIDIA GPU of compute capability 7.5 or newer, NVIDIA’s driver, cuBLAS and NVRTC, as the first GPU program says, and the gpu runtime flag chooses the card on a machine with several (runtime flags). A device tensor does not cross into a routine’s parameter, and a program that places a tensor on the GPU links in the single-task flavor only, which is a known limitation.

A debug build waits for every kernel and every product the card runs, so a failing one is named where it failed. An optimizing build queues them and lets the program go on, and a copy back to the processor waits for their results; a failure is then reported at the next act that waits, with a note that says so.

{
    OnTheCard - one sentence of matrix arithmetic, computed on the processor and on the GPU.

    'a' and 'b' live in the processor's memory and 'da' and 'db' in the card's. 'ToDevice' and 'ToHost' are the only
    roads between the two, and the same expression over the device's operands runs on the device. Every input is a
    multiple of one eighth, so the two lines agree to the last digit.
}
program OnTheCard;

imp
    WriteLn, ToDevice, ToHost : std;

const
    Size = 2;
    Eighth : float32 = 0.125;
    Rate : float32 = 0.5;
    HostLine = "host   = %f %f %f %f";
    DeviceLine = "device = %f %f %f %f";

var
    a, b, c : matrix[2, 2] of float32;
    da, db, dc : matrix[2, 2] of float32 on gpu;
    row, column : int64;

begin
    for row := 0 to Size - 1 do
        for column := 0 to Size - 1 do
        begin
            a[row, column] := ((row * Size + column + 1) as float32) * Eighth;
            b[row, column] := ((Size - row - column) as float32) * Eighth;
        end;
    c := (a * b) * Rate + a;
    WriteLn(HostLine, c[0, 0], c[0, 1], c[1, 0], c[1, 1]);
    da := ToDevice(a);
    db := ToDevice(b);
    dc := (da * db) * Rate + da;
    c := ToHost(dc);
    WriteLn(DeviceLine, c[0, 0], c[0, 1], c[1, 0], c[1, 1]);
end.

The program computes one sentence twice, on the processor and on the card, and prints each result on a line of its own. Every input is a multiple of one eighth, so every product and every sum in it is exact.

The floating-point flags

Two words join an --optimize set of a program that computes on the GPU, and each changes the digits its products give, so neither belongs to a level:

A product’s digits are the program’s decision for every library it links, so a library’s build refuses both words: The optimization 'fp_tensor_float32' requires a program as output binary: linking aborted. The command line describes each optimization word.