.. _alg_gmres:

************************************
Generalized minimal residual (GMRES)
************************************

The Generalized Minimal Residual (GMRES) method is a Krylov subspace method for
solving general non-Hermitian linear systems :math:`A \mathbf{x} = \mathbf{b}`.
It builds an approximate solution within a Krylov subspace by minimizing the
residual norm over that subspace.

Mathematical formulation
========================

GMRES seeks to solve:

.. math::
    A \mathbf{x} = \mathbf{b}

by constructing an approximate solution :math:`\mathbf{x}_m` within the Krylov
subspace of dimension :math:`m` generated by the initial residual
:math:`\mathbf{r}_0 = \mathbf{b} - A\mathbf{x}_0`:

.. math::
    \mathcal{K}_m(A, \mathbf{r}_0) = \mathrm{span}\{\mathbf{r}_0, A\mathbf{r}_0, \ldots, A^{m-1}\mathbf{r}_0\}

The method finds the vector :math:`\mathbf{x}_m` that minimizes the residual norm
:math:`\|\mathbf{b} - A\mathbf{x}_m\|` over this subspace.

Arnoldi iteration
-----------------

An orthonormal basis :math:`\{\mathbf{v}_1, \ldots, \mathbf{v}_m\}` for
:math:`\mathcal{K}_m` is built via the Arnoldi iteration with modified Gram-Schmidt
orthogonalization:

1. Start with :math:`\mathbf{v}_1 = \mathbf{r}_0 / \|\mathbf{r}_0\|`
2. For :math:`j = 1, \ldots, m`:

   a. Compute :math:`\mathbf{w} = A \mathbf{v}_j`
   b. Orthogonalize against previous vectors: :math:`h_{ij} = \langle \mathbf{v}_i, \mathbf{w} \rangle`
   c. :math:`\mathbf{w} \leftarrow \mathbf{w} - \sum_{i=1}^{j} h_{ij} \mathbf{v}_i`
   d. :math:`h_{j+1,j} = \|\mathbf{w}\|`
   e. :math:`\mathbf{v}_{j+1} = \mathbf{w} / h_{j+1,j}`

This produces an upper Hessenberg matrix :math:`H_m` of size :math:`(m+1) \times m`.

Least squares solution
----------------------

Once the basis is built, the residual minimization reduces to a small
:math:`(m+1) \times m` least-squares problem:

.. math::
    \mathbf{y} = \mathrm{argmin}_\mathbf{y} \|\|\mathbf{r}_0\| \mathbf{e}_1 - H_m \mathbf{y}\|

which is solved using standard linear algebra routines. The approximate solution
is then constructed as:

.. math::
    \mathbf{x}_m = \mathbf{x}_0 + \sum_{j=1}^{m} y_j \mathbf{v}_j

MPS implementation
==================

In the MPS formulation, each Arnoldi iteration involves:

- One MPO-MPS product (:math:`A \mathbf{v}_j`)
- Multiple MPS scalar products for orthogonalization
- Linear combinations of MPS
- Simplification steps to control bond dimension growth

Since each iteration increases the bond dimension through MPO-MPS contractions and
linear combinations, simplification steps are essential to keep the MPS representation
tractable.

Restarted GMRES
---------------

The implementation supports restarted GMRES, where upon reaching the maximum subspace
dimension ``nvectors``, the current approximation :math:`\mathbf{x}_m` becomes the
initial guess for a fresh Krylov subspace. This bounds memory usage while still
allowing convergence for difficult problems.

The algorithm detects ill-conditioning in the Krylov basis (when :math:`h_{j+1,j}`
becomes very small) and can restart early if needed.

When to use GMRES
=================

GMRES is best suited for:

- **General non-Hermitian systems**: No symmetry requirements on :math:`A`
- **Problems requiring high accuracy**: GMRES minimizes the residual optimally over the Krylov subspace
- **Moderately ill-conditioned systems**: The Krylov subspace naturally adapts to the problem structure

Compared to BiCGSTAB:

- GMRES guarantees monotonic residual decrease (within each restart cycle)
- GMRES requires storing :math:`m` basis vectors, while BiCGSTAB uses constant storage
- For very large problems, restarted GMRES may converge slower than BiCGSTAB

Example
=======

.. code-block:: python

    import numpy as np
    from seemps.state import random_uniform_mps
    from seemps.operators import MPO
    from seemps.solve import gmres_solve

    # Create your operator and right-hand side
    n = 10
    # ... construct MPO A and MPS b ...

    # Solve A @ x = b with restarted GMRES
    x, residual = gmres_solve(
        A, b,
        nvectors=10,      # Krylov subspace size
        max_restarts=5,   # Maximum number of restarts
        tolerance=1e-8
    )
    print(f"Final residual: {residual}")

.. autosummary::

    ~seemps.solve.gmres_solve

See also
========

- :doc:`cgs` - Conjugate gradient method for Hermitian positive-definite systems
- :doc:`bicgs` - Biconjugate gradient stabilized method
- :doc:`dmrg_solve` - DMRG-based solver with local tensor optimization
- :doc:`arnoldi` - Arnoldi method for eigenvalue problems (related Krylov technique)