Data fitting questions on multiple y-functions and finding the abscissa

[Towie__Ewan]

Hi,

I expect there is a very short answer to these questions, but let me ask anyway. I have looked through the documentation and the answers aren't completely clear to me.

I see that the Data Fitting library supports cubic spline interpolation with multiple y-functions. Am I right in thinking that any request for an interpolated/integrated value will return a list of results corresponding to each y-function?

In other words, is it possible to request interpolation/integration for only one function of the multiple stored?

Also, is there functionality to locate the abscissa from a y-value? This usually involves locating multiple roots of a cubic polynomial and is a pain to compute. If not, is there any plans to implement such functionality?

Thanks!

[Ruqiu_C_Intel · Accepted solution]

Hi Towie, Ewan,

Thanks for your questions. Basically we don't support only one request interpolation/integration function to generate multiple output. But users can call the Data Fitting function several times to create multiple tasks.

There is no functionality to locate the abscissa from a y-value currently, while we will think about this feature if customers highly recommend to use.

[Towie__Ewan]

Thanks for the answers, even if they aren't really what I wanted to hear!

I think it would be a good addition to the Data Fitting library to have functionality to locate the abscissa. It is a common request when working with cubic splines, and the methods to compute can be complicated and/or inaccurate.

[Pavel_D_Intel1]

Hi Ewan,
glad to see new questions from you.
I apologize but the previous answer was not completely correct. Let me to clarify the situation.

Intel MKL DataFitting interpolation/integration routines support work with the multiple functions simultaneously.
Please find a small sample of interpolation routine for the multiple functions:

#include "mkl.h"
#include <math.h>

#define N                          4 // number of breakpoints
#define NY                        3 // number of functions
#define NSITE                   1 // number of sites for spline-based interpolation
#define NDORDER            1 // size of array describing derivative orders
#define NSCOEFF              NY*(N-1)*DF_PP_CUBIC // total number of spline coefficients

int main()
{
    DFTaskPtr task;                                          // Data Fitting task descriptor
    MKL_INT dorder = 1;                                 // only spline value will be computed  
    double x = {0, 0.5, 1.5, 2};                     // breakpoints                
    double y[NY*N];                                         // function values
    double scoeff[NSCOEFF];                         // array of spline coefficients
    double site[NSITE];                                   // array of interpolation sites
    double r[NY*NSITE];                                 // interpolation results

    /***** Generate function y = sin(2 * Pi * freq * x) *****/
    int yi, xi;
    double freq = 1.7;

    for( yi=0; yi<NY; yi++ )
        for (xi = 0; xi < N; xi++)
            y[yi*N+xi] = sin( 2.0 * 3.14 * freq * x[xi] );

    /***** Create Data Fitting task *****/
    dfdNewTask1D( &task, N, x, DF_NON_UNIFORM_PARTITION, NY, y, DF_NO_HINT );

    /***** Edit task parameters for natural cubic spline construction *****/
    dfdEditPPSpline1D( task, DF_PP_CUBIC, DF_PP_NATURAL, DF_BC_FREE_END, 0, DF_NO_IC, 0,
                                 scoeff, DF_NO_HINT );

    /***** Construct natural cubic spline using STD method *****/
    dfdConstruct1D( task, DF_PP_SPLINE, DF_METHOD_STD );

    /***** Interpolate *****/
    // r[0] – interpolation result for the first y-function
    // …
    // r[3] - interpolation result for the forth y-function

    dfdInterpolate1D( task, DF_INTERP, DF_METHOD_PP, NSITE, site, DF_NON_UNIFORM_PARTITION, NDORDER,
                                    dorder, 0, r, DF_MATRIX_STORAGE_ROWS, 0 );

    // …

    /***** Delete Data Fitting task *****/
    dfDeleteTask( &task );
    return 0;
}

Speaking about your questions:

• In other words, is it possible to request interpolation/integration for only one function of the multiple stored?

Interpolation/integration routines perform for all the functions from the multiple stored. Selection of functions subset is not supported in Intel MKL DataFitting now. Can such feature be useful for your use cases?

• Also, is there functionality to locate the abscissa from a y-value? This usually involves locating multiple roots of a cubic polynomial and is a pain to compute. If not, is there any plans to implement such functionality?

As Ruqio said – Intel MKL DataFitting doesn’t support such functionality for now. Could you please provide an information about the algorithms which you are interested in (any references/papers/links)? It would be really useful for us in terms of further analysis.

Best regards,
Pavel.

[Towie__Ewan]

Hi Pavel,

Before I respond to your question, I actually have another question!

Can you get the interpolated first and second derivative for a given x-point? I know that the derivatives are computed, but not sure if we can access them.

Dyakov, Pavel (Intel) wrote:

1. Interpolation/integration routines perform for all the functions from the multiple stored. Selection of functions subset is not supported in Intel MKL DataFitting now. Can such feature be useful for your use cases?
2. As Ruqio said – Intel MKL DataFitting doesn’t support such functionality for now. Could you please provide an information about the algorithms which you are interested in (any references/papers/links)? It would be really useful for us in terms of further analysis.

See responses below:

1. Yes, being able to independently request an interpolation/integration for a given y-function is something that we use quite often. An example would be capturing a 2D surface with certain rotational symmetries. In this case we use a set of 1D splines (with the same x-points) at certain rotational angles to represent the surface.

