# [Scipy optimization tutorial] VI. root seeking, fixed point problem and its acceleration algorithm

Posted by abriggs on Sun, 27 Feb 2022 02:36:58 +0100

Refer to the official website: Scipy.

# Fixed point problem

The problem closely related to finding the zero point of a function is to find the fixed point of a function. The fixed point of a function refers to the point returned when the function is evaluated: g ( x ) = x g(x)=x g(x)=x. Obviously, this is f ( x ) = g ( x ) − x f(x)=g(x)-x The zero point problem of f(x)=g(x) − X. fixed_point provides a simple iterative method, which uses Aitkens sequence acceleration to estimate g ( x ) g(x) The fixed point of g(x).

```from scipy import optimize
import numpy as np
def func(x, c1, c2):
return np.sqrt(c1/(x+c2))
c1 = np.array([10,12.])
c2 = np.array([3, 5.])
ans = optimize.fixed_point(func, [1.2, 1.3], args=(c1,c2))
print(ans)
```

The result is:

```[1.4920333  1.37228132]
```

# Root finding of univariate equation

For example:

```import numpy as np
from scipy.optimize import root
def func(x):
return x + 2 * np.cos(x)
sol = root(func, 0.3)
print(sol.x)
```

The results are:

```[-1.02986653]
```

# Root finding of multivariable equation (method = 'lm')

For example:

```import numpy as np
from scipy.optimize import root
def func2(x):
f = [x[0] * np.cos(x[1]) - 4,
x[1]*x[0] - x[1] - 5]
df = np.array([[np.cos(x[1]), -x[0] * np.sin(x[1])],
[x[1], x[0] - 1]])
return f, df
sol = root(func2, [1, 1], jac=True, method='lm')
print(sol.x)
```

The results are:

```[6.50409711 0.90841421]
```

# Large scale root problem

hybr and lm cannot solve large-scale problems because they have to calculate n × N Jacobian matrix, with the increase of variables, the amount of calculation will be greater and greater.
For example, consider the following problems: we need to solve the following differential equations:

The limiting conditions are: the upper edge meets P ( x , 1 ) = 1 P(x,1)=1 P(x,1)=1, satisfied on the boundary P = 0 P=0 P=0. This can be achieved by approximating the continuous function p with the values on the grid (personal understanding is the finite difference method for fluid) P n , m ≈ P ( n h , m h ) P_{n,m}≈P(nh,mh) Pn,m​≈P(nh,mh). Then the derivative and integral can be approximated; For example:

In this way, large solvers can be used to solve: such as krylov, broyden2 and anderson. These methods use the so-called inexact Newton method, which does not accurately calculate the Jacobian matrix, but forms an approximate value.

```import numpy as np
from scipy.optimize import root
from numpy import cosh, zeros_like, mgrid, zeros

# parameters
nx, ny = 75, 75
hx, hy = 1./(nx-1), 1./(ny-1)

P_left, P_right = 0, 0
P_top, P_bottom = 1, 0

def residual(P):
d2x = zeros_like(P)
d2y = zeros_like(P)

d2x[1:-1] = (P[2:]   - 2*P[1:-1] + P[:-2]) / hx/hx
d2x[0]    = (P[1]    - 2*P[0]    + P_left)/hx/hx
d2x[-1]   = (P_right - 2*P[-1]   + P[-2])/hx/hx

d2y[:,1:-1] = (P[:,2:] - 2*P[:,1:-1] + P[:,:-2])/hy/hy
d2y[:,0]    = (P[:,1]  - 2*P[:,0]    + P_bottom)/hy/hy
d2y[:,-1]   = (P_top   - 2*P[:,-1]   + P[:,-2])/hy/hy

return d2x + d2y + 5*cosh(P).mean()**2

# solve
guess = zeros((nx, ny), float)
sol = root(residual, guess, method='krylov', options={'disp': True})
#sol = root(residual, guess, method='broyden2', options={'disp': True, 'max_rank': 50})
#sol = root(residual, guess, method='anderson', options={'disp': True, 'M': 10})
print('Residual: %g' % abs(residual(sol.x)).max())

# visualize
import matplotlib.pyplot as plt
x, y = mgrid[0:1:(nx*1j), 0:1:(ny*1j)]
plt.colorbar()
plt.show()
```

The result is:

