Matrix must be positive definite. Fit.PolynomailWeighted


(Eduard) #1

Hi everyone.
I’m using the method Fit.PolynomialWeighted() to interpolate a set of points where each points has its own weight. The data I’m using is the following:

xs;ys;ws
37747;0.36083865622578;0.99909942360515
37748;0.322406927332496;0.999267756924382
37749;0.323716492503407;0.999413696669535
37750;0.345244140364006;0.999538837307928
37753;0.377662559799022;0.999805429695319
37754;0.383248256140257;0.999863344417809
37755;0.353248523398455;0.999908450059564
37756;0.366464441296737;0.999942346907922
37757;0.364219472432317;0.99996663567195
37760;0.308242242830799;0.999997864660776
37761;0.291765774915146;0.999999733082431
37762;0.289119918753508;1
37763;0.335422401582168;0.999999733082431
37764;0.337680733356496;0.999997864660776
37768;0.396196915840393;0.999942346907922
37769;0.414330384584547;0.999908450059564
37770;0.440267792714541;0.999863344417809
37771;0.45184007269423;0.999805429695319
37774;0.433880321778871;0.999538837307928
37775;0.456944704278804;0.999413696669535
37776;0.503487719485795;0.999267756924382
37777;0.513429724456798;0.99909942360515

I’m calling the method as follows:

double[] polynomial = Fit.PolynomialWeighted(
x: xs,
y: ys,
w: ws,
order: 2);

The exception I’m getting is the following:
Additional information: Matrix must be positive definite.

What am I doing wrong?

thank you for your time! :slight_smile:


(Peter Vanderwaart) #2

You can get this message if either the X or W values are all zero. It happened to me (perils of cut and paste) when I tried to reproduce your result.

Edit: Expanding on my original comment. Here is my program which fits a straight line and an unweighted polynomial for comparison:

 Console.WriteLine("X,Y,W");
            for (int i = 0; i < X.Length; i++)
            {
                Console.WriteLine(X[i].ToString()+", "+Y[i].ToString()+", "+W[i].ToString());
            }
                var L = MathNet.Numerics.Fit.Line(X, Y);
            Console.WriteLine("Linear Fit Coefficients: " + L.ToString());
            var P = MathNet.Numerics.Fit.Polynomial(X, Y, 2);
            Console.WriteLine("Quadradic Fit Coefficients: " + P[0].ToString() + ", " + P[1].ToString() + ", " + P[2].ToString());

            Console.WriteLine("");
            Console.WriteLine("original problem");

            double[] WLS = new double[3];

            try
            {
                var R = MathNet.Numerics.Fit.PolynomialWeighted(X, Y, W, 2);
                Console.WriteLine("Weighted Fit Coefficients: " + R[0].ToString() + ", " + R[1].ToString() + ", " + R[2].ToString());
                WLS = R;
            }
            catch (Exception eWeighted)
            {

                Console.WriteLine(eWeighted.Message);
            }

            Console.WriteLine("Line fit: X, Y,Linear,Quadratic,Weighted");
            for (int i = 0; i < X.Length; i++)
            {
                Console.Write(X[i].ToString() + ", ");
                Console.Write(Y[i].ToString() + ", ");
                double l_fit = L.Item1 + L.Item2 * X[i];
                Console.Write(l_fit.ToString() + ", ");
                double q_fit = P[0] + P[1] * X[i] + P[2] * X[i] * X[i];
                Console.Write(q_fit.ToString() + ", ");
                double w_fit = WLS[0] + WLS[1] * X[i] + WLS[2] * X[i] * X[i];
                Console.WriteLine(w_fit.ToString());
            }

This results in the following output:

