xsample() : An R Function for Sampling Linear Inverse Problems

An R function is implemented that uses Markov chain Monte Carlo (MCMC) algorithms to uniformly sample the feasible region of constrained linear problems. Two existing hit-and-run sampling algorithms are implemented, together with a new algorithm where an MCMC step reﬂects on the inequality constraints. The new algorithm is more robust compared to the hit-and-run methods, at a small cost of increased calculation time.


Introduction
In linear programming and system theory, a linear model is conventionally written in matrix notation as 1 Ax = b + , with x a vector of unknowns, and an error vector.Additional equality and inequality constraints can be present, leading to a general formulation: This kind of problems are usually overdetermined, meaning that there is no solution for which = 0.They can then be solved with quadratic programming (Lawson and Hanson 1995) techniques, in which case a norm of the error term = Ax − b is minimized, for example In the absence of inequality conditions, it is straightforward to construct a series of samples x for which Ax − b = has the proposed distribution.However, when x is subject to inequality constraints (Gx ≥ h), cannot be normally distributed.
Instead, a truncated normal distribution is proposed for x: where the weight matrix W = diag(s −1 ).This formulation penalizes samples x when ||Ax − b|| increases, and leads to a normal distribution of Ax − b ∼ N (0, s) when there are no constraints.
Equation ( 1) is overdetermined when there is no exact fit Ax = b.then represents a model error term rather than uncertainties in the data: Here the model standard deviation σ is a scalar parameter that is estimated together with the other parameters x (Gelman, Carlin, Stern, and Rubin 2004).In the absence of inequality constraints, the mean estimate of σ equals the standard deviation of the residuals of a weighted linear regression.
The R (R Development Core Team 2008) function xsample() is currently part of the limSolve package (Soetaert, Van den Meersche, and van Oevelen 2009), available under the GPL (General Public License) from the Comprehensive R Archive Network at http://CRAN.R-project.
org/.limSolve contains several tools for linear inverse modeling.Function xsample() takes the matrices A, E, G and the vectors b, f , h as input, together with a vector of standard deviations for b and a number of technical input parameters.In the next sections, the function and contained algorithms are explained, and some examples are provided.

Method
The xsample() function aims to produce a sample set of vectors x that fulfill a number of equality constraints, and are confined by a number of inequality constraints.They are either uniformly distributed within their feasible region, or their distribution depends on the value of linear combinations Ax.This is done in two steps: (1) eliminate the equality constraints Ex = f and (2) perform a random walk on the reduced problem.

Step 1: Eliminate equality constraints
The elements x (i) of x are not linearly independent; they are coupled through the equations in Ex = f .They are first linearly transformed to a vector q for which all elements q (i) are linearly independent.If solutions exist for the equations in (2) and a vector x 0 is a particular solution of Ex = f , then all solutions x can be written as: Z is an orthonormal matrix, obtained from the QR-decomposition or singular value decomposition of E (Press, Teukolsky, Vetterling, and Flannery 1992), and serves as a basis for the null space of E: Z Z = I and EZ = 0.
There are no equality constraints for the elements in q.Thus, the problem is reduced to: In xsample(), a particular solution x 0 of Ex = f can either be provided as one of the input parameters or be calculated by xsample() as a particular solution using the Least Squares with Equalities and Inequalities (LSEI) algorithm (Haskell and Hanson 1981), available in the limSolve package as lsei().
Because p meets the inequality constraints Gp ≥ h, there is already one trivial solution of q: the null vector 0. From this point, new points are sequentially sampled.
We want to know which distribution of q is necessary to obtain the targeted distribution of the sample set x.If a vector x(q) is a function of q, the PDF (probability density function) of q is a product of the PDF of x and the Jacobian determinant: In this case, as Z is orthonormal, the Jacobian is || ∂x ∂q || = |Z| = 1.Therefore p(x) = p(q).This means that if q is sampled uniformly, then x is too.

