A Geometrical Approach to Anisotropic Covariance Synthesis. R. James Purser

1. A Geometrical Approach to Anisotropic Covariance Synthesis.R. James Purser S.A.I.C. National Centers for Environmental Prediction, Washington D.C. Tuesday October 19th, 2004

2. The Focus of this Talk: To describe and discuss some fairly specific geometrical features of the algorithms we are developing at NCEP in order to be able to implement, in an efficient way, geographically adaptive statistical assimilation techniques for our operational atmospheric (and oceanic?) models. Acknowledgment: Much of the work described here benefitted from the collaboration with, and feedback from, Dave Parrish and Wan-Shu Wu. Also, the continuing support and encouragement of Geoff DiMego, Steve Lord and John Derber is greatly appreciated.

3. Topics to be discussed • Analysis, Covariances and Filters. • “Aspect” tensors and their geometry. • Basic “Triad” and “Hexad” algorithms. • Inhomogeneous statistics and the snags they cause. • Galois Fields and “Chromatic” filter sequences. • Exploiting geometry to improve triad and hexad methods. • Some applications of symmetry in the “Blended Hexad” • Prospects for extension to four dimensions. • Conclusion.

4. Statistical Objective Analysis • Construction of an increment field involves the convolution of data forcing terms with a spatial covariance. • Many repetitions of this operation occur in the iterative solution algorithm. • There is a pressing need to accommodate anisotropic, inhomogeneous covariances.

5. Properties of a covariance: • Positive definite • Self-adjoint (symmetric) • Smooth • Computationally efficient to apply • Locally isotropic • Spatially homogeneous • Gaussian?

6. A Gaussian isotropic homogeneous covariance is efficient, smooth, symmetric and self-adjoint. But VERY limited! Several methods that approximate such a covariance are available. Direct evaluation by rational functions (Gaspari and Cohn) Simulated diffusion (Derber and Rosati) Cartesian product of one-dimensional (recursive?) filters (Purser, Wu and Parrish) Other methods: Multigrid; Spectral, etc.

7. The spectral methods are very efficient when the covariance is homogeneous, and allow precise control over the shape of the covariance. But strictly limited to homogeneousfunctions. The recursive filter is an efficient line-smoother. The only shape it makes sense to simulate is the Gaussian because it is only this special shape that, by sequential transverse applications, fills out a three-dimensional form free of the imprint of the underlying grid orientation – contours are ellipsoidal. ANY Gaussian in D dimensions can be synthesized by D line smoothers at (in general) carefully selected oblique orientations. However, we work on a grid where the set of practical smoothing orientations is limited. Can we still synthesize any Gaussian by operations confined to the lines of the given grid?

8. Triads and Hexads The affirmative answer to this question led to the “Triad” filtering algorithm in two dimensions and the “Hexad” algorithm in three dimensions. In general we claim that, in D dimensions, any homogeneous Gaussian may be synthesized on a uniform grid by appropriate Gaussian line smoothers applied sequentially along (D*(D+1))/2 generalized (possibly oblique) grid directions. Moreover, it is possible to show that in D=2 or D=3 (but not for D > 3) dimensions, the combination of orientations and line smoothing “spreads” is UNIQUE in these minimal algorithms.

9. We adopt the principle that adequate approximations to the SHAPE of the covariance may be obtained by the superposition of a small number of quasi-Gaussian components (with different scales), possibly augmented by the operation of Laplacians. The basic unit, or building block, is then the Gaussian, and the effort to achieve computational efficiency is focused on the construction of the most efficient ways to generate adequate approximations to these Gaussian filters. We

14. Geodesics in a slice through the space of two-dimensional aspect tensors using the “natural” Riemann metric (which provides a “distance” invariant to arbitrary linear transformations of the physical space). The domain of proper aspect tensors corresponds to the cone whose boundaries are shown here by the two heavy diagonal lines.

15. The Basic Triad Method The crucial property of aspect tensors that we exploit in the filtering algorithms for both 2 and 3 dimensions is: ASPECT TENSORS OF HOMOGENEOUS SMOOTHING FILTERS ADD. Note: The aspect tensor of any line filter is of rank ONE. Therefore, given any “triad” of generalized grid lines, a given aspect tensor resolves uniquely, by “linear projection”, into three components associated with these three grid lines. BUT ONLY ONE MINIMAL TRIAD GIVES +ve NUMBERS

