Generators for Discrete Polynomial L₁ Approximation Problems^*

1. Introduction

Recent years have seen increased interest in least absolute deviation (L₁) data fitting, either as an alternative to or in conjunction with the usual least squares approach [11, 19].¹ Estimation in the L₁ norm has certain desirable statistical properties (such as “robustness” [14, 16] when the underlying error distribution is long-tailed) and it can be carried out using reasonably efficient mathematical programming procedures.

While a number of algorithms for discrete L₁ approximation have recently been advanced [1, 3, 5, 6, 21, 22], computational comparisons of specific implementations (“codes”) for these approaches have not only been limited [3, 6, 22) but have often produced conflicting evidence. Resolution of these conflicts awaits the establishment of reliable and comprehensive methodology for testing such mathematical software.

As a first step in developing a sound evaluation methodology for comparing mathematical programming codes, four L₁ codes (representing a range of solution techniques and implementation strategies) were evaluated on particular classes of test problems [13]. The first group of test problems were “hand-picked” problems [20] – those created with a specified structure in mind or arising from actual applications. The second group consisted of pseudo-randomly generated problems, produced by means of a test problem generator [17] and representative of a range of general L₁ data fitting problems; such problems could be constructed with a number of controllable characteristics such as degeneracy, rank loss and optimal solution.

The two generators POLY1 and POLY2 described in this paper have subsequently been utilized to produce other types of problem classes for which L₁ approximation is appropriate. In particular, these are problems in which it is desired to obtain the best L₁ polynomial approximation to a specific function (POLY1) or to a set of discrete observations (POLY2). Polynomial approximation problems are frequently encountered in practice (e.g., among our hand-picked problems), and also enjoy an extensive theoretical basis [7, 18, 23]. Moreover, such problems are known to admit a range of “ill-conditioning” and numerical difficulties [10], notably evidenced in our previous testing efforts using 27 hand-picked problems [13] and in recent results obtained using the first of these two generators [8].

The use of generated test problems in code evaluation offers several advantages. Generators provide a virtually inexhaustible source of test problems, and a source that can be made relatively transportable from one computer to another. Since generators are able to provide test problems with controllable characteristics (such as size and structure), fairly well-defined classes of problems can be efficiently created. Moreover, through the invocation of pseudo-random number generators, test problems can be “randomly” selected from such problem classes and valid statistical conclusions can be drawn about code performance with respect to these classes. Another desirable property of test problem generators, and one that is useful in judging the accuracy or correctness of a code, involves the ability to specify in advance the solution to a generated problem.

Two generators are presented here for obtaining polynomial approximation test problems for L₁ curve-fitting. Both generators permit automatic generation of test problems with stipulated characteristics, and one of them (POLY2) allows the specification of known and unique solutions to the generated problems. The theoretical design of both generators is discussed in section 2, and the following two sections describe computer implementations of the theoretical design. The Appendices contain FORTRAN listings for both generators, written to be machine-independent.²

2. Design of the Generators

The discrete linear L₁ approximation problem can be formulated as follows. Given a set of n observations on a single dependent variable y and each of m + 1 independent variables z₀, …, z_m, find parameters β₀, …, β_mthat minimize

The L₁ polynomial approximation problem, with which we will be concerned, is a restricted case of the above formulation: namely, minimize by choice of β = (β₀, …, βm)

In the above, y_i refers to the i-th observation on variable y and x_i refers to the i-th observation on (the single independent) variable x. Equivalently, it is required to find a polynomial of degree m which best fits the given data, in the sense of minimizing the sum of absolute values of the residuals

A vector β* which yields the minimum value ϕ* of ϕ(β) is termed a solution vector, with optimum objective function value ϕ*. Unlike the case of L₁ polynomial approximation over a continuous interval [23, p. 38], the optimum solution vector β* to a discrete L₁ polynomial approximation problem need not be unique. It is known [2, 4, 12] that such a solution vector can always be found which interpolates at least m + 1 of the data points; that is, at least m + 1 residuals e_i are identically zero at the solution. A given problem exhibits degeneracy if more than m + 1 residuals equal zero at the optimum.