Markov chain Monte Carlo (MCMC)
What's left to do, is to properly sample q.This can be done numerically using an MCMC random walk.Especially for high-dimensional problems, this is more efficient than a gridbased approach.The Metropolis algorithm (Roberts 1996) produces a series of samples whose distribution approaches an underlying target distribution.In xsample(), new samples q 2 are drawn randomly from a jump distribution with PDF j(.|q 1 ) that only depends on the previously accepted point q 1 .The new sample point q 2 is either accepted or rejected based on the following criterion: if r ≤ p(q 2 ) p(q 1 ) accept q 2 else keep q 1 (8) with 0 < r ≤ 1 and p(•) the PDF of the target distribution.The only prerequisite for the sample distribution to converge to the target distribution with PDF p(•), is that the jump distribution from which a new sample is drawn, is symmetrical in the following sense: the probability to jump from q 1 to q 2 , j(q 2 |q 1 ), has to be the same as the probability to jump from q 2 to q 1 , j(q 1 |q 2 ).Three different jump distributions are implemented and are discussed further below.
In absence of matrix A and vector b, the target distribution of q is uniform and thus: else p(q 2 ) = 0 ⇒ reject q 2 If A and b are present, combining equations (4), ( 6) and ( 7): else p(q) = 0 The expression for fixed standard deviations is easily obtained from (3) by setting σ = 1 and W = diag(s −1 ).Otherwise, σ is estimated from fitting of the unconstrained model Ax − b ∼ N (0, σ).

Sampling the feasible region
New samples in the MCMC are taken from a symmetric jump distribution.A major challenge is to only sample points that fulfill the inequality constraints.Three algorithms that ensure this, are discussed in the next paragraphs.As a consequence, the sample set of vectors q and the derived sample set of vector x, has a distribution that is bounded by the inequality constraints.
In a euclidean space, every inequality constraint defines a boundary of the feasible subspace.Each boundary can be considered a multidimensional plane (a hyperplane).One side of the hyperplane is the feasible range, where the inequality is fulfilled.The other side of the hyperplane is non-feasible.The hyperplanes are defined by the following set of equations: Three jump algorithms for selecting new points q 2 were implemented: Two hit-and-run algorithms (Smith 1984): the random directions and coordinates directions algorithms and a novel mirror algorithm that uses the inequality bounds as reflective planes.All three algorithms produce sample points that fulfill all inequality constraints, and they fulfill the symmetry prerequisite for the Metropolis algorithm.

Random directions algorithm (rda)
The random directions algorithm (Smith 1984) consists of two steps: first a random direction is selected by drawing and normalizing a randomly distributed vector.Starting point and direction define a line in solution space.Then the intersections of this line with the hyperplanes defined by the inequality constraints are determined.A new point is then sampled uniformly along the line segment that fulfills all inequalities.

Coordinates directions algorithm (cda)
The only difference with the random directions algorithm, is that the coordinates directions algorithm (Smith 1984) starts with selecting a direction along one of the coordinate axes.This leads to a simpler formulation of the algorithm.

The mirror algorithm
The mirror algorithm was inspired by the reflections in mirrors and uses the inequality constraints as reflecting planes.New samples are taken from a normal jump distribution with q 1 as average and a fixed standard deviation, called the jump length.With an increasing number of inequality constraints, more and more samples from an unmodified normal distribution will be situated outside of the feasible region and have to be rejected based on criterion (8).While this is a correct approach and the sample distribution will also converge to the targeted distribution, it is inefficient because many points are rejected.We propose an alternative sampling routine that uses the inequalities to ensure that every newly sampled point is situated in the feasible region.
If q 1 is a point for which the inequality constraints are fulfilled, a new point q 2 can be sampled in the following way: first q 2−0 is sampled from a normal distribution in the unrestricted space, ignoring all inequality constraints: with η drawn from a normal distribution with mean 0 and a fixed standard deviation.If q 2−0 is in the feasible range (all inequalities are met), q 2−0 is accepted as a sample point q 2 and evaluated in the Metropolis algorithm (8).
If some inequalities are violated (Figure 1), then the new point q 2−0 is mirrored consecutively in the hyperplanes representing the unmet inequalities: the line segment q 1 → q 2−0 crosses these hyperplanes.For each hyperplane, a scalar α (i) can be calculated for which with η = q 2−0 − q 1 .The hyperplane with the smallest non-negative α (i) , call it α (s) , is the hyperplane that is crossed first by the line segment.q 2−0 is mirrored around this hyperplane.
If the new point (q 2−1 in Figure 1) still has unmet inequalities, a new set of α (i) 's is calculated from the line segment between the new point and the intersection of the previous line segment and the first hyperplane, i.e., q 1 +α (s) η. q 2−1 is again reflected in the hyperplane with smallest non-negative α (i) .This is repeated until all inequalities are met.The resulting point q 2 is in the feasible subspace and is accepted as a new sample point.
In most cases, the directional algorithms and the mirror algorithm converge to the same distributional result.However, we found that especially in high-dimensional problems, the mirror algorithm is still able to move away from the initial particular solution when the directional algorithms fail to do so.One possible explanation for this can be found in the initialisation of the MCMC with LSEI.LSEI often returns a solution in a corner of the feasible region, at the intersection of inequality constraints.In some circomstances, the line segment used by a random directions algorithm has then length zero and the algorithm fails to move away from the initial point.
In the mirror algorithm, η is drawn from a normal distribution with zero mean and a set of fixed standard deviations, which we call the jump lengths of the Markov chain.These jump lengths have a significant influence on the efficiency of the mirror algorithm, as they define the distance covered within the solution space in one iteration, but also the number of reflections in the solution boundaries.They can be set manually with the parameter jmp in xsample().When sampling the feasible region uniformly, a suitable jump length is often in the same order of magnitude as the ranges of the unknowns.When the default parameter setting jmp = NULL is used, a jump length is calculated internally, which gives quick and suitable results in most cases.Sometimes, these internally calculated jump lengths are too large, and the calculation time is too long.One can then turn to manually setting small jump lengths and gradually increasing them, until all elements in x are properly sampled.This can be checked by looking at the trace of the elements x (i) , which need to have an obviously random pattern, as illustrated in Figure 2.
Note that the hit-and-run algorithms rda and cda only work if G and H define a bounded feasible region.In an open or half open space, these algorithms will generate error messages because they draw from a uniform distribution confined by this feasible region.The mirror algorithm is not affected by this problem because new samples are drawn from a normal distribution instead of a uniform distribution.

