Root finding (previously: "Constant class")

[egpbos]

Following up on the discussion in #254, it might be nice to have a Constant class next to the Variable and Parameter classes. Consider the following example:

import numpy as np
import symfit as sf

ZERO = np.zeros(20)
X = np.linspace(0, 1, 20)
Y = 1 + X * 0.5 + np.random.normal(0, 0.05, 20)

a, b = sf.parameters('a, b')
x, y, zero = sf.variables('x, y, zero')

model = {zero: a + b*x - y}

fit = sf.Fit(model, x=X, y=Y, zero=ZERO)
res = fit.execute()

This will give a divide-by-zero warning:

.../symfit/core/fit_results.py:277: RuntimeWarning: divide by zero encountered in double_scalars
  return 1 - SS_res/SS_tot

The reason is that the SS_tot term in the coefficient of determination is proportional to the variance of the data of the variable that is being fit. In this case, that variable is zero. However, this isn't conceptually a variable at all, and indeed, the variance in the "data" is zero, because it is not actually data, but rather a work-around for the fact that constants are not present.

Note that in the example of this issue (y = a + bx), the model could have easily been written without the zero, but there are other cases where this doesn't work (see e.g. #254) and a constant term is really necessary.

If we could define zero to be a Constant, and still use it in place of Variables, this would increase the modeling/parameterization flexibility. Also, it would remove the current need to create mock data arrays for such values, rather only having to pass in actual data. And, finally, it would make the warning go away :)

Would this be feasible? How much of the existing codebase would be touched by this?

[tBuLi]

I'm not sure I would use Constants for this. I don't know actually how ROOT defines Constants so it may be a solution, but the way I see it we turned a fitting problem into a root-finding problem. So I think it would be nice if we can somehow indicate this change in the nature of the problem.

Currently, with the solution I proposed for #254, we are finding the roots by finding the least squares solution to f_i = 0. Although for all intents and purposes this is fine, it might not be the best way to go about root finding. So I think it would be better to have a separate objective function and/or minimizer for root finding.

We should be careful not to end up in a recursive jungle here though, those root finders might use the least squares solution and a normal minimizer. But even if that is the case, a syntax like

a, b, c, d = sf.parameters('a, b, c, d')
x, y, z, f = sf.variables('x, y, z, f')
plane_model = {f: a * x + b * y + c * z - d}

fit = Fit(plane_model, x=xdata, y=ydata, z=zdata)

might be enough for symfit to infer that root finding is desired.

[egpbos]

That indeed sounds like a superior alternative, especially given that you can just use constant values in models already, so indeed the only special case is zero on one side of the equation.

Perhaps it would be slightly more expressive to add an explicit sf.zero singleton object, so that you could do Fit(plane_model, x=xdata, y=ydata, z=zdata, f=sf.zero). Or to replace Fit with FindRoot. But, yes, it would be nice if root finding were treated as such!

[tBuLi]

I think this might already exist since we have everything sympy has: sf.Zero. So that might indeed be a good solution. Another option is to use a normal int 0, but that is less clear. To use a symbolical Zero is a lot better.

Concerning the internals: we would still have to turn this into the exact same problem, e.g. solve f_i = 0. It will just give a convenient API for doing so.

[egpbos]

It will just give a convenient API for doing so.

Indeed, that would be a great benefit I think, it would make it consistent with the general great ease of use that the rest SymFit already has.

Concerning the internals: we would still have to turn this into the exact same problem, e.g. solve f_i = 0.

Except for the implementation then? Indeed, as a first working version, you could just do it like in #254, but you'd like to also be able to use the SciPy root finders, maybe. I'm not an expert on the differences between all the algorithms, though. As you say, they may overlap.

[tBuLi]

Yeah I agree, some of the root finders might offer benefits over the #254 implementation, it is definitely worth investigating which ones, and wrapping them separately.

And this would definitely be an API improvement, since this kind of problem is quite common.