# How to fit this linear regression with constraints?

I would like to fit the following model:

$$Y=\beta_0+\beta_1(\sum_{i=1}^kw_iX_i)+\beta_2(\sum_{i=1}^kw_iX_i)^2+\epsilon$$ where $$\beta_0, \beta_1, \beta_2, w_1,...,w_k$$ are the parameters, and $$\epsilon$$ is some normal noise. It does not look like something I have encountered before, and I know it's different from simply including all the second order and interaction terms, since the coefficients are related/fixed in a specific way through the sharing of the weights $$w_i$$'s. It seems that this is some sort of linear regression with constraints that relate the coefficients. Could someone point me in the right direction how such model can be fitted?

• It doesn't look like linear regression. – kevin012 Jun 5 at 1:30

The model is overparametrised: you don't need $$\beta_1$$, which can be set to anything convenient, like 1.

One thing I thought of was to fit iteratively. Start out with some guess at $$w$$ and $$\beta_2$$. Then compute $$Z=(\sum_i \hat{w}_iX_i)^2$$ and fit the linear model

Y~ X1+X2+...+X_k + Z

The coefficients of the $$X$$s are the new $$\hat{w}_i$$, and the coefficient of $$Z$$ is $$\hat \beta_2$$. And then recompute Z, iterate and hope it converges. Sadly, it doesn't.

But if $$k$$ isn't too large, it's easy to just compute the residual sum of squares as a function of the parameters and run it through a general purpose optimiser. In R I'd use minqa::newuoa, but there are lots of alternatives.

> X<-matrix(rnorm(50*100),ncol=5)
> w<-1:5
> Y<- (X%*%w)+2*(X%*%w)^2+rnorm(100)
>
>
+   beta2<-theta
+   w<-theta[-1]
+   mu<- (X%*%w)+beta2*(X%*%w)^2
+   sum((Y-mu)^2)
+ }
>
parameter estimates: 1.99478699135839, 1.00032043499982, 2.00140284432351, 3.00312315850919, 4.00284240744153, 5.00537517104468
objective: 1047.51402563294
number of function evaluations: 1689


Then use the bootstrap to get standard error estimates.

With $$k=50$$ it doesn't work (without tuning -- I'm sure it would work if the optimiser defaults were changed or the starting values were better)

• You are correct that the model is over-parameterized (+1). It is worth mentioning that the model could be cast as a nonlinear least squares problem. (Rather, the MLE solution of the model, assuming Gaussian $\epsilon$). So a general nonlinear optimizer may be overkill (e.g. I think nls would be applicable in R). – GeoMatt22 Jun 5 at 19:43
• Good point. nls should be more reliable for large $k$. – Thomas Lumley Jun 5 at 22:26

If you write out the expression, you get a polynomial in terms of $$X_1,X_2,..,X_k$$, including their interactions, where the new "coefficients" are all function of $$\beta$$s and $$w$$s and twos. For k=2, you get a polynomial that has 5 coefficients (or 6 including the intercept) with 4 unknowns:

\begin{align*} Y &= \beta_0+(\beta_1w_1)X_1+(\beta_1w_2)X_2+(\beta_2w_1^2)X_1^2 + (\beta_2 w_2^2)X_2^2+(2\beta_2 w_1w_2)X_1X_2 +\varepsilon \\ &= \alpha_0+\alpha_1X_1+\alpha_2X_2+\alpha_3X_1^2 + \alpha_4X_2^2+\alpha_5X_1X_2 +\varepsilon \end{align*}

If you fit this regression, you will get the new $$\alpha$$ coefficients, which gives you a system of non-linear equations:

\begin{align*} \alpha_0 &= \beta_0 \\ \alpha_1 &= \beta_1w_1 \\ \alpha_2 &= \beta_1w_2 \\ \alpha_3 & =\beta_2w_1^2\\ \alpha_4 &= \beta_2 w_2^2 \\ \alpha_5 &= 2\beta_2 w_1w_2 \end{align*}

In principle, that system of equations should be solvable numerically, at least sometimes. It should remain solvable with $$k>3$$ since you don't have the curse of dimensionality since each new variable adds only one parameters but multiple new equations that help pin it down.

Here's a toy $$k=2$$ simulation example using Stata where I ignore the intercept equation since it is trivial:

. clear

. set obs 1000
number of observations (_N) was 0, now 1,000

. set seed 10011979

. gen b0 = 1

. gen b1 = 2

. gen b2 = 3

. gen w1 = 4

. gen w2 = 5

. gen x1  = rnormal(0,1)

. gen x2  = rnormal(10,2)

. gen eps = rnormal()

. gen y = b0 + b1*(w1*x1 + w2*x2) + b2*(w1*x1 + w2*x2)^2 + eps

. reg y (c.x1 c.x2)##(c.x1 c.x2)

