Wrapper function to do polynomial fits with gsl

I have two arrays of double, say x and y containing some physical data. I would like to find the best-fitting polynomial function at a given order:

$y=pol_n(x)=c_0 + c_1*x + c_2*x^2 + ... + c_n * x^n$

and return the coefficients ($c$) in an std::vector. There are no constraints on the function interface. The GNU Scientific Library, gsl, contains some handy functions to properly do the maths, however those functions are in plain old C. I have encapsulated them in this piece of C++ code:

std::vector<double> gsl_polynomial_fit(const double * const data_x,
const double * const data_y,
const int n, const int order,
double & chisq) {
gsl_vector *y, *c;
gsl_matrix *X, *cov;
y = gsl_vector_alloc (n);
c = gsl_vector_alloc (order+1);
X   = gsl_matrix_alloc (n, order+1);
cov = gsl_matrix_alloc (order+1, order+1);

for (int i = 0; i < n; i++) {
for (int j = 0; j < order+1; j++) {
gsl_matrix_set (X, i, j, pow(data_x[i],j));
}
gsl_vector_set (y, i, data_y[i]);
}

gsl_multifit_linear_workspace * work = gsl_multifit_linear_alloc (n, order+1);
gsl_multifit_linear (X, y, c, cov, &chisq, work);
gsl_multifit_linear_free (work);

std::vector<double> vc;
for (int i = 0; i < order+1; i++) {
vc.push_back(gsl_vector_get(c,i));
}

gsl_vector_free (y);
gsl_vector_free (c);
gsl_matrix_free (X);
gsl_matrix_free (cov);

return vc;
}


Although it is pretty much self-contained it doesn't look very safe. With speed and code clarity in mind, what could be done to improve it?

One issue is numerical stability. Although the GSL uses stable routines, the problem of fitting a polynomial in the monomial basis is so ill-conditioned (especially for large ranges in x and large (eg 50) degree), that you risk getting significant errors. Unless you really, really need the coefficients in this basis then I think you should use another basis, for example Legendre or Chebycheff polynomials. This will mean that any code that uses the coefficients will need to know how to do so, but the evaluation is in fact almost as efficient as evaluation in the polynomial basis.

• Many thanks for the comment. I'm actually interested in the monomial coefficients, but at quite low order (<5) and the fit should always stay well below order 10. But I will keep this in mind! Do you have any reference for further reading? Dec 2 '14 at 9:49
• I like Numerical Linear Algebra by Treferthen and Bau. Dec 2 '14 at 10:13

std::vector<double> gsl_polynomial_fit(const double * const data_x,
const double * const data_y,
const int n, const int order,
double & chisq) {


Why are you passing arrays as pointers around.
Avoid passing pointers at all costs as they do not convey ownership semantics. This makes re-using the code hard.

You should probably use a namespace:

std::vector<double> gsl_polynomial_fit();


Prefixing the function name to prevent collisions is real old school C. We have better constructs for that. Use a namespace it prevents clashes.

namespace gsl
{
std::vector<double>  polynomial_fit();
}


It also makes the code more readable (as I can see the name of the function).

Also I am betting gsl_ is reserved for the library you are using. By adding this prefix to your own code you are asking for trouble (the whole reason the gsl_ library is adding the prefix is to avoid collisions with user code. If you add it to your functions then you void the whole point of the prefixing processes).

Declare one variable per line:

 gsl_vector *y, *c;
gsl_matrix *X, *cov;


Use RAII to do resource allocation and release:
Any code that does an allocate block like this.

  y = gsl_vector_alloc (n);
c = gsl_vector_alloc (order+1);
X   = gsl_matrix_alloc (n, order+1);
cov = gsl_matrix_alloc (order+1, order+1);


Followed by code and finally a deallocate block:

  gsl_vector_free (y);
gsl_vector_free (c);
gsl_matrix_free (X);
gsl_matrix_free (cov);


Is just screaming out to be converted into C++ (Rather than C).

The following can be made much simpler by using an iterator:

  std::vector<double> vc;
for (int i = 0; i < order+1; i++) {
vc.push_back(gsl_vector_get(c,i));
}

return vc;


It would be much easier to read as:

 return std::vector<double>(gsl::vector::iterator(gsl_vector_get, c, 0), gsl::vector::iterator(gsl_vector_get, c, order+1));


Based on comment:

You don't need a specific get when calling C functions.

 class My_ gsl_vector
{
gsl_vector*  data;
public:
My_ gsl_vector()
: data(gsl_vector_alloc())
{
// Presumably the alloc can fail with NULL being returned.
// You should check if NULL was returned and throw an exception
// on failure.
}
~My_ gsl_vector()
{
gsl_vector_free(data);
}
// If you add a conversion operator.
// You can call the normal C routines with your object.
// This automatically retrieves the internal pointer at call point.
operator gsl_vector*() const {return data;}
};

int main()
{
My_ gsl_vector    vec1;

gsl_vector_set (vec1, 1, 5); // vec1 auto converted to gsl_vector*
}


When you wrap a C object but still have to call lots of C code (because the interface is too extensive to wrap every function call in the interface (you may do it over time but not all today)). Then you can add a conversion operator that returns the C object when it us used in a context where you need the original C object. This allows you to pass your C++ object as a parameter to the C interface and it auto converts.

• While I agree that pointers convey little to no ownership, I'm not completely against passing a "const const" array pointer to a function. With a const pointer to a const array, it is quite clear, I think, that you are only supposed to read from it. Dec 1 '14 at 20:30
• Thanks for the comments. I tried to apply RAII but I had to add .get() everywhere to extract the naked pointers to pass to gsl. This seemed a bit ugly, but probably it's worth it. Could you please explain or give a reference for gsl::vector::iterator(gsl_vector_get,0). I do not understand how you use the function gsl_vector_get here. Did you mean std::vector:: instead of gsl::vector::? Dec 2 '14 at 10:02
• Rather than add get() add a conversion operator to your class. See above. Dec 2 '14 at 15:45
• No. I want to pass the function as a parameter to the iterator object. The idea being that the iterator keeps the parameters internally. When you call the operator* it uses the function gsl_vector_get to retrieve the appropriate object from the underlying vector and when you call operator++() it increments the index (third parameter). Creating your own iterator may be a bit advanced at this point. Dec 2 '14 at 15:49