diff --git a/src/quasi-newton.jl b/src/quasi-newton.jl
index 6b6ecb0..974eb98 100644
--- a/src/quasi-newton.jl
+++ b/src/quasi-newton.jl
@@ -4,6 +4,7 @@ export AbstractDiagonalQNModel,
   LSR1Model,
   DiagonalPSBModel,
   DiagonalAndreiModel,
+  DiagonalBFGSModel,
   SpectralGradientModel,
   get_model,
   get_op
@@ -130,6 +131,23 @@ function DiagonalAndreiModel(
   return DiagonalQNModel{T, S, typeof(nlp), typeof(nlp.meta), typeof(op)}(nlp.meta, nlp, op)
 end
 
+"""
+    DiagonalBFGSModel(nlp; d0 = fill!(S(undef, nlp.meta.nvar), 1.0))
+
+Construct a `DiagonalBFGSModel` from another type of nlp, in which the Hessian is approximated
+via a diagonal BFGS quasi-Newton operator.
+`d0` is the initial approximation of the diagonal of the Hessian, and by default a vector of ones.
+See the
+[`DiagonalBFGS operator documentation`](https://juliasmoothoptimizers.github.io/LinearOperators.jl/stable/reference/#LinearOperators.DiagonalBFGS).
+"""
+function DiagonalBFGSModel(
+  nlp::AbstractNLPModel{T, S};
+  d0::S = fill!(S(undef, nlp.meta.nvar), one(T)),
+) where {T, S}
+  op = DiagonalBFGS(d0)
+  return DiagonalQNModel{T, S, typeof(nlp), typeof(nlp.meta), typeof(op)}(nlp.meta, nlp, op)
+end
+
 """
     SpectralGradientModel(nlp; σ = 1.0)
 
diff --git a/test/nlp/quasi-newton.jl b/test/nlp/quasi-newton.jl
index 02f9982..efdb376 100644
--- a/test/nlp/quasi-newton.jl
+++ b/test/nlp/quasi-newton.jl
@@ -14,6 +14,7 @@
         (LBFGSModel, LBFGSOperator),
         (DiagonalPSBModel, DiagonalPSB),
         (DiagonalAndreiModel, DiagonalAndrei),
+        (DiagonalBFGSModel, DiagonalBFGS),
         (SpectralGradientModel, SpectralGradient),
       ],
       T in [Float64, Float32],
@@ -24,7 +25,7 @@
       m = nlp.meta.ncon
 
       s, y = randn(T, n), randn(T, n)
-      if QNO ∈ (DiagonalPSB, DiagonalAndrei)
+      if QNO ∈ (DiagonalPSB, DiagonalAndrei, DiagonalBFGS)
         B = QNO(ones(T, n))
       elseif QNO == SpectralGradient
         B = QNO(one(T), n)
@@ -153,7 +154,14 @@
   end
 
   @testset "Show" begin
-    for QNM ∈ [LSR1Model, LBFGSModel, DiagonalPSBModel, DiagonalAndreiModel, SpectralGradientModel]
+    for QNM ∈ [
+      LSR1Model,
+      LBFGSModel,
+      DiagonalPSBModel,
+      DiagonalAndreiModel,
+      DiagonalBFGSModel,
+      SpectralGradientModel,
+    ]
       nlp = QNM(SimpleNLPModel())
       io = IOBuffer()
       show(io, nlp)