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


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

Get the plugin now

View by Category
About This Presentation

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


'Aspect' tensors and their geometry. Basic 'Triad' and 'Hexad' algorithms. ... Exploiting geometry to improve triad and hexad methods. ... – PowerPoint PPT presentation

Number of Views:68
Avg rating:3.0/5.0
Slides: 66
Provided by: wd287


Write a Comment
User Comments (0)
Transcript and Presenter's Notes

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

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
  • Galois Fields and Chromatic filter sequences.
  • Exploiting geometry to improve triad and hexad
  • Some applications of symmetry in the Blended
  • 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
  • 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
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.
(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
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
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
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.
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
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
  • Create the replacement pair by extending outwards
    by 2 the twin centers of the old quadrilaterals
    that did not contain the old disqualified
  • 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
(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
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
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.
  • 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!
About PowerShow.com