9.3.8 Selecting Realizations for Interpolation or Least Squares

9.3.8  Selecting Realizations for Interpolation or Least Squares

Realizations 1r[k] for interpolation or least squares should be selected with care. A poor choice may distort results. The simple approach of generating pseudorandom realizations based upon the distribution of 1R is inadvisable, as it will cause most realizations to cluster near 1|0μ 0E(1R). If a quadratic remapping is to be a good approximation over a large region of values for 1R, we need to apply interpolation or least squares to a dispersed set of realizations.

A general approach to specifying realizations is to let the point 1|0μ be one realization 1r[l], and select other realizations that are dispersed in some manner on the ellipsoid4


where 1|0Σ is the conditional covariance matrix of 1R and q is a constant. If 1R is joint-normal, this ellipsoid defines a level curve (surface) of its distribution—that is, a curve on which the probability density of 1R is constant. In practice, q is set equal to 1, 2 or some value in-between. In a sense, it reflects the number of standard deviations from 1|0μ at which realizations are placed, but this is precise only in one-dimension.

We may disperse l – 1 realizations on ellipsoid [9.30] by selecting points on the unit sphere centered at 0 and projecting these onto the ellipsoid. There are various ways this might be done. The following approach provides a good dispersion of realizations as long as 1|0Σ is not multicollinear.

Let 1|0σ be a diagonal matrix with diagonal elements equal to the conditional standard deviations of 1R. Let 1|0ρ be the conditional correlation matrix of 1R. Given points pk on the unit sphere centered at the origin 0, define corresponding realizations 1r[k] as


This leaves the question of how to select l – 1 points pk on the unit sphere centered at 0. We might randomly disperse them. Generate l – 1 pseudorandom vectors vk ~ Un((–1,1)n). Discard and generate replacements for any vectors that equal the zero vector 0 or have norm greater than 1. Then set


[for all k. This approach has little to recommend itself other than the fact that it is easy. Distortions may result from points randomly clustering in certain regions and not in others.

If a quadratic remapping is to be constructed with interpolation, we may directly select points pk based upon the coefficients ci, j , b i , and a to be determined. For each coefficient c i,i , select a point pk whose components are all 0’s except the i th component, which is 1. For each coefficient ci, j  for which i ≠ j, select a point pk whose components are all 0’s except the i th and j th components, which are . For each coefficient bi, select a point pk whose components are all 0’s except the i th component, which is –1. This procedure will yield precisely l – 1 points unless your quadratic form has a = 0. In that case, the approach will yield l points. To reduce this to l – 1 points, discard two points, pk and pj, and replace them with the single point


Obviously, discarded points pk and pj should be selected so that their replacement point is not the same as one of the other points already selected.

To construct a quadratic remapping of form


we would select points pk as indicated in Exhibit 9.18.

Exhibit 9.18: An arrangement of points pk on the unit sphere suitable for use in interpolation of a quadratic remapping of form [9.34].

Considering this configuration of points, we may wonder if better results might be obtained with a symmetrical configuration, such as that in Exhibit 9.19. Such symmetrical configurations tend to work best if a remapping is to be constructed using least squares and the number of points l exceeds the number m of coefficients to be selected by a reasonable margin.

Exhibit 9.19: A symmetrical alternative to the configuration of Exhibit 9.18.

Ignoring trivialities,5 a perfectly symmetrical configuration of l – 1 points on a unit sphere in n dimensions is not always well defined. In three dimensions, such an arrangement is possible with 4, 6, 8, 12, or 20 points. These symmetrical configurations are achieved by inscribing one of the five Platonic solids within the sphere and placing a point at each of the solid’s vertices.

Exhibit 9.20: The Platonic solids.

In three dimensions, five points cannot be distributed with such symmetry. Perhaps the most symmetrical configuration is the one illustrated in Exhibit 9.21. Here, points at the north and south poles of the sphere are symmetrical to each other, but are configured differently from those on the equator.

Exhibit 9.21: A nontrivial perfectly symmetrical configuration does not exist for five points on the surface of a sphere in three dimensions. Displayed here is a less-than-perfect solution. Points at the north and south poles are symmetrical to each other, but are configured differently from three points on the equator.

In higher dimensions, the situation is similar. Certain numbers of points allow for perfectly symmetrical configurations, but for most numbers of points, any nontrivial configuration affords less-than-perfect symmetry.

Accordingly, we don’t seek perfect symmetry, but only an arrangement of points that is as uniform as possible. A convenient solution is to distribute the points in the same manner in which l – 1 electrons would distribute themselves on the surface of a sphere based upon the mutual repulsive forces between them. This concept is defined6 in three dimensions, and the mathematics generalizes to higher dimensions. Treated as electrons, l – 1 points pk distribute themselves to minimize the sum


The configuration of Exhibit 9.21 is such a “minimum energy” configuration. Minimum energy configurations for other values of n and l – 1 can be obtained by computer simulation. First generate a set {} of l – 1 pseudorandom points on the surface of the sphere. A procedure for doing so was described earlier in this section. Next, shift the points around iteratively based upon applicable electrostatic forces to obtain subsequent point sets {}, {}, {}, … Continue until the points stop moving discernibly—perhaps until max(||||) < α1 for some suitable value α1.

At each iteration, the new point set {} is obtained from the current one {} as follows. For each point , calculate a vector-valued “force” ,


where α2 is a suitable scaling factor. Shift each point to obtain the subsequent point:


The algorithm may converge slowly; especially for large n and l. Values α1 = .001 and α2 = 1/(20l) work well for most value-at-risk applications.7

If n and l remain the same each time a value-at-risk measure is used, the above algorithm only needs to be run once. The resulting point set {pk} can be stored for reuse each time the value-at-risk measure is applied. Of course, the points will have to be projected onto different realizations 1r[k] on ellipsoid [9.30] each time because 1|0Σ will change.

With least squares, it may be advantageous to weight the realization 1r[l] = 1|0μ more heavily than others. Our discussion of least squares in Section 2.9 does not mention nonuniform weightings of points. However, this is easily accomplished by considering additional realizations of r[k] and setting them all equal to 1|0μ.