mathnet / mathnet-numerics

Math.NET Numerics
http://numerics.mathdotnet.com
MIT License
3.44k stars 891 forks source link

Solving systems of linear equations gives bad results #1012

Open Corcharelli opened 1 year ago

Corcharelli commented 1 year ago

Hello, I'm testing the library to find solutions of systems of linear equations, and the result I get is always wrong. I give a couple of examples: The following system has no solution: x-y = 0 x + 3y = 4 3x + 5y = 9

However, I get as a solution x = 1.2083333333333335 y = 1.0416666666666667

which clearly violates the first equation x-y = 0. I have implemented it as follows: var A = Matrix.Build.DenseOfArray(new double[,] { { 1.0, -1.0 }, { 1.0, 3.0 }, { 3.0, 5.0 } }); var b = Vector.Build.Dense(new double[] { 0.0, 4.0, 9.0 }); var x = A.Solve(b);

It seems that the Solve method always returns approximate values, never the exact solution, because I tried one of the examples in the documentation, which is the following:

var A = Matrix.Build.DenseOfArray(new double[,] { { 3, 2, -1 }, { 2, -2, 4 }, { -1, 0.5, -1 } }); var b = Vector.Build.Dense(new double[] { 1, -2, 0 }); var x = A.Solve(b);

And instead of getting the expected result x = -1, y = -2, z = -2, I get x = 1 y = -1.9999999999999996 z = -1.9999999999999992

Does anyone know what's the problem?

jkalias commented 1 year ago

Hello, I'm testing the library to find solutions of systems of linear equations, and the result I get is always wrong. I give a couple of examples: The following system has no solution: x-y = 0 x + 3y = 4 3x + 5y = 9

However, I get as a solution x = 1.2083333333333335 y = 1.0416666666666667

which clearly violates the first equation x-y = 0. I have implemented it as follows: var A = Matrix.Build.DenseOfArray(new double[,] { { 1.0, -1.0 }, { 1.0, 3.0 }, { 3.0, 5.0 } }); var b = Vector.Build.Dense(new double[] { 0.0, 4.0, 9.0 }); var x = A.Solve(b);

It seems that the Solve method always returns approximate values, never the exact solution,

Hello. The solver is using a QR decomposition internally to solve the system. If the system is overdetermined as in your case, Q will be a 3x3 matrix, whereas R a 3x2 matrix, and the process of "solving" the system is essentially finding the x which minimizes ||Ax - b||. So, yes, not all equations can be satisfied. (see this link for more info).

I tried one of the examples in the documentation, which is the following:

var A = Matrix.Build.DenseOfArray(new double[,] { { 3, 2, -1 }, { 2, -2, 4 }, { -1, 0.5, -1 } }); var b = Vector.Build.Dense(new double[] { 1, -2, 0 }); var x = A.Solve(b);

And instead of getting the expected result x = -1, y = -2, z = -2, I get x = 1 y = -1.9999999999999996 z = -1.9999999999999992

I can't tell how you reached the conclusion that the expected solution is "x = -1, y = -2, z = -2", but it clearly isn't. Take the first equation for example:

3x + 2y - 1z = 3(-1)+2(-2)-1(-2) = -3-4+2 = -4 which is not equal to b[0] = 1

whereas for the Numerics solution = 3+2(-2)-1(-2) = 3-4+2 = 1