Friday, April 12, 2019

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:

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

You can exempt an index from the repetition count with braces, so a(i)bi 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)bici, which would be a quadratic form with diagonal matrix A as kernel.

I use a special notation Ij or Ijk 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 apq, multiply by bq, and put the result in yp. 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 Ts of measured temperatures, or anomalies, and so the result can be represented as a weighted sum, wsTs. 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 Ingsxs. The sum has to be divided by the count, which in matrix terms is the product IngsIs. The result then has to be multiplied by the grid areas ag and added. Overall

Integral over earth = ag (ItIn(g)t)-1 Ings xs

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 xsa = areas of triangles
Loess order 0 Ip (ItW(p)t)-1 Wps xsW 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.


  1. Off-topic: I see you engaged with Dr. Spencer over at WUWT (e.g., Will Upload Whatever Trash) on his very simple CO2 model... it is amazing how threatened the contrarian community is by a) the rise of CO2, b) the attribution of that rise to humans, and then c) the concept that a portion of our emissions are effectively permanent (for the timescales we care about).

    While (a) seems to have died down (fewer and fewer citations to Ernst Beck), and (b) is hopefully fighting a losing battle, it seems like (c) is still going strong.

    And you'd think that Spencer would have learned some humility after he used a simple model to justify questioning (b) 10 years ago ( and But, no, he's just updated to a slightly less outrageous theory.

    I like simple models. They are a good way to experiment with assumptions, and build some physical intuition. But if you are using your simple model to challenge the consensus, you should at least put a little bit of effort into understand what that consensus was built on. And it seems like Spencer, Curry, Salby, Essenhigh, etc. etc. assume that all that carbon cycle experts do is to see if they can fit observed atmospheric concentrations. They usually ignore (or it they don't ignore, they completely misunderstand) the fact that the carbon cycle modelers have ocean data, isotope data, geographically resolved data, and understanding of physical laws to constrain their work.

    But when it comes down to it, there's a real simple question for Spencer and his conclusion that "if CO2 emissions were stabilized and kept constant at 2018 levels, the atmospheric CO2 concentration would eventually stabilize at close to 500 ppm, even with continued emissions": where's the infinite sink, and what's its timescale? Carbon cycle experts know that conversion into carbonate from weathering and sedimentation can take up very large quantities of CO2... but it operates on a timescale of 10s of thousands of years. Does Spencer know this? If so, does he assume that scientists have their kinetics wrong? That they've missed another mechanism? That there are naiads who spirit the excess carbon away to sacrifice on the altar of their Lord Poseidon?

    Just some thoughts...


    1. I didn't see much in Spencer's model at all. It was just tracking concentration in terms of emissions, with a formula that gave sinks as proportional to the difference from a notional equilibrium. I prefer my explanation, which interprets that airborne fraction (the other part) in terms of exponential rise. But the end result isn't that much different. There are sinks, and they are more active when CO2 is higher. That isn't controversial.

    2. Another OT comment. Has the NCEP/NCAR data stopped updating? If so, is it on your end or theirs? The last date I see there is 8 April.

    3. My end. The general download program had stopped. I'm tinkering with a lot of stuff as a fnalise TempLS V4; should be out next week. Will fix.

      The last table on the page records when files come in.