I have a data set, composed of points $(x_i, y_i)$ for $i=1,N$. I also have a known function $F$, which maps these points $x_i$ to $y_i$ as such $F(x_i, a(x_i),b(x_i)) = y_i$, where $a(x_i)$ and $b(x_i)$ are unknown functions of $x_i$. What I'd like to do is run an optimization (ideally using Python packages) in order to determine the functions $a$ and $b$ such that the error is minimized and also such that the functions $a$ and $b$ are smooth.

One, possibly naive way of doing this is to bin my data with respect to the $x_i$ variable, and then solve the optimization problem $F(x_i, a,b) = y_i$ for the single values $a,b$ for all points within the bin. Then, using this discretized solution for $a,b$, perform a numerical fit which results in a smooth curve.

However, I'm wondering if this is the best way to go about it. Is there a method that I could use which would determine the $a,b$ values simultaneously while trying to create a smooth curve from them? Or even if there were, would the effect be the same as if I perform the naive binned-optimization as I've described above?


Your problem is ill posed in the sense that if a solution exists, it is not unique. To see this, let ($a_i^*$, $b_i^*$), denote a solution of the equation $F(x_i, a_i,b_i)=y_i$, $i\in\{1,\ldots,N\}$ (or the minimizers of an appropriate cost function if equality cannot be achieved). Then, any pair of functions $\{a^*(x)$, $b^*(x)\}$ such that $a^*(x_i)=a_i^*$, $b^*(x_i)=b_i^*$, for all $i\in\{1,\ldots,N\}$, is a solution to your problem. Clearly, the set of all possible solutions is uncountable, even with the additional requirement that they are 'smooth'.

In order to avoid this ambiguity in the solution, one approach is to impose additional conditions. For example, you could restrict the solution functions to belong to a finite dimensional (and 'smooth') space (e.g., the space of polynomials of order $n$). In that case, the (difficult) functional optimization problem is converted to the (much simpler) problem of determination of a finite set of variables.

  • $\begingroup$ Yes, thank you for the clarity in what I'm looking for. Certainly I have no constraints on any aspect of the function outside the range of x_i, so I cannot have a unique solution. While I'm not yet sure what sort of basis set I should use to represent the functions a and b, that is a topic I know how to explore. Thanks again! $\endgroup$ – gammapoint Nov 22 '16 at 21:43

You could discretize the unknown functions $a$ and $b$, then formulate this as a least squares optimization problem with smoothing regularization. Specifically:

$$\min_{\mathbf{a}, \mathbf{b}} \frac{1}{2}\sum_{i=1}^N\left( f_i - F(x_i, \mathbf{m}_i^T \mathbf{a}, \mathbf{m}_i^T \mathbf{b}) \right)^2 + \frac{\alpha}{2} \mathbf{a}^T R \mathbf{a} + \frac{\alpha}{2} \mathbf{b}^T R \mathbf{b}$$


  • $\mathbf{a}$ and $\mathbf{b}$ are coefficients of $a$ and $b$ with respect to some discrete set of basis functions $\{\phi_j\}_{j=1}^m$ (say, a piecewise linear basis of "hat" functions), $$a(x) = \sum_{j=1}^m \mathbf{a}_j \phi_j (x)$$
  • $f_i$ is the $i$'th piece of data you are trying to fit,
  • $\mathbf{m}_i^T := \begin{bmatrix}\phi(x_1) & \phi(x_2) & \dots & \phi(x_N)\end{bmatrix}$ is the "measure at $x_i$ vector" such that $$\mathbf{m}_i^T \mathbf{a} = a(x_i)$$
  • R is a discrete form of the Laplacian with Neumann boundary conditions, plus some small multiple $\tau$ of the identity: $$R := \Delta_N + \tau I.$$

This could then be solved with the Gauss-Newton method, or L-BFGS, or something like that, (it is straightforward enough to figure out the derivative of the objective function by hand). In Python you could do the discretization with FEniCS, and the optimization with Scipy, though it will involve knowing a little about finite elements and a little about optimization.

  • $\begingroup$ Thank you for your suggestion, Nick. Modifying the objective function to explicitly have a smoothing term is something I thought about vaguely, but didn't know exactly the best way to approach that. Your suggestion is useful in that regard, though I think I may go with a naive binning as I've stated in my question or with a particular basis set for the a and b functions as Stelios suggested. If those don't work for my application as needed, I'll revisit your thoughts. $\endgroup$ – gammapoint Nov 22 '16 at 21:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.