# 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. Jun 5, 2020 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). Jun 5, 2020 at 19:43
• Good point. nls should be more reliable for large $k$. Jun 5, 2020 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. Jun 5, 2020 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. Jun 5, 2020 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? Jun 5, 2020 at 12:55
• I am not sure you expansion in the comments matches your question. In any case, I tried to clarify above. Jun 5, 2020 at 19:05