JOSS paper: changes across different Julia versions

[MohamedLaghdafHABIBOULLAH]

The example in the paper is not reproduced in Julia 1.12 and produces the following error:

julia> include("example.jl")
[ Info:  outer  inner     f(x)     h(x)  √(ξ1/ν)        ρ        σ      ‖x‖      ‖s‖      ‖B‖ R2N
[ Info:      0      1  7.6e+03  7.8e+02  1.5e+02  4.5e-01  7.4e-04  2.8e+01  1.5e+02  1.0e+00 =
[ Info:      1    200  2.9e+03  6.5e+02  1.1e+02 -7.2e-03  7.4e-04  1.4e+02  3.8e+04  1.6e+00 ↗
[ Info:      2    200  2.9e+03  6.5e+02  1.1e+02 -9.2e-03  2.2e-03  1.4e+02  2.3e+04  1.6e+00 ↗
[ Info:      3    200  2.9e+03  6.5e+02  1.1e+02 -1.4e-02  6.7e-03  1.4e+02  1.0e+04  1.6e+00 ↗
[ Info:      4    200  2.9e+03  6.5e+02  1.1e+02 -2.4e-02  2.0e-02  1.4e+02  4.0e+03  1.6e+00 ↗
[ Info:      5    199  2.9e+03  6.5e+02  1.1e+02 -4.0e-02  6.0e-02  1.4e+02  1.3e+03  1.6e+00 ↗
[ Info:      6     82  2.9e+03  6.5e+02  1.1e+02 -3.6e-02  1.8e-01  1.4e+02  4.4e+02  1.6e+00 ↗
[ Info:      7     23  2.9e+03  6.5e+02  1.1e+02  1.0e-01  5.4e-01  1.4e+02  1.5e+02  1.6e+00 =
[ Info:      8      8  1.1e+03  1.2e+03  3.3e+01  1.6e+00  5.4e-01  2.4e+02  2.3e+01  1.0e+00 ↘
[ Info:      9      9  8.6e+02  8.2e+02  1.8e+01  1.5e+00  1.8e-01  2.3e+02  1.9e+01  1.0e+00 ↘
[ Info:     10     21  6.1e+02  7.8e+02  1.3e+01  1.1e+00  6.0e-02  2.3e+02  3.7e+01  1.0e+00 ↘
[ Info:     11     25  3.8e+02  7.1e+02  3.8e+00  1.5e+00  2.0e-02  2.2e+02  1.5e+01  9.8e-01 ↘
[ Info:     12     24  3.6e+02  6.7e+02  3.3e+00  1.7e+00  6.7e-03  2.2e+02  2.2e+01  4.6e-01 ↘
[ Info:     13    200  3.4e+02  6.3e+02  2.0e+00  1.0e+00  2.2e-03  2.3e+02  4.2e-01  4.8e+00 ↘
[ Info:     14    200  3.4e+02  6.3e+02  2.0e+00  1.0e+00  7.4e-04  2.3e+02  4.1e-01  4.9e+00 ↘
[ Info:     15    200  3.4e+02  6.3e+02  2.1e+00  1.0e+00  2.5e-04  2.3e+02  1.6e+00  1.3e+00 ↘
[ Info:     16    200  3.4e+02  6.3e+02  2.9e+00  1.2e+00  8.2e-05  2.3e+02  1.1e+01  2.9e-01 ↘
[ Info:     17    200  3.4e+02  6.0e+02  2.4e+00  1.1e+00  2.7e-05  2.2e+02  8.1e+00  3.3e-01 ↘
[ Info:     18    134  3.3e+02  5.8e+02  1.8e+00  1.0e+00  9.1e-06  2.2e+02  9.4e+01  1.7e+00 ↘
[ Info:     19    200  3.3e+02  4.8e+02  2.5e+00  8.5e-01  3.0e-06  2.3e+02  1.4e+00  1.8e+00 =
[ Info:     20     20  3.3e+02  4.8e+02  2.1e+00  1.2e+00  3.0e-06  2.3e+02  8.6e+00  6.1e-01 ↘
...
[ Info:    140     88  1.9e+02  2.2e+02  4.5e-02  1.6e+00  2.6e-07  9.5e+02  2.2e+01  2.7e-02 ↘
mks = 1.731127610666441e7
ψ(s) = 9688.435470120465
hkn = 9688.435470120465
ν1 = 65329.83168707116
norm s1 = 46370.51666522656
norm Bk = 1.5209597991248702e-5
ERROR: LoadError: R2N: failed to compute a step: ξ = -1.731105604181186e7

If we print some parameters, we observe that $$B_k$$ at this iteration and $$\sigma_k$$ are very small, and consequently $$\nu_1$$ and $$s_1$$ are large.

When displaying the sub-iteration in question, we obtain:

[ Info:    140     88  1.9e+02  2.2e+02  4.5e-02  1.6e+00  2.6e-07  9.5e+02  2.2e+01  2.7e-02 ↘
### Sous-itérations à l'itération de l'erreur
[ Info:   iter     f(x)     h(x)   √(ξ/ν)        ρ        σ      ‖x‖      ‖s‖ R2
[ Info:      0  1.7e+07  9.7e+03  7.5e+02 -5.3e+02  1.5e-05  4.6e+04  4.9e+07 ↗
[ Info:      1  1.7e+07  9.7e+03  7.5e+02 -1.8e+02  4.6e-05  4.6e+04  1.6e+07 ↗
[ Info:      2  1.7e+07  9.7e+03  7.5e+02 -5.8e+01  1.4e-04  4.6e+04  5.5e+06 ↗
[ Info:      3  1.7e+07  9.7e+03  7.5e+02 -1.9e+01  4.1e-04  4.6e+04  1.8e+06 ↗
[ Info:      4  1.7e+07  9.7e+03  7.5e+02 -5.5e+00  1.2e-03  4.6e+04  6.1e+05 ↗
[ Info:      5  1.7e+07  9.7e+03  7.5e+02 -1.2e+00  3.7e-03  4.6e+04  2.0e+05 ↗
[ Info:      6  1.7e+07  9.7e+03  7.5e+02  2.7e-01  1.1e-02  4.6e+04  6.7e+04 =
[ Info:      7  3.4e+06  8.1e+03  3.6e+02  3.6e-01  1.1e-02  2.2e+04  3.2e+04 =
[ Info:      8 -6.6e+05  7.5e+03  2.6e+02  1.1e+00  1.1e-02  1.7e+04  2.3e+04 ↘
[ Info:      9 -7.3e+06  9.3e+03  4.7e+02  2.8e+00  3.7e-03  3.3e+04  1.3e+05 ↘
[ Info:     10 -1.7e+08  2.2e+04  2.2e+03  6.6e+00  1.2e-03  1.6e+05  1.8e+06 ↘
...
[ Info:    196 -1.2e+306  6.2e+78 1.6e+152  7.4e+00  6.2e+13 1.3e+154 2.4e+138 ↘
[ Info:    197 -1.2e+306  6.2e+78 1.8e+152     -Inf  2.1e+13 1.3e+154 8.6e+138 ↗
[ Info:    198 -1.2e+306  6.2e+78 1.6e+152  7.0e+00  6.2e+13 1.3e+154 2.7e+138 ↘
[ Info:    199 -1.2e+306  6.2e+78 1.9e+152     -Inf  2.1e+13 1.3e+154 9.1e+138 ↗
ERROR: LoadError: R2N: failed to compute a step: ξ = -1.731105604181186e7

It appears that the subproblem is unbounded. But, the iteration should be automatically classified as unsuccessful, since $$h$$ is prox-bounded.
I believe that R2N requires improved robustness when $$B_k$$ (SR1) is ill-conditioned.

cc @MaxenceGollier

[MaxenceGollier]

@MohamedLaghdafHABIBOULLAH, I will write it here instead of zulip so others can read.

I think what we should first do is to add more statuses to the set_status! function, in particular for your case, we could add an unbounded case. This would also improve JSO compliance, we should write for each one of

julia> SolverCore.show_statuses()

STATUSES:
  :acceptable     => solved to within acceptable tolerances
  :exception      => unhandled exception
  :first_order    => first-order stationary
  :infeasible     => problem may be infeasible
  :max_eval       => maximum number of function evaluations
  :max_iter       => maximum iteration
  :max_time       => maximum elapsed time
  :neg_pred       => negative predicted reduction
  :not_desc       => not a descent direction
  :small_residual => small residual
  :small_step     => step too small
  :stalled        => stalled
  :unbounded      => objective function may be unbounded from below
  :unknown        => unknown
  :user           => user-requested stop

This will also eliminate throwing errors whenever $$\xi$$ is negative as we could just set the status :not_desc and i have use cases for :small_step and :stalled.
I will open a separate issue + PR for this.

The question is then, what should R2N do in each case, should it just reject the step ? I think in cases like :small_step this wouldn't be right.

Another needed improvement is the lack of robustness observed for approximating the spectral norm.
For the LSR1, I think it would be more efficient and robust to use the fact that

$$\|B_k\|_2 = \|B_0 + \sum_{j=1}^k u_j u_j^T\|_2 \leq \|B_0\|_2 + \sum_{j=1}^k \|u_j\|_2^2.$$

[MohamedLaghdafHABIBOULLAH]

Thank you, very useful.
It remains unclear how to compute $$u_i$$ for LSR1.
For LBFGS, we have easier way, at each update

$$\|B_{new}\| \le \|B_{old}\| + \|y\|^2/|s^T y|$$

Thus,

$$\|B_{k}\| \le \|B_{k,0}\| + \sum_{j=1}^m \|y_j\|^2 / |s_{j}^T y_j|$$

[MohamedLaghdafHABIBOULLAH]

The simplest approach is to reject the step. However, if $$\sigma_k$$ is already sufficiently small, rejecting the step may hinder progress: the algorithm has already defined this region as "flat" (very successful) and would need to penalize it. A more efficient alternative is to increase the norm (or estimated norm) of $$B_k$$ without modifying $$\sigma_k$$.

[d1Lab]

Both suggestions are good but I am afraid that the norm bound can be big which may slow the convergence. A suggestion, to refine the estimation of the norm, could be to compute the prox. Then, if xi is negative, increase the number of iterations in the power iteration method til you get a positive value of xi. That way you could guarantee that the estimation of the norm is more accurate.

[MohamedLaghdafHABIBOULLAH]

Here are the steps to run and reproduce the experiments on macOS using Julia 1.11.7.

Running the experiments (branch paper) — RegularizedOptimization.jl

1) Clone the repository on the paper branch

git clone --branch paper --single-branch https://github.com/JuliaSmoothOptimizers/RegularizedOptimization.jl.git
cd RegularizedOptimization.jl

Alternative (if you already cloned the repository):

git fetch origin
git checkout paper

---

2) Start Julia-1.11 and activate the examples environment

From the RegularizedOptimization.jl folder, start Julia:

juliaup add 1.11.7
julia +1.11.7

Then in the Julia REPL:

cd("paper/examples/")
using Pkg
Pkg.activate(".")
Pkg.instantiate()

---

3) Run the example

Still in the Julia REPL:

include("example.jl")

With MacOs you should be able to recover the example of the paper

[ Info:  outer  inner     f(x)     h(x)  √(ξ1/ν)        ρ        σ      ‖x‖      ‖s‖      ‖B‖ R2N 
[ Info:      0      1  7.6e+03  7.8e+02  1.5e+02  4.5e-01  7.4e-04  2.8e+01  1.5e+02  1.0e+00 =
...
[ Info:    179    200  1.9e+02  2.0e+02  2.1e-02 -8.6e-01  6.3e-05  1.0e+03  2.3e+02  1.5e-01 ↗
[ Info:    180    200  1.9e+02  2.0e+02  2.1e-02 -8.6e-01  1.9e-04  1.0e+03  1.9e+02  1.5e-01 ↗
[ Info:    181    200  1.9e+02  2.0e+02  2.1e-02 -6.8e-01  5.6e-04  1.0e+03  8.1e+01  1.5e-01 ↗
[ Info:    182    200  1.9e+02  2.0e+02  2.1e-02  4.9e-01  1.7e-03  1.0e+03  1.2e+01  1.5e-01 =
[ Info:    183     24  1.9e+02  2.1e+02  1.9e-02  1.1e+00  1.7e-03  1.0e+03  2.6e+00  1.9e-01 ↘
[ Info:    184      0  1.9e+02  2.1e+02  7.5e-03  1.1e+00  5.6e-04  1.0e+03  2.6e+00  1.3e-01                
[ Info: R2N: terminating with √(ξ1/ν) = 0.0074767797420856276
"Execution stats: first-order stationary"

4) Run Benchmarks

julia> include("Benchmark.jl")
main (generic function with 1 method)

julia> main()

=== SVM: solver comparison ===

→ TR
  status          = first_order
  time (s)        = 4.3677
  # f eval        = 347
  # ∇f eval       = 291
  # prox calls    = 4037
  final objective = 179.837

→ R2N
  status          = first_order
  time (s)        = 1.8588
  # f eval        = 185
  # ∇f eval       = 101
  # prox calls    = 27932
  final objective = 192.4927

→ LM
  status          = first_order
  time (s)        = 16.9891
  # f eval        = 6
  # ∇f eval       = 2876
  # prox calls    = 1001
  final objective = 201.1862

→ LMTR
  status          = first_order
  time (s)        = 10.8845
  # f eval        = 11
  # ∇f eval       = 1614
  # prox calls    = 432
  final objective = 188.2737

=== Full Benchmark Table ===
(name = "TR", status = "first_order", time = 4.3677, fe = 347, ge = 291, prox = 4037, obj = 179.837)
(name = "R2N", status = "first_order", time = 1.8588, fe = 185, ge = 101, prox = 27932, obj = 192.4927)
(name = "LM", status = "first_order", time = 16.9891, fe = 6, ge = 2876, prox = 1001, obj = 201.1862)
(name = "LMTR", status = "first_order", time = 10.8845, fe = 11, ge = 1614, prox = 432, obj = 188.2737)

You can check the behavior across different operating systems here:
https://github.com/MaxenceGollier/RegularizedOptimization.jl/actions/runs/21079257881/job/60628720557?pr=3

As shown in the logs, the results vary across operating systems and Julia versions. This kind of variability is expected in numerical experiments (differences in floating-point behavior and underlying linear algebra libraries), and it is not specific to Julia.
Note: On macOS with Julia 1.12, we observed a specific issue, a fix is proposed in this PR.