Yesterday I wrote about using the first derivative to create a curve fit where the peeks/valleys (the stationary points) were at the specified point locations. I could have released the demo with this, but I had the idea I could improve upon it.

 

Derivative on each point

Point 1:	Point 2:	Point 3:	Point 4:
Point 5:	Point 6:	Point 7:	Point 8:

In this demo, each point can use have the first and/or second derivative forced to zero at the point's location. Setting the first derivative to zero will cause the point to occur on a stationary point in the polynomial—that is, a point where the polynomial is at a local minimum or maximum. Setting the second derivative to zero will cause the point to have a local change in curvature.

This is done by allowing each point to specify additional restrictions to the overall system of equations. With no derivative, each point adds one term to the polynomial, and one equation to the system. Each derivative adds an other term, and an other equation. So a plot involving 3 points, the first with no derivative, the second with the first derivative, and the third with the second derivative will result in a polynomial with (1 + 2 + 3) = 5 terms and 5 equations.

Part of the solving takes advantage of this relationship:

Where is the m^th derivative of the function (this is Lagrange's notation). Using this, it is easy to make a function that can construct an equation (or a row in a matrix). The order in which the rows in the matrix isn't important, the construction loop for the matrix simply loops for the number of points, and then the number of derivatives for each point.

What this means is that we could have any number of points, with any number of derivatives. I limited the derivative number to 2, but I really didn't have a reason to even use the 2^nd derivative—it could easily go higher. Also, I have the derivatives being forced to zero. We could also specify a value for this derivative, although I haven't thought out how that could be useful yet.

As of this writing, I'm considering releasing a library with these functions in it. So for now, no source code—it's a mess and needs to be cleaned up.

One item I had read about polynomial curve fitting was using the derivative to make the points at the peeks of the function. We touched on this idea in my linear algebra class, but like most of the class, it was a vague mention. However, using my notes I was able to come up with what was being attempted.

The system of equations we are trying to solve for polynomial curve fitting is of this form:

Where n is the number of known points for which the polynomial is trying to pass through. As a consequence the resulting polynomial is always of degree n. That is, 4 known points will result in a polynomial .

If we desire the points to be at peeks/valleys in the curve, we need to go back to first semester calculus. Stationary points will occur when the derivative of the function is equal to zero. So we have a second equation:

This will double the degree polynomial as we will have twice the number of equations to solve. The result is the following matrix:

It looks a bit messy, but it isn't all that bad. The top half of the matrix is much like solving for an n-point polynomial, except there are twice the number of terms. The bottom half is the derivative of the top half.

This demo can show the effects of a linear approximation on trigonometric functions. The initial setup shows sine from 0 to 2 π made up of 20 points. The linear approximation is constructed using a 19^th degree polynomial from the data points. The coefficients are listed in the background.

Sine

Cosine

Sine centered

Cosine centered

 

Try removing points using the “Delete point” button. Notice how the right end stays around 0.5 until around 10 points. It has subtle variations which continue to get worse.

Now hit the “Sine centered->20 points” button. This will produce what looks to be a sine wave spanning 2 π. Now try “Sine centered->30 points.” Notice how deformed the sine wave has become before and after the points. Try the 10 point version.

One last thing to try is moving just slightly one point close to either end. The graph will not change much. Now try moving a point near the center slightly. The graph will change wildly near the ends, especially when a large number of points are used.

These are all effects of the linear approximation being used to represent the function.

This article a follow-up to my May 10, 11, and 12 articles. This demo will accurately find a polynomial crossing through up to 40 points. It uses PHP's BCMath Arbitrary Precision Mathematics functions to overcome the limitations of the floating point precision.

On May 12, I wrote that I had created this implementation, but that it was slow. What I found was that rending the polynomial plot was the slowest part of the system. After looking into this, I changed from using the “power-of” function to accumulating the power using multiplication. This noticeably improved the speed and I decided to add a few more fetchers.

I have allowed up to 40 points to be placed in this implementation, which takes the web server around 3 and a half seconds to render. The coefficients are display in the background of the plot to 30 digits, although the actual precision used might actually be higher.

The title of this article is a tribute to the demo group Triton. Early in their 1993 demo Crystal Dreams 2 is the introduction “now more linear algebra.”

Just some follow-ups to yesterday's story:

There are some interesting results for situations that are impossible. For example, two points at different y locations, but the same x location. Polynomials are functions, and thus this isn't allowed. Let's analysis one simple scenario.

Now expand the augmented matrix and the row-equivalent matrix for this scenario:

This system is inconsistent as row two is saying is 0 = 0.1—a clearly false statement. How the implementation deals with inconsistencies is by simply ignoring them. The resulting equation from this system is y = -2 x² + 1.8 x. There should be an x⁰ term for 3 points, but there isn't. One of the points is simply ignored—in this case, y₁.

The other variant of this is when two points are identical.

In this case, the matrix is consistent, but there are an infinite number of solutions. No points are ignored, but the polynomial has been reduced from a 2^nd order to 1^st order. Reading the last column, the equation y = 0 x² + 1 x or y = x. There is no x² coefficient.

This is similar to what happens if 3 separate points are defined, but for a line rather then a parabola.

This leads to our last interesting scenario, which is a special case of two identical x vales.

What happened? Earlier I claimed that a function could not have two different y values for a single x value or it wouldn't be a function. The simple answer: it's wrong—these coefficients do not pass through (0.5,0.6) and (0.5,0.5).

Our resulting equations is y = 66666675.05555565 x² - 53333338.23333337 x + 10000001. Using an arbitrary-precision arithmetic math system with x = 0.5, we get around 0.647222227499999. This is not .6 or .5 as expected. However, the coefficients are large enough that the function looks almost vertical around x=0.5. At x=0.4999, y≈-1332 and at x=0.5001, y≈1335. It is this transition that makes it look as though the functions is passing through both (0.5,0.6) and (0.5,0.5) when in reality, it does not pass through either.

This is a limitation of floating point. If we set the precision higher, to say 25 decimal places, y = 1000000000000000000000001-5333333333333333333333338.2333333333333333333333337 x + 6666666666666666666666675.0555555555555555555555565 x².

Remember that there is no solution, and the reason a solution is being found is because of rounding errors. However, the error looks close to the impossible situation we desire. This result is like the small-angle approximation where sines, cosines and tangents of small values are treated as the value rather then the function. An other method of doing something similar is to use the following values:

This will produce an almost identical graph, but the system has an exact solution. We engineers like this fact of math. I call it “close enough.”

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.

Now that I am in linear algebra, I decided to come back to this problem and use elementary row operations to reduce an augmented matrix in order to obtain the polynomial coefficients.

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 2^nd 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.