Variant of the Secant Method
12-09-2020, 04:33 AM (This post was last modified: 12-09-2020 09:13 PM by ttw.)
Post: #1
 ttw Member Posts: 244 Joined: Jun 2014
Variant of the Secant Method
The order of the secant method for finding zeros of a real function can be improved by using spline-like interpolants. This gives a method which has order just under 2 and uses only 1 function evaluation per step.

https://arxiv.org/pdf/2012.04248.pdf
12-09-2020, 05:40 AM
Post: #2
 Valentin Albillo Senior Member Posts: 777 Joined: Feb 2015
RE: Variant of the Second Method
(12-09-2020 04:33 AM)ttw Wrote:  The order of the secant method for finding zeros of a real function can be improved by using spline-like interpolants. This gives a method which has order just under 2 and uses only 1 function evaluation per step.

https://arxiv.org/pdf/2012.04248.pdf

Interesting. Perhaps you might consider implementing it yourself for some HP calc and telling us how it compares in terms of speed and memory versus plain old Newton's method {two evaluations of f(x) per iteration), which can be implemented in as few as 34 steps even in the simple HP-10C model (fewer in RPN models with subroutine capability like the HP-11C).

Anyway, thanks for sharing and take care.
V.

All My Articles & other Materials here:  Valentin Albillo's HP Collection

12-09-2020, 04:38 PM (This post was last modified: 12-09-2020 11:21 PM by Albert Chan.)
Post: #3
 Albert Chan Senior Member Posts: 1,676 Joined: Jul 2018
RE: Variant of the Second Method
Let's try solving x, for sinc(x) = 1/6.5

lua> function f(x) return sin(x)/x - 1/6.5 end
lua> function secant(x1,y1,x2,y2) return x1-y1/(y2-y1)*(x2-x1) end
lua> x1, x2 = 2.711, 2.712
lua> y1, y2 = f(x1) , f(x2)
lua> x3 = secant(x1,y1, x2,y2)
lua> y3 = f(x3)
lua> x3
2.7113129474671775

-- Inverse Interpolation, assumed f(x) is strictly increasing or descreasing around root
-- Fit x = g(y), then evaluate g(0)

lua> x12 = (x2-x1)/(y2-y1)
lua> x23 = (x3-x2)/(y3-y2)
lua> x123 = (x23-x12)/(y3-y1)
lua> x3 + (0-y3)*(x23 + (0-y2)*x123)
2.711312910356982

-- Aikten interpolation, as 1 liner: inner secant = x4, outer secant = x5
-- y1 x1
-- y2 x2 x3
-- y3 x3 x4 x5

lua> secant(secant(x3,y3, x1,y1),y3, x3,y2)
2.711312910356982

-- improved secant, interpolate for slope

lua> y23 = (y3-y2)/(x3-x2)
lua> y12 = (y2-y1)/(x2-x1)
lua> y123 = (y23-y12)/(x3-x1)
lua> slope = y23 + (x3-x2)*y123
lua> x3 - y3/slope
2.7113129103569826

-- Newton's method, actual slope

lua> x, y = x3, y3
lua> slope = (x*cos(x) - sin(x))/x^2
lua> x - y/slope
2.711312910356983

Comparing apples to apples, input is float(1/6.5) = 0x1.3b13b13b13b14p-3

5542891849071380 2 55 X^Y ÷ XEQ "ASINC"

2.71131291035698331469084928874754

All 3 methods gives excellent results, doubling accuracy of x3
12-09-2020, 11:51 PM
Post: #4
 Albert Chan Senior Member Posts: 1,676 Joined: Jul 2018
RE: Variant of the Secant Method
I extend previous post, with another round of iteration, using all points.
(First iteration matched previous post, done in double precision)

inverse-interpolation:
2.71131291035698220988289207111279
2.71131291035698331469084928874683

improved secant:
2.71131291035698244759230088334763
2.71131291035698331469084928874750

newton's method:
2.71131291035698307703744887554087
2.71131291035698331469084928874753
12-10-2020, 01:54 PM (This post was last modified: 12-10-2020 02:05 PM by Namir.)
Post: #5
 Namir Senior Member Posts: 813 Joined: Dec 2013
RE: Variant of the Secant Method
I implemented the Generalized Secant algorithm in Excel VBA. I used the same order as in the article's example (divided differences for first and second derivatives). The algorithm did better than Newton's method in most cases. The Generalized Secant algorithm needs two initial guesses, while Newton's method needs one. It made a difference in some test if I gave Newton's method the second guess OR averaged the two guesses. The closer of these two values to the neighboring root gave better results.

