# Regression equation as constraint function in an optimization problem?

I have an idea and I want to check if it is already a thing.

I have an optimization problem currently solvable using sequential least squares. If it matters, the objective function and constraints are determined according to this paper but I think my idea is more general.

Edited to add this paragraph about the original problem. I am fitting parameters of a dynamical model (e.g. a system of ODEs). The paper linked above describes a method of formulating an optimization problem for this purpose. We have a model of a system, an observation of one state variable of the system, and unknown parameters and other state variables. The method proceeds roughly as follows. First, we restate our system of ODEs as a difference equation (e.g. using RK4 or Euler's method). We add to this a synchronization term that drives the model toward the observed trajectory at every time step. So if our original difference equation is $$\mathbf{y}_{n+1} = \mathbf{y}_n + \mathbf{F}(\mathbf{y}_n, \mathbf{r})$$ where $\mathbf{y}_n$ is the system state at step $n$ and $\mathbf{r}$ is the unknown parameter vector, adding the synchronization term we have $$\mathbf{y}_{n+1} = \mathbf{y}_n + \mathbf{F}(\mathbf{y}_n, \mathbf{r}) + u_n(x_n - y_{1,n})$$ where $u_n$ is a time-varying synchronization constant and $x_n$ is the observation of the first state variable.

We then formulate the optimization as follows. The objective function is $$\frac{1}{N}\sum^n ((x_n - y_{1,n})^2 + u_n^2),$$ in other words we are trying to choose the parameters that minimize the necessary amount of synchronization as well as the squared difference between the observed and model trajectories.

The difference equation is enforced as constraint functions, and the optimization solution vector is the concatenation of $\mathbf{y}_n$, $\mathbf{r}$, and $u_n$ for all $n$. So we are asking the solver to find the parameter values, synchronization constants, and the associated unobserved states that correspond to a simulation of the system with the minimal necessary synchronization to match the observed trajectory.

This method is called Dynamical Parameter Estimation by Abarbanel. I am now trying to add regression-like relationships among the parameter values $\mathbf{r}$. The rest of the original post follows.

Let's say my solution-space vector has 20 elements. I hypothesize that 10 of these are linear functions of 10 known values with unknown intercept $\beta_0$ and coefficient $\beta_1$. E.g., a simple regression. I see two ways to modify the optimization problem to include this assumption. Either

• replace the 10 values in question with $\beta_0$ and $\beta_1$, leaving a 12-D solution space. Of course, modify the constraint and objective functions accordingly. Or,
• add $\beta_0$ and $\beta_1$ to the solution vector, resulting in a 22-D solution-space, but reduce the degrees of freedom by adding the regression equations as 10 constraint functions.

So far so good, I think. So my first question is whether this is already a thing that someone has done.

My second question is extending this from simple regression to mixed models. Imagine now that those 10 values are sampled from a normal distribution with unknown mean and variance. Intuitively this is an equivalent reduction in the overall degrees of freedom as the previous example. Instead of $\beta_0$ and $\beta_1$ we now have $\mu$ and $\sigma^2$.

However, it doesn't appear to be so simple, because in order to allow the values to vary we must also include the residuals in the solution vector. The best solution I can think of is to add a penalty term to the objective function. In particular, the log likelihood of the normal distribution computed from the 10 values in question. Thus, in addition to the previous objective, the solver would try to maximize the likelihood that the 10 values come from a normal distribution. This method would not actually reduce the overall degrees of freedom but you could think of the added penalty in the objective function as a "soft" constraint.

A remaining question is, if this is all feasible, how to do hypothesis testing, i.e., how to decide if the addition of a fixed or random effect improves the resulting solution. But for now I'll stick to the main question: does this already exist? And if not, does it sound feasible? I'm also a bit over my head if it's not apparent so literature suggestions are welcome as well.

• I didn't really follow what you're doing. But for hypothesis testing, how about bootstrapping the whole procedure? – Mark L. Stone Nov 7 '16 at 4:52
• And if it helps clear it up, what I am trying to do could be summed up as "introducing a mixed-effect regression among the parameters in a pre-existing optimization problem." – ontologist Nov 7 '16 at 5:07
• I don't have access to the paper, and don't understand what you're trying to do. It's not that I can't make guesses about that, however. I have no idea what your original optimization problem is. Saying that it is solvable using sequential least squares doesn't tell me what it is, just that you have a particular method (whose merits I have no clue about, absent knowing what problem you're trying to solve) to solve it. Your exposition is at far too vague/wishy washy a level. – Mark L. Stone Nov 7 '16 at 12:33
• Fair enough. I added more information to the original post. – ontologist Nov 7 '16 at 15:25
• Please explain how the "regression-like relationships among the parameter values $r$" relate to the original model, which already has some difference equation constraints involving $r$. Will the regression "constraints" replace or augment the existing constraints involving $r$? – Mark L. Stone Nov 7 '16 at 17:35

Q1. I have no idea what's been done on this problem. I don't even have access to the paper you linked.

Q2. It sounds like you want to combine maximum likelihood and best fit of the rest of the model. In such case, the critical matter is to determine how to combine these 2 "sub-objectives" into one overall objective function, or otherwise trad them off.

One way of doing this is to add a monotonically decreasing function, $f$, of likelihood, $L$, to a (minimized) "best fit" objective, $g$. That function $f(L)$ could be $-m * log(L)$, or it could be $f(L) = -m * L$, in both cases for some non-negative value $m$ which you will have to choose. Or it could be something else. Unlike the situation in which only likelihood is maximized, i.e., negative likelihood is minimized, the choice of function $f$ can and generally will affect the optimal solution. You will also need to include any constraints in the parameters being estimated, such as non-negativity constraints on estimated variance, or semidefinite constraint on covariance matrix (if estimating correlated Multivariate Normal), i.e., constraint that covariance matrix being estimated is symmetric positive seni-definite (yeah, I'm throwing you in the water on that, but there are nonlinear objective semidefinite optimization solvers to handle such problems).

Alternatively, you can map out an efficient frontier by either:

a. For each of several values of lower bound on Likelihood, minimize $g$ subject to all constraints plus constraint consisting of specified lower bound on likelihood

or

b. For each of several values of upper bound on loss function from fit, minimize negative log likelihood (or negative likelihood) subject to all constraints plus constraint consisting of upper bound on loss function from fit

Q3. For hypothesis testing, bootstrap the whole procedure. Yes, that can be a lot of computing.