Polynomial fitting

[ivan-pi]

Motivation

Functions for low-order polynomial fitting using (weighted) least squares regression are common in many scientific libraries.

Here's an interface for the non-weighted case, inspired by the numpy version cited below:

  function polyfit(x, y, deg, rcond, rank, singular_values) result(p)
    real, intent(in) :: x(:), y(:)
    integer, intent(in) :: deg
    real, intent(in), optional :: rcond
    integer, intent(out), optional :: rank
    real, intent(out), allocatable, optional :: singular_values(:)

Prior Art

Prior art in popular scripting languages:

• numpy.polyfit
• polyfit (MATLAB)
• Polynomials.polyfit (Julia)
• lm (R) (in R one would use a combination of lm for linear regression and poly as demonstrated here)

Prior art in Fortran (in no particular order):

• POLFIT - A 704 program for polynomial least squares fitting (Fortran II)
• LSQFT - A nonlinear least squares data fitting subroutine suitable for minicomputers
• SLATEC: dpolft - Fit discrete data in a least squares sense by polynomials (related routines include dpcoef and dp1vlu)
• IMSL Fortran Library: RCURV (also see poly_regression in the C interface)
• NAG E02 Curve and Surface Fitting (see routines e02adf and e02aef)
• polyfit - Toolbox for performing and manipulating polynomial fits to data with statistical uncertainties.
• Least-Squares Polynomial Curve-Fitting Utilizing Orthogonal Polynomials (thesis by Robert E. Knight, contains Fortran code)
• FITLOS: A FORTRAN program for fitting low-order polynomial splices by the method of least squares
• An Evaluation of Linear Least Squares Computer Programs (article by Roy H. Wampler from National Bureau of Standard, lists over 20 codes, many in Fortran)
• spfit (part of MATH77/mathc90 library, available on Netlib)
• Polyfit program users' guide : a statistical polynomial curvefitting program
• curvefit
• E201 Least Squares Polynomial Fit (part of CERNLIB)
• Numerical Algorithms in Fortran by Gisela Engeln-Müllges & Frank Uhlig (see Chapter 8 on linear and nonlinear approximation)
• LS_POLY routine by Jean-Pierre Moreau (part of the LSQPLY.F90 file)
• PROGRESS: a program for robust regression
• IBM Scientific Subroutine Package (some Fortran source files are available here)

Prior art in other programming languages

• GSL: Linear Least-Squares Fitting (can be used via the Fortran interface fgsl, the relevant routine is called multifit which may be used for any linear model)
• alglib (C++)

Additional Information

This will probably require LAPACK for the factorization of the Vandermonde matrix either using a rank-revealing QR or SVD factorization.

