Use Math.NET's Fit.Polynomial method on functions of multiple parameters

Question

I previously used Math.NET Numerics library's Fit.Polynomial method to fit a cubic polynomial on a set of data that could be modeled as a function of one parameter y=f(x).
Now I would like to similarly find a 2 or 3 order polynomial that fits data that could be modeled as a function depending on multiple parameters y=f(x1, x2, x3, x4).

Is there already a built-in function in Math.NET that can compute that polynomial?
If not, do you see how I could manipulate my data in order to submit it to Fit.Polynomial?

Rachellerachis

What would your f look like exactly in case of, say, 2 variables and order 2? – Abridge

To keep it simple, I would say something like: f(x1, x2) = a*x1*x1 + b*x1 + c*x2*x2 + d*x2 + e – Rachellerachis

Maybe I can use the Fit.MultiDim method with an approach similar to what you describe in your post about linear regression: using x1*x1 and x2*x2 as separate parameters. However I am afraid this will produce results less accurate than using Fit.Polynomial (it was the case when I tried to do Cubic Function Fitting using your trick to leverage Linear Regression). Do you see any better method? – Rachellerachis

Answer 1

The Fit class is just a facade that is good enough in most scenarios, but you can always use the algorithms directly to get exactly what you need.

Fit.Polynomial: Polynomial curve fitting with high orders is a bit problematic numerically, so specialized algorithms and routines to tune/refine parameters at the end have been developed. However, Math.NET Numerics just uses a QR decomposition for now (although it is planned to replace the implementation at some point):

public static double[] Polynomial(double[] x, double[] y, int order)
{
    var design = Matrix<double>.Build.Dense(x.Length, order + 1, (i, j) => Math.Pow(x[i], j));
    return MultipleRegression.QR(design, Vector<double>.Build.Dense(y)).ToArray();
}

Fit.MultiDim on the other hand uses normal equations by default, which is much faster but less numerically robust than the QR decomposition. That's why you've seen reduced accuracy with this method.

public static double[] MultiDim(double[][] x, double[] y)
{
    return MultipleRegression.NormalEquations(x, y);
}

In your case I'd try to use the MultipleRegression class directly, with either QR (if good enough) or Svd (if even more robustness is needed; much slower (consider to use native provider if too slow)):

var x1 = new double[] { ... };
var x2 = new double[] { ... };
var y = new double[] { ... };

var design = Matrix<double>.Build.DenseOfRowArrays(
    Generate.Map2(x1,x2,(x1, x2) => new double[] { x1*x1, x1, x2*x2, x2, 1d }));
double[] p = MultipleRegression.QR(design, Vector<double>.Build.Dense(y)).ToArray();

(Using Math.NET Numerics v3.0.0-alpha7)

Selfwinding

Thank you for the suggestions, I will try them asap. – Rachellerachis

I managed to get results using MultipleRegression.QR, but I want more precision (no problem if the process is slow). I tried to use MultipleRegression.Svd instead, but I block on an exception on line 2212 of this file (revision 857751979..). Are you aware of a revision where this problem does not happen? – Rachellerachis

What kind of exception? – Abridge

The error is Index was outside the bounds of the array. It appears that the size of array double[] u passed is 0. It seems to be created here, but is not set afterwards, and passed two lines later to Control.LinearAlgebraProvider.SingularValueDecomposition via u.Values. – Rachellerachis

Answer 2

RosettaCode proposes this solution for polynomial regression(using Math.Net):

    public static double[] Polyfit(double[] x, double[] y, int degree)
    {
        // Vandermonde matrix
        var v = new DenseMatrix(x.Length, degree + 1);
        for (int i = 0; i < v.RowCount; i++)
            for (int j = 0; j <= degree; j++) v[i, j] = Math.Pow(x[i], j);
        var yv = new DenseVector(y).ToColumnMatrix();
        QR<double> qr = v.QR();
        // Math.Net doesn't have an "economy" QR, so:
        // cut R short to square upper triangle, then recompute Q
        var r = qr.R.SubMatrix(0, degree + 1, 0, degree + 1);
        var q = v.Multiply(r.Inverse());
        var p = r.Inverse().Multiply(q.TransposeThisAndMultiply(yv));
        return p.Column(0).ToArray();
    }

Astrobiology

This code requires a small update, initialize the QR like this: QR<double> qr = v.QR();. Now it supports Generic Types. Also updated in RossetaCode – Alewife

flip QR<double> qr = v.QR(); to var qr = v.QR(); this does away with a using statement :D – Roz

Answer 3

Notice that the x in the linear model can also be a vector x=[x1 x2 ⋯ xk] and the arbitrary functions fi(x) can accept vectors instead of scalars.
Here is something nearly about what you want.

Beghtol