PJFNK Matlab

This repository is for an implementation of the Preconditioned Jacobian-Free Newton-Krylov method in Matlab.

PJFNK MATLAB

This repository is for an implementation of the Preconditioned Jacobian-Free Newton-Krylov method in Matlab, as well as methods for standard full Newton method and formulation of full and sparse Numerical Jacobian matrices. With this method, users can choose a Krylov method from any of the following:

• 'gmres' - Generalized Minimum RESidual

  • (default) Most flexible method, but may require restarts to be practical

• 'bicgstab' - BiConjugate Gradients Stabilized

  • Enhanced BiCG method with a GMRES(1) step for stabilization

• 'bicgstabl'

  • Enhanced BiCGStab method with a GMRES(2) step for stabilization

• 'pcg'

  • Optimal method for Symmetric Positive Definite Jacobians

• 'minres'

  • Jacobians must be Symmetric, but does not need to be Positive Definite

• 'symmlq'

  • Jacobians must be Symmetric, but does not need to be Positive Definite

• 'cgs'

  • Similar to BiCG, but does not require transposition of matrix

• 'tfqmr'

  • Optimal for Symmetric Jacobians that are indefinite

More information on Krylov methods here

Requirements

• MATLAB
• [Optional] Optimization Toolbox

NOTE: The Optimization Toolbox in MATLAB is not required to run these codes, but is used for validation checks in some of the unit tests.

Basic Usage

To use this function, you must provide (at a minimum) a:

• Function Handle to evaluate a set of equations (F = @(x) ...)

• An initial guess vector (x0)

The size of the vector (x0) and the size of the evaluated function handle (F(x)) must both be Nx1.

% Example 1 - Basics

F = @(x) [2*x(1)*x(2) - x(2)^2; 
            -x(1)+2*x(2)-exp(-x(3)); 
            -x(2)+2*x(3)^3];

x0 = [1; 1; 1];

[x,stats,opts] = JacobianFreeNewtonKrylov(F,x0);

The function will return the solution to the system (x) and 2 struct() objects containing: (i) the status information of the solve (stats) and (ii) a set of options used in the solve (opts). Several default options are provided.

You can also provide a struct() to the input arguments of the function to change the behavior of the solve as needed.

% Example 2 - Solver Options
options = struct();
options.maxiter = 20;  % Set the maximum number of non-linear steps
options.ftol = 1e-6;   % Set the solver tolerance
options.krylov_solver = 'gmres';        % Pick from a selection of Krylov methods
options.linesearch = 'backtracking';    % Pick from a selection of linesearch methods

% [Optional] Provide a function handle for a known Jacobian
options.jacfun = J = @(x)    [2*x(2), 2*x(1)-2*x(2), 0;
                                -1, 2, exp(-x(3));
                                 0, -1, 6*x(3)^2];

This is a Jacobian-Free method, so you technically never need to supply the Jacobian. However, you can actually choose to use the full (or sparse) Jacobian (if you desire).

options.use_matrix = true;   % Specifies to use the Jacobian for the solve

By default, the method will set use_matrix to false, and just use the JacobianOperator function in conjunction with a Krylov solver. However, if you tell the solver to use the matrix form, it will do that instead.

For the Krylov solvers, you can also control how those methods are solved by adding to the options data structure.

% Example 3 - Setting Krylov options for 'gmres'
options.krylov_opts = struct();
options.krylov_opts.restart = 50;  % Iterations before restart
options.krylov_opts.tol = 1e-6;    % Linear solver tolerance

NOTE: Each Krylov solver will have its own set of options. The options available (including preconditioners) follows exactly with the Matlab standard set of options for each iterative solver (see here).

Lastly, there is also a full NewtonMethod implementation. You can use this method if you do not want to use Krylov solvers underneath the non-linear problem. In this case, the search direction is resolved using the Matlab mldivide () function. This method also has almost all the same options available to it as the JacobianFreeNewtonKrylov method.

% Example 4 - Using full Newton Method
F = @(x) [2*x(1) - x(2); 
            -x(1)+2*x(2)-x(3); 
            -x(2)+2*x(3) + 1];

x0 = [1; 0; 0];

% Exact Jacobian
J = @(x) [2, -1, 0;
            -1, 2, -1;
            0, -1, 2];

% Store the Jacobian function handle in the data structure 
options = struct();
options.jacfunc = J;

[x,stats,opts] = NewtonMethod(F,x0,options);

For more usage examples, see the tests folder of this repository.

---

Functions

• JacobianOperator

Performs an action of a Jacobian on any vector v without forming the Jacobian explicitly.

% Pseudo Code
Jv = (F(x+epsilon*v) - F(x))/epsilon;

where F(x) evaluates the non-linear system at state x.

where v is the vector being multiplied by the Jacobian.

where epsilon is a small perturbation value.

---

• NumericalJacobianMatrix

Approximates the full Jacobian with finite-differences approach. Can return either a 'dense' or 'sparse' Jacobian via user request. Default is to return a 'dense' Jacobian.

% Pseudo Code
for i=1:N
    dx = x;
    dx(i,1) =  x(i,1) + epsilon;
    J(:, i) = (F(dx) - F(x))/epsilon;
end

where F(x) evaluates the non-linear system at state x.

where epsilon is a small perturbation value.

---

• NewtonMethod

Performs a basic Newton's method with direct inversion of the Jacobian matrix. The user may provide an exact Jacobian function as a function handle. If they do not, then the method will use the NumericalJacobianMatrix function to approximate the Jacobian.