```0:  |F(x)| = 40.1231; step 1
1:  |F(x)| = 16.821; step 1
2:  |F(x)| = 5.69025; step 1
3:  |F(x)| = 1.41135; step 1
4:  |F(x)| = 0.0386739; step 1
5:  |F(x)| = 0.00171956; step 1
6:  |F(x)| = 0.000132096; step 1
7:  |F(x)| = 6.2608e-06; step 1
8:  |F(x)| = 3.69188e-07; step 1
Residual: 3.69188e-07
```

# Acceleration method

For krylov method, most of its time is used to calculate the inverse of Jacobian matrix. If you have an approximation of the inverse matrix M ≈ J − 1 M≈J^{-1} M ≈ J − 1, then you can use it to preprocess the linear inversion problem.
We're not going to solve it J s = y Js=y Js=y, but solve it M J s = M y MJs=My MJs=My. Because the matrix is "closer" to the identity matrix than before, the Krylov method should be easier to deal with this equation.
Matrix M can be passed to Krylov method as an option options['jac_options']['inner_M '].
For the problem in the previous section, we notice that the function to be solved consists of two parts: the first part is Laplace operator and the second is integral. In fact, we can easily calculate the Jacobian coefficient corresponding to the Laplacian part:

J 1 J_1 J1. You can use SciPy sparse. linalg. Splu (or approximate with scipy.sparse.linalg.splu). So we use M ≈ J 1 − 1 M≈J^{-1}_1 M ≈ J1 − 1 ≈ is pretreated to accelerate the calculation.

```import numpy as np
from scipy.optimize import root
from scipy.sparse import spdiags, kron
from scipy.sparse.linalg import spilu, LinearOperator
from numpy import cosh, zeros_like, mgrid, zeros, eye

# parameters
nx, ny = 75, 75
hx, hy = 1./(nx-1), 1./(ny-1)

P_left, P_right = 0, 0
P_top, P_bottom = 1, 0

def get_preconditioner():
"""Compute the preconditioner M"""
diags_x = zeros((3, nx))
diags_x[0,:] = 1/hx/hx
diags_x[1,:] = -2/hx/hx
diags_x[2,:] = 1/hx/hx
Lx = spdiags(diags_x, [-1,0,1], nx, nx)

diags_y = zeros((3, ny))
diags_y[0,:] = 1/hy/hy
diags_y[1,:] = -2/hy/hy
diags_y[2,:] = 1/hy/hy
Ly = spdiags(diags_y, [-1,0,1], ny, ny)

J1 = kron(Lx, eye(ny)) + kron(eye(nx), Ly)

# Now we have the matrix `J_1`. We need to find its inverse `M` --
# however, since an approximate inverse is enough, we can use
# the *incomplete LU* decomposition

J1_ilu = spilu(J1)

# This returns an object with a method .solve() that evaluates
# the corresponding matrix-vector product. We need to wrap it into
# a LinearOperator before it can be passed to the Krylov methods:

M = LinearOperator(shape=(nx*ny, nx*ny), matvec=J1_ilu.solve)
return M

def solve(preconditioning=True):
"""Compute the solution"""
count = [0]

def residual(P):
count[0] += 1

d2x = zeros_like(P)
d2y = zeros_like(P)

d2x[1:-1] = (P[2:]   - 2*P[1:-1] + P[:-2])/hx/hx
d2x[0]    = (P[1]    - 2*P[0]    + P_left)/hx/hx
d2x[-1]   = (P_right - 2*P[-1]   + P[-2])/hx/hx

d2y[:,1:-1] = (P[:,2:] - 2*P[:,1:-1] + P[:,:-2])/hy/hy
d2y[:,0]    = (P[:,1]  - 2*P[:,0]    + P_bottom)/hy/hy
d2y[:,-1]   = (P_top   - 2*P[:,-1]   + P[:,-2])/hy/hy

return d2x + d2y + 5*cosh(P).mean()**2

# preconditioner
if preconditioning:
M = get_preconditioner()
else:
M = None

# solve
guess = zeros((nx, ny), float)

sol = root(residual, guess, method='krylov',
options={'disp': True,
'jac_options': {'inner_M': M}})
print('Residual', abs(residual(sol.x)).max())
print('Evaluations', count[0])

return sol.x

def main():
sol = solve(preconditioning=True)

# visualize
import matplotlib.pyplot as plt
x, y = mgrid[0:1:(nx*1j), 0:1:(ny*1j)]
plt.clf()
plt.pcolor(x, y, sol)
plt.clim(0, 1)
plt.colorbar()
plt.show()

if __name__ == "__main__":
main()
```

Good preprocessing may be crucial, and it can even determine whether the problem can be solved in practice.

Topics: Python Algorithm scipy