CR: Update iterate on first iteration for negative curvature with line search and when zero curvature is detected

[farhadrclass]

This PR fixes an issue in the Conjugate Residual (CR) implementation (in src/cr.jl) where, on the first iteration (iter == 0), if the computed curvature is non-positive (i.e. either zero or negative), the algorithm may return without updating the iterate—even though we have
$$b = -\nabla f(x_k) \neq 0.$$

To address this, when line search is enabled and a non-positive curvature (zero or negative) is detected at iteration 0, we now explicitly update the iterate by performing

x .= b

This change guarantees that the method behaves consistently by setting x appropriately when curvature information is either zero or negative during the first iteration.

---

Changes Made:

1. First Iteration Check (iter == 0):

   • When either zero or negative curvature is detected (i.e. s <= 0), and linesearch is enabled, update x to b before returning.

2. Zero Curvature Check (general case):

   • (Existing handling) When a zero curvature is detected later in the method (e.g. at a different check), we continue to update x to b if linesearch is enabled.

[dpo]

A unit test would be good.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 93.55%. Comparing base (9536ef7) to head (0526df6).
Report is 45 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #958      +/-   ##
==========================================
- Coverage   94.68%   93.55%   -1.13%     
==========================================
  Files          45       47       +2     
  Lines        8027     9177    +1150     
==========================================
+ Hits         7600     8586     +986     
- Misses        427      591     +164     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[farhadrclass]

@dpo added couple tests and added some function to create the needed systems.
I think in future we need to expand the tests to CG solver as well since it is missing them

[dpo]

LGTM. How about you @amontoison ?

[amontoison]

I think we need to replace iter == 0 && kcopy!(n, x, b) by iter == 0 && kcopy!(n, x, p) because p is the initial residual b - Ax0.

We could have a wrong solution if linesearch is combined with warm-start (x0 != 0).

[farhadrclass]

@amontoison Thanks for pointing it out. However, I don’t think we should do that.
If we replace b with p, we’ll end up setting x to zero because, in line 170, we have:

if ρ == 0
    stats.niter = 0
    stats.solved, stats.inconsistent = true, false
    stats.timer = start_time |> ktimer
    stats.status = "p is a zero-curvature direction"
    history && push!(ArNorms, zero(T))
    solver.warm_start = false
    linesearch && kcopy!(n, x, p)  # x ← p
    return solver
end

I assume that when we use a warm start, the probability of p being zero is lower, so the chance of triggering this condition would be reduced.

The original code was written this way. If we set p to x instead of b to x, my R2N solver would encounter the same issue as before (i.e., x == 0).

What do you think?

[amontoison]

The original code was written this way. If we set p to x instead of b to x, my R2N solver would encounter the same issue as before (i.e., x == 0).

p and b have the same content at iter == 0, I don't understand how it is possible to have the old issue.

[amontoison]

@farhadrclass This is not a question of probability; I just want to ensure that x contains a relevant solution for the user.

It’s quite complex to test all the options together (like linesearch, preconditioner, warm-start, etc.).
I should frame the problem like this:
If x0 = 0, we want to solve A x = r0 where r0 = b - A x0 = b, and we set x = b = p if we encounter negative curvature.
If x0 ≠ 0, we solve a different system in Krylov.jl, where the system is shifted.
We solve A Δx = r0 where r0 = b - A x0, and at the end, we compute x = Δx + x0.
If we encounter negative curvature on the shifted system, we should set Δx = r0 from what I understand from your previous messages.
Therefore, Krylov.jl needs to return x = x0 + r0. Since at iter == 0, p = r0, we need something like this:

linesearch && kcopy!(n, x, p)  # x ← p
warm_start && kaxpy!(n, one(FC), Δx, x)

@dpo I think I made a mistake in all the Krylov methods where I implemented warm start...

For historical reasons, we checked whether the initial residual was zero or not.
In the case of a warm start, if the residual is zero, then Δx = 0, but we still need to set x = x0.
I only do x += x0 at the end of the code and not in specific corner cases where a solution is returned earlier...

It should be fixed by #961.
If we forget the line warm_start && kaxpy!(n, one(FC), Δx, x), I still think that we should replace

linesearch && kcopy!(n, x, b)  # x ← b

by

linesearch && kcopy!(n, x, p)  # x ← p

[dpo]

I think the linesearch and warmstart options are incompatible. In an inexact Newton method, we want to solve $Bs = -g$ where $s$ is a step, i.e., something that starts at zero, so that we can later update $x_{k+1} = x_k + s$ in the optimization method.

There is something special about the rhs $-g$ here: it's a descent direction. If we start from $s_0 = 0$ and encounter negative curvature at the first CG or CR iteration (along the direction $-g$, as it happens), we want to return $s = -g$.

If we supplied a "better" initial guess $s_0$ and solved instead $B \Delta s = r$, where $r = B s_0 + g$, and encountered negative curvature at the first iteration, we still want to return $-g$, because we don't know whether $r$ is descent direction.

If we encounter negative curvature at iteration $k \geq 1$, the code would currently return $s_0 + \Delta s_{k-1}$, but I'm really not sure what that iterate would represent from the point of view of the optimization method.

At least for the time being, I think we should make sure that linesearch and warmstart are not true at the same time.

The same goes for radius and warmstart.

Maybe linesearch is not a good name for that option because what I'm describing is specific to optimization. If you're solving nonlinear equations, there is no natural concept of curvature.

[amontoison]

Thanks for the explanation @dpo!

You can merge the PR and I will take care to return errors if some options like linesearch / radius and warmstart are used together.

[dpo]

Thanks @farhadrclass !

[farhadrclass]

thanks @amontoison for help and working on the follow-up
@dpo thank you, do we need to update the release number (version number)?

[dpo]

Sure. This is a bug fix.

[farhadrclass]

Can I trigger the V0.9.10?