Skip to content
Gijs Molenaar edited this page Feb 12, 2014 · 3 revisions

Tensor Math

A Result can contain a number of VellSets, which we used to call "planes". It has now been extended with an optional dims attribute, which is a vector of integer tensor dimensions. The constituent VellSets (the total number of which must correspond to the product of all dimensions) are then treated as a tensor value of the given shape, in row-major (i.e. "C" as opposed to "FORTRAN") order. Note that each tensor element is a VellSet, and thus may be represented by an n-dimensional array, indicating variability along axes such as time and freq, and may contain a number of spids and perturbed values. Each VellSet in the tensor may have its own axes of variability and perturbed values, completely independently of the other elements.

If no dims attribute is present, the Result is considered to be a scalar (given one VellSet) or a vector (with multiple VellSets).

Function-derived nodes (i.e. as all the math operations) are applied on an element-by-element basis. If any child results are tensors, their tensor dimensions must match. An operation on a tensor and a scalar causes the scalar to be applied to every element of the tensor.

An empty VellSet (i.e. no value attribute) is treated as a "missing" value. Functions propagate missing values into their output (i.e. any function with a "missing" argument produces a "missing" output), also on an element-by-element basis. Matrix multiplication, on the other hand, treats missing values as zeroes.

Tensor Operations

I don't want to breed ooddles of similar nodes when a few simple ones can accomplish the same thing. Various aliases can be implmeneted in the Python tree generator though. Towards this end the following basic kernel set is implemented:

Creating Tensors

MeqComposer (done) : extend MeqComposer with an optional dims field, which is simply a vector of tensor dimensions. For dims=[k,l,m], it will expect k*l*m planes in the child results, and produce a tensor result accordingly. Can also be used as a "reshaper" (i.e. feed in a 4-vector, get out a 2x2 matrix).

MeqFunction (done) : functions are applied vellset-by-vellset. Arguments should be either tensors of equal dimensions, or tensor/scalar combinations. The dims of a function result are inherited from its children.

MeqConstant (done) : allow vector/array constants, which return their value as a tensor.

Null/unity elements (done) : there should be a special check for null or unity elements, to speed up calculations. Nulls can also be possibly represented by empty Vells and empty VellSets. This optimization is now only implemented for addition and multiplication.

MeqSpigot/MeqSink (done) : operate with 2x2 coherency matrices (see Hamaker's paper IV btw), i.e., always generate a 2x2 matrix regardless of how many correlations there are in the measurement set. This is not wasteful since null matrix elements are ignored (see optimization above). Note that this is merely the default behaviour. The Spigot is flexible enough to produce a Result of any shape; it's init record tells it how to map columns of the measurement set to elements of the tensor. Thus it can be easily configured to produce scalar or vector results as needed. Likewise, the Sink is also configured with a mapping from tensor elements to MS columns.

MeqCondeq/MeqSolver (done) : condeqs should support tensor results on both sides with congruent tensors, or a scalar/tensor combination, to yield as many equations as there are elements in the tensor. Solver should deal with multiple equations from a single condeq. One thing to think about: perhaps a missing VellSet in a tensor should indicate missing data, and thus not yield any equations (even if it's not missing on the predict side?)

Element manipulation

MeqSelector (done) : can specify selection as index=[i,j,...] to extract a single VellSet. Any element of the selection can be set to <0 (<=0 in Glish) to select a "slice". E.g., index=[0,2] (in Glish) selects the second column of a matrix, returning a 1xn tensor. NB: If no index field is specified, a MeqSelector is a 'no-operation' node, i.e. it just passes on its chld(ren) with a minimum of overhead. This could be used to attach a step-child etc.

MeqPaster (done) : expects two children; "pastes" the second child at the specified position index=[i,j,...] of the first child. First child must obviously be a tensor. Just like Selector, slices may be specified by an index of <0. The tensor dimensions of the second child must of course match the shape of the selected slice.

Vector/matrix operations

MeqMatrixMultiply (done) : it's a completely different operation from MeqMultiply so it should be a different node class (regular MeqMultiply is element-by-element).

MeqTranspose (done) : transposes a matrix (or a vector). Set initrec.conj=T (default is false) for a conjugate transpose.

MeqMatrixInvert22 (done) : inversion of 2x2 matrices

MeqMatrixInvert : generalized matrix inversion (very powerful! Now who wants to write it? :-? )

MeqDirectProduct, MeqDotProduct, MeqCrossProduct : trivial in definition, but in practice requires a lot of housekeeping code. Defer implementation for now.


JEN's original picture

(This is for historical reference only, actual implementation follows the outline above).

Create/fill matrices

NB: At this moment, a child may be defined as a record or a string (node-name). It should be simple to extend this to numeric values. Internally, this would be treated as the equivalent of a MeqConstant, and that is what one would get out when extracting such an element.

Extract individual elements (or submatrices)

    • MeqSelector extract a specific element [i,j,k] * (or even sub-matrix [ii,jj]?)

Some special matrices for convenience (and consistency)

Vector/matrix operations

Mostly extensions of existing MeqMath node classes (unary and binary)::

    • MeqMultiply scalar/vector/matrix multiplication * depending on children

Clone this wiki locally