Novel Optimization Techniques for Parameter Estimation

In a parameter estimation problem, one has a function that maps inputs and a parameter set to a desired output. Then, for a given dataset one is interested in finding the best fitting (in some sense) parameter set. How one defines ‘best fitting’ is a choice that the modeler makes based on their understanding of the dataset and its generation. In particular, assume we have a function $f$ that takes inputs $x$ and parameter set $\theta$ and generates a prediction of a desired output, $y=f(x;\theta)$. If we have a dataset of paired inputs and outputs, $(x_{i},y_{i})_{i=1}^{n}$, and assume a simple statistical model generating the outputs, $y_{i}=f(x_{i};\theta)+Z_{i}$, where $Z_{i}$ are independent and identically distributed random variables with probability density function $\phi$ that are used to represent potential observation noise, we can then write the negative log-likelihood of the parameter set given the observed data as

To complete the parameter estimation problem, we must minimize the function $L\left(\theta|(x_{i},y_{i})_{i=1}^{n}\right)$ over $\theta$; however, as mentioned before, the components of $\theta$ might have a physical interpretation and therefore be constrained. For example, they might be constrained to be non-negative, or bounded, or, as stated in the introduction, they might be proportions that sum to one. As a result, we have reduced the parameter estimation problem to a constrained optimization problem. Note that when writing the negative log-likelihood we will often drop the explicit dependence on the observed data.

Our focus in this work is creating algorithms to numerically solve a broad class of constrained optimization problems that arise from parameter estimation problems of the form

where the objective function $L:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is possibly non-convex and the domain is the intersection of a linear and box constraint. Observe, all the examples given for constraints on the parameters can be represented by the conditions in (1), and we can actually further generalize our optimization framework by noting the box constraint, $l\leq\theta\leq u$, can be rewritten as a conic constraint. Letting $\theta_{1}:=\theta-l$ and $\theta_{2}:=u-\theta$ we see the condition $\theta_{1},\;\theta_{2}\geq 0$ and $\theta_{1}+\theta_{2}=u-l$ is equivalent to the box constraint. So, we can replace the box constraint with an additional linear equality constraint and restrictions to the non-negative orient, i.e., $\mathbb{R}^{n}_{+}:=\{x=(x_{1},\ldots,x_{n})\;|\>x_{i}\geq 0,\;i=1,\ldots,n\}$.

The set $\mathbb{R}^{n}_{+}$ is an example of a pointed convex cone, that is, a set $\mathcal{K}\subseteq\mathbb{R}^{n}$ such that for all $x,y\in\mathcal{K}$ and $t\geq 0$: $tx\in\mathcal{K}$, $x+y\in\mathcal{K}$, and $-\mathcal{K}\cap\mathcal{K}=\{0\}$. If additionally the span of $\mathcal{K}$ is the entire ambient space we say the cone is solid. Thus, due to this equivalence, we shall instead develop our algorithm to solve the more general class of minimization problems

where $\mathcal{K}\subseteq\mathbb{R}^{n}$ is a pointed convex solid cone and $L:\mathcal{K}\rightarrow\mathbb{R}$ is smooth but possibly non-convex. For a first-order method, one typically expects to find a first-order solution, which is defined as $\nabla L(\theta)=0$ without constraints. However, like the examples in Sections 4 and 5, functions may have many saddle points which satisfy first-order conditions while failing to be local minimums, and second-order methods are ideal to avoid converging to such saddle points. In the optimization literature, one of the major recent developments is the introduction of second-order methods which converge to second-order stationary points [CRS19, OW21, ROW20], but relatively few papers in the literature deal with developing second-order methods for non-convex constrained models. Independent of our work, in a very recent arXiv paper, Dvurechensky and Staudigl [DS24] proposed a similar approach to ours to solve non-convex constrained optimization models. The main difference between our work and theirs lies in the fact that we use the Dikin ellipsoid and a cubic regularization of the objective in defining our subproblems, while they use only a cubic regularization of the objective to form their subproblems.