Using the xsample() function in R
The default input for the xsample() function in R is: with the following arguments: A and B: a numeric matrix and a vector that contain the coefficients of the equations used to penalize the samples according to a normal distribution: Ax = b + .
E and F: a numeric matrix and a vector that contain the coefficients of the equality constraints, Ex = f .
G and H: a numeric matrix and a vector containing the coefficients of the inequality constraints, Gx ≥ h.
sdB: a vector with fixed standard deviations of b.
iter: the number of iterations.
outputlength: the number of samples of x in the output; this value is smaller than or equal to iter.
burninlength: the number of initial iterations that are not included in the output.
type: the algorithm used to sample new points: one of: "mirror", (mirroring algorithm), "rda" (random directions algorithm) or "cda" (coordinates directions algorithm) jmp: jump length of the transformed variables q: x = x 0 + Zq (only if type == "mirror"); if not provided, a reasonable jump length is calculated internally.
tol: tolerance for equality and inequality constraints.
x0: initial (particular) solution x 0 for the equality and inequality constraints.
fulloutput: if TRUE, the transformed variables q are included in the output.
test: if TRUE, xsample() will test for hidden equalities (i.e.equalities that are implemented as a set of two inequalities).This may be necessary for large problems, but slows down execution a bit.
The returned value of xsample() is a list containing: X: a numerical matrix with outputlength rows containing the samples.
acceptedratio: ratio of acceptance (i.e. the ratio of the accepted runs in the MCMC over the total number of iterations).
Q: only returned if fulloutput = TRUE: a numerical matrix with outputlength rows containing the transformed samples q.
p: only returned if fulloutput = TRUE: probability vector for all samples (one value for each row of X) jmp: the jump length used for the random walk.Can be used to check the automatically generated jump length.

Applications of xsample() and interpretations of the results
xsample() is primarily aimed at underdetermined systems.It samples their solution space and in addition, attempts to account for error in the data vector b.However, it also works if the input is overdetermined.A model error can then be estimated from the difference between model and data as well.In this section, we illustrate two underdetermined scenarios for which xsample() may be used: underdetermined problems with A = B = NULL, underdetermined problems with A and B not NULL.We also briefly discuss the use of xsample() for overdetermined problems.

