forked from pymor/pymor
-
Notifications
You must be signed in to change notification settings - Fork 0
/
newton.py
139 lines (118 loc) · 5.13 KB
/
newton.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
# This file is part of the pyMOR project (http://www.pymor.org).
# Copyright 2013-2020 pyMOR developers and contributors. All rights reserved.
# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
import numpy as np
from pymor.core.defaults import defaults
from pymor.core.exceptions import InversionError, NewtonError
from pymor.core.logger import getLogger
@defaults('miniter', 'maxiter', 'rtol', 'atol', 'relax', 'stagnation_window', 'stagnation_threshold')
def newton(operator, rhs, initial_guess=None, mu=None, error_norm=None, least_squares=False,
miniter=0, maxiter=100, rtol=0., atol=0., relax=1.,
stagnation_window=3, stagnation_threshold=np.inf,
return_stages=False, return_residuals=False):
"""Basic Newton algorithm.
This method solves the nonlinear equation ::
A(U, mu) = V
for `U` using the Newton method.
Parameters
----------
operator
The |Operator| `A`. `A` must implement the
:meth:`~pymor.operators.interface.Operator.jacobian` interface method.
rhs
|VectorArray| of length 1 containing the vector `V`.
initial_guess
If not `None`, a |VectorArray| of length 1 containing an initial guess for the
solution `U`.
mu
The |parameter values| for which to solve the equation.
error_norm
The norm with which the norm of the residual is computed. If `None`, the
Euclidean norm is used.
least_squares
If `True`, use a least squares linear solver (e.g. for residual minimization).
miniter
Minimum amount of iterations to perform.
maxiter
Fail if the iteration count reaches this value without converging.
rtol
Finish when the residual norm has been reduced by this factor relative to the
norm of the initial residual.
atol
Finish when the residual norm is below this threshold.
relax
Relaxation factor for Newton updates.
stagnation_window
Finish when the residual norm has not been reduced by a factor of
`stagnation_threshold` during the last `stagnation_window` iterations.
stagnation_threshold
See `stagnation_window`.
return_stages
If `True`, return a |VectorArray| of the intermediate approximations of `U`
after each iteration.
return_residuals
If `True`, return a |VectorArray| of all residual vectors which have been computed
during the Newton iterations.
Returns
-------
U
|VectorArray| of length 1 containing the computed solution
data
Dict containing the following fields:
:error_sequence: |NumPy array| containing the residual norms after each iteration.
:stages: See `return_stages`.
:residuals: See `return_residuals`.
Raises
------
NewtonError
Raised if the Netwon algorithm failed to converge.
"""
logger = getLogger('pymor.algorithms.newton')
data = {}
if initial_guess is None:
initial_guess = operator.source.zeros()
if return_stages:
data['stages'] = operator.source.empty()
if return_residuals:
data['residuals'] = operator.range.empty()
U = initial_guess.copy()
residual = rhs - operator.apply(U, mu=mu)
err = residual.l2_norm()[0] if error_norm is None else error_norm(residual)[0]
logger.info(f' Initial Residual: {err:5e}')
iteration = 0
error_sequence = [err]
while True:
if iteration >= miniter:
if err <= atol:
logger.info(f'Absolute limit of {atol} reached. Stopping.')
break
if err/error_sequence[0] <= rtol:
logger.info(f'Prescribed total reduction of {rtol} reached. Stopping.')
break
if (len(error_sequence) >= stagnation_window + 1
and err/max(error_sequence[-stagnation_window - 1:]) >= stagnation_threshold):
logger.info(f'Error is stagnating (threshold: {stagnation_threshold:5e}, window: {stagnation_window}). '
f'Stopping.')
break
if iteration >= maxiter:
raise NewtonError('Failed to converge')
if iteration > 0 and return_stages:
data['stages'].append(U)
if return_residuals:
data['residuals'].append(residual)
iteration += 1
jacobian = operator.jacobian(U, mu=mu)
try:
correction = jacobian.apply_inverse(residual, least_squares=least_squares)
except InversionError:
raise NewtonError('Could not invert jacobian')
U.axpy(relax, correction)
residual = rhs - operator.apply(U, mu=mu)
err = residual.l2_norm()[0] if error_norm is None else error_norm(residual)[0]
logger.info(f'Iteration {iteration:2}: Residual: {err:5e}, '
f'Reduction: {err / error_sequence[-1]:5e}, Total Reduction: {err / error_sequence[0]:5e}')
error_sequence.append(err)
if not np.isfinite(err):
raise NewtonError('Failed to converge')
data['error_sequence'] = np.array(error_sequence)
return U, data