Some math used in TempLS - index notation and sparse matrices.

Index notation for arrays, and the summation convention

Index notation for arrays became popular with the tensor calculus, partly because it also elegantly embraced the concepts of contravariance and covariance. It is sometimes called the Einstein notation. I'm not using the tensor aspect here, so no superscripts. Arrays, including vectors, are simply represented by writing down their coefficients with subscripts, with the understanding that those vary over a known range. There is no attempt to visualise them geometrically. The coefficients, like arrays, can be added if the subscripts match in range, and be multiplied by scalars or pother subscripted arrays.

But the core of the system is the summation convention, also named after Einstein. It expresses the idea of inner product, which is what distinguishes matrices from just vectors. If an index is repeated in a term, it implies summation over the range of that index. Some familiar cases:

aibj	is the outer product of two vectors, functioning as a 2-index array, or matrix, but
aibi	is the inner product, with a single number (scalar) result. i is referred to as a dummy variable, because it won't be referenced again.
aijbj	is a matrix right multiplied by a vector. The result is a vector, indexed with i
aijbi	would be left multiplied. But that terminology isn't needed; the indices say it all.
ajj	is the trace of matrix A
aijbjk	is the matrix product A*B, indices i and k.
aijbjkck	is the product A*B*c, a vector, index i

You can exempt an index from the repetition count with braces, so a_(i)b_i is not summed. The multiplication is real, so it is as if A was a diagonal matrix. It often appears as a kernel, as in a_(i)b_ic_i, which would be a quadratic form with diagonal matrix A as kernel.

I use a special notation I_j or I_jk for an array consisting entirely of 1's.

I have used this notation at least since introducing V2 of TempLS, but I'll be making more extensive use of it with V4. A feature is that once you have written down an index version of what you want to do, it maps directly onto R coding.

Sparse matrices

Often the indices have a small range, as over space dimensions. But sometimes the range is large. Most integration methods in TempLS involve mapping from one large set, like stations, to another, like grid cells. At some point there will be a matrix involving those two indices. In V4 we might have 30,000 stations, and 10,000 cells. You would not want to enter such a matrix into RAM.

The saving grace is that most nodes have no interaction with most cells. As with partial differential equations in maths, the relations are local. Most terms are zero, a sparse matrix. So it is possible to enter the non-zero terms, but you also have to say where they are in the matrix. I just list the N nonzero terms in a vector and in a Nx2 matrix, I list first the row and then the column of each term.

The key issue again is matrix-vector multiplication y=A*b. This is what can take a set of numbers over stations to one over grids. In a language like C, with this structure, the algorithm is simple. You go through that list, take each a_pq, multiply by b_q, and put the result in y_p. But that is a lot of operations, and in an interpreted language like R, they carry too much overhead. It would be possible to stretch out the terms of b to match a, so the multiplication could be between vectors. But they can't be added into the result vector in one operation, because several numbers would have to be added into one location. I use an algorithm that sorts the indices of a into subsets that do correspond to a unique location, so the multiplication can be reduced to a small number of vector operations.

Use in TempLS

I have written. for example here and very recently here about ways in which spatial integration is done in TempLS, and why. Integration is a linear operation on a set T_s of measured temperatures, or anomalies, and so the result can be represented as a weighted sum, w_sT_s. For each method, the first task of TempLS is to calculate those weights.

I'll go through the basic arithmetic of simple gridding from this viewpoint. To integrate directly, you work out which cells of the grid the stations belong to, calculated an average for each cell, and then multiply those by the areas of each cell and add. It's fairly simple to program, but I'll describe it in sparse matrix terms, because more complex cases have a similar structure.

The key matrix being formed is the incidence matrix In over indices g for grid (rows) and s for stations (columns). This is zero, except where station s is in cell g, when it has a 1. So the action of summing station data x in cells is the product In_gsx_s. The sum has to be divided by the count, which in matrix terms is the product In_gsI_s. The result then has to be multiplied by the grid areas a_g and added. Overall

Integral over earth = a_g (I_tIn_(g)t)^-1 In_gs x_s

The result has no indices, so all must be paired. Note the appearance of a bracketed (g) in the denominator. The repetition of separate station indices s and t indicates two separate summations.

Forming the index set for the sparse matrix is usually simple. Here the right column just numbers the stations, starting from 1. The left column just contains the cell number that that station is in.

That was reasoned starting from known x. To calculate the weights w, x must be removed, and the linear algebra started from the other end, with matrices transposed. Transposing a sparse matrix is trivial - just swap the columns of the index.

I'll go through the algebra for more complex cases in a future post, probably in conjunction with the code. But here is a table showing how the basic structure repeats. These are the various integration algorithms in index notation.

No weighting	(ItIt)-1Isxs
Simple grid	ag (ItIn(g)t)-1 Ings xs
Triangle mesh	ag Ings xs	a = areas of triangles
Loess order 0	Ip (ItW(p)t)-1 Wps xs	W is a function (with cutoff) of distance between station s and nearby node p.
Loess order 1	ypj (ztjztkW(p)t)-1 Wpszsk xs	zsk are the station of station s; ypj are the coordinates of node p.