mrpro.algorithms.optimizers.cg

mrpro.algorithms.optimizers.cg(operator: Callable[[Tensor], tuple[Tensor]], right_hand_side: tuple[Tensor] | Tensor, *, initial_value: tuple[Tensor] | Tensor | None = None, preconditioner_inverse: Callable[[Tensor], tuple[Tensor]] | None = None, max_iterations: int = 128, tolerance: float = 0.0001, callback: Callable[[CGStatus], bool | None] | None = None) → tuple[Tensor][source]
mrpro.algorithms.optimizers.cg(operator: Callable[[Unpack[tuple[Tensor, ...]]], tuple[Tensor, ...]], right_hand_side: Sequence[Tensor], *, initial_value: Sequence[Tensor] | None = None, preconditioner_inverse: Callable[[Unpack[tuple[Tensor, ...]]], tuple[Tensor, ...]] | None = None, max_iterations: int = 128, tolerance: float = 0.0001, callback: Callable[[CGStatus], bool | None] | None = None) → tuple[Tensor, ...]

(Preconditioned) Conjugate Gradient for solving \(Hx=b\).

This algorithm solves systems of the form \(H x = b\), where \(H\) is a self-adjoint positive semidefinite linear operator and \(b\) is the right-hand side.

The method performs the following steps:

  1. Initialize the residual \(r_0 = b - Hx_0\) (with \(x_0\) as the initial guess).

  2. Set the search direction \(p_0 = r_0\).

  3. For each iteration \(k = 0, 1, 2, ...\):

    • Compute \(\alpha_k = \frac{r_k^T r_k}{p_k^T H p_k}\).

    • Update the solution: \(x_{k+1} = x_k + \alpha_k p_k\).

    • Compute the new residual: \(r_{k+1} = r_k - \alpha_k H p_k\).

    • Update the search direction: \(\beta_k = \frac{r_{k+1}^T r_{k+1}}{r_k^T r_k}\),

      then \(p_{k+1} = r_{k+1} + \beta_k p_k\).

The operator can be either a LinearOperator or a LinearOperatorMatrix. In both cases, this implementation does not verify if the operator is self-adjoint and positive semidefinite. It will silently return the wrong result if the assumptions are not met. It can solve a batch of \(N\) systems jointly if right_hand_side has a batch dimension and the operator interprets the batch dimension as \(H\) being block-diagonal with blocks \(H_i\) resulting in \(b = [b_1, ..., b_N] ^T\).

If preconditioner_inverse is provided, it solves \(M^{-1}Hx = M^{-1}b\) implicitly, where preconditioner_inverse(r) computes \(M^{-1}r\).

See [Hestenes1952], [Nocedal2006], and [WikipediaCG] for more information.

Parameters:
Returns:

An approximate solution of the linear system \(Hx=b\).

References

[Hestenes1952]

Hestenes, M. R., & Stiefel, E. (1952). Methods of conjugate gradients for solving linear systems. Journal of Research of the National Bureau of Standards , 49(6), 409-436

[Nocedal2006]

Nocedal, J. (2006). Numerical Optimization (2nd ed.). Springer.