# Presentation: 376175

## 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

## 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.

## 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.

## 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.

## Properties of a covariance:

- Positive definite
- Self-adjoint (symmetric)
- Smooth
- Computationally efficient to apply
- Locally isotropic
- Spatially homogeneous
- Gaussian?

## 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.

## 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 _**homogeneous**_** **functions.
- 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?**

## 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.

## We 

- 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.

## Slide 10

## Slide 11

## Slide 12

## Slide 13

## 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.

## 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**

## 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.**

## 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.**

## 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***)

## 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!**

## A Poincare representation of the hyperbolic “plane” of normalized

- aspect vectors and the honeycomb of triads by which they are sorted

## 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.**

## 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).**

## 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!

## 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.

## 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)].**

## Slide 26

## Slide 27

## Slide 28

## Slide 29

## Slide 30

## 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)

## 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.

## 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:

## Slide 34

## Slide 35

## Slide 36

## Slide 37

## Slide 38

## Slide 39

## Slide 40

## Slide 41

## 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.

## 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).**

## 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!

## 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}.**

## 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.

## 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. 

- 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.

## Models representing the 16 hexads surrounding a 3-fold degeneracy. Six hexad directions are 

- 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.

## Slide 49

## 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.)

## Not surprisingly, or coincidentally, the 16 conjoined hexads

- around a co-dimension-3 junction share a set of 13 grid
- orientations that exhibit all 13 “colors” of the Galois Field,
- GF(27).
- Recall that GF(27) involved a repeating 3*3*3 pattern in
- the grid of physical orientations. Each “cluster” of 16 hexads
- of a given junction have their 13 different grid generators
- (together with their negatives) occupying the exposed points
- of a particular 3*3*3 parallelepiped within the larger grid
- of all possible generators.

## The “Star-Burst”

- We need to know: how large a kernel can we have at a given
- point before we “contaminate” our local cluster’s 13
- orientations with additional contributions found beyond
- the** next **junction?

- The “Star-Burst” configuration helps to formulate the answer.
- It depicts, in the most symmetrical way, the 16 interlocking
- “clusters” that contain the same hexad (in the center), as an
- arrangement of generators (and their negatives), represented
- as physical grid points. The parallelepipeds denoting clusters
- are found to fall into two categories {4 + 12}.

## Any given 3-fold degeneracy and its surrounding “cluster” of 16 hexads can be associated with a particular 3*3*3 parallelepiped of generalized lattice directions. But each given hexad belongs to 16 such clusters. In this geometric model, it is the vertices that associate with grid line directions (“primal” form) and the lattice here is of the “close packing” (fcc) type. This makes the union of the 16 cluster-parallelepipeds take on the most symmetrical possible configuration.

- The same 13 “colors” , in various arrangements, suffice to properly color each cluster (no hexad has the same color twice; no cluster lacks a color.)

## Every target aspect vector belongs to a unique hexad,  so we may take that hexad to be the center of the star-burst of the 16 relevant clusters. 

- To proceed, we must endow the affine space of aspect vectors with a suitable Euclidean metric. A good choice is to impose throughout this 6-D space the metric that, at the target aspect point, matches the intrinsic Riemannian metric we described earlier.
- But at the same time, it is convenient to simplify the geometry of hexad junctions by confining attention to the 5-D “tangent plane” which, at the target point, is tangent to the “surface” of constant aspect tensor determinant (and orthogonal to the ray that points to the apex of the aspect
- cone). Every hexad intersects this tangent plane in the shape of a “simplex”***** and the adoption of our Euclidean “tangent metric”  incurs only minimal distortion of the ideal Riemannian distance in the relevant immediate vicinity of our target point.
*** A “simplex” in d dimensions is the convex hull of d+1 non-coplanar**
**points.**

## The smoothing kernel at this target point can now be represented

- by an isotropic density inside a 5-D sphere centered on the target.
- The question: How large can this sphere be? now has the
- following well-defined answer in terms of the 5-D Euclidean
- geometry of the representation of the neighboring hexads.
- For each of the 16 clusters belonging to the star-burst centered
- on the hexad containing our target, we record the** MINIMUM **
- Euclidean distance to the boundary of this cluster. (This
- 4-D polyhedral boundary consists of 48 simplex panels.)
- Then we find the **MAXIMUM** of our 16 minima and note the
- particular cluster to which it belongs.
- This distance dictates how large a radius for our smoothing
- kernel we are allowed; the winning cluster dictates the 13
- smoothing directions we must use at this aspect tensor.

## In practice, we don’t want to have to do this calculation

- from scratch for every point in the analysis grid.
- Instead, we should tabulate results of this and other
- geometrical calculations in a regular table spanning
- a representative portion of the 5-dimensional section
- through a representative hexad region that contains all
- six of the hexad’s aspect space basis vectors.
- As previously noted, the 5-dimensional image of a hexad
- in any 5-D planar section has the shape of a SIMPLEX.
- Five dimensions is a BIG space. However, the hexad
- simplex possesses a symmetry group of order 24
- (from the octahedral group symmetry of the cuboctahedron
- of hexad generators). We can exploit this symmetry.