A number of algorithms for solving discrete L₁ approximation problems make essential use of the above “interpolation” or “extreme point” result. Namely, these algorithms examine only basic solutions, where precisely m + 1 residuals are zero, and move in a systematic way from one basic solution to another. Any such basic solution is defined by a set of m + 1 indices (or rows) i ϵ {1, 2, …, n} where e_i = 0. For convenience, these indices are said to correspond to active constraints, while the remaining indices correspond to inactive constraints. Thus, in the absence of degeneracy the set N = {l, 2, …, n} can for any basic solution be partitioned as N = N⁰ ∪ N⁺ ∪ N⁻, where N⁰ indicates the active constraints, N⁺ the inactive constraints with positive residuals, and N⁻ the inactive constraints with negative residuals.

The two generators for polynomial approximation problems discussed here will, for specified n and m, produce the following test problem data:

³ and

The basic difference between the generators resides in the type of problem structure they simulate.

The first generator, POLY1, is intended to model the situation where one wishes to approximate, over a discrete set, a continuous function f(x) by a polynomial of specified degree. Accordingly, the values y_i, in (2) are given by y_i, = f(x_i). Controllable features of this generator include: the function f(x), the real interval I over which f(x) is to be approximated, the number and distribution of observation points x_i ϵ I, and the degree m of the approximating polynomial. A more detailed description of this generator is provided in section 3. It should be noted that POLY1 does not produce a known solution to the test problem created; therefore, other means are necessary in order to verify whether or not a given L₁ code has actually “solved” a generated problem. For this reason, we have also devised a computer routine to check optimality of codeproduced solutions by examination of the Kuhn-Tucker conditions for an associated linear program [4, 17].

The second generator, POLY2, models the situation where an L₁ fit is required to a given discrete set of data (y_i, x_i), i = 1, …, n. Unlike the case for POLY1, the y_i are not assumed to be generated via some known functional relation. Controllable features of this second generator include: the real interval of approximation I, the number and distribution of observation points x_i ϵ I, the degree of the approximating polynomial, and the statistical distribution of the residuals e_i Moreover, one is also able to specify in advance a solution vector β* which will be guaranteed to be the unique L₁ solution to the generated problem.

The approach used for generating a data set (X, y) with known optimal solution β* will now be outlined. The matrix X in (1) is partitioned as

corresponding to indices i ϵ N⁰, N⁺, and N⁻ respectively. Note that X⁰ is a square (m + 1) × (m + 1) nonsingular matrix, with the proviso that the x_i’s for i ϵN⁰ are distinct. By a result proved in [17], a sufficient condition for the optimality and uniqueness of β* in (X, y) is that

and that

where

It will now be shown how values x_i can be found that satisfy (3).

First, we note that for j = 0, equation (3) requires ∣N⁺∣ = ∣N ⁻∣. For notational convenience, let r = ∣N⁺∣ = ∣N⁻∣ and redefine the x_i’s (i ϵ N⁺) as v₁, …, v_r; similarly, let the x_i’s (i ϵ N⁻) be denoted as w₁ …, w _r. Also, write

if

Then, following the development in [15], it is straightforward to verify that for any d

Relation (6) can be used recursively to produce values x_i that satisfy (3). For example, since

equation (6) gives

and so forth. Accordingly, the consecutive integers 0, 1, 2, …, 2^q+1 — 1 can be divided into two disjoint subsets {v₁, …, v_r} and {w₁, …, w_r}, with r = 2^q, such that [v₁, …, v_r]_q = [w₁, …, w_r]_q. The first subset consists of all integers l, 0 ≤ l ≤ 2^q+1 − 1 wilth an even number of ones in their binary expansion, and the second subset consists of all such integers with an odd number of ones. Dividing each v_p and w_p by 2^q+1 −1 produces a set of numbers equally spaced on [0, 1], which in turn can be scaled and translated to yield a set of numbers $\left\{v_p^{\prime },w_p^{\prime }:p=1,\dots ,r\right\}$ equally spaced on an arbitrary interval I = [a, b]. Moreover, relation (5) still holds for the scaled and translated $v_p^{\prime },w_p^{\prime }$. Finally, these values can be used in (3) as the required x_i’s for i ϵN⁺ and i ϵN⁻, respectively.

More generally the recursive process based on (6) can be initiated, using any positive integer d₀, by means of