16. Minimal Triads In a minimal triad, the determinant of the matrix formed by any two of the lines’ generators (shortest grid steps) has a determinant of +1 or -1. The aspect vectors of a unit-spread filter in each of these directions, in units measured by the generators, form a BASIS for any aspect vector.

17. Geometry of Minimal Triads Each potential triad basis has its aspect vector on the surface of the aspect cone The triads of such vectors subdivide the interior of the aspect cone in a sort of (triangular) “honeycomb”. The interior of each triad’s region of aspect space are the vectors expressible as positive linear combinations of the 3 basis vectors.

18. Remarks It is possible to make a physical 3-dimensional model of the convex shell formed by these basis vectors on the aspect cone but almost all the facets away from the apex become highly elongated. An easier 3-dimensional depiction, equally informative, is obtained by taking the Legendre-Fenchel DUAL. (E.g., See R. T. Rockafellar, Convex Analysis)

19. Remarks (continued) In any three dimensional model of the basis-triads, different triads seem to have different shapes and sizes. However, when projected in line with the cone’s apex onto the unit-hyperboloid (of aspect tensors that have unit-determinant), and given the metric appropriate to a hyperbolic geometry, ALL TRIADS BECOME EQUAL!

20. A Poincare representation of the hyperbolic “plane” of normalized aspect vectors and the honeycomb of triads by which they are sorted

21. The Basic Triad Algorithm 1. Take any triad 2. Resolve the aspect tensor into this triad’s spread components. 3. If all three spreads are nonnegative, STOP. 4. If one spread is negative, exchange that generator by the unique alternative that makes up a new valid triad, and repeat 2, etc.

22. The Basic Hexad Algorithm The analog, in 3 physical dimensions, of the Triad algorithm in 2 physical dimensions is the “Hexad” algorithm. Six lines are configured as the diameters of a minimal skewed grid-imbedded cuboctahedron (a figure with 8 triangles alternating with 6 parallelograms). Alternatively, by Legendre-Fenchel duality, the six orientations may be characterized by the normals of a minimal skewed (dual-)grid-embedded rhombic dodecahedron (12-hedron).

23. Hexad transition rules • Remove the offending vertex-pair of the cuboctahedron. • Create the replacement pair by extending outwards by *2 the twin centers of the old quadrilaterals that did not contain the old disqualified vertices. • A new cuboctahedron is formed!

24. Snags when triads/hexads vary across the domain 1. The set of orientations proliferates – how do we order them to ensure that one line’s filtering can never derange the intermediate computations of another that crosses at a shared grid point? 2. Dislocations occur along triad/hexad boundaries because the “spread” coefficients drop to, or rise from, zero too fast at such junctions.

25. Chromatization via Galois Fields Every generalized grid line may be systematically assigned one of a finite palette of “colors” according to the non-null elements of an appropriate Finite, or “Galois” Field. Each triad or hexad never contains the same color twice. (I.e., lines of the same color never cross in the analysis grid). Perform the filter sequence by color and spread the work-load without fear of data conflict. Simplest chromatic triad scheme has 3 colors [GF(4)]; simplest chromatic hexad scheme has 7 colors [GF(8)].

31. Gaussian covariances may be synthesized in two and three dimensions by the “triad” And “hexad” procedures – sequential applications of 1-D gaussian filters of carefully chosen widths along carefully chosen generalized sets of 3 or 6 grid lines. The “chromatic” Triad and hexad procedures assign “colors” to each line direction so that each triad Comprises three distinctly colored direction and each hexad comprises six distinctly Colored directions. Color coding allows line filtering to proceed sequentially and without Conflicts, even as triads/hexad change from region to region. 3-color triads (Poincare representation) 4-color triads (Poincare representation)

32. A solution to the second difficulty – the dislocations, involves mixing of “blending” the aspect vectors in a symmetrical distribution (“kernel”) centered on the target aspect point. Linearity of aspect tensors implies that, by a properly weighted attribution to the triads/hexads overlapped by the kernel, and associating to each involved triad/hexad the overlapping portion’s centroid as that triad or hexad’s effective target, the correct aspect tensor emerges from the filter sequence that this prescription dictates. For example, in the triad case, we may choose a kernel confined to a disc small enough at each point in a projection of aspect space never to overlap more than two triads. In each triad it overlaps, the centroid of the overlap portion may be found. The basic triad algorithm is run with both these secondary targets. Their line-spreads are aggregated.

