Iterative vector valued non-linear equation solver using Newton's method

This embedded spreadsheet solves three non-linear equations in three unknowns using up to 10 iterations of the vector version of Newton' method. This animation shows the scaler version of Newton's method at work for y = 0. To try it out simply enter a 0 into the yellow reset cell (C1) at the top. Enter a 1 back in this same cell to reset the spreadsheet. Use is made of the circular reference facility in Excel to do the iterations: when reset (yellow cell C1) is not zero then x.2 is a function of x.1 and x.1 is a copy of initial guess x.0 (a user set constant), but when reset = 0 then x.1 is instead a copy of x.2, thus the circular reference. 10 was the predefined maximum number of iterations selected. Instructions are on the "Solutions" sheet. The other two sheets define the non-linear vector valued function (f) and some constants respectively. To see it in action, set cell C1 (reset) to 0 and the estimate should appear in x.1 (cells C11:C14) after up to 10 (unseen) iterations. Compare with the true answer (x_true) in cells B5:B7. To reset it to the initial guess x0 (A11:A13), set reset  back to a non-zero value. Try changing any of the orange or yellow cells on the Solutions sheet or the NonLinFunct sheet. Chances are you can break it since Newton's method is famous for converging quickly near the solution, but diverging if too far away (should the equations being solved are too non-linear). If after iterating (setting reset to 0), you see that "err" (C16:C18), representing y - f(x.1), is non-zero, you can attempt to improve the estimate by entering any change to any of the unused white cells (which causes the sheet to automatically recompute).

Even though you can't see the iterations themselves (since they are all completed before Excel redraws the spreadsheet), you can still get a sense for the iterations required to converge by examining the default sheet when reset = 1. Refresh the page and then compare the orange "x_true" 3x1 vector (the true solution) with the green "x.2 (next x.1)" 3x1 vector below and to the right of x_true. If f() were perfectly linear  then x.2 would equal x_true, which it clearly does not: thus we know that more than one iteration will be required to get us close with a non-linear f(). Now while still in reset (reset = 1) try setting the non-linear component of f() to zero by setting the orange cell C2 ("NLmult: eps.B") to 0. Then you should see x.2 be equal (or very close) to x_true, indicating that just one iteration is required. Enjoy!



You might ask why I don't use more sophisticated built in Excel tools like "Goal Seek" or "Solver," perhaps in a Visual Basic macro, or even just a specialized macro on it's own. The answer is that those tools are not available in the stripped down online version of Excel available for free on Microsoft's OneDrive website (which is where I created this). You can, however, initiate an off-line spreadsheet using the full version of Excel, set the check box to allow circular references (in the options) along with a change threshold to stop iterations (which I made very small) and a maximum number of iterations to perform (100 by default). Once uploaded to the OneDrive site, those setting cannot be changed, but the sheet still obeys them. The alternative for an online spreadsheet is to fix the number of iterations by breaking out each iteration into its own set of cells (like I did here), which gets ugly fast.

Also note that calculating a matrix inverse is not the way to go in general. I had a professor once in grad school (whose specialty was finding numerically stable methods of solving matrix equations, especially those important in optimal feedback control problems, such as the matrix Ricatti equation) tell us that he'd never seen a legitimate reason to calculate a matrix inverse, and if we ever thought we found one to please contact him! Basically you always want the inverse multiplied by something else, i.e. you don't want inverse(A) = A^-1, you want (A^-1)*x, for example. The downside to calculating the inverse is two-fold:
  1. It's inefficient and unnecessary
  2. It can be numerically unstable
Of these two, 2 is by far more important for an application like this (since efficiency is not a big concern). Why can it be numerically unstable? Because matrix A might be close to singular (i.e., non-invertable: for example the 1x1 matrix [0] is singular), and when that happens, some of its inverse elements get very large, and then you're probably relying on differences of large number to be accurate, which is always a bad idea (e.g. your machine probably isn't using enough digits of precision to accurately return the answer "1" when computing (1e20 + 1) - 1e20). In Matlab, for example, you would avoid calculating the inverse when solving a problem like this for x:

y = A*x

By entering this:

x = A\y

instead of this:

x = inverse(A)*y

A\y is equivalent to inverting A (when A is non-singular) and multiplying it by y, but in a numerically stable way. Either a QR decomposition is used or a pseudo-inverse is used. Essentially it's solving a least squares problem, but in a case like this, it is a fully determined system (square non-singular A) rather than the usual overdetermined system (A has more rows than linearly independent columns) associated with typical least squares (linear regression) problems, like fitting a curve to a data set. A pseudo-inverse is even well defined for an under-determined system (more columns than rows). In either case it's a better idea than inverting A. Excel has a built in way to do this with their linear regression functions, but they are not convenient to use. Plus I have no idea what's going on inside them: you might hope they'd use a good method, but they might actually be forming a Gramian matrix, which is a really bad idea, since you're essentially squaring the elements and further ruining your chances for numerical stability. Bizarrely enough, the most inconvenient thing about Excel's linear regression tools is the way the answer is returned: as a row vector (not that bad because of the built in TRANSPOSE() function) with the elements in reverse order. You'd think that wouldn't be a problem, but to reverse the order of an array in Excel automatically is NOT easy!... once you figure out the arcane and frustrating way to use Excel functions to do that, I guarantee you'll forget them again in five minutes... even doing it once using contradictory online help pages for inspiration, it took me over an hour! (the worst part was getting the unhelpful "VALUE!" message with no explanation of why... and once I figured out why I was near homicidal)... absolutely stupid! You'd think that would be a built in function. But then again, you might also expect MLEFTDIVIDE() might be a built in function to mimic Matlab's left divide operation (\). Of course you can code up your own macro (again, not available online) or just reverse the elements by hand. Anyway, for this example it was easier just to use MINVERSE() and hope for the best: it's by far the clearest way to go, and that's what I was shooting for here. Once some other bugs are fixed in Excel having to do with using other matrix manipulations without having to populate cells with the intermediate answer, and using the new (in 2013) MUNIT() function to create an identity matrix, I think I can eliminate almost all the intermediate results above in the light purple cells (and certainly the ones in the gray "junk" cells, which are there ONLY because of bugs in Excel).

If you mess up the spreadsheet and want to start over again, just reload the page. You can download a copy of the spreadsheet using the green "X" Excel symbol on the right hand side of the black bar across the bottom of the spreadsheet. You can also see the full spreadsheet (including the NonLinFunct sheet defining the example non-linear vector valued function f() I'm using here) and all the cell comments by clicking on the "View full-size workbook" button in the extreme lower right corner.

No comments :

Post a Comment