In the following demo, try dragging the black dots around. The image below will be recalculated to find a polynomial function that crosses each of these points. More or less points can be added using their respective buttons. Up to 75 points can be added.
I wrote an article back in June of 2009 about finding coefficients for a polynomial that will pass through n-points. Turns out this is covered in the second week of my linear algebra class, so not much of an accomplishment. However, I was using Cramer's Rule to solve the system of equations, which took (n+1)(2n!+n) multiplications, where n is the number of points. By the time n is 10 there are 79,833,710 multiplications to preform. I knew a better method to solve the system of equation must exist, but at the time I didn't really know much about it.
To start, take for instance the following system of equations, it's matrix equivalent, and it's form when placed in an augmented matrix:
This system can be solved by using row operations, and reducing the left side down to the identity matrix.
The last column holds the result.
This is probably the easiest system to do on paper when having to solve a system of linear equations. It also can be done in a spreadsheet, and I made just such a spreadsheet for assisting my class homework. One problem with the spreadsheet is that it can not row-swap. So sometimes the initial matrix has to have rows swapped so the spreadsheet can solve it. This isn't a problem when using the spreadsheet, but isn't great when trying to implement software to do row reduction. We could swap the rows in software, but that slows down everything.
The key to putting this system together was doing out-of-order row-reduction. We do not swap rows during the reduction. The outer loop goes through all but the last column (the column that has been augmented) of the matrix. This loop will first search for a row that has a term in this column (i.e. not zero).
Take for instance, the three matrices above. Normally, rows would be swapped in a matrix like B—row 1 and 3. The outer loop's first iteration will find row 1 for matrix A, row 2 for matrix B, and nothing for matrix C. If there is a result from the search, the row is normalized by dividing each term by the first term.
Above are the results after first term normalization. Matrix A was already normalized, matrix B row 2 is now normalized, and there wasn't a row found in matrix C. Once normalization, the active row is subtracted from all rows that have a leading term and have not been normalized.
The active row is marked as complete, and the loop continues onto the next column. The active row will be 2 for matrix A, 1 for matrix B and 1 for matrix C. Now the normalization happens again:
For subtraction, row 1 in matrix A is complete, so subtraction only needs to be done to row 3. For matrix B, row 2 is complete, so only row 3 remains. And in matrix C, rows 2 and 3 need subtraction because the initial pass did not reduce anything.
This process continues one more time for the 2nd to last column.
Notice that matrix A and C have multiple solutions due to the rows of zeros. This is perfectly alright. During this process, the order of the rows is saved. For matrix A, the order is [1,2,3]; matrix B [2, 1, 3] and matrix C [1,2,3]. This is needed to know which coefficient is represented by which row. One could find this be looking at all the proceeding zeros in the row, but that is time consuming—easier just to store the order in an array.
The last step is the second loop, where back substitution is preformed. The outer loop goes through the columns in reverse, and selects the row representing this column. So on the first iteration, this is row 3 for matrix all the matrices. Once selected, this row is subtracted from all remaining rows.
The next selected row will be 2 for matrix A and C, and 1 for matrix B.
Each of these matrices is now reduced. The resulting coefficients would be
Notice that the results for B have row 2 and row 1 swamped. This was accounted for in the order array, which for B was [2,1,3]. This entire process can be scaled to a matrix of any size.
Here is the updated function for finding the coefficients:
Looking at the loop in software, we can obtain the a set of summations that represent the complexity of the calculation. Using summation identities, it can then be reduced. This will give an loop maximum, as the loops may end early depending on the data in the equations being solved. I could step through the summation reductions, but it takes forever—I cheated and used software.
What we get is an equation that is the maximum number of operations needed to reduce any matrix of n rows and columns. For 10 points, that turns out to be 750 multiplications—100,000 times less then using Cramer's Rule. That is worst-case—the actual number of multiplications may be less.
The increased performance allows 50 points (85,750 multiplications) to be calculated with around the same complexity as 7 points in the old system. Even at 75 points (286,750 multiplications), the performance is good enough that the web server can handle it.
The draw back is that the floating point precision is good to around 16-points. Around this point the accuracy of the calculations becomes so poor due to accumulated rounding errors that the resulting polynomial begins to miss their targets. Nonetheless, I have allowed up to 75 points to be added to the demo despite the fact they will often be missed by the resulting polynomial.
In addition, there is more optimization that can be done to improve the speed at which the solutions can be derived. At some point, I plan on looking into implementing LU decomposition. This should further improve speeds. However, it would be necessary to use an arbitrary-precision arithmetic to make this worthwhile, since the precision of floating point is a much more limiting factor then the speed at which the matrices can be solved.