RAFF . jl : Robust Algebraic Fitting Function in Julia

Let f : R → R be a function whose mathematical description is not available. This function can be, for example, a black-box, a proprietary computer program or an experiment. Suppose that a dataset S = {(x1, y1), . . . , (xm, ym)} is available, where yi is an approximation of f(xi) (from an experimental procedure, numerical approximation, etc.) and we want to approximate f by a known model φ. Model φ can be defined as φ(x, θ), where x are the n independent variables of f and θ represents some parameters of φ. RAFF.jl (Robust Algebraic Fitting Function) is a Julia package developed to find parameters θ for φ in order to adjust it to the observed values S of the unknown function f . Following Liu & Wang (2008) and Keleş (2018), in general, the adjustment can be related to

1. Classical least squares (algebraic fit): which considers the sum of deviations of type |ϕ(x i , θ) − y i | 2 , also known as regression; 2. Orthogonal least squares (geometric fit): which considers the sum of deviations of type min x ∥(x, ϕ(x, θ)) − (x i , y i )∥ 2 (orthogonal projection on the curve to be adjusted).
RAFF.jl was developed to solve a generalization of the first case.
Linear and nonlinear regression is essentially the adjustment of mathematical functions to data and is a problem that appears in many areas of science.When data comes from real experiments, non-expected errors may cause the appearance of outliers, which might be responsible for causing the regression calculated by sum of deviations to result in misleading approximations.Regression is strongly connected to Statistics but practical methods to detect outliers are not very common.Motulsky & Brown (2006), for example, develop a method for outlier detection based on the assumption that the error follows a Lorentzian distribution around the function and use nonlinear regression based on least squares.RAFF.jl provides automatic detection of outliers using a voting system.It is an optimization-based package, based on algorithms for Lower Order-Value Optimization (LOVO) which were introduced by Andreani, Dunder, & Martıńez (2005) and revisited by Andreani, Martıńez, Martıńez, & Yano (2009).
Recently, Martıńez (2012) performed a complete review about LOVO problems considering theoretical aspects of algorithms to solve it and potential applications.

Background
To elucidate the essence of how RAFF.jl works, let us detail some aspects related to the LOVO problem and its resolution.Let us consider m functions F i : R n → R, i = 1, ..., m.
Given θ ∈ R n , we can sort the set {F i (θ), i = 1, ..., m} in ascending order: Considering a value 1 ≤ p ≤ m, we can define the LOVO function as and the LOVO problem as min θ∈R n S p (θ).
Assuming that F i , i = 1, ..., m are continuous functions we have that S p is a continuous function, but assuming that F i 's are differentiable functions we cannot conclude that S p is differentiable.This can be seen by reformulating the LOVO problem as follows.Denoting C = {C 1 , ..., C r } as the set of all combinations of {1, ..., m} taken p at time, we can define for each i ∈ {1, ..., r} the following function It can be observed that, for a given θ, f min (θ) is obtained by a combination C j which contains the smallest sum of p elements of the set {F i (θ), i = 1, ..., m}.Therefore f min (θ) = S p (θ) and, consequently, the LOVO function is non differentiable.The LOVO problem description can become more clear by considering an example.In this sense, let us consider the dataset given by x y -0.5 0.119447 0.0 0.3 0.5 0.203551 0.75 0.423998 and the model defined by ϕ(x, θ) = θ(sin(x) + cos(x)).Naturally, we have m = 4 and let us consider p = 3.The F i 's can assume diferent forms.To leave the example closest to our approach, let's consider F i 's as: Since m = 4 and p = 3, we have 4 possible subsets with 3 elements each from set {1, 2, 3, 4}: Thus, associated to each C i , i = 1, ..., 4, we can define function f i as follows and consequently, As previously pointed out, this function is continuous but it is not differentiable as illustrated in Figure 2. Andreani al. (2009) introduced line search methods and handled the possible singularities in a clever way, using the following approximation for ∇f min (θ) where i ∈ I(θ) = {k ∈ {1, ..., r}; f k (θ) = f min (θ)}.This approach can naturally be extended for second order derivatives.
An important point for practical purposes is when we consider the LOVO problem with p = m and F i (θ) = (ϕ(x i , θ) − y i ) 2 .In this case, the LOVO problem coincides with classical least squares and, consequently, it can be seen as a generalization of the least squares problem.When p < m and F i (θ) = (ϕ(x i , θ) − y i ) 2 , the solution θ provides a model ϕ(x, θ) free from influence of the m − p points with the highest deviation.The number p can be interpreted as the number of trusted points, that is, m − p possible outliers were identified.
One of the most usual ways to solve the problem of nonlinear least squares is by using the Levenberg-Marquardt method (Moré, 1978).This method is a first-order method, where derivatives of the model ϕ with respect to θ are used to compute the gradient of the objective function in the associated least squares problem.The reason for the wide use of Levenberg-Marquardt method is, in general, associated with quadratic convergence properties even using only first-order derivatives.In this direction, it is relevant to ask about Levenberg-Marquardtbased methods to solve LOVO problems in the context of adjustment functions.
RAFF.jl implements a Levenberg-Marquardt algorithm in the context of LOVO problems, i.e., it solves the problem of minimizing f min (θ), where F i (θ) = (ϕ(x i , θ) − y i ) 2 , for i = 1, . . ., m.In this sense, first-order derivatives are necessary and the same strategy of Andreani et al. (2009) is used.It uses first-order derivatives of the model ϕ with respect to θ to approximate the gradient of f min (θ), which is a non differentiable function.Moreover, LOVO problems have the limitation that the number p of possible trusted points needs to be given by the user.RAFF.jl solves this limitation by implementing a voting system.In this voting system, several LOVO subproblems are solved with different values for p, the number of possible trusted points.Each solution of a LOVO subproblem is associated to a vector parameter θ.The vector parameters are compared against each other using the Euclidean distance, where small distances (using a threshold) are considered the same solution.The parameter θ * which most occurs among them is declared as the solution.