In this paper, we use the definitions for first-order (FOSP) and second-order stationary points (SOSP) as introduced in [HL23]. The authors of [HL23] also introduced a so-called Newton-conjugate gradient (Newton-CG) based barrier method to find an $(\epsilon,\sqrt{\epsilon})$-SOSP for non-convex conically constrained models with $O(\epsilon^{-3/2})$ Cholesky factorizations and $\tilde{O}(\epsilon^{-3/2}\min\{n,\epsilon^{-1/4}\})$ other fundamental operations. In comparison, the authors of [DS24] apply a cubic regularized model to get a $O(\epsilon^{-3/2})$ worst-case iteration bound. CRNAS also obtains a $O(\epsilon^{-3/2})$ worst-case iteration bound to $\epsilon$-SOSP for (2); CRNAS is simple to implement and is well suited for the situation when the Hessian is easily obtainable. Additionally, our approach avoids the need to solve a possibly ill-conditioned Newton’s equation; we only require the solution to a simple subproblem which we discuss in Appendix A.2. Next, we introduce the notion of self-concordant barrier functions for conic constraints, which play an important role in our analysis, and the corresponding notion of $\epsilon$-FOSP and $\epsilon$-SOSP arising from them.

Logarithmic homogeneous self-concordant barrier functions [NN94] play an important role in interior-point methods. In this paper, we assume the cone $\mathcal{K}$ is convex pointed and solid and has a logarithmic homogeneous self-concordant barrier function $B:{\rm int}(\mathcal{K})\to\mathbb{R}$, that is, $B$ is convex, three-times continuously differentiable over ${\rm int}(\mathcal{K})$, $B(\theta_{k})\rightarrow\infty$ for all sequences $\{\theta_{k}\}_{k\in\mathbb{N}}\subseteq{\rm int}(\mathcal{K})$ which converge to a point on the boundary of $\mathcal{K}$, and for any $\theta\in{\rm int}(\mathcal{K})$, $\tau>0$, and direction $h\in\mathcal{R}^{n}$, the following properties are satisfied:

where $D$ is a positive constant. As an example, $B(\theta)=-\sum_{i=1}^{n}\log(\theta_{i})$ is a logarithmic homogeneous self-concordant barrier function for the cone $\mathbb{R}^{n}_{+}$ with $D=1$.

An important aspect of barrier functions comes from the fact they can define induced local norms. These norms are standard in the conic optimization literature and play a crucial role in our algorithm development and analysis. Following convention, we let

for all $v\in\mathbb{R}^{n}$ and $C\in\mathbb{R}^{n\times n}$. Note, unless otherwise stated, $\|\cdot\|$ denotes the regular Euclidean norm. Following the pattern of [HL23], a feasible solution $\hat{\theta}\in\mbox{\rm int}(\mathcal{K})$ with $A\hat{\theta}=b$ is said to be an $\epsilon$-approximate first-order stationary point ($\epsilon$-FOSP) for (2) if

which means that the angle between $(\nabla^{2}B(\hat{\theta}))^{-1/2}\nabla L(\hat{\theta})$ and the orthogonal complement subspace of $(\nabla^{2}B(\hat{\theta}))^{-1/2}\mbox{\rm Range}(A^{\top})$, namely $(\nabla^{2}B(\hat{\theta}))^{1/2}\mbox{\rm Null}(A)$, is of order $\epsilon$. The latter statement can be restated as in terms of the local norms as

for all $Ad=0$. Moreover, a solution $\hat{\theta}$ is called an $\epsilon$-approximate second-order stationary point ($\epsilon$-SOSP) if in addition to being an $\epsilon$-FOSP the point $\hat{\theta}$ also satisfies

for all $Ad=0$; see Remark 1 (ii) of [HL23] for further details. Next, we introduce the two methods which inspired our algorithm design.

The cubic regularized Newton’s method was first proposed by Nesterov and Polyak to solve smooth unconstrained optimization problems [NP06]. Their approach adds a cubic regularization term to the second-order Taylor expansion of the objective function about the current iterate, $\theta^{k}$, and solves this subproblem to compute the next iterate, i.e.,

If the objective function has a gradient Lipschitz Hessian with constant $L_{H}\geq 0$ and $M\geq L_{H}$, then the above subproblem amounts to minimizing a cubic upper bound of the objective function about $\theta^{k}$. Hence, under mild assumptions, cubic regularized Newton monotonically decreases the value of the objective function.

This method boasts multiple desirable qualities. The algorithm converges globally to second-order stationary points and computes a point satisfying the approximate first-order unconstrained optimally conditions, i.e., $\|\nabla L(\theta)\|\leq\epsilon$, within $O(\epsilon^{-3/2})$ iterations, which bests gradient descent, and the method converges quadratically near strict local minimums [NP06]. Furthermore, though (4) is generally a non-convex problem, the subproblem is actually equivalent to minimizing a convex function in one variable. Thus, cubic regularized Newton is an implementable, globally convergent, and locally fast converging method to second-order stationary points. For these reasons, we seek to incorporate aspects of this procedure into our design.

