## 2 Materials and Methods  
### A New Approach: Continued Fractions with Splines  
In previous contributions [8; 9], a memetic algorithm was employed to find the approximations. Here, we present another method to fit continued fraction representations by iteratively fitting splines. Splines provide a regression technique that involves fitting piecewise polynomial functions to the given data [15]. The domain is partitioned into intervals at locations known as \"knots\". Then, a polynomial model of degree \(n\) is separately fitted for each interval while generally enforcing boundary conditions, including the continuity of the function and continuity of the first \((n-1)\)-order derivatives at each of the knots. Splines can be represented as a linear combination of basis functions, of which the standard is the B-spline basis. Thus, fitting a spline model is equivalent to fitting a linear model of basis functions. We refer to [PERSON] et al. [16] for the particular definition of the B-spline basis. First, when all the functions \(b_{i}(\mathbf{x})=1\), we have a _simple continued fraction_ representation for all \(i\), which we can write as  
\[f(\mathbf{x})=g_{0}(\mathbf{x})+\frac{1}{g_{1}(\mathbf{x})+\frac{1}{g_{2}( \mathbf{x})+\frac{1}{g_{3}(\mathbf{x})+...}}}. \tag{5}\]  
Note that for a term \(g_{i}(\mathbf{x})\) we say that it is at a \"depth\" of \(i\). Finding the best values for the coefficients in the set of functions \(\{g_{i}(\mathbf{x})\}\) can be addressed as a nonlinear optimization problem as in [8; 9]. However, despite the great performance of that approach we aim to introduce a faster variant that can scale well to larger datasets such as this one. Towards that end, and thinking about the scalability, we fit the model iteratively by depth as follows: we first consider only the first term \(g_{0}(\mathbf{x})\) at depth 0, ignoring all other terms. We fit a model for the first term using the predictors \(\mathbf{x}\) and target \(f(x)\). Next, we consider only the first and second depths, with the terms \(g_{0}(\mathbf{x})\) and \(g_{1}(\mathbf{x})\), ignoring the rest. We then fit \(g_{1}(\mathbf{x})\) using the previously fit model for \(g_{0}(\mathbf{x})\). For example, truncating the expansion at depth 1, we have  
\[g_{1}(\mathbf{x})=\frac{1}{f(\mathbf{x})-g_{0}(\mathbf{x})}. \tag{6}\]  
Thus, we fit \(g_{1}(\mathbf{x})\) using the predictors \(\mathbf{x}\) and the target \((f(\mathbf{x})-g_{0}(\mathbf{x}))^{-1}\). We label this target as \(y^{(1)}\). We repeat this process, fitting a new model by truncating at the next depth by using the models fit from previous depths and iterations. We have that at depth \(i>0\), the target \(y^{(i)}\) for the model \(g_{i}(\mathbf{x})\) is \((\epsilon_{i-1})^{-1}(\mathbf{x})\), where \(\epsilon_{i-1}(\mathbf{x})\) is the residual of the previous depth's model, \(y^{(i-1)}-g_{i-1}(\mathbf{x})\). One notable characteristic of this approach is that if any model \(g_{i}(\mathbf{x})\), \(i>0\) evaluates to 0, then we will have a pole in the continued fraction, which is often spurious. To remedy this, we modify the structure of the fraction such that each fitted \(g_{i}(\mathbf{x})\), \(i>0\) is encouraged to be strictly positive on the domain of the training data. To do this, we add a constant \(C_{i}\) to \(\epsilon_{i}\) when calculating the target \(y^{(i+1)}\), where \(C_{i}=|\min_{x}\epsilon_{i}|\). Thus, the targets \(y^{(i)}\) for \(i>0\) are all non-negative, encouraging each \(g_{i}(\mathbf{x})\), \(i>0\), to be strictly positive. For example, for \(g_{1}(\mathbf{x})\), we would have that the target \(y^{(1)}=(f(\mathbf{x})-g_{0}(\mathbf{x})+C_{1})^{-1}\). Of course, we must then subtract \(C_{i}\) from \(g_{i-1}(\mathbf{x})\) in the final continued fraction model. We have found that data normalization often results in a better fit using this approach. It is sufficient to simply divide the targets uniformly by a constant when training and multiply by the same constant for prediction. We denote this constant parameter norm.