mljs / regression-polynomial

Polynomial Regression
MIT License
15 stars 3 forks source link

"LU matrix is singular" Error for well-behaving data (LuDecomposition.solve in ml-matrix) #10

Open saunus opened 1 year ago

saunus commented 1 year ago

Hi,

for the data as shown in the example below, I get the following error if I try to fit well-behaved data with fit-order 2:

Error: LU matrix is singular
    at LuDecomposition.solve (/home/damner/code/node_modules/ml-matrix/matrix.js:3201:13)
    at solve (/home/damner/code/node_modules/ml-matrix/matrix.js:3977:43)
    at regress (/home/damner/code/node_modules/ml-regression-polynomial/lib/index.js:174:45)
    at new PolynomialRegression (/home/damner/code/node_modules/ml-regression-polynomial/lib/index.js:54:28)

The x-data is monotonous, the y-data looks like a noisy parabola. Matlab has no issue with fitting.

Any ideas, why this is failing?

Example:

const mlr = require('ml-regression-polynomial');
x = [404.26220703125, 404.27349853515625, 404.28448486328125, 404.29541015625, 404.3040771484375, 404.3128356933594, 404.3204040527344, 404.33154296875, 404.3409423828125, 404.3497314453125, 404.3594665527344, 404.36865234375, 404.37451171875, 404.3879699707031, 404.3951416015625, 404.40545654296875, 404.4149475097656, 404.4233703613281, 404.431884765625, 404.4431457519531, 404.4515686035156, 404.4609069824219, 404.4708251953125, 404.479248046875, 404.48883056640625, 404.4971923828125, 404.50799560546875, 404.5179443359375, 404.5263977050781, 404.5354919433594, 404.5458068847656, 404.5533447265625];
y = [0.3631210884030956, 0.38180362489729697, 0.391453923949123, 0.4037110236952629, 0.4163452490432181, 0.4288353165530643, 0.4430856652809317, 0.45358586770200177, 0.4562565407929964, 0.4650388402896439, 0.47573011091885675, 0.4837785283897874, 0.4905844679839888, 0.49022240374688425, 0.49039250930602896, 0.49193229796429716, 0.490495250991277, 0.504985674313608, 0.49096439100589284, 0.4781324538477614, 0.48063366656261153, 0.46976621057011786, 0.4587271348758588, 0.44576288402868225, 0.43918945480010974, 0.42895113530071644, 0.41687014727334876, 0.3938250656311999, 0.38222703735756586, 0.3651908364214456, 0.34985174702375216, 0.34842848544450433];

const regression = new mlr.PolynomialRegression(x, y, 2);

Thanks for your help.

ghost commented 1 year ago

Hi @saunus and sorry about that.

I don't have Matlab installed, but does this result look correct ? (At least the last term is negative)

PolynomialRegression {
  degree: 2,
  powers: [ 0, 1, 2 ],
  coefficients: [
    0.0004257708983743669,
    0.08609266575876973,
    -0.00021019886491302486
  ]
}

Also r.predict(404.5) gives 0.43211824353537054

I have tested old and newer versions of the software, and I think we just need to pass true to the solve function.

(The mathematical reason could be that by rounding errors those rows end up being multiple of each other, that would make the system of equations unsolvable without SVD, probably.)

So I may send a PR with this fix.

saunus commented 1 year ago

Thanks for digging into that!

Unfortunately, solving the singular matrix did not yield the correct result, but suppresses the quadratic term: image

Correct result would be: [-1126496.33992227, -6.88810931148790, 5571.1518772892]

I think the current implementation of polynomial regression is not suited for fitting data with offsets >> variation. The question is, whether data normalization should be part of a fitting function, or it should be the user's responsibility.

Thanks for your help

ghost commented 1 year ago

I'd assume this is left in user's hands at least for now, since it would require quite a few changes.

Maybe @targos has a different idea, otherwise we may update the docs and close the issue.

ghost commented 1 year ago

I've been doing some research on this issue, and it seems to be that calculating XtX which is what the current algorithm does, increases the error in the result.

I have replaced it by a QR decomposition, and got this results:

PolynomialRegression {
  degree: 2,
  powers: [ 0, 1, 2 ],
  coefficients: [ -1126496.3326794403, 5571.151841472426, -6.888109267208169 ]
}

Which although they are in a different order, are very close to what is reported:

Correct result would be: [-1126496.33992227, -6.88810931148790, 5571.1518772892]

I could try to find some backing for the claim, if there is will to change (only a few lines) of the algorithm.

Basically, now it does not use the transpose but does this:

  const qrF = new QrDecomposition(F);
  console.log(qrF, F, y);//temporary
  return {
    coefficients: qrF.isFullRank()
      ? qrF.solve(Y).to1DArray()
      : solve(F, Y).to1DArray(),
    degree: Math.max(...powers),
    powers,
  };

