Loading...

PPT – A Geometrical Approach to Anisotropic Covariance Synthesis. R. James Purser PowerPoint presentation | free to download - id: 4679c-ZDc1Z

The Adobe Flash plugin is needed to view this content

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(D1))/2 generalized (possibly oblique)

grid directions. Moreover, it is possible to

show that in D2 or D3 (but not for D gt 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 filters.

We

(No Transcript)

(No Transcript)

(No Transcript)

(No Transcript)

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 triads

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 cones 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 triads 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 lines

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

(No Transcript)

(No Transcript)

(No Transcript)

(No Transcript)

(No Transcript)

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 portions centroid as that triad or

hexads 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

(No Transcript)

(No Transcript)

(No Transcript)

(No Transcript)

(No Transcript)

(No Transcript)

(No Transcript)

(No Transcript)

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 associat

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

s 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! L

ikewise, 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, 412.

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

(No Transcript)

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

ll that GF(27) involved a repeating 333 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 333 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 clusters 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 333 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 d1 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 dont 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 hexads 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 d1 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 dltd 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. (Su

ch profiles describe the forms of the simplest

self-similar expanding shocks of the wave

equation in d dimensions.)

(No Transcript)

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 Greens 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 d5 dimensions, the tables grid must be diced

and sliced. Having diced, there remain (d!)/2

60 ways to slice! The different ways fall into

6024121284 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 481260 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!