The affine scaling method [LV90, Bar86, VMF86] is an algorithm for solving linear programming problems by rescaling the variables and constraining the next iterate to lie within a ball contained inside the cone $\mathbb{R}^{n}_{+}$ to get a better iterate at a reduced computational cost. Consider the linear program

In each iteration, we scale the problem based on the current point $\theta^{k}>0$. Let $\theta^{\prime}=D_{k}^{-1}(\theta-\theta^{k})=D_{k}^{-1}\theta-\textbf{1}_{n}$, where $D_{k}=\mbox{\rm Diag}(\theta^{k})$ is a diagonal matrix with diagonal elements $\theta^{k}$ and $\textbf{1}_{n}$ is the vector of all ones in $\mathbb{R}^{n}$. Then, the non-negativity constraint is equivalent to $\theta^{\prime}\geq-\textbf{1}_{n}$. We desire to have the next iterate stay inside the interior of $\mathbb{R}^{n}_{+}$, so we replace $\theta^{\prime}\geq-\textbf{1}_{n}$ with a ball constraint, that is, we instead enforce $\|\theta^{\prime}\|\leq 1-\alpha$, where $0<\alpha<1$. This then ensures the next iterate cannot lie on the boundary of the cone with the proximity to the boundary dictated by $\alpha$, i.e., $\alpha$ close to zero can yield new iterates near the boundary. With this new constraint and change-of-variable, we form the following subproblem

The main benefit of (5) is it has a closed formed solution. So, the affine scaling method proceeds by forming and solving (5) and using the variable transformation to obtain the next iterate. Utilizing the notation introduced in Section 2.2, we can then concisely write the affine scaling update as

where $B(\theta)=-\sum_{i=1}^{n}\log(\theta_{i})$ is the barrier function used to define the norm in the constraint.

Our new method seeks to combine the ideas of cubic regularized Newton and affine scaling to develop a procedure which solves (2). The crux of our approach is to utilize the affine scaling method to handle the constraints, leveraging the concept of homogeneous self-concordant barrier functions to deal with the general conic constraints, but replace the linear objective in (6) with the local cubic approximation of the objective function in (4). Thus, our method generates new iterates by solving the following subproblem

where $M$ is a positive number. Since our method brings together these two different procedures we coined it cubic regularized Newton with affine scaling (CRNAS). A precise description of CRNAS is provided below.

Cubic Regularized Newton with Affine Scaling (CRNAS) Step 0: Provide an interior point $\theta^{0}$, i.e., $A\theta^{0}=b$ and $\theta^{0}\in\mbox{\rm int}(\mathcal{K})$; choose the constants $\eta>0$, $M>0$, and $\alpha\in(0,1)$; set $k=0$ Step 1: Solve the following subproblem and proceed to Step 2: $\displaystyle\theta^{k+1}$ $\displaystyle=$ $\displaystyle\arg\min_{\theta:A\theta=b,\|\theta-\theta^{k}\|_{\theta^{k}}\leq 1-\alpha}\left(\left\langle\nabla L\left(\theta^{k}\right),\theta-\theta^{k}\right\rangle+\frac{1}{2}\nabla^{2}L\left(\theta^{k}\right)[\theta-\theta^{k}]^{2}+\frac{M}{6}\left\|\theta-\theta^{k}\right\|_{\theta^{k}}^{3}\right)$ Step 2: If $\|\theta^{k+1}-\theta^{k}\|_{\theta^{k}}<\eta$, let $K=k+1$ and stop. Otherwise, go back to Step 1 with $k=k+1$

Though (7) appears to be more difficult to solve than (4), this subproblem forming the backbone of CRNAS is solvable. A theoretical analysis and practical approach to solving (7) is detailed in Appendix A.2. In the next section, we provide an overview of the convergence guarantees we have for CRNAS; the technical proofs of the results are left to Appendix A.1.

We first declare CRNAS is well-defined in the sense all of the iterates produced by the algorithm remain inside the cone $\mathcal{K}$. The following lemma guarantees this follows from the constraints in the subproblem.

For our convergence theory to hold, we assume the following scaled Lipschitz smoothness of the second-order derivative of the objective:

Therefore, in view of the discussions on the optimality conditions in Subsection 2.2, Theorem 2.2 implies that in at most $O(\epsilon^{-3/2})$ iterations CRNAS is guaranteed to find an $\epsilon$-SOSP. So, for the highly nonlinear and non-convex optimization problems coming from parameter estimation problems, CRNAS is able to avoid the many sub-optimal first-order stationary points leading to a proposed solution of often superior quality compared to first-order methods.

