From 841138db9fef1c8d065d33032177f444a43926f0 Mon Sep 17 00:00:00 2001
From: "RAYNAUD Paul (raynaudp)" <raynaudp@grenoble-inp.lan>
Date: Fri, 24 Jun 2022 16:17:26 -0300
Subject: [PATCH] fix s and y collinear by cancelling data.scaling if they are
 collinear

---
 src/lsr1.jl       |  9 +++++++++
 test/test_lsr1.jl | 10 ++++++++++
 2 files changed, 19 insertions(+)

diff --git a/src/lsr1.jl b/src/lsr1.jl
index 6afe1d82..aad0124c 100644
--- a/src/lsr1.jl
+++ b/src/lsr1.jl
@@ -136,6 +136,13 @@ function push!(op::LSR1Operator, s::AbstractVector, y::AbstractVector)
   ϵ = eps(eltype(op))
   well_defined = abs(dot(ymBs, s)) ≥ ϵ + ϵ * norm(ymBs) * sNorm
 
+  # check if y and s are collinear: make a componentwise division and check if the resulting components get approximately the same factor.
+  and(bool1, bool2) = bool1 && bool2 
+  s_y = (s ./ y)
+  collinear = (length(s) == 1) || mapreduce(i-> i ≈ s_y[1], and, s_y) # ≈ to avoid numerical issues
+  collinear && (tmp_scaling = data.scaling) # save data.scaling
+  collinear && (data.scaling = false) # if the scaling is applied then all the components of Matrix(op_updated) are NaN
+
   sufficient_curvature = true
   scaling_condition = true
   if data.scaling
@@ -179,6 +186,8 @@ function push!(op::LSR1Operator, s::AbstractVector, y::AbstractVector)
     end
   end
 
+  collinear && (data.scaling = tmp_scaling) # set data.scaling to its original value
+
   return op
 end
 
diff --git a/test/test_lsr1.jl b/test/test_lsr1.jl
index 2fb9a93a..a274bb94 100644
--- a/test/test_lsr1.jl
+++ b/test/test_lsr1.jl
@@ -110,6 +110,16 @@ function test_lsr1()
     resid = vals[end]
     @test norm(resid) ≤ sqrt(eps(eltype(B)))
   end
+
+  @testset "LSR1 (s,y) collinear" begin
+    n = 10
+    lsr1 = LSR1Operator(n)
+    s = rand(n)
+    alpha = rand(1) * 1000
+    y = s .* alpha
+    push!(lsr1, s, y)
+    Matrix(lsr1) != I
+  end
 end
 
 test_lsr1()