The algorithm is very interesting!!!

Thanks!

Namir
12-10-2020, 02:56 PM (This post was last modified: 12-11-2020 01:18 PM by Namir.)
Post: #6
 Namir Senior Member Posts: 813 Joined: Dec 2013
RE: Variant of the Secant Method
I was able to implement in Excel VBA a version of the Generalized Secant algorithm using divided differences for the first, second and third derivatives.

Here is the VBA listing for the version that uses two initial guesses and divided differences for the first and second derivatives.
Code:
 Cell Mapping =========== A1 : "X0" A2 : value for X0 A3 : "X1" A4 : value for X1 A5 : "Tolerance" A6 : value for tolerance A7 : "F(x)" A8: String containing expression for f(x)=0, such as exp(x)-3*x^2. B1 : "X" C1 : "F(X)"  (for Generalized Secant method) D1 : "X" E1 : "F(X)"  (for Newton's method, used for comparison) VBA Listing =========== Function Fx(ByVal sFx As String, ByVal X As Double) As Double   sFx = Replace(sFx, "EXP", "!")   sFx = Replace(sFx, "X", "(" & CStr(X) & ")")   sFx = Replace(sFx, "!", "EXP")   Fx = Evaluate(sFx) End Function Sub genSecant2()   Dim sFx As String, X0 As Double, X1 As Double, X2 As Double, X3 As Double   Dim Toler As Double, Df1 As Double, Df2 As Double   Dim Row As Integer, h As Double   Dim Fx0 As Double, Fx1 As Double, Fx2 As Double, Diff As Double      Range("B2:X1000").Clear      X0 = [A2].Value   X1 = [A4].Value   Toler = [A6].Value   sFx = [A8].Value   sFx = Replace(sFx, " ", "")   sFx = UCase(sFx)      Fx0 = Fx(sFx, X0)   Fx1 = Fx(sFx, X1)   Df1 = Fx0 / (X0 - X1) + Fx1 / (X1 - X0)   X2 = X1 - Fx1 / Df1   Fx2 = Fx(sFx, X2)   Diff = X2 - X1      Row = 2   Cells(Row, 2) = X2   Cells(Row, 3) = Fx(sFx, X2)   Row = Row + 1   Do Until Abs(Diff) < Toler     Df1 = Fx2 / (X2 - X1) + Fx1 / (X1 - X2)     Df2 = Fx2 / (X2 - X1) / (X2 - X0) + Fx1 / (X1 - X2) / (X1 - X0) + Fx0 / (X0 - X2) / (X0 - X1)     Diff = Fx2 / (Df1 + Df2 * (X2 - X1))     X3 = X2 - Diff     X0 = X1     Fx0 = Fx1     X1 = X2     Fx1 = Fx2     X2 = X3     Fx2 = Fx(sFx, X3)     Cells(Row, 2) = X3     Cells(Row, 3) = Fx2     Row = Row + 1   Loop      ' Newton's method   X1 = [A4].Value   Fx1 = Fx(sFx, X1)   Row = 2   Diff = 2 * Toler   Do Until Abs(Diff) < Toler     h = 0.01 * (1 + Abs(X1))     Diff = h * Fx1 / (Fx(sFx, X1 + h) - Fx1)     X1 = X1 - Diff     Fx1 = Fx(sFx, X1)     Cells(Row, 4) = X1     Cells(Row, 5) = Fx1     Row = Row + 1   Loop    End Sub Output ====== Columns B and C show the refined guesses and their function values for the Generalized Secant Algorithm. Columns D and E show the refined guesses and their function values for Newton's method.

Here is the version of the VBA code that uses three initial guesses and the divided differences for the first three derivatives:

Code:
 Cell Mapping =========== A1 : "X0" A2 : value for X0 A3 : "X1" A4 : value for X1 A5 : "X2" A6 : value for X2 A7 : "Tolerance" A8 : value for tolerance A9 : "F(x)" A10: String containing expression for f(x)=0 B1 : "X" C1 : "F(X)"  (for Generalized Secant method) D1 : "X" E1 : "F(X)"  (for Newton's method, used for comparison) VBA Listing =========== Function Fx(ByVal sFx As String, ByVal X As Double) As Double   sFx = Replace(sFx, "EXP", "!")   sFx = Replace(sFx, "X", "(" & CStr(X) & ")")   sFx = Replace(sFx, "!", "EXP")   Fx = Evaluate(sFx) End Function Sub genSecant3()   Dim sFx As String, X0 As Double, X1 As Double, X2 As Double, X3 As Double, X4 As Double   Dim Toler As Double, Df1 As Double, Df2 As Double, Df3 As Double   Dim Row As Integer, h As Double   Dim Fx0 As Double, Fx1 As Double, Fx2 As Double, Fx3 As Double, Diff As Double      Range("B2:X1000").Clear      X0 = [A2].Value   X1 = [A4].Value   X2 = [A6].Value   Toler = [A8].Value   sFx = [A10].Value   sFx = Replace(sFx, " ", "")   sFx = UCase(sFx)      Fx0 = Fx(sFx, X0)   Fx1 = Fx(sFx, X1)   Fx2 = Fx(sFx, X2)   Df1 = Fx0 / (X0 - X1) + Fx1 / (X1 - X0)   Df2 = Fx2 / (X2 - X1) / (X2 - X0) + Fx1 / (X1 - X2) / (X1 - X0) + Fx0 / (X0 - X2) / (X0 - X1)   X3 = X2 - Fx2 / (Df1 + Df2 * (X2 - X1))   Fx3 = Fx(sFx, X3)   Diff = X3 - X2      Row = 2   Cells(Row, 2) = X3   Cells(Row, 3) = Fx3   Row = Row + 1   Do Until Abs(Diff) < Toler Or Row > 100     Df1 = Fx3 / (X3 - X2) + Fx2 / (X2 - X3)     Df2 = Fx3 / (X3 - X1) / (X3 - X2) + Fx2 / (X2 - X1) / (X2 - X3) + Fx1 / (X1 - X2) / (X1 - X3)     Df3 = Fx3 / (X3 - X2) / (X3 - X1) / (X3 - X0) + _           Fx2 / (X2 - X3) / (X2 - X1) / (X2 - X0) + _           Fx1 / (X1 - X3) / (X1 - X2) / (X1 - X0) + _           Fx0 / (X0 - X3) / (X0 - X2) / (X0 - X1)     Diff = Fx3 / (Df1 + Df2 * (X3 - X2) + Df3 * (X3 - X2) * (X3 - X1))     X4 = X3 - Diff     X0 = X1     Fx0 = Fx1     X1 = X2     Fx1 = Fx2     X2 = X3     Fx2 = Fx3     X3 = X4     Fx3 = Fx(sFx, X4)     Cells(Row, 2) = X4     Cells(Row, 3) = Fx3     Row = Row + 1   Loop      ' Newton's method   X1 = [A4].Value   Fx1 = Fx(sFx, X1)   Row = 2   Diff = 2 * Toler   Do Until Abs(Diff) < Toler     h = 0.01 * (1 + Abs(X1))     Diff = h * Fx1 / (Fx(sFx, X1 + h) - Fx1)     X1 = X1 - Diff     Fx1 = Fx(sFx, X1)     Cells(Row, 4) = X1     Cells(Row, 5) = Fx1     Row = Row + 1   Loop    End Sub Output ====== Columns B and C show the refined guesses and their function values for the Generalized Secant Algorithm. Columns D and E show the refined guesses and their function values for Newton's method.

Use different worksheets for each of the above methods as the input values in column A are not the same.

I wrote the above code to make it easy to ready. A programmer can optimize the above code by using arrays, matrices, and loops.
12-11-2020, 04:34 PM (This post was last modified: 02-07-2021 11:57 AM by Albert Chan.)
Post: #7
 Albert Chan Senior Member Posts: 1,676 Joined: Jul 2018
RE: Variant of the Secant Method
Quadratic interpolation is simply 2 linear interpolations, followed by linear interpolation.

x3 y3
x2 y2 y23
x1 y1 y13 y123

Similarly, cubic interpolation is 2 quadratic interpolations, followed by linear interpolation.

x4 y4
x3 y3
x2 y2 y234
x1 y1 y134 y1234

Using this recursive nature of polynomial interpolation, I get this neat Lua code

Code:
function interp(x1,y1,x2,y2,x3,y3,x4,y4)    -- estimate y = p(0)     if x3 then         y2 = interp(x2,y2,x3,y3,x4,y4)         y1 = interp(x1,y1,x3,y3,x4,y4)     end     local y, p = y2, x2/(x2-x1)     if p > 0.5 then y, p = y1, x1/(x2-x1) end     return y - p*(y2-y1) end

Extending interp(...) is trivial, by adding (x5,y5), ...
The down-side is each bump up in polynomial degree, doubled the work.
For high degree fit, divided-difference style is better, because it can cache previous results.

OTTH, assuming we have convergence, old points are very, very far away.
Fitting those probably may not help much ...

---

For slope of improved secant, we can also consider this an interpolation problem.
Example, for a cubic fit, we generate divided-difference table:

x4 y4
x3 y3 y34
x2 y2 y24 y234
x1 y1 y14 y123 y1234

p(x) = y4 + (x-x4)*(y34 + (x-x3)*(y234 + (x-x2)*y1234)) = y4 + (x-x4)*M

p'(x) = M + (x-x4)*M' ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → p'(x4) = M

Basing all values relative to (x4, y4). Example X3 = x3-x4, Y3 = y3-y4 ...
M is just interpolation of individual slopes.

X3 ﻿ ﻿ ﻿ ﻿ ﻿ Y3/X3
X2 ﻿ ﻿ ﻿ ﻿ ﻿ Y2/X2
X1 ﻿ ﻿ ﻿ ﻿ ﻿ Y1/X1 ﻿ ﻿ ﻿ ﻿ ﻿ M

Again, extending this is simply adding more slopes to fit

Code:
function secant(x1,y1,x2,y2,x3,y3,x4,y4)    -- estimate root of p(x)     if not x3 then return interp(y1,x1,y2,x2) end     if x4 then x3,y3,x4,y4 = x4,y4,x3-x4,(y3-y4)/(x3-x4) end     x1, y1, x2, y2 = x1-x3, y1-y3, x2-x3, y2-y3     return x3 - y3 / interp(x1,y1/x1, x2,y2/x2, x4,y4)     end

Example:

lua> function f(x) return sin(x)/x - 1/6.5 end
lua> x1, x2 = 2.711, 2.712
lua> y1, y2 = f(x1), f(x2)
lua> x3 = secant(x1,y1, x2,y2)
lua> y3 = f(x3)

lua> interp(y1,x1, y2,x2, y3,x3) -- inverse interpolate, x = g(0)
2.711312910356982
lua> secant(x1,y1, x2,y2, x3,y3) -- improved secant, x = x3 - y3/slope
2.7113129103569826
12-11-2020, 04:46 PM (This post was last modified: 02-07-2021 11:59 AM by Albert Chan.)
Post: #8
 Albert Chan Senior Member Posts: 1,676 Joined: Jul 2018
RE: Variant of the Secant Method
Let's compare inverse-interpolation vs improved-secant (code from previous post)

Code:
function test(f, x0, x1, n)   -- n = 1 (linear) to 3 (cubic)     local y0, y1 = f(x0), f(x1)     local t0 = {y0,x0, y1,x1} -- inverse interpolate     local t1 = {x0,y0, x1,y1} -- improved secant     n = (n or 3)*2 + 1     for i = 4,18,2 do         -- i = table length         x0 = interp(unpack(t0, math.max(1,i-n)))         x1 = secant(unpack(t1, math.max(1,i-n)))         t0[i+1], t0[i+2] = f(x0), x0         t1[i+1], t1[i+2] = x1, f(x1)         print(i/2, x0, x1)     end end

lua> function f(x) return x^3 - 8 end

lua> test(f, 5, 4, 2) -- https://arxiv.org/pdf/2012.04248.pdf, Table 2
2 ﻿ ﻿ ﻿ ﻿ ﻿ 3.081967213114754 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 3.081967213114754
3 ﻿ ﻿ ﻿ ﻿ ﻿ 2.3945608933295457 ﻿ ﻿ ﻿ ﻿ ﻿ 2.2862188297178117
4 ﻿ ﻿ ﻿ ﻿ ﻿ 2.0980163342039635 ﻿ ﻿ ﻿ ﻿ ﻿ 2.010344209437878
5 ﻿ ﻿ ﻿ ﻿ ﻿ 2.0088865345029276 ﻿ ﻿ ﻿ ﻿ ﻿ 1.99979593345267
6 ﻿ ﻿ ﻿ ﻿ ﻿ 2.0001129292461193 ﻿ ﻿ ﻿ ﻿ ﻿ 2.0000000722313933
7 ﻿ ﻿ ﻿ ﻿ ﻿ 2.0000000388749792 ﻿ ﻿ ﻿ ﻿ ﻿ 2.000000000000015
8 ﻿ ﻿ ﻿ ﻿ ﻿ 2.0000000000000164 ﻿ ﻿ ﻿ ﻿ ﻿ 2
9 ﻿ ﻿ ﻿ ﻿ ﻿ 2 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 2

Last column is improved secant, fitting data up to quadratic (slope by linear interpolation)
Middle column is inverse-interpolation for x, at y = 0, also fitting up to quadratic.

Improved secant uses the last point as the "base", then does correction.
Inverse-interpolation is not as good, probably due to putting equal weight to its points.

Let's try again, but with up to cubic fit.

lua> test(f, 5, 4, 3)
2 ﻿ ﻿ ﻿ ﻿ ﻿ 3.081967213114754 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 3.081967213114754
3 ﻿ ﻿ ﻿ ﻿ ﻿ 2.3945608933295457 ﻿ ﻿ ﻿ ﻿ ﻿ 2.2862188297178117
4 ﻿ ﻿ ﻿ ﻿ ﻿ 2.082744738020821 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 2.034337291023909
5 ﻿ ﻿ ﻿ ﻿ ﻿ 2.0058425366743506 ﻿ ﻿ ﻿ ﻿ ﻿ 2.000576313421517
6 ﻿ ﻿ ﻿ ﻿ ﻿ 2.0000383915262443 ﻿ ﻿ ﻿ ﻿ ﻿ 2.000000166004798
7 ﻿ ﻿ ﻿ ﻿ ﻿ 2.0000000023200664 ﻿ ﻿ ﻿ ﻿ ﻿ 2.0000000000000138
8 ﻿ ﻿ ﻿ ﻿ ﻿ 1.9999999999999998 ﻿ ﻿ ﻿ ﻿ ﻿ 2
9 ﻿ ﻿ ﻿ ﻿ ﻿ 2 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 2

Based on above, anything more than a cubic fit probably not worth the trouble.
For example, lets see how Newton's method does (actual slope)

lua> function df(x) return 3*x*x end
lua> x = 4
lua> for i=2, 9 do x = x - f(x)/df(x); print(i, x) end
2 ﻿ ﻿ ﻿ ﻿ ﻿ 2.833333333333333
3 ﻿ ﻿ ﻿ ﻿ ﻿ 2.221068819684737
4 ﻿ ﻿ ﻿ ﻿ ﻿ 2.0212735368091126
5 ﻿ ﻿ ﻿ ﻿ ﻿ 2.0002231146078984
6 ﻿ ﻿ ﻿ ﻿ ﻿ 2.0000000248863623
7 ﻿ ﻿ ﻿ ﻿ ﻿ 2.0000000000000004
8 ﻿ ﻿ ﻿ ﻿ ﻿ 2
9 ﻿ ﻿ ﻿ ﻿ ﻿ 2
12-12-2020, 04:22 AM
Post: #9
 Namir Senior Member Posts: 813 Joined: Dec 2013
RE: Variant of the Secant Method
I went back to my Excel VBA code and added statements to implement Halley's method. The Generalized Secant method did better than Newton's method but ws behind Halley's method in the number of steps. Of course results depend on the function tested and the initial guess.

Namir
12-12-2020, 02:41 PM
Post: #10
 ttw Member Posts: 244 Joined: Jun 2014
RE: Variant of the Secant Method
Methods with high convergence can have problems numerically; they also tend to have a smaller domain of convergence. Multi-point methods seem to be a bit more stable than high order methods.
12-12-2020, 04:27 PM
Post: #11
 Albert Chan Senior Member Posts: 1,676 Joined: Jul 2018
RE: Variant of the Secant Method
(12-12-2020 02:41 PM)ttw Wrote:  Methods with high convergence can have problems numerically; they also tend to have a smaller domain of convergence.

This might not always true.
In fact, recently we come across opposite case, LambertW function implementation.

Solving f(x) = x*e^x - a, with Newton's method is very unstable.

Higher convergence method, say Halley, it becomes "safe"

P.S. the next post, I had suggested another way, to suppress the exponential.
This version, Newton's method is fast, and safe.

>>> from mpmath import *
>>> a = -1e99
>>> findroot(lambda x: x + ln(x/a), 200+3j) # bad guess
mpc(real='222.55067066150301', imag='3.1275404175517207')
>>> lambertw(a)
mpc(real='222.55067066150301', imag='3.127540417551721')
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)