A = NULL and B = NULL: The mink diet composition
The most straightforward application of xsample() is to uniformly sample underdetermined systems with a set of equality constraints (E,F) and inequality constraints (G,H), and no input A or B. The first example is a low-dimensional, underdetermined problem estimating the relative proportion of food sources of a predator based on stable isotope data.
The mink data set, from Ben-David, Hanley, Klein, and Schell ( 1997) is a small data set used in the program Isosource, which estimates the contributions of food sources to the diet of a consumer, based on stable isotope compositions (Phillips and Gregg 2003).Isosource searches all possible compositions with a grid-based approach and retains those compositions that are sufficiently close to the equality constraints imposed by the isotopic data.Our method obtains similar results, but obtains exact solutions instead of approximations.It also searches the solution space in a more direct and efficient way.
Stable isotope ratios of carbon (δ 13 C) and nitrogen (δ 15 N ) in a predator, the Canadian Mink, q q q q q q q −24 −22 −20 −18 −16 −14 and seven candidate food sources are provided (Figure 3).The data are given without error estimate, and are considered exact.As the isotopic composition of the predator is a simple weighted average of the compositions of its food sources, and the fractions sum to one, there are three equality constraints.Also, all fractions have to be positive, which gives seven inequality constraints.In matrix notation this can be written as: Filling in the numbers, this results into: 14.7 9.7 11.2 12.0 10.1 14.0 14.9 −14.2 −18.5 −15.3 −16.9 −24.6 −18.7 −21.3 The problem has three equations with seven unknowns and is thus underdetermined.The first two equations from (15) are implemented in limSolve as dataset Minkdiet.A set of samples that are uniformly distributed over the solution space defined by (15), are generated in R as follows: R> xs <-xsample(E = rbind(Minkdiet$Prey, rep(1,7)), + F = c(Minkdiet$Mink,1), G = diag( 7), H = rep(0,7), iter = 3000) R> pairs(xs$X) The results using the mirror algorithm are shown in Figure 4.With only 3000 jumps, the solution space is sampled thoroughly and the conclusion is that fish and crab are the most probable food sources for the mink.These results match the results presented by Phillips and Gregg (2003), but are obtained in a more straightforward manner.
When the only information we have on a set of parameters are a number of linear relationships and parameter ranges, as in this example, then to our best knowledge, every parameter set that fulfills these relationships has equal probability of occurring.The samples produced by xsample() are distributed uniformly over the feasible region of these parameters.They can be used to estimate value ranges, quantiles and means of the expected parameter values.
In the diagonal of Figure 4, histograms of the parameter values are shown.While the probability of each parameter set is equal, this is obviously not the case for the marginal distributions of the individual parameters. is considered exactly 3.5, then the xsample() input becomes: and the solution is trivial: the estimates of x (1) and x (2) each have a uniform distribution between 0 and 3.5, while their sum is constant.The result is shown in the left panel of Figure 5. Now consider a normal distribution for the sum of x: u = x (1) + x (2) ∼ N (3.5, 5).As argued in the introduction, this is not possible because u ≥ 0, and expression (3) is used instead: The probability outside the feasible region ( x (1) < 0 or x (2) < 0), is still assumed zero.
The xsample() input then becomes: The results are illustrated in the right panel of Figure 5.Although p(x) is highest when x (1) + x (2) = 3.5, this is not the case for the distribution of the sum itself.Because of the shape of the feasible region, there are more feasible vectors x for larger values of u = x (1) +x (2) (Figure 5), causing the distribution of u to be skewed towards these larger values.
Generally, the shape of the feasible region will influence the distribution of Ax in a complex way which makes it difficult to produce a sample set x for which Ax has a predictable distribution.

Applying xsample() on overdetermined systems
In case the problem is overdetermined, a set of samples can also be generated from bootstrapping a constrained linear regression, e.g., using orlm() from the R package ic.infer (Groemping 2008).The results will however be different: samples that are located outside of the feasible region in the unconstrained model, are mapped on the borders of the feasible region.This is consistent with formulations of the underlying distribution of the parameters in a constrained linear regression (Shapiro 1988).
In contrast, xsample() will ignore the points outside of the feasible region and thus lead to a truncated normal distribution.The distribution of the samples can therefore not be interpreted as a probability distribution of the parameter estimate in a constrained linear regression.One should always be aware of these differences.

Concluding remarks
The xsample() algorithm successfully and efficiently produces a sample set that is uniformly distributed over the feasible region of an underdetermined linear problem.The sample set can also be non-uniformly distributed around an extra set of linear equations, and has a truncated normal distribution when the problem is overdetermined.
Quick convergence of the sample set to the target distribution (3) or ( 4) is one of the advantages of the algorithm.Convergence can be tested using tools for output analysis and diagnosis of MCMC, such as the R package coda (Plummer, Best, Cowles, and Vines 2009).
Calculation times are in most cases in the order of seconds to minutes on a modern pc.This makes the function an attractive alternative to linear or quadratic programming, offering a more complete answer to underdetermined linear problems and complementing the more traditional approaches.

Figure 2 :
Figure 2: A good random walk of a parameter x i , using xsample() with 1000 iterations.

Figure 3 :
Figure 3: Isotopic composition of Canadian mink and seven candidate food sources.

Figure 4 :Figure 5 :
Figure 4: 3000 MCMC samples of the solution space of the mink diet problem: The lower half of the figure contains pairwise scatter plots of the fractions of the different food sources contributing to the diet.In the diagonal are density plots of the marginal distributions of these fractions.