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
@stevengj You are right, I forgot to handle the complex numbers with my example of custom workspace.
It is not what I use internally by default in Krylov.jl.
I will fix that, the idea was just to illustrate what to define if you don’t have a classic AbstractVector where the operations norm, dot, etc… are not defined.
The argument n is here to directly call BLAS / rocBLAS / cuBLAS if you want.
@Veenty We need more details to understand what you want as inner product.
If you just want to minimize the residual norm in a different metric than the L2 norm, using a left preconditioner (argument M in Krylov.gmres) could be what you want.
Can’t you get that from the size of the arrays? Why would you need to pass it as a separate argument?
It depends on your custom workspace.
You can also imagine that it is a Ptr{Cvoid} or CuPtr{Cvoid}, and you would need to use unsafe_wrap on it with the length.
Also, why would I want to recompute the length of the vector if it is already available?
You can just ignore this argument if you don’t need it.
In any case, these functions kdot, knorm, kaxpy, etc., are intended for advanced users and specific applications.
You should define your own versions only if the default dispatch can’t handle your type:
Almost nobody asked for this feature before, but I needed to make this custom API available for Oceananigans.jl, where they use “Fields” (3D offset arrays with specific halo regions) as vectors for Krylov workspaces.
This custom API could also be useful for doing distributed Krylov solvers with a custom “MPI” vector type but I don’t want to diverge from the initial topic.
Because (a) passing/storing redundant data is a notorious magnet for bugs and (b) it’s not “computed”, it’s usually stored in any high-level container data structure (e.g. in a Vector).
@stevengj I agree with your points.
For (a), I have some consistency checks at the beginning for the lengths of the vectors. I just compute properties of the vectors outside the call to these methods.
We could discuss the design, but most of these things were developed back when Julia 0.6 (or earlier) was current, and they haven’t evolved much since then.
These functions were initially created because Julia didn’t dispatch to BLAS automatically, so we handled it manually with these “k” routines.
For (b), I remember that the length of a CuVector wasn’t stored in the structure and was recomputed at every length call. I just checked, and it seems that was changed at the end of 2022. However, the GPU support in Krylov.jl was added way earlier, back when CuArrays.jl was still around.
So, quite a long time ago.
@Veenty Please explain what you are trying to achieve. 
1 Like
This is good. But why is this topic put under optimization category?
I edited the category, thanks.
I’m solving an integral equation problem and wanted to see if constructing the Krylov spaces in the natural dot product of the space would give better results. It all boils down into doing the naive dot product with a diagonal matrix in the middle.
I did try the suggestion of @stevengj of using that matrix instead as a left preconditioner, but it was giving me pretty bad results.
Thanks for your explanation @Veenty.
Did you pass the diagonal matrix D as left preconditioner or sqrt.(D) ?
You need sqrt.(D) for recovering the norm that you want.
If you already provide the correct matrix, I suggest to try to plug a custom kdot and knorm to quickly test if it helps.