mr2.algorithms.optimizers.bicg

mr2.algorithms.optimizers.bicg(operator: LinearOperator | LinearOperatorMatrix | Callable[[Unpack[tuple[Tensor, ...]]], tuple[Tensor, ...]], right_hand_side: Sequence[Tensor], *, initial_value: Sequence[Tensor] | None = None, preconditioner_inverse: LinearOperator | LinearOperatorMatrix | Callable[[Unpack[tuple[Tensor, ...]]], tuple[Tensor, ...]] | None = None, max_iterations: int = 128, tolerance: float = 1e-6, callback: Callable[[BiCGStatus], bool | None] | None = None) → tuple[Tensor, ...][source]

(Preconditioned) Bi-Conjugate Gradient Stabilized method for \(Hx=b\).

This algorithm solves systems of the form \(H x = b\), where \(H\) is a linear operator (potentially non-symmetric) and \(b\) is the right-hand side. The method can solve multiple systems jointly if the operator and tensors handle batch dimensions.

The method performs the following steps (simplified):

1. Initialize \(x_0\), compute residual \(r_0 = b - Hx_0\). Set shadow residual \(\hat{r}_0 = r_0\).

2. Initialize \(\rho_0 = \alpha = \omega = 1\), \(p_0 = v_0 = 0\).

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

   1. Compute inner product: \(\rho_k = \hat{r}_0^H r_{k-1}\).

   2. Compute step direction modifier: \(\beta_{k-1} = (\rho_k / \rho_{k-1}) (\alpha_{k-1} / \omega_{k-1})\).

   3. Update search direction: \(p_k = r_{k-1} + \beta_{k-1} (p_{k-1} - \omega_{k-1} v_{k-1})\).

   4. Apply preconditioner: \(y_k = M^{-1} p_k\).

   5. Compute operator applied to search direction: \(v_k = H y_k\).

   6. Compute step size \(\alpha_k = \rho_k / (\hat{r}_0^H v_k)\).

   7. Compute intermediate residual: \(s_k = r_{k-1} - \alpha_k v_k\).

   8. Check convergence: If \(\|s_k\|_2\) is small, update \(x_k = x_{k-1} + \alpha_k y_k\) and stop.

   9. Apply preconditioner to intermediate residual: \(z_k = M^{-1} s_k\).

   10. Compute operator applied to preconditioned intermediate residual: \(t_k = H z_k\).

   11. Compute stabilization factor: \(\omega_k = (t_k^H s_k) / (t_k^H t_k)\).

   12. Update solution: \(x_k = x_{k-1} + \alpha_k y_k + \omega_k z_k\).

   13. Update residual: \(r_k = s_k - \omega_k t_k\).

   14. Update \(\rho_{k-1} = \rho_k\) for the next iteration.

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

See [VanDerVorst1992] and [WikipediaBiCGSTAB] for more details.

Parameters:

• operator (LinearOperator | LinearOperatorMatrix | Callable[[Unpack[tuple[Tensor, ...]]], tuple[Tensor, ...]]) – Linear operator \(H\). Does not need to be self-adjoint.

• right_hand_side (Sequence[Tensor]) – Right-hand-side \(b\). Must be a sequence of tensors.

• initial_value (Sequence[Tensor] | None, default: None) – Initial guess \(x_0\). If None, defaults to zeros. Must be a sequence of tensors.

• preconditioner_inverse (LinearOperator | LinearOperatorMatrix | Callable[[Unpack[tuple[Tensor, ...]]], tuple[Tensor, ...]] | None, default: None) – Preconditioner \(M^{-1}\). If None, no preconditioning is applied (identity operator).

• max_iterations (int, default: 128) – Maximum number of iterations.

• tolerance (float, default: 1e-6) – Tolerance for the L2 norm of the residual \(\|r_k\|_2\).

• callback (Callable[[BiCGStatus], bool | None] | None, default: None) – Function called at the end of each iteration. See BiCGSTABStatus. If the callback returns False, the iteration stops.

Returns:

An approximate solution \(x\) of the linear system \(Hx=b\) as a tuple of tensors.

References

[VanDerVorst1992]

Van der Vorst, H. A. (1992). Bi-CGSTAB: A Fast and Smoothly Converging Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems. SIAM Journal on Scientific and Statistical Computing, 13(2), 631-644.

[WikipediaBiCGSTAB]

Wikipedia: Bi-conjugate gradient stabilized method https://en.wikipedia.org/wiki/Bi-conjugate_gradient_stabilized_method