From c53cf06ccc5bbeed2e5fa756c44f968ce1d2ea9f Mon Sep 17 00:00:00 2001
From: MohamedLaghdafHABIBOULLAH <mohamed-laghdaf-2.habiboullah@polymtl.ca>
Date: Thu, 12 Sep 2024 12:49:27 -0400
Subject: [PATCH 1/3] Add Diagonal BFGS

---
 src/DiagonalHessianApproximation.jl | 58 ++++++++++++++++++++++++++++-
 1 file changed, 57 insertions(+), 1 deletion(-)

diff --git a/src/DiagonalHessianApproximation.jl b/src/DiagonalHessianApproximation.jl
index 15bdf2ee..4e6de932 100644
--- a/src/DiagonalHessianApproximation.jl
+++ b/src/DiagonalHessianApproximation.jl
@@ -1,4 +1,4 @@
-export DiagonalPSB, DiagonalAndrei, SpectralGradient
+export DiagonalPSB, DiagonalAndrei, SpectralGradient, DiagonalBFGS
 
 """
     DiagonalPSB(d)
@@ -203,3 +203,59 @@ function push!(
   B.d[1] = dot(s, y) / dot(s, s)
   return B
 end
+
+"""
+    DiagonalBFGS(d)
+
+A diagonal approximation of the BFGS update inspired by 
+Marnissi, Y., Chouzenoux, E., Benazza-Benyahia, A., & Pesquet, J. C. (2020). 
+Majorize–minimize adapted Metropolis–Hastings algorithm.
+https://ieeexplore.ieee.org/abstract/document/9050537. 
+
+# Arguments
+
+- `d::AbstractVector`: initial diagonal approximation.
+"""
+mutable struct DiagonalBFGS{T <: Real, I <: Integer, V <: AbstractVector{T}, F} <:
+               AbstractDiagonalQuasiNewtonOperator{T}
+  d::V # Diagonal of the operator
+  nrow::I
+  ncol::I
+  symmetric::Bool
+  hermitian::Bool
+  prod!::F
+  tprod!::F
+  ctprod!::F
+  nprod::I
+  ntprod::I
+  nctprod::I
+  args5::Bool
+  use_prod5!::Bool # true for 5-args mul! and for composite operators created with operators that use the 3-args mul!
+  allocated5::Bool # true for 5-args mul!, false for 3-args mul! until the vectors are allocated
+end
+
+@doc (@doc DiagonalBFGS) function DiagonalBFGS(d::AbstractVector{T}) where {T <: Real}
+  prod = (res, v, α, β) -> mulSquareOpDiagonal!(res, d, v, α, β)
+  n = length(d)
+  DiagonalBFGS(d, n, n, true, true, prod, prod, prod, 0, 0, 0, true, true, true)
+end
+
+# update function
+# s = x_{k+1} - x_k
+# y = ∇f(x_{k+1}) - ∇f(x_k)
+function push!(
+  B::DiagonalBFGS{T, I, V, F},
+  s0::V,
+  y0::V,
+) where {T <: Real, I <: Integer, V <: AbstractVector{T}, F}
+  s0Norm = norm(s0, 2)
+  if s0Norm == 0
+    error("Cannot update DiagonalQN operator with s=0")
+  end
+  # sᵀBs = sᵀy can be scaled by ||s||² without changing the update
+  s = (si / s0Norm for si ∈ s0)
+  y = (yi / s0Norm for yi ∈ y0)
+  sT_y = dot(s, y)
+  B.d .= sum(abs.(y)) / sT_y .* abs.(y) 
+  return B
+end

From 1b49400e77de6319eaf3321d10679fc99ccb25f9 Mon Sep 17 00:00:00 2001
From: MohamedLaghdafHABIBOULLAH <mohamed-laghdaf-2.habiboullah@polymtl.ca>
Date: Mon, 4 Nov 2024 14:01:11 -0500
Subject: [PATCH 2/3] add unit test for DiagonalBFGS

---
 test/test_diag.jl | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/test/test_diag.jl b/test/test_diag.jl
index 93490390..3c02886b 100644
--- a/test/test_diag.jl
+++ b/test/test_diag.jl
@@ -120,6 +120,10 @@ end
   mul!(u, C, v)
   @test (@allocated mul!(u, C, v)) == 0
   @test (@wrappedallocs push!(C, u, v)) == 0
+  D = DiagonalBFGS(d)
+  mul!(u, D, v)
+  @test (@allocated mul!(u, D, v)) == 0
+  @test (@wrappedallocs push!(D, u, v)) == 0
 end
 
 @testset "reset" begin

From 0012c31c8e20fe9fb52d50b187cd3ec65270eb4e Mon Sep 17 00:00:00 2001
From: MohamedLaghdafHABIBOULLAH <mohamed-laghdaf-2.habiboullah@polymtl.ca>
Date: Thu, 20 Feb 2025 21:37:00 -0500
Subject: [PATCH 3/3] Remove allocation from DBFGS update

---
 src/DiagonalHessianApproximation.jl | 17 ++++++++---------
 1 file changed, 8 insertions(+), 9 deletions(-)

diff --git a/src/DiagonalHessianApproximation.jl b/src/DiagonalHessianApproximation.jl
index 4e6de932..40a59930 100644
--- a/src/DiagonalHessianApproximation.jl
+++ b/src/DiagonalHessianApproximation.jl
@@ -245,17 +245,16 @@ end
 # y = ∇f(x_{k+1}) - ∇f(x_k)
 function push!(
   B::DiagonalBFGS{T, I, V, F},
-  s0::V,
-  y0::V,
+  s::V,
+  y::V,
 ) where {T <: Real, I <: Integer, V <: AbstractVector{T}, F}
-  s0Norm = norm(s0, 2)
-  if s0Norm == 0
+  sNorm = norm(s, 2)
+  if sNorm == 0
     error("Cannot update DiagonalQN operator with s=0")
   end
-  # sᵀBs = sᵀy can be scaled by ||s||² without changing the update
-  s = (si / s0Norm for si ∈ s0)
-  y = (yi / s0Norm for yi ∈ y0)
-  sT_y = dot(s, y)
-  B.d .= sum(abs.(y)) / sT_y .* abs.(y) 
+  sNorm2 = sNorm^2
+  sT_y = dot(s, y) / sNorm2
+  B.d .= abs.(y)
+  B.d .*= sum(B.d) / sT_y
   return B
 end