Is it possible to define a custom inner product for gmres in Krylov.jl ?
1 Like
Yes you can implement your own kdot
1 Like
The example implementations of kdot and knorm seem wrong? The element type is FloatOrComplex, but there are no complex conjugations?
That is, kdot should use res += conj(_x[i,j]) * _y[i,j] and knorm should use res += abs2(_x[i,j]), no? (Putting aside the risk of spurious overflow in an un-rescaled norm algorithm.)
Also, in knorm, the result should be initialized to res = zero(real(T)) in order to give a real-type result for complex vectors.
(And why do kdot and knorm have an n argument that is not used? Also, aren’t the loops in the wrong order for locality?)
1 Like
Another alternative might be rescale your problem to implicitly use the desired norm. It depends on what, precisely, you want to do, and why you want to use a different norm or inner product.
Suppose you want to work with a weighted norm \Vert x \Vert_W = \sqrt{x^* W x} = \sqrt{(Rx)^* (Rx)} = \Vert Rx \Vert_2, where W = R^* R is a Hermitian positive-definite weight matrix and R is e.g. a Cholesky factor (or perhaps the matrix square root in cases where it is easy like a diagonal W).
Now, suppose you are solving Ax = b and you want GMRES to minimize the residual \Vert Ax-b\Vert_W = \Vert RAx-Rb\Vert_2 in your W norm, instead of in the usual Euclidean norm, over x in the Krylov space at each step. By inspection, however, that is equivalent nearly equivalent to applying GMRES to the equations RAx = Rb (except that the Krylov space changes). That is, you pass Rb for your right-hand-side and RA as your matrix (perhaps doing the multiplication lazily as a composition, e.g. via LinearMaps.jl), much like a left preconditioner.
(If A is Hermitian/real-symmetric and you want to take advantage of Lanczos, then RA breaks the symmetry, but you can restore the symmetry by a further change of variables x = R^* y, much like a right preconditioner.)
Note that this is not equivalent to changing all your inner products: the basis for the Krylov space is still orthonormalized in the Euclidean inner product. But I’m not sure why you should care about that?
PS. If your norm is not L2-like, e.g. you want to use an L1 norm or some such radical change, then it is a very different algorithm than GMRES and you probably need to re-think and re-implement it from scratch.
1 Like