Gauss elimination code (Julia) is wrong

[VikingScientist]

Basically

• Don't check for singularity as Gauss Elimination can be done on singular matrices
• Pivots are not neccessarily on the diagonal
• Cover the general case and don't assume triangular matrices

This issue affects PR #200 (contains more details on the bugs) and PR #195

[jiegillet]

Also the back substitution is more complex since the echelon form doesn't necessarily have non-zero elements on the diagonal. See wikipedia.
The current implementation is correct if we only ever care about non-singular matrices, but the algorithm is more general than that.

[jiegillet]

#183 (which is already merged) is also wrong

[leios]

Hey guys, I mentioned this in #200, but this is actually the correct location to talk about this. Here's my current understanding.

Standard implementations of Gaussian Elimination do not usually work on non-singular matrices. Most implementations are used specifically to solve a system of equations (as outlined in the chapter), so pivots must be on the diagonal and the output must be triangular, otherwise the solution is either over or underdetermined. This is why we check all the time whether the matrix is is singular. If we realize no solution exists, we quit to save some compute time. Because of this, the current implementation in the AAA is not necessarily incorrect; however, you guys are right in that the algorithm can be used in the more general case.

This puts us in an odd position, though. We can do a number of things:

1. Scrap the current chapter and re-write for the general case (as it seems you guys are suggesting?). I don't like this solution because the case we have should be close to the correct implementation for things later. We obviously want to stop the code as soon as no solution exists to save a bit of computational time.
2. Ignore the general case. Also obviously a bad idea, because then we are not able to find the rank of our matrix via gaussian elimination.
3. Combine the two cases. Because back-substitution doesn't make sense in the case of singular matrices, I suppose we could do the check after gaussian elimination has completed, but just before back substitution... but again, that's not a great solution either because then we cannot quit the computation early if no solution exists.
4. Add a few more sections for using Gaussian Elimination for certain cases and provide example code for each case.

I am leaning towards 4 right now, but would like to hear input... specifically because I am not sure exactly what implementations we need for what.

[jiegillet]

The best is 4. In general, almost nothing you state in the chapter is wrong so it doesn't need that much rewrite. Back substitution DOES make sense in the case of singular matrices, you just need to travel through the pivots, not the diagonal. Again, very little rewrite, I believe. And yes, putting a couple more examples of non-physical examples will help with intuition.

I'm very happy to try and add something to my PR tonight, I've already modified the chapter anyway.

[leios]

Yup, sounds good. Please feel free to do what you can when you have the time.

[jiegillet]

I've modified the chapter over in #193