Householder's method for systems of nonlinear equations?

[avachon100501]

I found an interesting link here: https://www.researchgate.net/publication/264833135_Modified_of_Householder_Iterative_Method_for_Solving_Nonlinear_Systems
It could be one generalization of Householder's method for systems of nonlinear equations. I tried to put it in as a reply in YouTube, but the comment kept deleting itself... what's going on?

[osveliz]

I took a look at the paper and it doesn't appear to be the Generalized Householder Method that you were expecting. Rather this appears to be based on "two-step methods" that use one method to make a prediction and correct that prediction using another method, all as part of one iteration to increase the order. I would have expected something like Householder-3 that when generalized would require taking 3rd-order partial derivatives which the paper does not do.

YouTube sometimes considers posts with links as spam and automatically removes them. I never even see them in the held-for-review section of the comments.

[osveliz]

It certainly looks close if not slightly off from what someone would expect with the slant that I suspect should instead be symmetric. To sanity check, try Traub's system version of z^3-1 applied to your generalized Householder and compare it the regular Householder fractal of z^3-1. I believe that they would look incredibly similar if not identical. You could do the same with Halley.

[avachon100510]

(Sorry, it's me, but I had lost access to my previous account) Back to the topic of the method of tangent hyperbolas, I had found a way to solve the dimension mismatch: Remember that a_i was the Newton step? What if each of the elements of the vector was arranged into a diagonal n x n matrix, like so:

X_{n+1} = X_n - (F'(X_n) + 0.5*F''(X_n)*a)^-1*F(X_n)

Where a = diag(-F'(X_n)^-1*F(X_n))

The next part here is going to be the third order Householder method:

X_{n+1} = X_n - (F' + 0.5*F''*a + (1/6)*F'''*a^2)^-1*F

Or if you use both the Halley and Newton corrections,

X_{n+1} = X_n - (F' + b*(0.5*F'' + (1/6)*F'''*a))^-1*F

Where b = diag(-(F'(X_n) + 0.5*F''(X_n)*a)^-1*F(X_n))

We can generalize this further with the following (I only use the Newton correction in this part for the sake of simplicity):

X_{n+1} = X_n - (F' + 0.5*F''*a + (1/6)*F'''*a^2 + ... + (1/k!)*F^(k+1)(a^k))^-1*F

We finally have one version of a generalized kth order Householder method, and further solves the dimension mismatch paradox for this problem.

In other developments, I had just managed to learn how to compute the square root of a matrix.
I know you might say it's going to be secret, but what version of generalized bisection will it be showcased, in the condition that I won't spoil it to anyone else.

[osveliz]

Good that you were able to derive this. Have you tried making fractals?

I've been delayed in making videos. Trying to finalized the code for a version of Generalized Bisection but not sure if I have a bug or if there is something wrong with the method (or both). Video will discuss multiple versions of GB but only focus on one. If I can't get this to work, then I'll have to switch to another one.

[avachon100510]

Yes, I did some fractals! Here's one for the system:
F(x, y) = (x^3 y^3 - 1, x^3 y - 1)
The roots are: X1* = (1, 1) and X2* = (-1, -1)

[avachon100510]

The same system, but from the x-i plane:

And from the y-i plane:

[avachon100510]

Here's the Ultra Fractal code used for the images. Again, I used Cramer's rule, because it's simpler, especially on a computer.

2DHouseholder4 {
init:
  z = pixel
  if (@planetype == "xy")
  x = real(z)
  y = imag(z)
  i = sqrt(-1)
  elseif (@planetype == "xi")
  x = real(z)
  y = flip(imag(z))
  i = sqrt(-1)
  elseif (@planetype == "yi")
  x = flip(real(z))
  y = imag(z)
  i = sqrt(-1)
  endif
loop:
  complex func f1 (complex x, complex y)
    return x^3*y^3 - 1
  endfunc
  complex func f2 (complex x, complex y)
    return x^3*y - 1
  endfunc

  complex func f11 (complex x, complex y)
    return 3*x^2*y^3
  endfunc
  complex func f12 (complex x, complex y)
    return 3*x^3*y^2
  endfunc
  complex func f21 (complex x, complex y)
    return 3*x^2*y
  endfunc
  complex func f22 (complex x, complex y)
    return x^3
  endfunc
  complex func f111 (complex x, complex y)
    return 6*x*y^3
  endfunc
  complex func f112 (complex x, complex y)
    return 6*x^3*y
  endfunc
  complex func f121 (complex x, complex y)
    return 6*x*y
  endfunc
  complex func f122 (complex x, complex y)
    return 0
  endfunc
  complex func f1111 (complex x, complex y)
    return 6*y^3
  endfunc
  complex func f1122 (complex x, complex y)
    return 6*x^3
  endfunc
  complex func f2211 (complex x, complex y)
    return 6*y
  endfunc
  complex func f2222 (complex x, complex y)
    return 0
  endfunc

  complex func detA2x2 (complex a11, complex a12, complex a21, complex a22)
    return a11*a22-a12*a21
  endfunc


  a11 = -detA2x2(f1(x,y), f12(x,y), f2(x,y), f22(x,y))/detA2x2(f11(x,y), f12(x,y), f21(x,y), f22(x,y)), a12 = 0
  a21 = 0, a22 = -detA2x2(f11(x,y), f1(x,y), f21(x,y), f2(x,y))/detA2x2(f11(x,y), f12(x,y), f21(x,y), f22(x,y))

  a112 = a11*a11 + a12*a21, a122 = a11*a12 + a12*a22
  a212 = a21*a11 + a22*a21, a222 = a21*a12 + a22*a22

  xold = x
  yold = y

  d11 = f11(x,y)+0.5*(f111(x,y)*a11 + f112(x,y)*a21)+(1/6)*(f1111(x,y)*a112 + f1122(x,y)*a212), d12 = f12(x,y)+0.5*(f111(x,y)*a12 + f112(x,y)*a22)+(1/6)*(f1111(x,y)*a122 + f1122(x,y)*a222)
  d21 = f21(x,y)+0.5*(f121(x,y)*a11 + f122(x,y)*a21)+(1/6)*(f2211(x,y)*a112 + f2222(x,y)*a212), d22 = f22(x,y)+0.5*(f121(x,y)*a12 + f122(x,y)*a22)+(1/6)*(f2211(x,y)*a122 + f2222(x,y)*a222)

  x = x - @relax * (detA2x2(f1(x,y), d12, f2(x,y), d22)/detA2x2(d11, d12, d21, d22))
  y = y - @relax * (detA2x2(d11, f1(x,y), d21, f2(x,y))/detA2x2(d11, d12, d21, d22))
bailout:
  sqrt((x-xold)^2 + (y-yold)^2) >= @bailout || sqrt(f1(x,y)^2 + f2(x,y)^2) >= @bailout
default:
  title = "2D Householder 4"
  helpfile = "Uf*.chm"
  helptopic = "Html\formulas\standard\newton.html"
$IFDEF VER50
  rating = recommended
$ENDIF
  maxiter = 100
  center = (0, 0)
  param planetype
    caption = "Plane"
    default = 0
    enum = "xy" "xi" "yi"
    hint = "Changes the plane type."
  endparam
  param bailout
    caption = "Bailout value"
    default = 0.00001
    min = 0
$IFDEF VER40
    exponential = true
$ENDIF
    hint = "This parameter defines how soon a convergent orbit bails out while \
            iterating. Smaller values give more precise results but usually \
            require more iterations."
  endparam

  param relax
    caption = "Relaxation"
    default = (1, 0)
    hint = "This can be used to slow down or speed up the convergence of \
            the formula."
  endparam
}

[avachon100510]

Okay, so I finally figured it out. Householder's method, for a system of nonlinear equations, is obtained from an [1/d]-Padé approximant for the function g(y), where y = f(x) and g(y) = 1/f(x). That is, the reciprocal of f(x). In Cuyt's original paper, this expansion is expressed as a j-linear operator of the form E_j*a_i, where a_n = -F'(x_i) F(x_i), i.e. the Newton correction step. The third order Householder method is obtained from the [1/2]-Padé approximant of the Taylor series 1/f(x). The j-linear operator E_j is obtained, using Lagrange's inversion theorem via reciprocal bijection as:
E_n a_i^n = x_i + sum((-1)^(n - 1)/n! F'(x_i)^(1 - n)) F^(n)(x_i) a_i^n.
Where E_0 a_i^0 = x_i.
In Cuyt's paper, Halley's method is described as the [1/1]-abstract rational approximant, so naturally, the [1/2]-ARA would give:

x_{i + 1} = (E_2 a_i^2)/(E_3 a_i^3).

Halley's method, on the other hand, is as follows:
x_{i + 1} = (E_1 a_i^1)/(E_2 a_i^2)

It turns out that obtaining the abstract Padé approximant and the abstract rational approximant is done so in the same exact way as each other: the same process of solving for the coefficients of the numerator polynomial and the coefficients of the denominator polynomial.

The expression would be very long, and here, Fréchet derivatives are used.
Note that we're going to use tensors for higher-order partial derivatives, if we choose not to use Fréchet derivatives.