a geometrical approach to anisotropic covariance synthesis r james purser l.
Skip this Video
Download Presentation
A Geometrical Approach to Anisotropic Covariance Synthesis. R. James Purser

Loading in 2 Seconds...

play fullscreen
1 / 65

A Geometrical Approach to Anisotropic Covariance Synthesis. R. James Purser - PowerPoint PPT Presentation

  • Uploaded on

A Geometrical Approach to Anisotropic Covariance Synthesis. R. James Purser. S.A.I.C. National Centers for Environmental Prediction, Washington D.C. Tuesday October 19 th , 2004. The Focus of this Talk:. To describe and discuss some fairly specific geometrical

I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.
Download Presentation

PowerPoint Slideshow about 'A Geometrical Approach to Anisotropic Covariance Synthesis. R. James Purser' - arnav

Download Now An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript
a geometrical approach to anisotropic covariance synthesis r james purser

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


National Centers for Environmental Prediction, Washington D.C.

Tuesday October 19th, 2004

the focus of this talk
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
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
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
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 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?

triads and hexads
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 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




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 Basic Triad Method

The crucial property of aspect tensors that we exploit in

the filtering algorithms for both 2 and 3 dimensions is:



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.


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


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


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



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
The Basic Triad Algorithm

1. Take any triad

2. Resolve the aspect tensor into this triad’s spread


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


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:


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


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

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

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.


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,


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



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


(Such profiles describe the forms of the simplest self-similar

expanding “shocks” of the wave equation in d dimensions.)


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

more about simplex tables
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


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


tables continued
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
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”.

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