diff --git a/src/DiagonalHessianApproximation.jl b/src/DiagonalHessianApproximation.jl
index 15bdf2e..40a5993 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,58 @@ 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},
+  s::V,
+  y::V,
+) where {T <: Real, I <: Integer, V <: AbstractVector{T}, F}
+  sNorm = norm(s, 2)
+  if sNorm == 0
+    error("Cannot update DiagonalQN operator with s=0")
+  end
+  sNorm2 = sNorm^2
+  sT_y = dot(s, y) / sNorm2
+  B.d .= abs.(y)
+  B.d .*= sum(B.d) / sT_y
+  return B
+end
diff --git a/test/test_diag.jl b/test/test_diag.jl
index 9349039..3c02886 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