Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.

MKL - General Polynomial Fitting function?

LWang77
Beginner
2,306 Views

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 ai of a polynomial:

P= a0 + a1Xi + a2Xi2 ... anXin

Where,

χ2 = sum( 1/σi2 * (Pi - Yi));

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)?

0 Kudos
4 Replies
mecej4
Honored Contributor III
2,306 Views

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.

0 Kudos
LWang77
Beginner
2,306 Views

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.

0 Kudos
mecej4
Honored Contributor III
2,306 Views

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.

0 Kudos
karlmarcox
Beginner
1,626 Views

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);

 

}
0 Kudos
Reply