mrpro.algorithms.optimizers.cg

mrpro.algorithms.optimizers.cg(operator: LinearOperator, right_hand_side: Tensor, initial_value: Tensor | None = None, max_iterations: int = 128, tolerance: float = 1e-4, callback: Callable[[CGStatus], bool | None] | None = None) → Tensor[source]

CG for solving a linear system \(Hx=b\).

This algorithm solves systems of the form \(H x = b\), where \(H\) is a self-adjoint linear operator and \(b\) is the right-hand side. The method can solve a batch of \(N\) systems jointly, thereby taking \(H\) as a block-diagonal with blocks \(H_i\) and \(b = [b_1, ..., b_N] ^T\).

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\).

This implementation assumes that \(H\) is self-adjoint and does not verify this condition.

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

Parameters:

  • operator (LinearOperator) – self-adjoint operator (named \(H\) above)

  • right_hand_side (Tensor) – right-hand-side of the system (named \(b\) above)

  • initial_value (Tensor | None, default: None) – initial value of the algorithm; if None, it will be set to right_hand_side

  • max_iterations (int, default: 128) – maximal number of iterations. Can be used for early stopping.

  • tolerance (float, default: 1e-4) – tolerance for the residual; if set to zero, the maximal number of iterations is the only stopping criterion used to stop the cg. If the condition number of \(H\) is large, a small residual may not imply a highly accurate solution.

  • callback (Callable[[CGStatus], bool | None] | None, default: None) – Function to be called at each iteration. This can be used to monitor the progress of the algorithm. If it returns False, the algorithm stops at that iteration.

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.