Newton-Krylov method - Revision history

Apwillis at 15:50, 1 February 2023

2023-02-01T15:50:54Z

← Older revision		Revision as of 08:50, 1 February 2023
Line 2:		Line 2:

	'''This page and the codes have been moved to https://github.com/apwillis1/JFNK-Hookstep '''		'''This page and the codes have been moved to https://github.com/apwillis1/JFNK-Hookstep '''

			Codes are in MATLAB (Octave compatible) and Fortran90.
			They can be integrated with codes in other languages via command-line calls.

Apwillis at 15:07, 1 February 2023

2023-02-01T15:07:59Z

← Older revision		Revision as of 08:07, 1 February 2023
Line 1:		Line 1:
	=== Jacobian-Free Newton-Krylov (JFNK) method ===		=== Jacobian-Free Newton-Krylov (JFNK) method ===

	This page and the codes ~~are now available at~~ https://github.com/apwillis1/JFNK-Hookstep		'''This page and the codes have been moved to https://github.com/apwillis1/JFNK-Hookstep '''

Apwillis at 15:06, 1 February 2023

2023-02-01T15:06:52Z

← Older revision		Revision as of 08:06, 1 February 2023
Line 1:		Line 1:
	=== Jacobian-Free Newton-Krylov (JFNK) method ===		=== Jacobian-Free Newton-Krylov (JFNK) method ===

	This page ~~has been moved to~~ https://github.com/apwillis1/JFNK-Hookstep		This page and the codes are now available at https://github.com/apwillis1/JFNK-Hookstep

Apwillis: Replaced content with "=== Jacobian-Free Newton-Krylov (JFNK) method === This page has been moved to https://github.com/apwillis1/JFNK-Hookstep"

2023-02-01T15:04:44Z

Replaced content with "=== Jacobian-Free Newton-Krylov (JFNK) method === This page has been moved to https://github.com/apwillis1/JFNK-Hookstep"

