Solving linear regression with weights and constraints

Question

I would like to solve a linear regression (in R) with weights $w$ and a constraint.

In other words, I would like to find $x$ that minimizes the sum of squares $$\sum_i w_i(b_i-Ax_i)^2$$

On top of that I have an external vector $d$, which I would like to use in a constraint, such that $d \cdot x \le 5$.

Is this something that would be possible to do in R with solve.QP or perhaps some other R script?

Edit: I am adding a bounty for a solution that doesn't require any other custom software except the cran packages. While rstan works perfectly unfortunately I am unable to install it on my production servers due to old versions of some libraries.

Answer 1

You're looking for the mgcv package. With the toy data we used before, it works just fine. (I'm uncertain why rstan is so confident in its results... I'm still looking into it.)

set.seed(1880)

N       <- 1500
d       <- c(1/2, 2/pi, 2/3)
x       <- c(2, 1, 3)
limit   <- 5
d%*%x <= limit

A       <- cbind(1, rnorm(N), rnorm(N))
b.hat   <- A%*%x
wgt     <- rexp(N)
b       <- rnorm(N, mean=b.hat, sd=wgt)

library(mgcv)

pin <- c(1.5, .75, 2.5)
Ain <- matrix(d, nrow=1)

M   <- list(y=b, w=wgt, X=A, p=pin, Ain=-Ain, bin=-limit, C=matrix(1, ncol=0, nrow=0))
pcls(M)

1.8844996 0.9421333 2.9770852

The inequality in this package is flipped the other direction by default. So we have to multiply both sides by $-1$.

Answer 2

Whenever I have a complicated model to fit, I usually just fit it directly in rstan because it's great at fitting highly constrained coefficients, and because it's easy to include penalties and transformations of variables. This is true even when I'm not explicitly fitting a Bayesian model.

This is what I've worked up for your particular problem.

library(rstan)

set.seed(1880)

N       <- 1500
d       <- c(1/2, 2/pi, 2/3)
x       <- c(2, 1, 3)
limit   <- 5
d%*%x <= limit
> TRUE
A       <- cbind(1, rnorm(N), rnorm(N))
b.hat   <- A%*%x
tau     <- 5
wgt     <- rexp(N)
Sigma   <- tau*wgt
b       <- rnorm(N, mean=b.hat, sd=Sigma)

constrained.reg <- "
    data{
        int<lower=1>        N;
        int<lower=1>        K;
        vector<lower=0>[N]  wgt;
        matrix[N,K]         A;
        vector[N]       b;
        vector[K]       d;
        real            limit; // s.t. d*x<=limit
    }
parameters{
    real<upper=limit>   c; // this is the largest possible value of x%*%d.
    simplex[K]      sim_x;
    real<lower=0>       tau;
}
transformed parameters {
    vector[K]   x;
    vector[N]   b_hat;
    vector[N]   Sigma;

    x       <- d .*sim_x /c;
    b_hat   <- A*x;
    Sigma   <- tau*wgt;
}
    model{
        b ~ normal(b_hat, Sigma);
        increment_log_prob(-2*log(tau)); // uniform prior on beta, noninformative prior on tau
    }
    generated quantities{
        vector[N]   resid;
        resid   <- (b_hat-b) ./Sigma;
    }
"
fake.data   <- list(N=N, A=A, K=3, b=b, wgt=wgt, d=d, limit=limit)

fit.test    <- stan(model_code=constrained.reg, data=fake.data, iter=10)

system.time(fit     <- stan(fit=fit.test, iter=1000, data=fake.data))
print(fit, c("x", "tau")); x

I realized that I was being dense and that we can enforce the inequality by sampling a value as large as the maximum permissible dot product result and then transforming appropriately.

     mean se_mean   sd 2.5%  25%  50%  75% 97.5% n_eff Rhat
x[1] 1.99       0 0.01 1.98 1.98 1.99 1.99  2.00  1645 1.00
x[2] 0.99       0 0.01 0.97 0.98 0.99 0.99  1.00   624 1.00
x[3] 3.00       0 0.01 2.98 2.99 3.00 3.01  3.02   945 1.00
tau  4.82       0 0.09 4.62 4.76 4.82 4.88  5.00   558 1.01

These results look fine to me.