Source |       SS           df       MS      Number of obs   =     1,000
-------------+----------------------------------   F(5, 994)       >  99999.00
Model |  1.1237e+10         5  2.2475e+09   Prob > F        =    0.0000
Residual |  1052.11816       994  1.05846897   R-squared       =    1.0000
Total |  1.1237e+10       999  11248523.6   Root MSE        =    1.0288

------------------------------------------------------------------------------
y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
x1 |   8.082131   .1573906    51.35   0.000     7.773275    8.390987
x2 |   9.852645    .110114    89.48   0.000     9.636562    10.06873
|
c.x1#c.x1 |    47.9813   .0233895  2051.40   0.000      47.9354     48.0272
|
c.x1#c.x2 |   119.9907   .0153233  7830.59   0.000     119.9606    120.0208
|
c.x2#c.x2 |   75.00664   .0053927  1.4e+04   0.000     74.99605    75.01722
|
_cons |    1.77947   .5532575     3.22   0.001      .693783    2.865156
------------------------------------------------------------------------------

.
. clear mata

. mata:
------------------------------------------------- mata (type end to exit) -----------------------------------------------------------------------------------------------------------------------------------------------
: void mysolver(todo, p, lnf, S, H)
>          {
>                  b1   = p
>                  b2   = p
>                  w1   = p
>                                  w2   = p
>                  lnf = (b1*w1 - 8.082131)^2\
>                        (b1*w2 - 9.852645)^2\
>                                            (b2*w1^2 - 47.9813)^2\
>                                            (b2*w2^2 - 75.00664)^2\
>                                            (2*b2*w1*w2 - 119.9907)^2
>                 }
note: argument todo unused
note: argument S unused
note: argument H unused

:
: S = optimize_init()

: optimize_init_evaluator(S, &mysolver())

: optimize_init_evaluatortype(S, "v0")

: optimize_init_params(S, (1,1,1,1))

: optimize_init_which(S,  "min" )

: optimize_init_tracelevel(S,"none")

: optimize_init_conv_ptol(S, 1e-16)

: optimize_init_conv_vtol(S, 1e-16)

: p = optimize(S)

: p
1             2             3             4
+---------------------------------------------------------+
1 |    2.1561597   3.521534782   3.691630188   4.614939185  |
+---------------------------------------------------------+

: end
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------


The solution is not very good (unless you squint and round to the nearest integer), since $$p = (2,3,4,5)$$ in the simulation. I am probably doing something wrong when when I solve the equations numerically. But even the intercept is pretty off with $$b_0 = 1.77947 \ne 1$$.

Code:

cls
clear
set obs 1000
set seed 10011979
gen b0 = 1
gen b1 = 2
gen b2 = 3
gen w1 = 4
gen w2 = 5
gen x1  = rnormal(0,1)
gen x2  = rnormal(10,2)
gen eps = rnormal()
gen y = b0 + b1*(w1*x1 + w2*x2) + b2*(w1*x1 + w2*x2)^2 + eps
reg y (c.x1 c.x2)##(c.x1 c.x2)

clear mata
mata:
void mysolver(todo, p, lnf, S, H)
{
b1   = p
b2   = p
w1   = p
w2   = p
lnf = (b1*w1 - 8.082131)^2\
(b1*w2 - 9.852645)^2\
(b2*w1^2 - 47.9813)^2\
(b2*w2^2 - 75.00664)^2\
(2*b2*w1*w2 - 119.9907)^2
}

S = optimize_init()
optimize_init_evaluator(S, &mysolver())
optimize_init_evaluatortype(S, "v0")
optimize_init_params(S, (1,1,1,1))
optimize_init_which(S,  "min" )
optimize_init_tracelevel(S,"none")
optimize_init_conv_ptol(S, 1e-16)
optimize_init_conv_vtol(S, 1e-16)
p = optimize(S)
p
end

• Unlike the usual constraints problem it's not a linear system of equations, because of the way the sum is inside the square in the second term. It could still work, but it's harder than usual. – Thomas Lumley Jun 5 at 5:53
• @ThomasLumley You're right that this system might not be identified. I tried a toy example above, and it did not quite solve it. – Dimitriy V. Masterov Jun 5 at 7:43
• I am a bit confused here. For the $k=2$ example, the original expression would be $Y=\beta_0+w_1X_1+w_2X_2+\beta_1(w_1^2X_1^2+w_2^2X_2^2+2w_1w_2X_1X_2)$. If I fit instead $Y=a_0+a_1X_1+a_2X_2+a_{12}X_1X_2+a_{11}X_1^2+a_{22}X_2^2$, doesn't this give me more unknowns? Once I match $\beta_0,w_1,w_2$ with $a_0,a_1, a_2$, it is not clear how the remaining can be matched since there are more unknowns in the latter regression. Could you explain? – Xiaohuolong Jun 5 at 12:55
• I am not sure you expansion in the comments matches your question. In any case, I tried to clarify above. – Dimitriy V. Masterov Jun 5 at 19:05