← Older revision		Revision as of 08:04, 1 February 2023
Line 1:		Line 1:
	~~{{latexPreamble}}~~		=== Jacobian-Free Newton-Krylov (JFNK) method ===
	Free ~~to download:~~

	~~'''CODE (MATLAB / GNU Octave): Template-Example: [[File:Newton_Lorenz_m.zip]] / [[File:Newton_Lorenz_m.tgz]]'''~~		This page has been moved to https://github.com/apwillis1/JFNK-Hookstep

	~~'''CODE (FORTRAN): Template-Example: [[File:Newton_Lorenz.f90]].'''~~

	~~'''Info on main subroutine:''' [[File:NewtonHook.f90]] (essentially same for FORTRAN/MATLAB versions)~~

	~~___ __~~

	~~'''(An Extended overview of the Newton-Krylov method is here (pdf,arxiv): [https://arxiv.org/pdf/1908.06730.pdf] See section 4 on using the code.)'''~~

	~~___NOTOC___ __~~



	~~The codes above implement the '''Jacobian-free Newton-Krylov (JFNK) method''' for solving~~
	~~\begin{equation}\label{eq:ROOTF}\vec{F}(\vec{x})=\vec{0},\end{equation}~~
	~~where $\vec{x}$ and $\vec{F}$ are $n$-vectors, '''supplemented with a Hookstep--Trust-region approach'''.~~

	This is a powerful method that can solve for $\vec{x}$ for a complicated nonlinear $\vec{F}(\vec{x})$. For example, to find an ''equilibrium'' solution or a ''periodic orbit'', let $\vec{F}(\vec{x}) = \vec{X}(\vec{x})-\vec{x}$, where $\vec{X}(\vec{x})$ is the result of time integration of an initial condition $\vec{x}$.

	~~==== Newton-Raphson method ====~~

	To find the roots $x$ of a function $f(x)$ in one dimension, given an initial guess $x_0$, the Newton-Raphson method generates improvements using the iteration $x_{i+1}=x_i-f(x_i)/f'(x_i)$, where the dash denotes the derivative. The iteration can be re-expressed as
	\[
	~~x_{i+1} = x_i + \delta x_i \,~~
	~~\quad\mbox{where}\quad~~
	~~f'(x_i)\,\delta x_i = -f(x_i) \, .~~
	\]
	~~The extension of Newton's method to an $n$-dimensional system is then~~
	~~\begin{equation} \label{eq:N-R}~~
	~~(a)\quad \vec{x}_{i+1} = \vec{x}_i + \vec{\delta x}_i \,~~
	~~\quad\mbox{where}\quad~~
	~~(b)\quad\left. \frac{\vec{\partial F}}{\vec{\partial x}}\right\|_{\vec{x}_i}~~
	~~\vec{\delta x}_i = -\vec{F}(\vec{x}_i) \, .~~
	~~\end{equation}~~
	~~In order to apply the update $(a)$, the linear system $(b)$ needs to be solved for the unknown $\vec{\delta x}_i$.~~

	~~==== Jacobian-Free Newton method ====~~

	The $n\times n$ Jacobian matrix $\left.\frac{\vec{\partial F}}{\vec{\partial x}}\right\|_{\vec{x}_i}$ is usually difficult to evaluate. However, the problem $(b)$ is in the form $A\,\vec{\delta x}=\vec{b}$, which can be solved using the '''Krylov-subspace method GMRES(m)'''. (See [[File:GMRESm.f90]].) The GMRES algorithm does not need to know the matrix $A$ itself, only the result of multiplying a vector by $A$. For a given starting vector $\vec{\delta x}_0$, e.g. $\vec{\delta x}_0=\vec{b}$, the method seeks solutions for $\vec{\delta x}$ in $\mathrm{span}\{\vec{x}_0,\,A\vec{x}_0,\,A^2\vec{x}_0,...\}$, but uses Gram-Schmidt orthogonalisation to improve the numerical suitability of this basis set. The set of orthogonalised vectors is called the Krylov-subspace, and m is the maximum number of vectors stored.

	~~Iterations of the GMRES algorithm for the problem $(b)$ involve calculating products of the Jacobian with given $\vec{\delta x}$, which may be approximated, e.g.~~
	~~\begin{equation}\label{eq:Japprox}~~
	~~(c)\qquad~~
	~~\left.\frac{\vec{\partial F}}{\vec{\partial x}}\right\|_{\vec{x}_i}\,\vec{\delta x} \approx~~
	~~\frac{1}{\epsilon}(\vec{F}(\vec{x}_i+\epsilon\,\vec{\delta x})-\vec{F}(\vec{x}_i))~~
	~~\end{equation}~~
	~~for some small value $\epsilon$. (Try $\epsilon$ such that $(\epsilon\|\|\vec{\delta x}\|\|)\,/\,\|\|\vec{x}_i\|\|=10^{-6}$.) The important point is that we do not need~~ to ~~know the Jacobian; '''only a routine for evaluating $\vec{F}(\vec{x})$ is required'''.~~

	Note that provided that each step of the Newton method, $\vec{\delta x}$, is essentially in the correct direction, the method is expected to converge. Therefore the tolerance specified in the accuracy of the solution for $\vec{\delta x}$ in each Newton step (calculated via the GMRES method) typically need not be so stringent as the tolerance placed on the Newton method itself for the solution $\vec{x}$. For example, we might seek a relative error for the Newton solution $\|\|\vec{F}(\vec{x})\|\|/\|\|\vec{x}\|\|=O(10^{-8})$, but a relative error for the GMRES steps $\|\|A\,\vec{\delta x}-\vec{b}\|\|/\|\|\vec{\delta x}\|\|=O(10^{-3})$ is likely to be sufficient for calculation of the steps $\vec{\delta x}$.

	~~==== Hookstep approach ====~~

	To improve the domain of convergence of the Newton method, it is commonplace to limit the size of the step taken. One approach is simply to take a 'damped' step in the direction of the solution to $(b)$, i.e. step by $\alpha\ \vec{\delta x}_i$, where $\alpha \in (0,1]$. In the '''hookstep approach''', we minimise subject to the condition that the magnitude of the Newton step is limited, $\|\|\vec{\delta x}_i\|\|<\delta$, where $\delta$ is the size of the '''trust region''':
	~~\begin{equation} \label{eq:dxerr}~~
	~~(d)\qquad~~
	~~\min_{\vec{\delta x}_i:\ \|\|\vec{\delta x}_i\|\|<\delta} \ \left\|\left\|\left.~~
	~~\,\vec{F}(\vec{x}_i) + \frac{\vec{\partial F}}{\vec{\partial x}}\right\|_{\vec{x}_i}~~
	~~\vec{\delta x}_i \right\|\right\|\ .~~
	~~\end{equation}~~
	Given the minimisation, the hookstep $\vec{\delta x}_i$ is expected to produce a better result than a simple damped step of the same size. It is also expected to perform much better in 'valleys', where it produces a bent/hooked step to a point along the valley,
	~~rather than jumping from one side of the valley to the other.~~

	The hookstep can be calculated with little extra work to the GMRES method, provided that the size of Krylov-subspace, m, is chosen sufficiently large to solve to the desired accuracy within m GMRES iterations.

	For a given $\vec{\delta x}_i$, the reduction in error predicted by the minimisation $(d)$ can be compared with the actual reduction in $\|\|\vec{F}(\vec{x})\|\|$. According to the accuracy of the prediction, the size of the trust region $\delta$ can be adjusted automatically.

	~~==== Preconditioning ====~~

	The GMRES implementations can be supplied with a preconditioner routine (see [[File:GMRESm.f90]]), but this can be avoided if the method is combined with time integration (see comment at [[File:Arnoldi.f]]).

	~~==== Further details ====~~

	~~Extended overview of method (arxiv;pdf) [~~https://~~arxiv~~.~~org~~/~~pdf~~/~~1908.06730.pdf] (see section 4 on use of code).~~

	~~Usage of NewtonHook subroutine [[File:NewtonHook.f90]] (essentially same for FORTRAN and MATLAB implementations).~~

	~~On GMRES and preconditioning, see [[File:GMRESm.f90]].~~

	~~Further details on the method, see e.g. [http://channelflow.org/dokuwiki/doku.php?id=docs:math:newton_krylov_hookstep Newton-Krylov~~-Hookstep~~] (channelflow.org).~~

Apwillis at 12:03, 3 March 2020

2020-03-03T12:03:36Z

← Older revision		Revision as of 05:03, 3 March 2020
Line 58:		Line 58:
	\begin{equation} \label{eq:dxerr}		\begin{equation} \label{eq:dxerr}
	(d)\qquad		(d)\qquad
	\min_{\vec{\delta x}_i:\ \|\|\vec{\delta x}_i\|\|<\delta} \ \left\|\left\|\left. \frac{\vec{\partial F}}{\vec{\partial x}}\right\|_{\vec{x}_i}		\min_{\vec{\delta x}_i:\ \|\|\vec{\delta x}_i\|\|<\delta} \ \left\|\left\|\left.
	\vec{\delta x}_i ~~+ \vec{F}(\vec{x}_i)~~ \right\|\right\|\ .		\,\vec{F}(\vec{x}_i) + \frac{\vec{\partial F}}{\vec{\partial x}}\right\|_{\vec{x}_i}
			\vec{\delta x}_i \right\|\right\|\ .
	\end{equation}		\end{equation}
	Given the minimisation, the hookstep $\vec{\delta x}_i$ is expected to produce a better result than a simple damped step of the same size. It is also expected to perform much better in 'valleys', where it produces a bent/hooked step to a point along the valley,		Given the minimisation, the hookstep $\vec{\delta x}_i$ is expected to produce a better result than a simple damped step of the same size. It is also expected to perform much better in 'valleys', where it produces a bent/hooked step to a point along the valley,

Apwillis at 11:58, 3 March 2020

2020-03-03T11:58:19Z

← Older revision		Revision as of 04:58, 3 March 2020
Line 66:		Line 66:
	The hookstep can be calculated with little extra work to the GMRES method, provided that the size of Krylov-subspace, m, is chosen sufficiently large to solve to the desired accuracy within m GMRES iterations.		The hookstep can be calculated with little extra work to the GMRES method, provided that the size of Krylov-subspace, m, is chosen sufficiently large to solve to the desired accuracy within m GMRES iterations.

	For a given $\vec{\delta x}_i$, the reduction in error predicted by the ~~linearisation~~ $(d)$ can be compared with the actual reduction in $\|\|\vec{F}(\vec{x})\|\|$. According to the accuracy of the prediction, the size of the trust region $\delta$ can be adjusted automatically.		For a given $\vec{\delta x}_i$, the reduction in error predicted by the minimisation $(d)$ can be compared with the actual reduction in $\|\|\vec{F}(\vec{x})\|\|$. According to the accuracy of the prediction, the size of the trust region $\delta$ can be adjusted automatically.

	==== Preconditioning ====		==== Preconditioning ====

Apwillis at 11:54, 3 March 2020

2020-03-03T11:54:52Z

← Older revision		Revision as of 04:54, 3 March 2020
Line 66:		Line 66:
	The hookstep can be calculated with little extra work to the GMRES method, provided that the size of Krylov-subspace, m, is chosen sufficiently large to solve to the desired accuracy within m GMRES iterations.		The hookstep can be calculated with little extra work to the GMRES method, provided that the size of Krylov-subspace, m, is chosen sufficiently large to solve to the desired accuracy within m GMRES iterations.

	For a given $\vec{\delta x}_i$, the reduction in error predicted by the linearisation $(d)$ can be compared with the actual reduction in $\|\|\vec{F}(\vec{x}~~_i+\vec{\delta x}_i~~)\|\|$. According to the accuracy of the prediction, the size of the trust region $\delta$ can be adjusted automatically.		For a given $\vec{\delta x}_i$, the reduction in error predicted by the linearisation $(d)$ can be compared with the actual reduction in $\|\|\vec{F}(\vec{x})\|\|$. According to the accuracy of the prediction, the size of the trust region $\delta$ can be adjusted automatically.

	==== Preconditioning ====		==== Preconditioning ====

Apwillis at 12:36, 25 February 2020

2020-02-25T12:36:48Z

Apwillis at 00:32, 31 January 2020

2020-01-31T00:32:02Z

← Older revision		Revision as of 17:32, 30 January 2020
Line 6:		Line 6:
	'''CODE (FORTRAN): Template-Example: [[File:Newton_Lorenz.f90]].'''		'''CODE (FORTRAN): Template-Example: [[File:Newton_Lorenz.f90]].'''

	'''Info on main subroutine:''' (essentially same for FORTRAN/MATLAB versions) ~~[[File:NewtonHook.f90]]~~		'''Info on main subroutine:''' [[File:NewtonHook.f90]] (essentially same for FORTRAN/MATLAB versions)

	___ __		___ __

Apwillis at 01:20, 29 January 2020

2020-01-29T01:20:51Z

← Older revision		Revision as of 18:20, 28 January 2020
Line 1:		Line 1:
	{{latexPreamble}}		{{latexPreamble}}
	Free to download:		Free to download:

			'''CODE (MATLAB / GNU Octave): Template-Example: [[File:Newton_Lorenz_m.zip]] / [[File:Newton_Lorenz_m.tgz]]'''

	'''CODE (FORTRAN): Template-Example: [[File:Newton_Lorenz.f90]].'''		'''CODE (FORTRAN): Template-Example: [[File:Newton_Lorenz.f90]].'''

	~~'''CODE (MATLAB / GNU Octave): Template-Example: [[File:Newton_Lorenz_m.zip]] / [[File:Newton_Lorenz_m.tgz]]'''~~

	'''Info on main subroutine:''' (essentially same for FORTRAN/MATLAB versions) [[File:NewtonHook.f90]]		'''Info on main subroutine:''' (essentially same for FORTRAN/MATLAB versions) [[File:NewtonHook.f90]]

← Older revision		Revision as of 05:36, 25 February 2020
Line 37:		Line 37:
	\vec{\delta x}_i = -\vec{F}(\vec{x}_i) \, .		\vec{\delta x}_i = -\vec{F}(\vec{x}_i) \, .
	\end{equation}		\end{equation}
	In order to apply the update ~~\ref{eq:N-R}~~(a), the linear system ~~\ref{eq:N-R}~~(b) needs to be solved for the unknown $\vec{\delta x}_i$.		In order to apply the update $(a)$, the linear system $(b)$ needs to be solved for the unknown $\vec{\delta x}_i$.

	==== Jacobian-Free Newton method ====		==== Jacobian-Free Newton method ====

	The $n\times n$ Jacobian matrix $\left.\frac{\vec{\partial F}}{\vec{\partial x}}\right\|_{\vec{x}_i}$ is usually difficult to evaluate. However, the problem ~~\ref{eq:N-R}~~(b) is in the form $A\,\vec{\delta x}=\vec{b}$, which can be solved using the '''Krylov-subspace method GMRES(m)'''. (See [[File:GMRESm.f90]].) The GMRES algorithm does not need to know the matrix $A$ itself, only the result of multiplying a vector by $A$. For a given starting vector $\vec{\delta x}_0$, e.g. $\vec{\delta x}_0=\vec{b}$, the method seeks solutions for $\vec{\delta x}$ in $\mathrm{span}\{\vec{x}_0,\,A\vec{x}_0,\,A^2\vec{x}_0,...\}$, but uses Gram-Schmidt orthogonalisation to improve the numerical suitability of this basis set. The set of orthogonalised vectors is called the Krylov-subspace, and m is the maximum number of vectors stored.		The $n\times n$ Jacobian matrix $\left.\frac{\vec{\partial F}}{\vec{\partial x}}\right\|_{\vec{x}_i}$ is usually difficult to evaluate. However, the problem $(b)$ is in the form $A\,\vec{\delta x}=\vec{b}$, which can be solved using the '''Krylov-subspace method GMRES(m)'''. (See [[File:GMRESm.f90]].) The GMRES algorithm does not need to know the matrix $A$ itself, only the result of multiplying a vector by $A$. For a given starting vector $\vec{\delta x}_0$, e.g. $\vec{\delta x}_0=\vec{b}$, the method seeks solutions for $\vec{\delta x}$ in $\mathrm{span}\{\vec{x}_0,\,A\vec{x}_0,\,A^2\vec{x}_0,...\}$, but uses Gram-Schmidt orthogonalisation to improve the numerical suitability of this basis set. The set of orthogonalised vectors is called the Krylov-subspace, and m is the maximum number of vectors stored.

	Iterations of the GMRES algorithm for the problem ~~\ref{eq:N-R}~~(b) involve calculating products of the Jacobian with given $\vec{\delta x}$, which may be approximated, e.g.		Iterations of the GMRES algorithm for the problem $(b)$ involve calculating products of the Jacobian with given $\vec{\delta x}$, which may be approximated, e.g.
	\begin{equation}\label{eq:Japprox}		\begin{equation}\label{eq:Japprox}
			(c)\qquad
	\left.\frac{\vec{\partial F}}{\vec{\partial x}}\right\|_{\vec{x}_i}\,\vec{\delta x} \approx		\left.\frac{\vec{\partial F}}{\vec{\partial x}}\right\|_{\vec{x}_i}\,\vec{\delta x} \approx
	\frac{1}{\epsilon}(\vec{F}(\vec{x}_i+\epsilon\,\vec{\delta x})-\vec{F}(\vec{x}_i))		\frac{1}{\epsilon}(\vec{F}(\vec{x}_i+\epsilon\,\vec{\delta x})-\vec{F}(\vec{x}_i))
Line 54:		Line 55:
	==== Hookstep approach ====		==== Hookstep approach ====

	To improve the domain of convergence of the Newton method, it is commonplace to limit the size of the step taken. One approach is simply to take a 'damped' step in the direction of the solution to ~~\ref{eq:N-R}~~(b), i.e. step by $\alpha\ \vec{\delta x}_i$, where $\alpha \in (0,1]$. In the '''hookstep approach''', we minimise subject to the condition that the magnitude of the Newton step is limited, $\|\|\vec{\delta x}_i\|\|<\delta$, where $\delta$ is the size of the '''trust region''':		To improve the domain of convergence of the Newton method, it is commonplace to limit the size of the step taken. One approach is simply to take a 'damped' step in the direction of the solution to $(b)$, i.e. step by $\alpha\ \vec{\delta x}_i$, where $\alpha \in (0,1]$. In the '''hookstep approach''', we minimise subject to the condition that the magnitude of the Newton step is limited, $\|\|\vec{\delta x}_i\|\|<\delta$, where $\delta$ is the size of the '''trust region''':
	\begin{equation} \label{eq:dxerr}		\begin{equation} \label{eq:dxerr}
			(d)\qquad
	\min_{\vec{\delta x}_i:\ \|\|\vec{\delta x}_i\|\|<\delta} \ \left\|\left\|\left. \frac{\vec{\partial F}}{\vec{\partial x}}\right\|_{\vec{x}_i}		\min_{\vec{\delta x}_i:\ \|\|\vec{\delta x}_i\|\|<\delta} \ \left\|\left\|\left. \frac{\vec{\partial F}}{\vec{\partial x}}\right\|_{\vec{x}_i}
	\vec{\delta x}_i + \vec{F}(\vec{x}_i) \right\|\right\|\ .		\vec{\delta x}_i + \vec{F}(\vec{x}_i) \right\|\right\|\ .
Line 64:		Line 66:
	The hookstep can be calculated with little extra work to the GMRES method, provided that the size of Krylov-subspace, m, is chosen sufficiently large to solve to the desired accuracy within m GMRES iterations.		The hookstep can be calculated with little extra work to the GMRES method, provided that the size of Krylov-subspace, m, is chosen sufficiently large to solve to the desired accuracy within m GMRES iterations.

	For a given $\vec{\delta x}_i$, the reduction in error predicted by the linearisation ~~\ref{eq:dxerr}~~ can be compared with the actual reduction in $\|\|\vec{F}(\vec{x}_i+\vec{\delta x}_i)\|\|$. According to the accuracy of the prediction, the size of the trust region $\delta$ can be adjusted automatically.		For a given $\vec{\delta x}_i$, the reduction in error predicted by the linearisation $(d)$ can be compared with the actual reduction in $\|\|\vec{F}(\vec{x}_i+\vec{\delta x}_i)\|\|$. According to the accuracy of the prediction, the size of the trust region $\delta$ can be adjusted automatically.

	==== Preconditioning ====		==== Preconditioning ====