Preconditioners#

Some optimization methods, like the conjugate gradient, benefit from using a preconditioner matrix.

Some optimization methods, like the conjugate gradient, can be accelerated by transforming the linear system they need to iteratively solve in such a way that improves the eigenvalue distribution of the matrix \(\mathbf{A}\). One way to achieve this is to perform a change of variables through a preconditioning matrix (see [NocedalWright1999], page 118).

We can pass a preconditioner to most of the minimizers of the new framework, and it offers some functions to generate appropriate preconditioners for a given objective function.

Preconditioning the conjugate gradient#

The minimize.conjugate_gradient function through the preconditioner argument.

For example, consider an objective function phi. We can use the get_jacobi_preconditioner() to obtain a suitable preconditioner matrix for it:

import inversion_ideas as ii

initial_model = np.zeros(n_params)
preconditioner = ii.get_jacobi_preconditioner(phi, initial_model)
preconditioner

<DIAgonal sparse array of dtype 'float64'
	with 5 stored elements (1 diagonals) and shape (5, 5)>

preconditioner.toarray()

array([[1.30998202e-05, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
        0.00000000e+00],
       [0.00000000e+00, 1.40400846e-05, 0.00000000e+00, 0.00000000e+00,
        0.00000000e+00],
       [0.00000000e+00, 0.00000000e+00, 1.44113431e-05, 0.00000000e+00,
        0.00000000e+00],
       [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.36213274e-05,
        0.00000000e+00],
       [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
        1.54970445e-05]])

Now let’s minimize it using conjugate_gradient() and use that matrix as a preconditioner:

inverted_model = ii.minimize.conjugate_gradient(
   phi, initial_model, preconditioner=preconditioner
)

Preconditioning the Gauss-Newton conjugate gradient#

When minimizing non-linear functions we’ll some non-linear minimizer like minimize.GaussNewtonConjugateGradient. When calling this object we can also pass a preconditioner object:

# Define the minimizer
minimizer = ii.minimize.GaussNewtonConjugateGradient()

# Run it using the `run()` method, and pass the preconditioner
inverted_model = minimizer.run(
   phi, initial_model, preconditioner=preconditioner
)

The preconditioner will be passed by the minimize.GaussNewtonConjugateGradient to the conjugate gradient algorithm that runs once per Gauss-Newton iteration.

The issue of this approach is that the linear system that gets passed to the conjugate gradient function is different for each iteration: the Hessian matrix depends on the current model in the Gauss-Newton step, and because the problem is non-linear, that Hessian will be different for each iteration. Trying to use the same preconditioner for each one of those systems might be unproductive. Instead, it’s better to calculate a new preconditioner for each Gauss-Newton step.

We can achieve that by using a dynamic preconditioner object, like the inversion_ideas.JacobiPreconditioner, and pass it to the minimizer as the preconditioner argument. The minimize.GaussNewtonConjugateGradient will run the update method of the passed precondtioner object before running the conjugate gradient algorithm, and it’ll pass it as argument to it.

For example:

# Define a dynamic Jacobi preconditioner
preconditioner = ii.JacobiPreconditioner(phi)
preconditioner

<5x5 JacobiPreconditioner with dtype=float64>

Note

These dynamic preconditioner are LinearOperator objects that under the hood they hold the preconditioning matrix after it gets updated or recomputed.

We can pass it to the minimizer when running it:

# Run it using the `run()` method, and pass the preconditioner
inverted_model = minimizer.run(
   phi, initial_model, preconditioner=preconditioner
)

References#

[NocedalWright1999]

Nocedal, J., & Wright, S. J. (1999). Numerical optimization. Springer.