33. In the triad case, a more direct method of smoothing the triad-basis polyhedral “shell” does essentially the same job. The geometry is sufficiently simple to allow the construction of such smooth “bridging functions”. Each target aspect tensor is centrally-projected onto the standardized smoothed “shell”. If the projection puts the result on a bridging function straddling two triads, the four vertices involved share the spread-weighting according to their respective barycentric weights assigned to this projected aspect point. A final multiplier, common to the four vertices (=line generators in physical space), inflates the projected aspect vector back to its original target value. This is the 4-color blended triad method illustrated:

42. In order to remove transition noise from the hexad procedure (the 3-D case) we must replace the simplest of the chromatic algorithms (seven colors) by a 13-color sequence. However, the generic transition singularities are considerably more complicated, since hexad weights may vanish three at a time and still leave a valid (positive definite) aspect tensor. The singularities cannot be described in less than three dimensions (out of the six dimensions of a general aspect tensor). Once their geometry is understood, it is possible to carry out the appropriate averaging of the line-weights associated with hexads within some prescribed “distance” (in aspect space). These can be tabulated.

43. Constructing a “Blended Hexad” In this MUCH more complicated case, the first order of business is to identify and characterize the geometry of the interior hexad-junctions of the six-dimensional aspect cone. The geometry is far too complicated to admit any simple “bridging function”. We adopt an explicit kernel smoothing strategy in this case. Thus, we must learn how many hexads can be simultaneously conjoined at the worst possible generic interior junction (since even in infinitesimal kernel will overlap all hexads involved).

44. Interior junctions occur when any hexad spread component vanishes. But as many as three hexad components may simultaneously vanish without leading to a degenerate aspect tensor – provided the generators of the three nonvanishing components that remain are not mutually coplanar. In the triad case, interior junctions had dimensionality-2, or co-dimensionality-1 (a spherically symmetrical kernel intersecting the interface could be replaced by its marginal projection into its central “axis” oriented to intersect that interface normally – it becomes essentially a 1-dimensional problem! Likewise, for the hexad, the essential dimensionality of an interior junction is really only three – not six!

45. Generic interior junctions among hexads in the 6-D aspect cone A careful analysis of the geometrical situation reveals that the form of the singularity or “junction” of highest co-dimension (=3) is a configuration of lines and planes arrayed about the central singularity in the form roughly resembling a strangely modified “radar reflector” separating precisely 16 distinct hexad regions, grouped asymmetrically, {4+12}.

46. A model showing the geometric form taken by a generic 3-fold degeneracy of the hexad decomposition projected from the six-dimensional space of aspect tensors to the three non-trivial dimensions that characterize the singularity. Each of the 16 corner regions correspond to a distinct hexad, but three of the contributing directions (1,2,3) are common to all of them. Another ten directions (4,..,13) become active arbitrarily close to the central singularity, each being invoked in any of the three or six corner regions that surround a particular one of the ten “rays” or edges extending out from the center.

47. Mercator representation of the generic transition singularity associated with three-fold degeneracies of the hexad decomposition of aspect tensors of anisotropic covariances in three dimensional lattices. Each hexad associated with the singularity occupies one of the 16 triangular-corner regions and comprises directions 1, 2, 3, together with the three numbered directions indicated by the labels at the relevant triangle’s three vertices.

48. Models representing the 16 hexads surrounding a 3-fold degeneracy. Six hexad directions are indicated by the vectors normal to the faces of these linearly-transformed rhombic-dodecahedra. (The “dual” representation.) The 13 directions, associated in this model with faces, edges and corners of a cube oriented with the lattice, are color-coded.

50. A close-up, showing in geometric form how a transition from one hexad to a neighbor results from the substitution of one of the six member directions with another. The left hexad, (1,2,3,6,8,12) loses “gold” (12) and gains “red” (7) to morph into the right hexad (1,2,3,6,7,8) here. (The fly appears to prefer the latter.)