MKL - General Polynomial Fitting function?

[LWang77]

Hi,

I’m looking for a general Polynomial Fitting function in MKL, which is similar to Polyfit function from Matlab, OptiVec and Armadillo.  This function uses two set of input data X and Y to determine the coefficients a_i of a polynomial:

P_i = a₀ + a₁X_i + a₂X_i² ... a_nX_iⁿ

Where,

χ² = sum( 1/σ_i² * (P_i - Y_i)² );

Is minimized.

I have searched MKL Data Fitting Library and found that all of routines are Spline-based.  In MKL LAPACK library, I found some routines the Least Squares ones.

I understand that,

1. as the order of the polynomial increases, Spline-based methods are preferable over polynomial interpolation because its interpolation error increases
2. when the number of points of X and Y > (1 + the degree) of the polynomial to be used, the Least Squares fit will be used

However, what I’m looking for is a general Polynomial Fitting function described above.  Is there such function in MKL?  Where can I find it in MKL (or IPP)?

[mecej4]

The Matlab polyfit routine is restricted to fitting with uniform weighting, i.e., σ_i is the same for all i. This capability is available in Lapack, in the subroutine ?GELS.

There is no ready-to-use-routine for the non-uniform weights case in MKL, as far as I know, but you should be able to put one together using the matrix manipulation routines that are available in BLAS/Lapack.

[LWang77]

Thank you for your advice, mecej4!

According to the MKL Reference, the function ?gels

Uses QR or LQ factorization to solve a overdetermined or underdetermined linear system with full rank matrix.

The uniform weighting is ok for our application.  In addition to the Least Squares fit, however, we are looking for Polynomial Fitting function that determine the coefficients of a polynomial for:

      Number of input points = (1 + degree) of the polynomial to be used

Is there such Polynomial Fitting function in MKL (or IPP)? In other words, we need to calculate the coefficients of a polynomial that neither overdetermined nor underdetermined.

[mecej4]

The interpolating polynomial can be obtained in the Lagrange form explicitly, without any matrix algebra. See https://en.wikipedia.org/wiki/Lagrange_polynomial .

If the coefficients are desired, rather than interpolants, you may use the Lapack routine ?GESV, which is available in MKL.

[karlmarcox]

Ok，let me fill and close this question.Here comes the code.I registerd this account just for posting this.

void cf1::pow_matrix(double *a,long n,double *res,int order){

    for(int i=0;i<=order;i++){

        vdPowx(n,a,i,res+i*n);

    }

}

void cf1::polyfitting(double *x,long n,int order0,double *y,double *y_){

    int order=order0+1;

    vector<double> vpx(n*(order));

    double *px=vpx.data();

    pow_matrix(x,n,px,order0);

    vector<double> vxtx(order*order);

    double *xtx=vxtx.data();

    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, order, order, n, 1.0, px, n, px, n, 0.0, xtx, order);  

    vector<double> vxty2beta(order);//This vector will be transferred from xty to beta,thanks for the anti-human design of the MKL-API!

    double *xty2beta=vxty2beta.data();

    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, order, 1, n, 1.0, px, n, y, n, 0.0, xty2beta, 1);

    vector<int> vipiv(order*order);

    int *ipiv=vipiv.data();

    LAPACKE_dgesv (LAPACK_ROW_MAJOR, order,1 , xtx , order , ipiv , xty2beta ,1 );

    cblas_dgemm(CblasRowMajor,CblasTrans,CblasNoTrans,n,1,order,1.0,px,n,xty2beta,1,0.0,y_,1);

}