and propagated using any d₁, d₂, …, d_q to produce 2^q+1 numbers satisfying (5), and therefore (3). We have chosen (see sect. 4) to generate the d_j’s from a uniform probability distribution. The resulting (scaled and translated) values for x_i are no longer equally-spaced on I, but tend to be approximately normally distributed within the (finite) interval I. With probability one, all such x_i values generated are distinct.

The above procedures form the theoretical basis for implementation of the second generator. Either equidistant x_i or nonequidistant x_i (corresponding to inactive constraints) can thus be generated. To ensure that (3) holds for j = 0, …, m, the number of observations n is required to be of the form

where k is an integer, k ≥ 0. Further details of POLY2 are described in section 4.

3. Computer Implementation of POLY1

This section describes a computerized version, in FORTRAN, for the first generator (POLY1) discussed in section 2. We will indicate those problem design characteristics available as user-specified input options to P0LY1, and then discuss how these inputs are utilized in creating test problem data.

A variety of characteristics are available to the user as controllable options and are specified through input parameters to POLY1. These options include:

1. The number of observations n.

2. The degree m of the approximating polynomial.

3. The specific function f(x) to be approximated. Ten such functions,⁴ representing a variety of shapes likely to be encountered in practice, have been incorporated within the program (see table 1 and fig. 1).

4. The interval of approximation I = [a, b].

5. A switch ICOL to indicate whether the approximating polynomial includes a constant term or not. (The latter case is useful in forcing the approximating polynomial to pass through the origin.)

6. The location of the n observation points x_i within I: either drawn from a uniform probability distribution or equally-spaced over I.

7. A perturbation factor EPS, possibly zero, used in perturbing equally-spaced observation points x_i.

8. Two initial starting seeds, utilized internally by a pseudo-random number generator.

Given these input parameters, generation of the data sets can begin. The first step involves selecting the n observation points x_i, from the desired interval, according to the distributional form specified. The interval I = [a, b] is determined by the user (if CSWTCH ≠ 0) through input parameters ENDL = a and ENDR = b, or (if CSWTCH = 0) by means of the default interval settings shown in table 1. The observation points are chosen according to a uniform probability distribution over I (if MSWTCH ≠ 0) or equally-spaced over I (if MSWTCH = 0). In the latter case, the user has the ability of perturbing the equally-spaced observation points x_i: namely, if I = [a, b] then

where $\Delta =\frac{b-a}{n-1}$ is the common subinterval length for equal spacing, the α_i are independent uniform random variates over [−1, 1], and EPS is a user-supplied input representing the (maximum) fraction of the subinterval length Δ used for perturbation. If EPS = 0.0, no perturbation is performed; otherwise, a reasonable range for EPS is 0 < EPS < 1/2.

Next, the values y_i — f(x_i) are calculated, as well as the successive powers $x_i^j$. If the user has specified a constant term in the approximating polynomial (ICOL = 1), then the above index j has range j = 0, 1, …, m. Otherwise (ICOL ≠ 1), a constant term is not included and j = 1, …, m. This stage of calculating successive powers $x_i^j$ is most sensitive to finite-precision arithmetic and to accumulated round-off error. Accordingly, such quantities are calculated in double-precision, using the recursion $x_i^{j+1}=x_i\left(x_i^j\right).$.

Finally, the problem data are returned in a single matrix PRBMAT. The first column of PRBMAT contains the vector y, and the remaining columns contain the design matrix X. In addition, the number of rows (NOBS) and the number of columns (NPAR) in matrix PRBMAT are returned to the calling routine, as well as a single-precision vector containing the observation points x_i. The final seeds available from the random number generator are also returned in ISEED and JSEED, thus allowing the creation of different test problems (having the same input specifications) in successive calls of POLY1 from a main routine.

All error messages are written out using unit IOUT = 6. Changing this specific value of IOUT once at the beginning of the program permits other print unit numbers to be accommodated. Appendix A gives a complete listing of POLY1, together with the subroutine UNIRAN [9], a well-tested uniform random number generator which it calls.

Figure 2 illustrates the type of test problem that can be created by POLY1. The n = 50 data points (y_i, x_i) are derived from function 7 of table 1 over I = [0, 5] and are shown here by the symbol “x”; in this case, the x_i’s were chosen to be equally-spaced (without perturbation) over I. By way of reference, the solid curve indicates the best L₁ fitting polynomial of degree m = 5 to these generated data.