% Pseudo Code
while (norm(F(x)) > ftol)
    % Perform a Newton Step
    %   J*s = -F --> s = -J\F
    x = x0 - J(x)\F(x);

    % Update x
    x0=x;
end

where F(x) evaluates the non-linear system at state x.

where J(x) evaluates the Jacobian matrix at state x.

where ftol is the error or residual tolerance for convergence.

---

• JacobianFreeNewtonKrylov

Performs a Jacobian-Free Newton-Krylov solve with associated user options. The user may choose any valid Krylov method to solve the Newton step.

% Pseudo Code
while (norm(F(x)) > ftol)
    % Perform a Newton Step
    %   J*s = -F --> s = -J\F
    x = x0 + KrylovSolve( @JacOp( @F(x0),x0,F), -F );

    % Update x
    x0=x;
end

where F(x) evaluates the non-linear system at state x.

where JacOp evaluates the Jacobian operator at state x.

where KrylovSolve evaluates the solution to the linear system.

where ftol is the error or residual tolerance for convergence.

---

Linesearch Methods

• StandardNewtonStep

Takes a standard Newton step without any linesearching.

% Invoked via the 'options' struct() for Newton methods
options.linesearch = 'none';

---

• BacktrackingLinesearch

Attempts to take a standard Newton step, then if the step fails to show sufficient reduction in the norm, it reduces the step size by the 'contraction_rate' factor. Reduction of step size stops when it reaches the 'min_step' tolerance, regardless of whether or not sufficient norm reduction was achieved.

Guarantees that the residual reduction is monotonic, but may result in step sizes too small to get a solution in reasonable timeframe.

% Invoked via the 'options' struct() for Newton methods
options.linesearch = 'backtracking';
options.linesearch_opts = struct('min_step',1e-3,'contraction_rate',0.5);

---

Preconditioning

Users may provide any custom preconditioning matrices or function handles to the JacobianFreeNewtonKrylov method. Those preconditioning functions must have the following format

% Custom preconditioner function used to help establish the search direction
%
%       Consider this function to be a Pseudo-Inverse for the following:
%           s = J(x)\b
%               where s = a vector to solve for, J(x) is the current 
%               Jacobian at the current non-linear state, and b is the 
%               known right-hand side vector
%
%   @param b = Nx1 vector representing the right-hand size of J*s = b
%   @param Jacfun = Function handle for the Jacobian (J = Jacfun(x))
%   @param x = Nx1 vector for the current non-linear state
%   @param options = An optional user defined struct() for any other information 
%               you may need to perform this action.
function s = custom_precon(b,Jacfun,x,options)
    ...
end

You provide this preconditioner to the JacobianFreeNewtonKrylov function through the options struct().

options.krylov_opts = struct();
options.krylov_opts.M1 = @(b,Jacfun,x,options) custom_precon(b,Jacfun,x,options);

NOTE: We use the same preconditioning syntax from the standard Matlab library of Iterative Krylov Methods (see here). This allows users to provide either a single preconditioner Matrix/Function [M1], or a pair of preconditioners [M1 and M2]. If you provide both preconditioners, note that they are resolve in the following order: M2(M1\b).

---

Automatic Matrix Reordering for JFNK

When using Krylov methods to iteratively solve the linear sub-problems, often times the condition number of the matrix will make solutions difficult to accomplish (even when preconditioning!). Thus, the JacobianFreeNewtonKrylov implementation here also allows the user to specify a reordering scheme.

% Example - Specifying a reordering scheme
options.use_matrix = true;    % Tells the method to use either the user 
                              % provided Jacobian function or a numerical 
                              % Jacobian matrix to apply ordering to

options.equilibrate = true;   % Tells the method that you will request an
                              % automatic reordering of the Jacobian

options.reordering_method = 'dissect';  % Tells the method to use the nested
                                        % dissection method for reordering

This option is only available to you if you specify to use_matrix for the solve. Otherwise, if you do not set use_matrix to true, the equilibrate option will just be ignored.

There are 3 supported reordering methods:

• 'dissect' = Nested dissection permutation

• 'amd' = Approximate minimum degree permutation

• 'symrcm' = Sparse reverse Cuthill-McKee ordering

More info on the equilibrate method

IMPORTANT NOTE: You MUST supply all preconditioners to the JacobianFreeNewtonKrylov function as function handles, AND you must use the newly calculated Jacobian from Jacfun(x) inside your preconditioning function to formulate the preconditioner result. This is because when doing reordering, the shape/form of the original problem has changed, and may have changed in an unknown way. Thus, using this option means that you may not know what your Jacobian looks like when it gets to your preconditioning function.

---

References

** JacobianOperator and JacobianFreeNewtonKrylov**

• Knoll D.A., Keyes D.E. (2003) "Jacobian-Free Newton-Krylov methods: a survey of approaches and applications", Journal of Computational Physics.

** BacktrackLinesearch **

• Armijo, Larry (1966). "Minimization of functions having Lipschitz continuous first partial derivatives". Pacific J. Math. 16 (1): 1–3.

** YBusFormation for Power Flow **

• Praviraj PG (2023). Admittance Bus (Y-Bus) Formation (https://www.mathworks.com/matlabcentral/fileexchange/19512-admittance-bus-y-bus-formation), MATLAB Central File Exchange. Retrieved March 19, 2023.

---

Citation

Ladshaw, A.P., "PJFNK MATLAB: A MATLAB implementation for Jacobian-Free Newton-Krlov solvers for non-linear systems," https://github.com/aladshaw3/pjfnk_matlab, Accessed (Month) (Day), (Year).