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?

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)

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

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.