Constrained Regression in R

Question

I'm using the R type-provider from F# to access some regression related R functionality. I would like to estimate a regression when there is a constraint on the regression coefficients, so that their weighted average is 0. The weights sum to 1. The below example is simplified as I have dozens of coefficients, with varying weights, I only show the R code below:

y1 <- runif(n = 50,min = 0.02,max=0.05)
y2 <- runif(n=50,min=0.01,max=0.03)
y <- c(x1,x2)
x1 <- c(rep(0,50),rep(1,50))
x2 <- c(rep(1,50),rep(0,50))
lm(y~x1+x2)

This gives the output of

> lm(y~x1+x2)

Call:
lm(formula = y ~ x1 + x2)

Coefficients:
(Intercept)           x1           x2  
    0.03468     -0.01460           NA

as expected. However I would like to place a constraint on x1 and x2, so their weighted average is (0.5 * x1 + 0.5 * x2) = 0. In that case the intercept becomes mean(y) = 0.02737966 and the x1 and x2 coefficients will show the offset from this value (-0.006 and +0.007 respectively). It seems the packages quadprog and mgcvare applicable however I wasn't able to apply the constraints.

@ZheyuanLi thanks for taking a look at this. All the covariates are numeric and they sum to 1, hence the need for the constraint. There are also some other unconstrained variables but I skipped them here. By the look of it PCLS in mgcv seemed somewhat relevant. Do you think this is better suited for cross validated? — s952163
@s952163 have you thought about forming an MLE and using constrained optimization? Nlopt works in F#. I can also take a look after work. — professor bigglesworth
@professorbigglesworth It is indeed a form of constrained optimization, hence I was looking into quadprog on R. Hmmm.... I will need to look into NLOPT. Thanks for the idea! — s952163

professor bigglesworth professor bigglesworth · Accepted Answer · 2016-07-23T05:41:57

Maybe not exactly an answer to your question, since it asks for doing the optimization in R. But maybe the following helps. It uses the NLopt library anyway which I think is what R uses? Let me know if you need help in formulating the MLE but for a linear model with gaussian assumptions and no endogeneity it should be straightforward enough.

Note that even though LN_COBYLA doesn't use user supplied gradients, the match with pattern in cFunc and oFunc ignores it. I tried with LD_LBFGS but that doesn't support AddEqualZeroConstraint().

[EDIT]

Adding complete example you can use as template. Its not idiomatic, and quite ugly, but illustrates the point. However, in this example, the constraints will cause this to degenerate. You need NLOptNet, MathNet.Numerics, Fsharp Charting. Maybe it helps other people looking to do constrained optimization in F#.

open System
open System.IO
open FSharp.Core.Operators
open MathNet.Numerics
open MathNet.Numerics.LinearAlgebra
open MathNet.Numerics.LinearAlgebra.Double
open MathNet.Numerics.Distributions
open DiffSharp.Numerical.Float64
open NLoptNet

let (.*) (m1 : Matrix<float>) (m2 : Matrix<float>) =
    m1.Multiply(m2)

let (.+) (m1 : Matrix<float>) (m2 : Matrix<float>) =
    m1.Add(m2)

let (.-) (m1 : Matrix<float>) (m2 : Matrix<float>) =
    m1.Subtract(m2)


let V = matrix [[1.;      0.5;    0.2]
                [0.5;      1.;    0.]
                [0.2;      0.;    1.]]
let dat = (DenseMatrix.init 200 3 ( fun i j -> Normal.Sample(0., 1.) )) .* V.Cholesky().Factor
let y = DenseMatrix.init 200 1 (fun i j -> 0.)
let x0 = DenseMatrix.init 200 1 (fun i j -> 0.)
let x1 = DenseMatrix.init 200 1 (fun i j -> 0.)
for i in 0 .. 199 do
    y.[i, 0] <- dat.[i, 0]
    x0.[i, 0] <- dat.[i, 1]
    x1.[i, 0] <- dat.[i, 2]

let ll (th : float array) =
    let t1 = x0.Multiply(th.[0]) .+ x1.Multiply(th.[1])
    let res = (y .- t1).PointwisePower(2.)
    res.ColumnAbsoluteSums().Sum() / 200.

let oFunc (th : float array) (gradvec : float array) =
    match gradvec with
    | null  -> ()
    | _     -> (grad ll th).CopyTo(gradvec, 0)

    ll th

let cFunc (th : float array) (gradvec : float array) =
    match gradvec with
    | null  -> ()
    | _     -> (grad ll th).CopyTo(gradvec, 0)

    th.[0] + th.[1]

let fitFunc () =
    let solver = new NLoptSolver(NLoptAlgorithm.LN_COBYLA, uint32(2), 1e-7, 100000)
    solver.SetLowerBounds([|-10.; -10.;|])
    solver.SetUpperBounds([|10.; 10.;|])
    //solver.AddEqualZeroConstraint(cFunc)
    solver.SetMinObjective(oFunc)
    let initialValues = [|1.; 2.;|]
    let objReached, finalScore = solver.Optimize(initialValues)
    objReached |> printfn "%A"
    let fittedParams = initialValues
    fittedParams |> printfn "%A"
    fittedParams

let fittedParams = fitFunc() |> DenseVector
let yh = DenseMatrix.init 200 1 (fun i j -> 0.)
for i in 0 .. 199 do
    yh.[i, 0] <- dat.[i, 1] * fittedParams.[0] + dat.[i, 2] * fittedParams.[1]


Chart.Combine([Chart.Line(y.Column(0), Name="y")
               Chart.Line(yh.Column(0), Name="yh")
               |> Chart.WithLegend(Title="Model", Enabled=true)] )
               |> Chart.Show

Constrained Regression in R

1 Answers