Functionality
RAFF.jl main methods expect as input a dataset of the observed data and a model function, whose parameters one intends to adjust.The model function is a regular Julia function with 2 arguments: θ represents the parameters of the model and x represents the arguments of function f .The following function is an example of a model representing the logistic function .
The observed data can be represented by the following In this example, the true function was given by f (x) = 1000 + 5000 1.0 + exp (−0.2x + 3) .
The observed data was generated as random normal perturbations around the graphic of f and is shown in Figure 1.The dots and triangles represent the observed data, where the red triangles were manually set to be the outliers.Using the least squares technique with the model above, the green function is found.When RAFF.jl is applied to the same problem, it correctly identifies the two outliers.The resulting function is depicted as the red one, very close to f .

Additional features
The user may also provide more information to RAFF.jl, such as an rough approximation to the expected number of trusted observations.Additional methods and options are also available to more advanced users, such as generation of random test data and multistart strategies.First-order derivatives of the model ϕ with respect to θ can also be provided, which results in a faster executing time.When they are not provided by the user, RAFF.jl uses Julia's ForwardDiff.jlpackage (Revels, Lubin, & Papamarkou, 2016).
RAFF.jl can be run in serial, parallel and distributed environments.Parallel and distributed methods use the native Distributed.jlpackage.The distributed version is a primary-worker implementation that does not use shared arrays, therefore, can be run both locally or on a cluster of computers.
This package is intended to be used by any experimental researcher with a little knowledge about mathematical modeling and fitting functions.

Installation and usage
RAFF.jl is an open-source software that can be downloaded from Github.It is a registered package and can be directly installed from Julia's package repository.The whole description for first time usage or its API is available at its documentation.

Figure 1 :
Figure 1: The red function represents the LOVO function.Observing the interval [0.2, 0.25] we can note a singular point even considering f1, f2, f3 and f4 as differentiable functions.

Figure 2 :
Figure 2: Points representing the logistic function.The red triangles are two outliers that should be ignored.The blue dashed function is the true one, while the green was obtained by traditional least squares techniques and the red one was obtained by RAFF.jl.