Integrating a preconditioner into Trunk's trust-region subsolver

[mpf]

I'm considering adding functionality to allow users to pass a preconditioner into Trunk's trust-region subsolver and would appreciate some guidance before I begin making changes.

Options for Passing the Preconditioner

There seem to be at least two options here:

1. Modify SolverCore.solve!: Introduce a new keyword argument M=I in SolverCore.solve! and pass it directly to Krylov.solve!. This approach is straightforward and would allow updating the preconditioner with the latest iterate just before the call to Krylov.solve!.

2. Use a Custom CgSolver: Pass in to trunk a custom CgSolver that includes the preconditioner. I'm not sure of the details here, particularly how the preconditioner would be updated at each iteration.

The first approach seems simpler and more direct.

Modifying Step-to-the-Boundary Calculation

As highlighted by @amontoison in this Krylov.jl issue, the step-to-the-boundary computation in to_boundary needs updating to solve the equation

‖x + σd‖_M = radius

where ‖v‖_M^2 = v'Mv. This modification requires changes to both the to_boundary function and the corresponding call in cg.jl:

σ = radius > 0 ? maximum(to_boundary(n, x, p, radius, dNorm2=pNorm²)) : α

Any advice or suggestions on these approaches would be greatly appreciated.

Thanks!

[amontoison]

Hi @mpf!

The best solution is 1., Krylov.jl already provides a way to change preconditioner at each call of a Krylov method. We just need to make it possible to provide a preconditioner from a solver here in JSOSolvers.jl.

For the function to_boundary, I wanted to limit the number of products with M. The simplest solution is to perform three products with M to compute ‖d‖_M, xᴴMd and ‖x‖_M at each iteration. Can we do better?

[mpf]

For the function to_boundary, I wanted to limit the number of products with M. The simplest solution is to perform three products with M to compute ‖d‖_M, xᴴMd and ‖x‖_M at each iteration. Can we do better?

Would it be possible to get away with only 2 products with M? Note

||x + σd‖²_M = σ² (d'Md) + σ (x'Md) + (x'Mx) = σ² (d'v) + σ (x'v) + (x'w),

where v:=Md and w:=Mx.

[amontoison]

I think it's optimal with 2 products with M. :+1: Michael, I will go to SIAM LA24 next week but I can add a preconditioner M in the function to_boundary when I will be back to Montréal.