Multiple equality constraints on parameters

I’ve been playing with the R package limSolve’s xsample function recently. For those not aware, it provides a MCMC approach to simulating from a distribution with equality and inequality constraints.

I was able to implement the same approach in Stan for equality constraints (most of the hard work is done outside of Stan). For some variable x, it expresses the problem in terms of x = x0 + Z * q, where x0 is the result of a least squares, equality-constrained problem and Z is the matrix that provides a basis for the null space (using MASS).

So far I haven’t been able to properly implement the inequality constraints in the general case (I was able to get something that produced what I expected where there are N independent inequality constraints, but with lots of divergences). I would be interested if anyone can get that to work.

The R code below creates constraints E and F that that restrict the sum of the variable x to be 1 and the sum of the first and second x’s to be 0.3.

rstan_options(auto_write = TRUE)

wd <- "C:/..." #location of where file is saved

N <- 10
E <- matrix(rep(1, N), 1, N)
E <- rbind(E, matrix(c(c(1, 1), rep(0, N - 2)), 1, N))
F <- matrix(1, 1, 1)
F <- rbind(F, matrix(0.3, 1, 1))

x0pre <- lsei(E = E, F = F)
x0 <- matrix(x0pre$X, N, 1)

tol <- sqrt(Machine$double.eps)
Z <- Null(t(E));
Z[abs(Z)<tol] <- 0

stan_data <- list(N = N, E = E, N_Z = dim(Z)[2], Z = Z, x0 = as.numeric(x0))

fit_eq <- stan(file = 'equality_constr.stan', data = stan_data, 
			   verbose = FALSE, iter = 500, chains = 4)

e_eq_x <- extract(fit_eq, pars = "x", include = TRUE, permuted=TRUE)
e_eq_x <- e_eq_x$x
rowSums(e_eq_x[, 1:2])

This is the accompanying stan code.

data {
    int<lower=0> N;
    int<lower=0, upper=N> N_Z;
    matrix[N, N_Z] Z;
    vector[N] x0;
parameters {
    vector<lower=-1, upper=1>[N_Z] q;
transformed parameters {
    vector[N] x;
    x = x0 + Z * q;
model {
1 Like

I don’t think I follow your logic fully - for me, I would have to see the math and it is hard for me to eyeball it from the code (especially as it uses packages I am not familiar with), but at least in linear programming, you can always transform inequality constraints to equality constraints (e.g. Linear programming - Wikipedia) so a similar technique could IMHO apply here.

Whenever you have F(x) < b, you can introduce a slack variable z > 0 and have an equality constraint F(x) + z = b. Would that help to convert ineqaulities to equalities in your case?

Best of luck with your model!

Hi @martinmodrak

I appreciate your reply.

The math is worked out in one of the limSolve vignettes [1]. Z is a matrix such that EZ=0 (and Z'Z=I), this means that EX=E(x_{0}+Z*q)=Ex_{0}. So if you pick x_{0} such that the equality Ex_{0}=F holds, then EX=F also does.

I’m not sure your solution would work with the equality constraints since (I believe) it would require re-computing the Z matrix to enforce the equality.

I spent a little more time playing with the limSolve package (xsample function) using actual data instead of the simple example I was working with and it turned out to work less well (my sdB variable turned out to be than what I had been testing with, which meant the result was closer to x_{0}, I also tested without inequalities to facilitate comparing to stan). I did something similar in stan using simplexes and got pretty much the same result.


So I’ve been trying to figure out a way to convert the constraints to equality via the slack variables, but I couldn’t make it work.

Some related posts are: Constraints during sampling and Additional constraint on unit simplex (which is yours, so I assume you’ve tried those approaches already :-) Also, the work on Quadratic optimizer seems relevant, but looking at the associated Github issue it looks like the implementation was never finished. But there was some discussion on constraints and convex polytope’s which might (and might not) be relevant.

A big problem for any way to parametrize a general convex polyhedron of dimension n constrained by linear inequalities for use in Stan is that you need a smooth bijective mapping between \mathbb{R}^n and points in the region. The smoothness seems to be a problem. If I for example choose to parametrize a convex polytope by choosing an interior point a unit n-dimensional vector indicating direction from the point and a number between 0 and 1 indicating position between the chosen point and a side of the polyhdreon / infinity, this mapping will not be smooth as the direction traverses individual sides.

EDIT: After reflection the idea in this post is probably misguided, hiding it for now.

Click to see a probably wrong direction

However, there are shapes that are easy to parametrize - e.g. a simplex (in the geometrical sense) can be easily parametrized as a convex combination of its vertices (which is both smooth and bijective), using Stan’s simplex type. Simplices with some sides “open” (i.e. unconstrained) would also be easy to parametrize.

A probably not very practical way to get around this would be to decompose the polyhedron into as few as possible segments that can be easily parametrized and then add a discrete random variable indicating which segment was chosen. This discrete variable could then be marginalized out and this could give a well-behaved parametrization. For non-trivial shapes in higher dimensions this would however introduce a lot of nuisance parameters and would be difficult to implement (and not guaranteed to actually work…)

So sorry, no good answers from me.

Tagging @bgoodri as I believe he has more insight than me into whether something like this is feasible or if there are other good alternatives.

@martinmodrak Thanks for thinking on the issue. That math is sounding like it’s getting beyond my skills.