X,Y,W
37747, 0.36083865622578, 0.99909942360515
37748, 0.322406927332496, 0.999267756924382
37749, 0.323716492503407, 0.999413696669535
37750, 0.345244140364006, 0.999538837307928
37753, 0.377662559799022, 0.999805429695319
37754, 0.383248256140257, 0.999863344417809
37755, 0.353248523398455, 0.999908450059564
37756, 0.366464441296737, 0.999942346907922
37757, 0.364219472432317, 0.99996663567195
37760, 0.308242242830799, 0.999997864660776
37761, 0.291765774915146, 0.999999733082431
37762, 0.289119918753508, 1
37763, 0.335422401582168, 0.999999733082431
37764, 0.337680733356496, 0.999997864660776
37768, 0.396196915840393, 0.999942346907922
37769, 0.414330384584547, 0.999908450059564
37770, 0.440267792714541, 0.999863344417809
37771, 0.45184007269423, 0.999805429695319
37774, 0.433880321778871, 0.999538837307928
37775, 0.456944704278804, 0.999413696669535
37776, 0.503487719485795, 0.999267756924382
37777, 0.513429724456798, 0.99909942360515
Linear Fit Coefficients: (-183.108936451217, 0.00485913033756305)
Quadradic Fit Coefficients: 562496.387141512, -29.7964390935051, 0.000394592892764656

original problem
Weighted Fit Coefficients: 2334594.96940846, -123.652422706614, 0.00163732075181406

Line fit: X, Y,Linear,Quadratic,Weighted
37747, 0.36083865622578, 0.308656400775305, 0.363185534370132, 0.534969519358128
37748, 0.322406927332496, 0.31351553111287, 0.356536880135536, 0.492076971102506
37749, 0.323716492503407, 0.318374661450434, 0.350677411654033, 0.452459063380957
37750, 0.345244140364006, 0.323233791787999, 0.345607128925622, 0.416115797590464
37753, 0.377662559799022, 0.337811182800664, 0.335131395841017, 0.326733850874007
37754, 0.383248256140257, 0.342670313138228, 0.33321785624139, 0.303489151410759
37755, 0.353248523398455, 0.347529443475793, 0.33209350251127, 0.283519093412906
37756, 0.366464441296737, 0.352388573813357, 0.331758334417827, 0.26682367734611
37757, 0.364219472432317, 0.357247704150922, 0.332212352310307, 0.253402902279049
37760, 0.308242242830799, 0.371825095163615, 0.338309520389885, 0.232788426801562
37761, 0.291765774915146, 0.37668422550118, 0.341920281643979, 0.232466217596084
37762, 0.289119918753508, 0.381543355838744, 0.34632022830192, 0.235418650321662
37763, 0.335422401582168, 0.386402486176308, 0.351509360829368, 0.241645724512637
37764, 0.337680733356496, 0.391261616513873, 0.357487679109909, 0.251147440634668
37768, 0.396196915840393, 0.410698137864131, 0.389292810345069, 0.321900717914104
37769, 0.414330384584547, 0.415557268201695, 0.3992170576239, 0.347775640897453
37770, 0.440267792714541, 0.42041639853926, 0.409930490655825, 0.376925205811858
37771, 0.45184007269423, 0.425275528876824, 0.421433109440841, 0.409349412191659
37774, 0.433880321778871, 0.439852919889489, 0.460676080547273, 0.526269879657775
37775, 0.456944704278804, 0.444712050227054, 0.475335442693904, 0.571792651899159
37776, 0.503487719485795, 0.449571180564618, 0.490783990360796, 0.620590065605938
37777, 0.513429724456798, 0.454430310902183, 0.507021724013612, 0.672662120778114

The output data can be plotted by copying into a text file, saving as .csv and opening it in Excel. Then Insert Scatter plot. Graph here: https://drive.google.com/open?id=1PZKqRZ5KbvyHZUbVtPjIVti_CowOMfIV

Series 1 is the Y, Series 2 is the linear fit, Series 3 is the unweighted polynomial fit and Series 4 is the Weighted Polynomial fit. I don’t think anyone would be the last as being the most representative of the first. I find that extremely odd, especially as the % difference between all the weighting factors is tiny.


(Eduard) #3

Hi PVanderwaart. So in conclusion this seems to be happening also to you and we don’t know why, doesn’t it?
It seems that this happens if almost all weights are approximately 1, doesn’t it?

Bye and thank you!


(Peter Vanderwaart) #4

The only explanation that seems even slightly plausible to me is that the program is interpreting the weights as serial autocorrelation. In that case, it would be trying to find the polynomial fit for which the error at each X[i] is most like the error at X[i-1]. There is a good explanation of using a covariance matrix in least squares in one of the books by Wonnacott and Wonnacott, but I’m not sure which one off the top of my head.