So at least seems to extend the range of usability. Currently all tests pass but one.

@saunus @targos

saunus commented 1 year ago

Hi,

I did not dare to ask for QRDecomposition in first place:

But maybe something to give you confidence: QR-decomposition is used by matlab internally for least-squares:


function [p,S,mu] = polyfit(x,y,n)

% Construct the Vandermonde matrix V = [x.^n ... x.^2 x ones(size(x))]
V(:,n+1) = ones(length(x),1,class(x));
for j = n:-1:1
    V(:,j) = x.*V(:,j+1);
end

% Solve least squares problem p = V\y to get polynomial coefficients p.
[p, rankV, QRfactor, perm] = matlab.internal.math.leastSquaresFit(V,y1);

where leastSquaresFit is

function [p, rankV, QRfactor, perm] = leastSquaresFit(V,y)

% Solve least squares problem p = V\y to get polynomial coefficients p.
[QRfactor, tau, perm, rankV, ~] = matlab.internal.decomposition.builtin.qrFactor(V, -2);

% use nonzero diagonal entries to determine rank for qrSolve.
rV = sum(abs(getDiag(QRfactor)) ~= 0);

% QR solve with rank = rV.
p = matlab.internal.decomposition.builtin.qrSolve(QRfactor, tau, perm, y, rV);

I can not contribute, on what downside the change to QRDecomposition might have. But I would definitely appreciate to have the option to use it in this package.

Thanks for your help

ghost commented 1 year ago

It makes sense. I think one approach could be to either do it when it fails, and do QR as a fallback, and another one would be use QR first.

I will try the second approach for now, and will add a benchmark to see how bad it gets.

For the moment, the results seem better (but I think it does not get so close as the results you shared.)

before1

after1

It will take a little while though, cause I don't want to introduce a more complex approach that we can't fix later on. Any extra data (just the arrays of numbers.) where it fails would be great.

saunus commented 1 year ago

Hi,

I extracted a full data set of 450 curves to fit: fitdata.txt

The file contains alternating x/y arrays.

Some of them throw a singularity error, some of them work, some of them give bad fit results albeit not throwing an error!. Unfortunately I can not spend the time on preparing each case.

The following wrapper solves the issue for me in all cases and does not slow down the regression noticeably. I only tested it on the 450 data sets with fit order 2 and on the same amount of data with fit order 1, but not for higher orders.

import PolynomialRegression from "ml-regression-polynomial";

export class NormalizedPolynomialRegression {
  constructor(x, y, order) {

    // Calculate normalization
    const offsX = -Math.min(...x);
    const factX = 1 / (Math.max(...x) - Math.min(...x));

    const offsY = -Math.min(...y);
    const factY = 1 / (Math.max(...y) - Math.min(...y));

    // Apply normalization
    const x1 = x.map((x) => (x + offsX) * factX);
    const y1 = y.map((y) => (y + offsY) * factY);

    // Regression on normalized data
    const regression = new PolynomialRegression(x1, y1, order);
    const normCoefficients = regression.coefficients;

    // Retrive coefficients from normalized coefficients
    let coeffs = Array(order + 1).fill(0);

    coeffs[0] = normCoefficients[0];

    for (let n = 1; n < order + 1; n++) {
      const powers = calculatePowersOfElements(n);

      const currentPreFact = normCoefficients[n] * factX ** n;

      for (let j = 0; j < powers.length; j++) {
        const currentExp = powers[j].powerA;
        const currentFact = powers[j].factor * offsX ** powers[j].powerB;
        coeffs[currentExp] = coeffs[currentExp] + currentPreFact * currentFact;
      }
    }

    coeffs = coeffs.map((coeff) => coeff / factY);
    coeffs[0] = coeffs[0] - offsY;

    this.coefficients = coeffs;
  }

  predict(x0) {
    let val = 0;
    for (let n = 0; n < this.coefficients.length; n++)
      val = val + this.coefficients[n] * x0 ** n;
    return val;
  }
}

function calculatePowersOfElements(n) {
  const powers = [];

  for (let i = 0; i <= n; i++) {
    const aPower = i;
    const bPower = n - i;

    powers.push({
      factor: binomialCoefficient(n, i),
      powerA: aPower,
      powerB: bPower,
    });
  }

  return powers;
}

function binomialCoefficient(n, k) {
  if (k < 0 || k > n) {
    return 0;
  }

  if (k === 0 || k === n) {
    return 1;
  }

  let result = 1;
  for (let i = 1; i <= k; i++) {
    result *= (n - i + 1) / i;
  }

  return result;
}

Thanks for helping out.