It is worth remarking one can devise a first-order version of CRNAS which uses a quadratic approximation of the objective function in the subproblem rather than a cubic approximation. A similar convergence result can be derived for this procedure as well. As situations arise where second-order derivatives are prohibitively expensive to compute, this version of CRNAS would be prudent to implement; therefore, for the sake of completeness, we include a description of our first-order version of CRNAS with associated convergence theory in Appendix A.3

The Hill equation is a fundamental model for describing dose-response relationships in pharmacological studies [GZKB⁺12]. For a given drug dose level $d$, the Hill equation specifies the fraction of viable cells as

where the parameters $b$, $E$, and $n$ represent the maximum drug effect, the half maximal effect dose (EC50), and the Hill coefficient respectively. This model has been well utilized in pharmacology, biochemistry, and various other fields over the past hundred years [Hil10].

An important question that arises in the study of the Hill equation is the accurate identification of model parameters based on observed data [GC15]. A particularly challenging aspect of dealing with the parameters of the Hill equation is the nonlinear term in the denominator, $E^{n}$. To mitigate the challenge of this term, we apply a variable transformation. In particular, we introduce a new parameter $\mathcal{E}=E^{n}$ to produce the modified Hill equation

Given the transformation $\mathcal{E}=E^{n}$ is a one-to-one function when $E$ is non-negative, we can determine the value of $E$ and $n$ from the estimated values of $\mathcal{E}$ and $n$. This parameter transformation is applied during the numerical experiments, while the original EC50 parameter $E$ is used throughout the manuscript for clarity.

The logistic growth model is popular for studying population growth with a carrying capacity. Recently, there has been significant interest in an extended version of the logistic growth model, which aims to represent cell populations consisting of heterogeneous subpopulations [JMS18]. These subpopulations can exhibit varying growth dynamics, such as differing growth rates and carrying capacities. For simplicity here we study a linear combination of logistic growth functions with varying parameter sets for each subpopulation to model a mixture of heterogeneous populations growing according to a logistic function. Note, we do not assume the populations are interacting via their carrying capacity as that would necessitate the numerical solution of a nonlinear system of differential equations which is outside the focus of this manuscript.

In particular, we start with the parameters of $S$ subpopulations: $\theta_{LG}(S)=(p_{i},\alpha_{i},\beta_{i}),i=1,\cdots,S$, and we then model the total cell count dynamics with respect to time as

where $F_{LG}(0)$ is the initial total cell count and $p_{i}$ is the initial proportion of each subpopulation, satisfying an equality constraint $\sum_{i=1}^{S}p_{i}=1$. Compared to the widely used three-parameter logistic growth model, we normalize the maximum carrying capacity for each subpopulation to 1, ensuring the identifiability of each parameter. We retain another two parameters, $(\alpha_{i},\beta_{i})$, to capture the growth behavior of type $i$ population at the inflection point. The parameter $\alpha_{i}$ encodes the maximum growth rate that the type $i$ population can achieve around the inflection point, while the parameter $\beta_{i}$ relates to the shift in time of the inflection point. From a modeling perspective, there are no restrictions on $\alpha_{i}$ and $\beta_{i}$. However, we assume $\alpha_{i}\geq 0$ for all $i=1,\cdots,S$ for a growing population and select a distinct $\beta_{i}$ for each subpopulation to ensure these parameters are identifiable. The detailed selection method is specified in Appendix A.6. For this experiment, we assume a zero observation noise. Thus, for a dataset observed at a given set of time points $\mathcal{T}$, $\mathcal{X}=\{x_{t};t\in\mathcal{T}\}$, we solve the least squares problem

to estimate the true parameter set $\theta^{*}_{LG}(S)$.

We compared the performance of all the algorithms in Table 1 for solving the least squares problem described in (13). Similar to the previous case study, we compared the methods based on the described metrics in Section 3. Experimental details are provided in Appendix A.6.

From Figure 7, we see again CRNAS boasted the smallest number of iterations required to obtain the optimal solution. Furthermore, the middle panel in Figure 7 shows CRNAS yielded the lowest best objective values among all five methods. Based on the advantages observed in both the number iterations and objective value, we see CRNAS outperformed four different implementations of state-of-the-art algorithms in solving this highly nonlinear and non-convex parameter estimation problem.