Any conjugate gradient method allowing a custom inner product?

Specific Domains Numerics
1

Is there any conjugate gradient (CG) implementation in Julia allowing one to specify the inner product (with a fixed symmetric and positive definite matrix M, say)?

In my use case, I know my system matrix A is symmetric and positive definite with respect to M, i.e., x' * M * A * x > 0 for all x != 0 and y' * M * (A * x) == (A * y)' * M * x for all x, y. Ideally, I would like to avoid scaling the system matrix and vectors myself for using a CG method. Just curious to know whether there’s such a flexible package out there.

2

I don’t know the answer to the question about specific implementations of CG, but since your property of being symmetric and positive definite with respect to M is equivalent to MA being symmetric positive definite, it seems reasonable to just solve MA x = Mb using the CG routine from IterativeSolvers.jl and form MA as a product of LinearMaps:

using IterativeSolvers
using LinearMaps
C = LinearMap(M) * LinearMap(A)
x=cg(C, M*b)

That’s pretty easy. In fact, since CG includes more inner products than matrix-vector multiplies, I think this actually involves fewer multiplies by M than if you tried to implement CG for a different inner product. If multiplying by M is nontrivial, this is probably faster than what you were hoping to do.

1 Like
3

Conjugate Gradients · IterativeSolvers.jl has the optional argument Pl to add a left precondition. (That should correspond to a change to the inner product as @mstewart explained)
As the documentation say it uses a different algorithm which is specialised for this setting.

4

It’s a tempting idea, and there may be some trick to make use of a preconditioner, but doing it in a direct way for this problem doesn’t seem like an option, given that PCG assumes that both the preconditioner and A are symmetric positive definite. That’s not the case here.

5

KrylovKit supports any vector type that have custom dot and norm methods, without assuming they are the standard ones. See eg Introduction · KrylovKit.jl

1 Like
6

Thank you all! I can only select one answer and have a hard time selecting one :sweat_smile:

1 Like