(Peter Vanderwaart) #5

The reference was:
Econometrics, Wonnacott & Wonnacott, John Wiley & Sons, Inc,.New York, 1970

Just out of curiosity, I wrote a routine based on Eq: 16-48 (pg 328 in my edition) which I include below. I was unable to reproduce the results from Fit.PolynomialWeighted. I still suspect the that routine is expecting the weights to represent something as yet unexplained. We could do with some help from the developers.

private void btnGLS_Click(object sender, EventArgs e)
{
    Matrix<double> X = Matrix<double>.Build.Dense(22, 3);
    Vector<double> Y = Vector<double>.Build.Dense(22);
    Matrix<double> U = Matrix<double>.Build.Dense(22, 22);
    Vector<double> x = Vector<double>.Build.Dense(22);

    using (StreamReader rdr = new StreamReader(@"w_ls data.csv"))
    {
        string line;
        int iRow = 0;
        while ((line = rdr.ReadLine()) != null)
        {

            int n = line.IndexOf(',');
            double d = double.Parse(line.Substring(0, n));
            X[iRow, 0] = 1;
            X[iRow, 1] = d;
            X[iRow, 2] = d * d;

            x[iRow] = d;

            line = line.Substring(n + 1);
            n = line.IndexOf(',');
            Y[iRow] = double.Parse(line.Substring(0, n));

            line = line.Substring(n + 1);
            d = double.Parse(line);
            for (int j = 0; j < 22; j++)
            {
                if (Math.Abs(iRow - j) <= 1) U[iRow, j] = Math.Pow(d, Math.Abs(iRow - j));
                else U[iRow, j] = 0.0;
            }

            iRow++;
        }

        rdr.Close();
    }

    var B = (X.Transpose() * U.Inverse() * X).Inverse() * X.Transpose() * U.Inverse() * Y;
    Console.WriteLine(B.ToString());
    for (int i = 0; i < 22; i++)
    {
        double y = B[0] + B[1] * x[i] + B[2] * x[i] * x[i];
        Console.WriteLine(x[i].ToString() + ", " + y.ToString() + ", " + Y[i].ToString());
    }
}

(Peter Vanderwaart) #6

One more comment, then maybe I can leave this alone.

For weighted least squares, the expression we want to minimize is (pardon my pseudo-code):

sum [ w[i] * ( c0+ c1x[i] + c2x[i]^2 - Y[i] )^2 ] for i = 1,…n

where c0, c1, c2 are the parameters to be optimized.

Now, put w[i] inside the parens (where it becomes sqrt(w[i]) and multiply across the expression. That gives n new equations that can be minimized with ordinary least squares. The answer I get to the original problem (and the answer I would have expected from PolynomialWeighted) is:

c0 = 901695
c1 = -47.7615
c2 = 0.000632464

Program:

            Matrix<double> X = Matrix<double>.Build.Dense(22, 3);
            Vector<double> Y = Vector<double>.Build.Dense(22);
            Matrix<double> U = Matrix<double>.Build.Dense(22, 22);
            Vector<double> x = Vector<double>.Build.Dense(22);

            using (StreamReader rdr = new StreamReader(@"c:\w-ls data.csv"))
            {
                string line;
                int iRow = 0;
                while ((line = rdr.ReadLine()) != null)
                {

                    int n = line.IndexOf(',');
                    x[iRow] = double.Parse(line.Substring(0, n));
                    line = line.Substring(n + 1);
                    n = line.IndexOf(',');
                    double y = double.Parse(line.Substring(0, n));
                    line = line.Substring(n + 1);
                    double w = Math.Sqrt(double.Parse(line));



                    X[iRow, 0] = w;
                    X[iRow, 1] = w * x[iRow];
                    X[iRow, 2] = w * x[iRow] * x[iRow];

                    Y[iRow] = w * y;

                    iRow++;
                }

                rdr.Close();
            }

            var B = (X.Transpose() * X).Inverse() * X.Transpose() * Y;
            Console.WriteLine(B.ToString());