4. Computer Implementation of POLY2

This section describes a computer implementation of the second generator (POLY2) for polynomial approximation problems. As indicated in section 2, an important feature of this generator is its ability to create data sets whose optimal L₁ solution vector can be specified in advance. Consequently, the correctness and accuracy of solutions produced by various L₁ codes can be readily assessed using the known solutions to such problems. The fact that these solutions are guaranteed to be unique is also advantageous, since in the case of nonunique solutions it becomes difficult to compare the efficiency and accuracy of codes that reach correct, but different, solutions.

A number of problem characteristics can be controlled by the user through appropriate input specifications. Options available in POLY2 include:

1. The degree m of the approximating polynomial.

2. The value k in (7), which together with m determines the number of observations n.

3. The (unique) solution vector BETA to the generated problem.

4. The interval of approximation I = [a, b].

5. The location of the observation points x_i corresponding to active constraints: either drawn from a uniform probability distribution or equally-spaced over I.

6. The location of the observation points x_i, corresponding to inactive constraints: either equally-spaced or unequally-spaced over I.

7. A perturbation factor EPS, possibly zero, used in perturbing equally-spaced x_i associated with active constraints.

8. The distribution of the residuals e_i associated with inactive constraints.

9. Two initial starting seeds, utilized internally by pseudo-random number generators.

The first step in developing a test problem is to select the solution vector BETA, if the user has not chosen to input a vector to POLY2. The solution vector is randomly selected, in this case, from a uniform probability distribution over an interval specified by the user (BSW = 1) or defined within the program (BSW ≠ 1).

Next, the observation points x_i are generated from the interval I according to the desired distribution. Interval I is either specified on input (CSW = 1) or assigned a default setting within the program (CSW ≠ 1). The x_i’s corresponding to active constraints are derived from a uniform probability distribution over I (MSW = 1) or are equally-spaced over I (MSW ≠ 1). In this latter case, the x_i can be perturbed using the factor EPS, exactly as described for POLY1. The observation points x_i corresponding to inactive constraints are generated so that equation (3) is satisfied. Such x_i are either equally-spaced (ISW = 1) or unequally-spaced (ISW ≠ 1) over the interval I; see section 2 for details of how these two methods of generation are achieved.

Also, the entries $x_i^j$ of the design matrix are successively calculated for j = 0, …, m using appropriate care to preserve numerical accuracy. At the same time, the residuals e_i, i ϵN⁺ or i ϵN⁻, are generated from either a uniform (RSW = 1) or a normal distribution (RSW ≠ l), are given the appropriate sign, and are then used to calculate y_i from (4a) – (4b). The optimum objective function value SUMRES is calculated in double-precision, based on these residuals. The input parameter RESID indirectly controls the magnitude of this objective function value. Namely, RESID acts as a scale factor with respect to the size of the residuals, which are either selected from a uniform distribution on [– RESID, RESID], whose standard deviation is thus RESID/√3, or from a normal distribution having mean 0 and standard deviation RESID.

Finally, the problem data are returned as before in a single matrix PRBMAT. The first column of PRBMAT contains the y vector, and the remaining columns contain the design matrix X. The number of rows (KKM1) and the number of columns (M2) are also returned, as well as BETA, SUMRES and LOC (a vector containing the m + 1 row locations of the active constraints). The final seeds available from the last call to the random number generator are returned in ISEED and JSEED.

All error messages are written out using unit IOUT, which can be easily modified to accommodate differing print unit assignments. Appendix B gives a complete listing of POLY2, together with the subroutines it calls: SORTM, SET1, SET2, EQUAL, NORRAN, SORTP and UNIRAN. The call structure of these subroutines is shown in figure 3; the latter three subroutines are fairly well-tested and are described further in [9].

Figures 4 and 5 illustrate the types of data sets that can be produced using POLY2. Both show as dots 134 data points (y_i, x_i) together with the best L₁ fitting polynomial of degree 5 (solid curve). In figure 4, the x_i are generated via the equally-spaced option (ISW = 1), while in figure 5 they are generated via the unequally-spaced option (ISW ≠ 1).