2. This can be quite a huge problem to capture as it depends on the order of the polynomial, as well as the magnitudes of the co-efficients obtained (where certain co-efficient can lead to numerical problems for the analytical methods). For computational efficiency, it may require selecting an appropriate method based on the cubic spline co-efficients.

Predominantly, the methods I have been currently trying to use are analytical and based on making some bounding of the problem using the co-efficients. For problems containing only real co-efficients there are several methods that work but can be problematic in certain cases. References for such approaches can be found on Wikipedia (https://en.wikipedia.org/wiki/Cubic_function), Wolfram Mathworld (http://mathworld.wolfram.com/CubicFormula.html), Numerical Recipes, and many other online resources if you search for cubic root finding. I strongly recommend reading J.F. Blinn's work on the topic, a 5-part series 'How to solve a cubic equation' originally published in IEEE Computer Graphics and Applications and now thankfully available at https://courses.cs.washington.edu/courses/cse590b/13au/lecture_notes/ (pdf files with cubic in the name). Bear in mind the topic of cubic root finding goes back some 1000+ years!

However, for some cases these analytical models are far from perfect and then we need to look at more numerical methods to solve the roots. Again there are many methods here and selecting one most suitable for the problem at hand, and being computationally efficient, probably requires some thought. In my opinion, at the moment, I would look at the Jenkins-Traub algorithm (https://en.wikipedia.org/wiki/Jenkins%E2%80%93Traub_algorithm) as an efficient and robust method. But there are several other methods employing matrices, such as the shifted QR method for finding the roots of an nth degree polynomial using LAPACK functionality.

As you can gather, I am not having fun computing the abscissa. I think I am reaching the end point of the analytical forms and moving to a numerical method may be essential. I am mostly struggling with cases where the floating-point accuracy caused by odd co-efficient magnitudes mean the bounding required for the analytical forms starts to fail. Anyway, this is beyond the scope of your question.

Thanks,

Ewan

[Pavel_D_Intel1]

Hi Ewan,

Answering your question:

You can compute derivatives for the interpolation site by setting ndorder and dorder parameters:

• ndorder - Maximal derivative order increased by one to be computed at the interpolation sites.
• dorder - Array of size ndorder that defines the order of the derivatives to be computed at the interpolation sites.

E.g. if you want to compute a value, the first derivative and the second derivative for the interpolation site you can set described parameters as follows:

• ndorder = 3;
• dorder = {1, 1, 1};

Note: the similar situation is presented in dfdcubicspline_interp.c example.

Speaking about previous questions:

1. Yes, being able to independently request an interpolation/integration for a given y-function is something that we use quite often. An example would be capturing a 2D surface with certain rotational symmetries. In this case we use a set of 1D splines (with the same x-points) at certain rotational angles to represent the surface.

As a workaround I propose to create several 1-dimensional tasks for the requested coordinates and after call the interpolation routine. But anyway we will think about APIs extension to support such use-cases.

One remaining question here – do you use the same type of boundary conditions for all functions (multiple stored) in your use cases?  

2. This can be quite a huge problem to capture as it depends on the order of the polynomial, as well as the magnitudes of the co-efficients obtained (where certain co-efficient can lead to numerical problems for the analytical methods). For computational efficiency, it may require selecting an appropriate method based on the cubic spline co-efficients.

...

Thank you for such detailed description of the problem. I will analyze the provided information and come back with the comments.

Best regards,
Pavel.

[Towie__Ewan]

No problem Pavel.

A TLDR summary of my answer to question 2 would simply be to implement a form of the Jenkins-Traub algorithm (https://en.wikipedia.org/wiki/Jenkins%E2%80%93Traub_algorithm).

If anyone wants to use an analytical cube root first of all as a starting guess, then it's relatively simple to implement themselves.

Thanks,

Ewan

[Gennady_F_Intel]

Ewan, regarding question #1 - we would encourage you to submit the Feature Request to the Intel Online Service Center and will proceed with this request in official ways.

[Towie__Ewan]

Gennady F. (Blackbelt) wrote:

Ewan, regarding question #1 - we would encourage you to submit the Feature Request to the Intel Online Service Center and will proceed with this request in official ways.

Done! If you want/need more information then let me know.

Thanks,

Ewan

[Pavel_D_Intel1]

Hello Ewan,

I apologize for the delay in response.

I have analyzed your feature requests in current forum thread and in Intel Online Service Center. Please find some analysis below:

1. Feature request #1 - function subset selection for Data Fitting interpolation/integration routines.

As was mentioned before now interpolation/integration routines perform computation for all the functions from the multiple stored.
Function subset selection is a good option to add - now I am working on API for this new feature.

Next possible use-cases types can be taken into account:

• Construct spline for all functions first, perform interpolation/integration/search for the selected functions subset.
• Construct spline for the selected function subset and after perform interpolation/integration/search for the same subset.

Could you please clarify: are you interesting only in the 1^st scenario? Or the 2^nd one can be also relevant for your work?
This information can be really helpful in terms of API development.

2. Feature request #2 – abscissa location from a y-value.

At the first – thanks a lot for the detailed information about this task and solving methods.

Intel MKL Data Fitting provides basic spline-based routines: spline construction/interpolation/integration.

Speaking about abscissa location from y-value – this task requires solving the system of non-linear/linear (depends on spline type) equations based on the specific group of the numerical methods that are out of scope of the Intel MKL Data Fitting component (probably for now).

So that I cannot expect adding this new feature in near future releases of Intel MKL Data Fitting.

Best regards,
Pavel