Add weighted GMRES(m)

[pescap]

I propose to implement the weighted GMRES(m), as this variant can lead to mesh independent convergence results for FEM/BEM simulations (see e.g. this article).

The weighted version requires only a few changes to the code as compared to its classic Euclidean version.
It would take H (a supposedly Hermitian positive definite) LinearOperator as an input.

def gmres(A, b, x0=None, tol=1e-5, restrt=None, maxiter=None, 
           M=None, H=None, callback=None, residuals=None, orthog='householder', 
           **kwargs): 

I think it could be interesting to add it to pyamg, as I expect the H^1(\Omega) Sobolev norm GMRES to outperform the Euclidean one (I am not so familiar with pyamg).

Let me know if you have further comments. I have already some work in progress and will try to send a PR soon.

Also refer to: scipy/scipy#15693 (comment) for its scipy variant (also work in progress).

[lukeolson]

Interesting. So is the approach to simply replace the Euclidean inner products with $H$ inner products?

[pescap]

I sent a PR for the H-norm _gmres_mgs. I attach below a simple code that tests it. I could add a similar proceeding to test_krylov.py later on.

So far, I could not implement the H-norm _gmres_householder as I am not totally understanding which parts have to be replaced in the reflections procedure. Also, it sounds like I have to access to the amg_core and modify:

pyamg/pyamg/amg_core/krylov.h

Line 50 in af81d67

T alpha = dot_prod(Bptr, z, n);

I think it would be judicious to create another issue (and related PR) later on for the householder case.

[pescap]

import numpy as np
import pyamg
import scipy.linalg

# Solving Ax=b with H-norm GMRES(k) is equivalent to solving
# (H^{1/2}AH^{-1/2})(H^{1/2}x)=H^{1/2}b with Euclidean norm GMRES(k)

N = 1000
A = np.random.random((N, N)) + 1j * np.random.random((N, N))
b = np.random.random(N)
H = np.random.random(N) * np.eye(N)
Hm = np.linalg.inv(H)

H12 = scipy.linalg.sqrtm(H)
Hm12 = scipy.linalg.sqrtm(Hm)

Adot = np.dot(np.dot(H12, A), Hm12)
bdot = np.dot(H12, b)

x0=None
tol=1e-5,
restrt=None
maxiter=None,
M=None
callback=None
residuals=None
reorth=False

xmgs, conv = pyamg.krylov.gmres(Adot, bdot, orthog = 'mgs')
xHmgs, convH = pyamg.krylov.gmres(A, b, H=H, orthog = 'mgs')

xhous, conv = pyamg.krylov.gmres(Adot, bdot)
#xHhous, convH = pyamg.krylov.gmres(A,b,H=I)

xhousdot = np.dot(Hm12, xhous)
xmgsdot = np.dot(Hm12, xmgs)

print(np.linalg.norm(np.dot(A, xmgsdot)-b))
print(np.linalg.norm(np.dot(A, xHmgs) - b))
print(np.linalg.norm(np.dot(A, xhousdot)-b))

[pescap]

According to Householder orthogonalization with a non-standard inner product (Meiyue Shao), it sounds like Householder reflections for H-inner products are of the form (see section 2):

I -  2*w*conjugate(w)^T H

Householder-based implementation sounds more complicated than what I was expecting. I think we are fine with the MGS-based H-GMRES. In the appropriate setting, H-GMRES(m) converges for all m\geq1, so I hope that mgs is sufficient.

[pescap]

Solving Ax=b with H-norm GMRES(k) is equivalent to solving
(H^{1/2}AH^{-1/2})(H^{1/2}x)=H^{1/2}b with Euclidean norm GMRES(k)

Refer to Some observations on weighted GMRES (Güttel, Stefan and Pestana, Jen) Algorithms 1&2 p. 13.