## Simplex Tables

- A simplex, in any number of dimensions, splits up very
- naturally into a uniform grid of mini-simplexes.
- Linear interpolation from a simplex tabulation to an
- interior target point needs d+1 source points (forming
- a mini-simplex surrounding the target).
- Kernel-radius limits need to be smoothed – this can
- also be done on the regular grid of a simplex table. A
- smoother may be constructed as a “diffusion” operation
- using a Laplacian based on the Riemannian metric.
- Only roughly 1/24 of the simplex need be tabulated!

## The form of the blending kernel

- The most natural family of radial density profiles having
- the “compact support” dictated by a hard upper limit for
- the radius is the family of symmetric “Beta” distributions.
- For the conventional applications to statistics, with a
- positive pair of parameters, these are positive densities,
- concentrated more to the center as the parameters increase.
- In orthogonal polynomial theory, the same functions form
- the positive density distributions that give rise to the
- ultraspherical (“Gegenbauer”) family of orthogonal
- polynomials (including Chebyshev and Legendre).

## Why are these distributions “natural”?

- Because, in special cases, they reduce to polynomials and,
- as “spherical” distributions in d dimensions, their “marginal”
- densities by parallel projection into d’<d dimensions gives
- another density distribution that belongs to the same family.
- Recurrence relations enable us to extend the dimensionality
- in** both **directions, including instances where the radial
- profile needs to be expressed in terms of  “generalized”
- distributions (Dirac delta functions, etc.) and are not always
- positive.
- (Such profiles describe the forms of the simplest self-similar
- expanding “shocks” of the wave equation in d dimensions.)

## Slide 60

## In order to match the blending that works so well for the

**triad algorithm, we may choose the Beta distribution profile**
**whose marginalization into the  three co-dimensions of**
**hexad junctions is the fortuitously simple spherical shell!**
**Thus, the kernel integrations for each point of our simplex table**
**are no longer six-dimensional, five-dimensional, or even**
**three-dimensional – we have reduced the integral to **
**two-dimensional fragments (up to 16 of them) of ordinary**
**spherical surfaces. Moreover, these integrals of area and**
**first-moment (needed to find weight and centroid of each**
**fragment intersecting each hexad belonging to the appropriate**
**“****cluster”) can be further reduced, by application of Green’s**
**Theorem, to **_**one**_**-dimensional integrals on smooth arcs that**
**form the boundaries of the spherical surface fragments!**
**Gauss-Legendre quadrature yields **_**extremely**_** ****accurate results.**

## More about Simplex Tables

- The simplex table covers the region occupied by 1/24 of one
- repesentative hexad. Each data entry supplies the 13 active
- generators for this particular tabulation point, together with
- the 13 corresponding normalized (proportionate) spread
- weights, all listed in the proper chromatic order. Symmetry
- principles (i.e., applied group theory) dictate how to map
- any aspect tensor into the table and extract the needed
- quantities.
- In order to interpolate linearly in a simplex table of d=5
- dimensions, the table’s grid must be diced and sliced. Having
- diced, there remain (d!)/2 = 60 ways to slice! The different
- ways fall into 60=24+12+12+8+4 equivalence classes with
- respect to operations of the octahedral group of hexad
- symmetries.

## Tables (continued)

- We would like to guarantee that, regardless of the “path”
- taken by the Hexad Algorithm to arrive at the proper hexad,
- the target is interpolated from the tabulated data in precisely
- the **same **way (the same 6 surrounding mini-simplex points
- of the table). An analysis involving the geometry, the
- symmetry, and with help from the Galois coloring, GF(8),
- it has been possible to formulate the new Hexad Algorithm
- and choose a dicing-slicing option that fully satisfies this
- requirement and that, moreover, simplifies the color assignment
- procedure in **ANY** chromatization scheme!

## Four dimensions?

- The Hexad method was discovered by dumb luck – an
- approach that unfortunately failed to find the Four-
- Dimensional equivalent!
- A more systematic approach is to map the aspect basis
- vectors of all the 4-D grid generators into 10-D aspect
- space and examine the geometry of the surface “shell”
- of the convex-hull of these aspect vectors.
- The shell comprises two kinds of “facet”. One kind is
- a 10-pointed simplex and the other is a 12-pointed
- object. Subdivision of the 12-pointed object into 48+12=60
- simplexes of two additional types supplies the necessary
- “trick” for a “Decad Algorithm”.

## Conclusion

**Triad and Hexad algorithms allow any Gaussian covariance components to be efficiently synthesized by successive line smoothers.**
**Galois Fields provide conflict-free sequencing.**
**Defects with inhomogeneous covariances synthesized by the basic 3-color Triad or 7-color Hexad method are eliminated (at the cost of more line operators) by adopting the 4-color Blended Triad or 13-color Blended Hexad methods.**
**Abstract algebraic methods have some surprising practical applications!**