A function to evaluate the polynomial (and it's derivatives) using Horner's scheme (or alternatively, in some orthogonal form), should be added in parallel.

An advanced interface could consider using derived types for polynomials, similar to the numpy Polynomial sub-package.

[ivan-pi]

Given how many times this has been re-implemented both in Fortran and other languages, I suspect that finding an interface which meets everyone's needs won't be easy. Hopefully we can reach an 80/20 compromise.

Here are some numbered items to help guide the discussion:

1. subroutine vs function
2. automatic vs assumed-shape arrays for the x and y arguments
3. weighted vs unweighted regression
4. manual vs automatic centering and scaling for better conditioning of the numerical problem
5. Taylor (1, x, x**2...) vs orthogonal polynomial (Chebyshev) form
6. statistics summary, goodness of fit (standard deviation, R^2, variance-covariance matrix...)
7. control over the factorization method (and inspection of the factorized matrix, condition number, singular values, rank)
8. recovery from input errors or an ill-conditioned (singular) matrix, etc.
9. a method to fix the intercept y = p(0).

[Beliavsky]

Regarding the interface, a function with intent(out) arguments should be a subroutine IMO.

Often the user will want to fit successively higher polynomial orders and choose the best order using an information criterion such as AIC. Ideally such fits would be done efficiently (there are special algorithms for successively adding predictors).

Least-squares regression by successively adding general basis functions (here they are powers) is sometimes done.

[ivan-pi]

How about this as a simple interface?

program polyfit_example
implicit none
interface
    subroutine polyfit(x,y,n,p) bind(c,name="c_arma_polyfit")
        use iso_c_binding, only: c_double, c_int
        real(c_double), intent(in), contiguous :: x(:)
        real(c_double), intent(in), contiguous :: y(:)
        integer(c_int), intent(in), value :: n
        real(c_double), intent(out) :: p(:)
    end subroutine
end interface

integer, parameter :: dp = kind(1.0d0)

real(dp) :: x(10), y(10)
real(dp) :: p(3)
integer :: i

x = [((i-1)*0.2_dp,i=1,10)]
y = x**2 - 2.0_dp*x + 3.0_dp

call polyfit(x,y,2,p)

print *, "coefficients = ", p
print *, all(abs(p - [1.0_dp,-2.0_dp,3.0_dp]) < 100*epsilon(1.0_dp))

end program

For the implementation I just wrote a little wrapper of the polyfit function in the Armadillo C++ library:

Click for details

// c_arma_polyfit.cpp
#include <cassert>

#include <ISO_Fortran_binding.h>
#include <armadillo>

template<class C>
C as_(CFI_cdesc_t *x);


template<>
auto as_<arma::Col<double>>(CFI_cdesc_t *x) -> arma::Col<double> {

    double* ptr = static_cast<double *>(x->base_addr);
    auto nx = static_cast<arma::uword>(x->dim[0].extent);

    arma::Col<double> x_{ ptr, nx, 
        false, /* copy_aux_mem */
        true   /* strict */
    };
    return x_;
}

extern "C"
void c_arma_polyfit(CFI_cdesc_t *fx, CFI_cdesc_t *fy, int n, CFI_cdesc_t *fp)
{
    const auto x = as_<arma::Col<double>>(fx);
    const auto y = as_<arma::Col<double>>(fy);
    assert(x.n_elem == y.n_elem);
    
    auto p = as_<arma::Col<double>>(fp);
    assert(p.n_elem == arma::uword(n+1));

    [[maybe_unused]] bool success = arma::polyfit(p,x,y,n);
}

# Makefile

FC=gfortran-12
CXX=g++-12

FCFLAGS=-Wall
CXXFLAGS=-Wall

ARMA_INC=-I/usr/local/Cellar/armadillo/12.0.1/include

ARMA_LDFLAGS=-L/usr/local/Cellar/armadillo/12.0.1/lib
ARMA_LDLIBS=-larmadillo -lstdc++

polyfit_example: polyfit_example.f90 c_arma_polyfit.o
	$(FC) $(FCFLAGS) $(ARMA_LDFLAGS) -o $@ $^ $(ARMA_LDLIBS)

c_arma_polyfit.o: c_arma_polyfit.cpp
	$(CXX) $(CXXFLAGS) $(ARMA_INC) -c $<

.PHONY: clean

clean:
	rm -f *.o
	rm -f polyfit_example

Here's an alternative implementation using Eigen:

Click for details

// c_eigen_polyfit.cpp
//
#include <cassert>
#include <Eigen/Dense>

#include <ISO_Fortran_binding.h>

template<class C>
Eigen::Map<C> as_(CFI_cdesc_t *x);

template<>
auto as_<Eigen::VectorXd>(CFI_cdesc_t *x) -> Eigen::Map<Eigen::VectorXd> {

    double* ptr = static_cast<double *>(x->base_addr);
    auto nx = static_cast<Eigen::Index>(x->dim[0].extent);

    Eigen::Map<Eigen::VectorXd> x_{ ptr, nx, 1};
    return x_;
}

// Implementation based upon
//    https://gist.github.com/ksjgh/4d5050d0e9afc5fdb5908734335138d0
//
Eigen::VectorXd polyfit(const Eigen::VectorXd& x, const Eigen::VectorXd& y, int n)
{
    assert(x.size() == y.size());
    assert(1 <= n && n <= x.size() - 1);

    Eigen::MatrixXd A(x.size(), n+1);

    // Calculate Vandermonde matrix of shape m-by-(n+1)
    //   
    //    x1^n x1^n-1 ... 1
    //    x2^n x2^n-1 ... 1
    //     :    :
    //    xm^n xm^n-1 ... 1
    //
    A.col(n).setOnes();
    for (int i = n; i > 0; i--) {
        A.col(i-1) = A.col(i).array() * x.array();
    }

    auto QR = A.householderQr();
    auto p = QR.solve(y);

    return p;
}

extern "C"
void c_eigen_polyfit(CFI_cdesc_t *fx, CFI_cdesc_t *fy, int n, CFI_cdesc_t *fp)
{

    const auto x = as_<Eigen::VectorXd>(fx);
    const auto y = as_<Eigen::VectorXd>(fy);
    assert(x.size() == y.size());

    auto p = as_<Eigen::VectorXd>(fp);
    assert(p.size() == Eigen::Index(n+1));

    p = polyfit(x,y,n);

}

To compile use: g++ -I/usr/local/include/eigen3 -c c_eigen_polyfit.cpp.

[ivan-pi]

Often the user will want to fit successively higher polynomial orders and choose the best order using an information criterion such as AIC. Ideally such fits would be done efficiently (there are special algorithms for successively adding predictors).

Least-squares regression by successively adding general basis functions (here they are powers) is sometimes done.

This sounds like something the qrupdate library could be used for. I believe this library is used in the Octave qrupdate routine. It was written by Jaroslav Hajek, a computing expert from the Aeronautical Research and Test Institute (VZLU) in Prague. A GitLab mirror for qrupdate can be found here.

[ivan-pi]

A prototype Fortran implementation I wrote 3 years ago can be found here: https://gist.github.com/ivan-pi/0fd517048c415ceca441eac626bf24c9

[ivan-pi]

cc @arjenmarkus, @nshaffer, @jvdp1

We need a few more voices to move this issue forward. Given how many polynomial fitting codes are out there, it looks like polynomial regression is quite popular among practitioners.

Do you prefer the simple interface or the complex one?

    ! simple
    subroutine polyfit(x,y,order,p)
        real(wp), intent(in) :: x(:), y(:)
        integer(ip), intent(in) :: order
        real(wp), intent(out) :: p(:)
    end subroutine

    ! complex
   subroutine polyfit(x, y, order, p, rcond, rank, singular_values)
      real(wp), intent(in) :: x(:), y(:)
      integer(ip), intent(in) :: order
      real(wp), intent(out) :: p(:)
      real(wp), intent(in), optional :: rcond
      integer(wp), intent(out), optional :: rank
      real(wp), intent(out), allocatable, optional :: singular_values(:)
   end subroutine

[arjenmarkus]

Let me read the mail thread (again) and I will get back to your question tomorrow. Op wo 1 mrt 2023 om 14:44 schreef Ivan Pribec ***@***.***>:

…

cc @arjenmarkus <https://github.com/arjenmarkus>, @nshaffer <https://github.com/nshaffer>, @jvdp1 <https://github.com/jvdp1> We need a few more voices to move this issue forward. Given how many polynomial fitting codes are out there, it looks like polynomial regression is quite popular among practitioners. Do you prefer the simple interface or the complex one? ! simple subroutine polyfit(x,y,n,p) real(wp), intent(in) :: x(:) real(wp), intent(in) :: y(:) integer(ip), intent(in) :: n real(wp), intent(out) :: p(:) end subroutine ! complex subroutine polyfit(x, y, deg, p, rcond, rank, singular_values) real(wp), intent(in) :: x(:), y(:) integer(ip), intent(in) :: deg real(wp), intent(out) :: p(:) real(wp), intent(in), optional :: rcond integer(wp), intent(out), optional :: rank real(wp), intent(out), allocatable, optional :: singular_values(:) end subroutine — Reply to this email directly, view it on GitHub <#601 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAN6YRZGBXTZDOB2D7LKORLWZ5HENANCNFSM5KDTWO2A> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[arjenmarkus]

I would opt for the extended interface: it provides just that more output that might otherwise be forgotten about. For instance, the possibility of getting back the condition number means that people are aware that such a thing as the condition number is of importance in fitting problems, even if they will not use it. By the way, I did not quite understand the argument "singular_values" - does this have to do with the SVD? Op wo 1 mrt 2023 om 16:17 schreef Arjen Markus ***@***.***>:

…

Let me read the mail thread (again) and I will get back to your question tomorrow. Op wo 1 mrt 2023 om 14:44 schreef Ivan Pribec ***@***.***>: > cc @arjenmarkus <https://github.com/arjenmarkus>, @nshaffer > <https://github.com/nshaffer>, @jvdp1 <https://github.com/jvdp1> > > We need a few more voices to move this issue forward. Given how many > polynomial fitting codes are out there, it looks like polynomial regression > is quite popular among practitioners. > > Do you prefer the simple interface or the complex one? > > ! simple > subroutine polyfit(x,y,n,p) > real(wp), intent(in) :: x(:) > real(wp), intent(in) :: y(:) > integer(ip), intent(in) :: n > real(wp), intent(out) :: p(:) > end subroutine > > ! complex > subroutine polyfit(x, y, deg, p, rcond, rank, singular_values) > real(wp), intent(in) :: x(:), y(:) > integer(ip), intent(in) :: deg > real(wp), intent(out) :: p(:) > real(wp), intent(in), optional :: rcond > integer(wp), intent(out), optional :: rank > real(wp), intent(out), allocatable, optional :: singular_values(:) > end subroutine > > — > Reply to this email directly, view it on GitHub > <#601 (comment)>, > or unsubscribe > <https://github.com/notifications/unsubscribe-auth/AAN6YRZGBXTZDOB2D7LKORLWZ5HENANCNFSM5KDTWO2A> > . > You are receiving this because you were mentioned.Message ID: > ***@***.***> >