diff --git a/bindings/AMPL/AMPLModel.cpp b/bindings/AMPL/AMPLModel.cpp
index 5882c4a6..3ee3a4c4 100644
--- a/bindings/AMPL/AMPLModel.cpp
+++ b/bindings/AMPL/AMPLModel.cpp
@@ -129,8 +129,8 @@ namespace uno {
       fint error_flag = 0;
       // prevent ASL to crash by catching all evaluation errors
       Jmp_buf err_jmp_uno;
-      asl->i.err_jmp_ = &err_jmp_uno;
-      asl->i.err_jmp1_ = &err_jmp_uno;
+      this->asl->i.err_jmp_ = &err_jmp_uno;
+      this->asl->i.err_jmp1_ = &err_jmp_uno;
       if (setjmp(err_jmp_uno.jb)) {
          error_flag = 1;
       }
@@ -199,68 +199,6 @@ namespace uno {
       }
    }
 
-   void AMPLModel::compute_lagrangian_hessian_sparsity() {
-      // compute the maximum number of nonzero elements, provided that all multipliers are non-zero
-      // int (*Sphset) (ASL*, SputInfo**, int nobj, int ow, int y, int uptri);
-      const int objective_number = -1;
-      const int upper_triangular = 1;
-      this->number_asl_hessian_nonzeros = static_cast<size_t>((*(this->asl)->p.Sphset)(this->asl, nullptr, objective_number, 1, 1, upper_triangular));
-      this->asl_hessian.reserve(this->number_asl_hessian_nonzeros);
-
-      // sparsity pattern
-      [[maybe_unused]] const fint* asl_column_start = this->asl->i.sputinfo_->hcolstarts;
-      // check that the column pointers are sorted in increasing order
-      assert(in_increasing_order(asl_column_start, this->number_variables + 1) && "AMPLModel::evaluate_lagrangian_hessian: column starts are not ordered");
-   }
-
-   size_t AMPLModel::number_objective_gradient_nonzeros() const {
-      return static_cast<size_t>(this->asl->i.nzo_);
-   }
-
-   size_t AMPLModel::number_jacobian_nonzeros() const {
-      return static_cast<size_t>(this->asl->i.nzc_);
-   }
-
-   size_t AMPLModel::number_hessian_nonzeros() const {
-      return this->number_asl_hessian_nonzeros;
-   }
-
-   const Collection<size_t>& AMPLModel::get_equality_constraints() const {
-      return this->equality_constraints_collection;
-   }
-
-   const Collection<size_t>& AMPLModel::get_inequality_constraints() const {
-      return this->inequality_constraints_collection;
-   }
-
-   const Collection<size_t>& AMPLModel::get_linear_constraints() const {
-      return this->linear_constraints_collection;
-   }
-
-   const SparseVector<size_t>& AMPLModel::get_slacks() const {
-      return this->slacks;
-   }
-
-   const Collection<size_t>& AMPLModel::get_single_lower_bounded_variables() const {
-      return this->single_lower_bounded_variables_collection;
-   }
-
-   const Collection<size_t>& AMPLModel::get_single_upper_bounded_variables() const {
-      return this->single_upper_bounded_variables_collection;
-   }
-
-   const Vector<size_t>& AMPLModel::get_fixed_variables() const {
-      return this->fixed_variables;
-   }
-
-   const Collection<size_t>& AMPLModel::get_lower_bounded_variables() const {
-      return this->lower_bounded_variables_collection;
-   }
-
-   const Collection<size_t>& AMPLModel::get_upper_bounded_variables() const {
-      return this->upper_bounded_variables_collection;
-   }
-
    void AMPLModel::evaluate_lagrangian_hessian(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
          SymmetricMatrix<size_t, double>& hessian) const {
       assert(hessian.capacity() >= this->number_asl_hessian_nonzeros);
@@ -297,6 +235,31 @@ namespace uno {
       this->asl->i.x_known = 0;
    }
 
+   void AMPLModel::compute_hessian_vector_product(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
+         Vector<double>& result) const {
+      fint error_flag = 0;
+      // prevent ASL to crash by catching all evaluation errors
+      Jmp_buf err_jmp_uno;
+      this->asl->i.err_jmp_ = &err_jmp_uno;
+      this->asl->i.err_jmp1_ = &err_jmp_uno;
+      if (setjmp(err_jmp_uno.jb)) {
+         error_flag = 1;
+      }
+
+      // scale by the objective sign
+      objective_multiplier *= this->objective_sign;
+      const int objective_number = -1;
+      // flip the signs of the multipliers: in AMPL, the Lagrangian is f + lambda.g, while Uno uses f - lambda.g
+      this->multipliers_with_flipped_sign = -multipliers;
+
+      // compute the Hessian-vector product
+      (this->asl->p.Hvcomp)(this->asl, result.data(), const_cast<double*>(x.data()), objective_number, &objective_multiplier,
+            const_cast<double*>(this->multipliers_with_flipped_sign.data()));
+      if (error_flag) {
+         throw HessianEvaluationError();
+      }
+   }
+
    double AMPLModel::variable_lower_bound(size_t variable_index) const {
       return this->variable_lower_bounds[variable_index];
    }
@@ -309,6 +272,30 @@ namespace uno {
       return this->variable_status[variable_index];
    }
 
+   const Collection<size_t>& AMPLModel::get_lower_bounded_variables() const {
+      return this->lower_bounded_variables_collection;
+   }
+
+   const Collection<size_t>& AMPLModel::get_upper_bounded_variables() const {
+      return this->upper_bounded_variables_collection;
+   }
+
+   const SparseVector<size_t>& AMPLModel::get_slacks() const {
+      return this->slacks;
+   }
+
+   const Collection<size_t>& AMPLModel::get_single_lower_bounded_variables() const {
+      return this->single_lower_bounded_variables_collection;
+   }
+
+   const Collection<size_t>& AMPLModel::get_single_upper_bounded_variables() const {
+      return this->single_upper_bounded_variables_collection;
+   }
+
+   const Vector<size_t>& AMPLModel::get_fixed_variables() const {
+      return this->fixed_variables;
+   }
+
    double AMPLModel::constraint_lower_bound(size_t constraint_index) const {
       return this->constraint_lower_bounds[constraint_index];
    }
@@ -325,6 +312,18 @@ namespace uno {
       return this->constraint_status[constraint_index];
    }
 
+   const Collection<size_t>& AMPLModel::get_equality_constraints() const {
+      return this->equality_constraints_collection;
+   }
+
+   const Collection<size_t>& AMPLModel::get_inequality_constraints() const {
+      return this->inequality_constraints_collection;
+   }
+
+   const Collection<size_t>& AMPLModel::get_linear_constraints() const {
+      return this->linear_constraints_collection;
+   }
+
    // initial primal point
    void AMPLModel::initial_primal_point(Vector<double>& x) const {
       assert(x.size() >= this->number_variables);
@@ -379,6 +378,18 @@ namespace uno {
       }
    }
 
+   size_t AMPLModel::number_objective_gradient_nonzeros() const {
+      return static_cast<size_t>(this->asl->i.nzo_);
+   }
+
+   size_t AMPLModel::number_jacobian_nonzeros() const {
+      return static_cast<size_t>(this->asl->i.nzc_);
+   }
+
+   size_t AMPLModel::number_hessian_nonzeros() const {
+      return this->number_asl_hessian_nonzeros;
+   }
+
    void AMPLModel::generate_constraints() {
       for (size_t constraint_index: Range(this->number_constraints)) {
          this->constraint_lower_bounds[constraint_index] = (this->asl->i.LUrhs_ != nullptr) ? this->asl->i.LUrhs_[2*constraint_index] : -INF<double>;
@@ -407,6 +418,20 @@ namespace uno {
       }
    }
 
+   void AMPLModel::compute_lagrangian_hessian_sparsity() {
+      // compute the maximum number of nonzero elements, provided that all multipliers are non-zero
+      // int (*Sphset) (ASL*, SputInfo**, int nobj, int ow, int y, int uptri);
+      const int objective_number = -1;
+      const int upper_triangular = 1;
+      this->number_asl_hessian_nonzeros = static_cast<size_t>((*(this->asl)->p.Sphset)(this->asl, nullptr, objective_number, 1, 1, upper_triangular));
+      this->asl_hessian.reserve(this->number_asl_hessian_nonzeros);
+
+      // sparsity pattern
+      [[maybe_unused]] const fint* asl_column_start = this->asl->i.sputinfo_->hcolstarts;
+      // check that the column pointers are sorted in increasing order
+      assert(in_increasing_order(asl_column_start, this->number_variables + 1) && "AMPLModel::evaluate_lagrangian_hessian: column starts are not ordered");
+   }
+
    void AMPLModel::determine_bounds_types(const std::vector<double>& lower_bounds, const std::vector<double>& upper_bounds, std::vector<BoundType>& status) {
       assert(lower_bounds.size() == status.size());
       assert(upper_bounds.size() == status.size());
diff --git a/bindings/AMPL/AMPLModel.hpp b/bindings/AMPL/AMPLModel.hpp
index f93fbd34..47d656df 100644
--- a/bindings/AMPL/AMPLModel.hpp
+++ b/bindings/AMPL/AMPLModel.hpp
@@ -38,6 +38,8 @@ namespace uno {
       void evaluate_constraint_jacobian(const Vector<double>& x, RectangularMatrix<double>& constraint_jacobian) const override;
       void evaluate_lagrangian_hessian(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
             SymmetricMatrix<size_t, double>& hessian) const override;
+      void compute_hessian_vector_product(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
+            Vector<double>& result) const override;
 
       [[nodiscard]] double variable_lower_bound(size_t variable_index) const override;
       [[nodiscard]] double variable_upper_bound(size_t variable_index) const override;
diff --git a/uno/model/BoundRelaxedModel.hpp b/uno/model/BoundRelaxedModel.hpp
index cea67123..d5861c96 100644
--- a/uno/model/BoundRelaxedModel.hpp
+++ b/uno/model/BoundRelaxedModel.hpp
@@ -32,6 +32,10 @@ namespace uno {
             SymmetricMatrix<size_t, double>& hessian) const override {
          this->model->evaluate_lagrangian_hessian(x, objective_multiplier, multipliers, hessian);
       }
+      void compute_hessian_vector_product(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
+            Vector<double>& result) const override {
+         this->model->compute_hessian_vector_product(x, objective_multiplier, multipliers, result);
+      }
 
       // only these two functions are redefined
       [[nodiscard]] double variable_lower_bound(size_t variable_index) const override;
diff --git a/uno/model/FixedBoundsConstraintsModel.cpp b/uno/model/FixedBoundsConstraintsModel.cpp
index f77c0c6b..6a07d08b 100644
--- a/uno/model/FixedBoundsConstraintsModel.cpp
+++ b/uno/model/FixedBoundsConstraintsModel.cpp
@@ -74,6 +74,11 @@ namespace uno {
       this->model->evaluate_lagrangian_hessian(x, objective_multiplier, multipliers, hessian);
    }
 
+   void FixedBoundsConstraintsModel::compute_hessian_vector_product(const Vector<double>& x, double objective_multiplier,
+         const Vector<double>& multipliers, Vector<double>& result) const {
+      this->model->compute_hessian_vector_product(x, objective_multiplier, multipliers, result);
+   }
+
    double FixedBoundsConstraintsModel::variable_lower_bound(size_t variable_index) const {
       if (this->model->variable_lower_bound(variable_index) == this->model->variable_upper_bound(variable_index)) {
       // remove bounds of fixed variables
diff --git a/uno/model/FixedBoundsConstraintsModel.hpp b/uno/model/FixedBoundsConstraintsModel.hpp
index 3f1d7943..ed0a42e6 100644
--- a/uno/model/FixedBoundsConstraintsModel.hpp
+++ b/uno/model/FixedBoundsConstraintsModel.hpp
@@ -29,6 +29,8 @@ namespace uno {
       void evaluate_constraint_jacobian(const Vector<double>& x, RectangularMatrix<double>& constraint_jacobian) const override;
       void evaluate_lagrangian_hessian(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
             SymmetricMatrix<size_t, double>& hessian) const override;
+      void compute_hessian_vector_product(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
+            Vector<double>& result) const override;
 
       [[nodiscard]] double variable_lower_bound(size_t variable_index) const override;
       [[nodiscard]] double variable_upper_bound(size_t variable_index) const override;
diff --git a/uno/model/HomogeneousEqualityConstrainedModel.cpp b/uno/model/HomogeneousEqualityConstrainedModel.cpp
index be7d178f..76722798 100644
--- a/uno/model/HomogeneousEqualityConstrainedModel.cpp
+++ b/uno/model/HomogeneousEqualityConstrainedModel.cpp
@@ -101,6 +101,11 @@ namespace uno {
       }
    }
 
+   void HomogeneousEqualityConstrainedModel::compute_hessian_vector_product(const Vector<double>& x, double objective_multiplier,
+         const Vector<double>& multipliers, Vector<double>& result) const {
+      this->model->compute_hessian_vector_product(x, objective_multiplier, multipliers, result);
+   }
+
    double HomogeneousEqualityConstrainedModel::variable_lower_bound(size_t variable_index) const {
       if (variable_index < this->model->number_variables) { // original variable
          return this->model->variable_lower_bound(variable_index);
diff --git a/uno/model/HomogeneousEqualityConstrainedModel.hpp b/uno/model/HomogeneousEqualityConstrainedModel.hpp
index 2d5c9ea1..ca1ebb90 100644
--- a/uno/model/HomogeneousEqualityConstrainedModel.hpp
+++ b/uno/model/HomogeneousEqualityConstrainedModel.hpp
@@ -26,6 +26,8 @@ namespace uno {
       void evaluate_constraint_jacobian(const Vector<double>& x, RectangularMatrix<double>& constraint_jacobian) const override;
       void evaluate_lagrangian_hessian(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
             SymmetricMatrix<size_t, double>& hessian) const override;
+      void compute_hessian_vector_product(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
+            Vector<double>& result) const override;
 
       [[nodiscard]] double variable_lower_bound(size_t variable_index) const override;
       [[nodiscard]] double variable_upper_bound(size_t variable_index) const override;
diff --git a/uno/model/Model.hpp b/uno/model/Model.hpp
index 389c53ea..b6372eb0 100644
--- a/uno/model/Model.hpp
+++ b/uno/model/Model.hpp
@@ -51,6 +51,8 @@ namespace uno {
       virtual void evaluate_constraint_jacobian(const Vector<double>& x, RectangularMatrix<double>& constraint_jacobian) const = 0;
       virtual void evaluate_lagrangian_hessian(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
             SymmetricMatrix<size_t, double>& hessian) const = 0;
+      virtual void compute_hessian_vector_product(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
+            Vector<double>& result) const = 0;
 
       // purely virtual functions
       [[nodiscard]] virtual double variable_lower_bound(size_t variable_index) const = 0;
diff --git a/uno/model/ScaledModel.cpp b/uno/model/ScaledModel.cpp
index 22b2880b..81d52d99 100644
--- a/uno/model/ScaledModel.cpp
+++ b/uno/model/ScaledModel.cpp
@@ -11,7 +11,8 @@ namespace uno {
          Model(original_model->name + " -> scaled", original_model->number_variables, original_model->number_constraints,
                original_model->objective_sign),
          model(std::move(original_model)),
-         scaling(this->model->number_constraints, options.get_double("function_scaling_threshold")) {
+         scaling(this->model->number_constraints, options.get_double("function_scaling_threshold")),
+         scaled_multipliers(this->number_constraints) {
       if (options.get_bool("scale_functions")) {
          // evaluate the gradients at the current point
          initial_iterate.evaluate_objective_gradient(*this->model);
@@ -66,12 +67,20 @@ namespace uno {
          SymmetricMatrix<size_t, double>& hessian) const {
       // scale the objective and constraint multipliers
       const double scaled_objective_multiplier = objective_multiplier*this->scaling.get_objective_scaling();
-      // TODO preallocate this vector
-      static Vector<double> scaled_multipliers(this->number_constraints);
       for (size_t constraint_index: Range(this->number_constraints)) {
-         scaled_multipliers[constraint_index] = this->scaling.get_constraint_scaling(constraint_index) * multipliers[constraint_index];
+         this->scaled_multipliers[constraint_index] = this->scaling.get_constraint_scaling(constraint_index) * multipliers[constraint_index];
       }
-      this->model->evaluate_lagrangian_hessian(x, scaled_objective_multiplier, scaled_multipliers, hessian);
+      this->model->evaluate_lagrangian_hessian(x, scaled_objective_multiplier, this->scaled_multipliers, hessian);
+   }
+
+   void ScaledModel::compute_hessian_vector_product(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
+         Vector<double>& result) const {
+      // scale the objective and constraint multipliers
+      const double scaled_objective_multiplier = objective_multiplier*this->scaling.get_objective_scaling();
+      for (size_t constraint_index: Range(this->number_constraints)) {
+         this->scaled_multipliers[constraint_index] = this->scaling.get_constraint_scaling(constraint_index) * multipliers[constraint_index];
+      }
+      this->model->compute_hessian_vector_product(x, scaled_objective_multiplier, this->scaled_multipliers, result);
    }
 
    double ScaledModel::variable_lower_bound(size_t variable_index) const {
diff --git a/uno/model/ScaledModel.hpp b/uno/model/ScaledModel.hpp
index 582c18a3..4648fed6 100644
--- a/uno/model/ScaledModel.hpp
+++ b/uno/model/ScaledModel.hpp
@@ -6,6 +6,7 @@
 
 #include <memory>
 #include "Model.hpp"
+#include "linear_algebra/Vector.hpp"
 #include "preprocessing/Scaling.hpp"
 
 namespace uno {
@@ -23,6 +24,8 @@ namespace uno {
       void evaluate_constraint_jacobian(const Vector<double>& x, RectangularMatrix<double>& constraint_jacobian) const override;
       void evaluate_lagrangian_hessian(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
             SymmetricMatrix<size_t, double>& hessian) const override;
+      void compute_hessian_vector_product(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
+            Vector<double>& result) const override;
 
       [[nodiscard]] double variable_lower_bound(size_t variable_index) const override;
       [[nodiscard]] double variable_upper_bound(size_t variable_index) const override;
@@ -53,6 +56,7 @@ namespace uno {
    private:
       const std::unique_ptr<Model> model{};
       Scaling scaling;
+      mutable Vector<double> scaled_multipliers{};
    };
 } // namespace
 
diff --git a/uno/optimization/EvaluationErrors.hpp b/uno/optimization/EvaluationErrors.hpp
index e83189cb..ee91d4bf 100644
--- a/uno/optimization/EvaluationErrors.hpp
+++ b/uno/optimization/EvaluationErrors.hpp
@@ -9,15 +9,21 @@ namespace uno {
       [[nodiscard]] const char* what() const noexcept override = 0;
    };
 
+   struct FunctionEvaluationError : EvaluationError {
+      [[nodiscard]] const char* what() const noexcept override {
+         return "A numerical error was encountered while evaluating a function\n";
+      }
+   };
+
    struct GradientEvaluationError : EvaluationError {
       [[nodiscard]] const char* what() const noexcept override {
          return "A numerical error was encountered while evaluating a gradient\n";
       }
    };
 
-   struct FunctionEvaluationError : EvaluationError {
+   struct HessianEvaluationError : EvaluationError {
       [[nodiscard]] const char* what() const noexcept override {
-         return "A numerical error was encountered while evaluating a function\n";
+         return "A numerical error was encountered while evaluating the Lagrangian Hessian\n";
       }
    };
 } // namespace
diff --git a/uno/reformulation/OptimalityProblem.cpp b/uno/reformulation/OptimalityProblem.cpp
index 29aa94bb..ceb98aa2 100644
--- a/uno/reformulation/OptimalityProblem.cpp
+++ b/uno/reformulation/OptimalityProblem.cpp
@@ -32,6 +32,10 @@ namespace uno {
       this->model.evaluate_lagrangian_hessian(x, this->get_objective_multiplier(), multipliers, hessian);
    }
 
+   void OptimalityProblem::compute_hessian_vector_product(const Vector<double>& x, const Vector<double>& multipliers, Vector<double>& result) const {
+      this->model.compute_hessian_vector_product(x, this->get_objective_multiplier(), multipliers, result);
+   }
+
    // Lagrangian gradient split in two parts: objective contribution and constraints' contribution
    void OptimalityProblem::evaluate_lagrangian_gradient(LagrangianGradient<double>& lagrangian_gradient, Iterate& iterate,
          const Multipliers& multipliers) const {
diff --git a/uno/reformulation/OptimalityProblem.hpp b/uno/reformulation/OptimalityProblem.hpp
index 9f488dae..600a34c2 100644
--- a/uno/reformulation/OptimalityProblem.hpp
+++ b/uno/reformulation/OptimalityProblem.hpp
@@ -16,6 +16,7 @@ namespace uno {
       void evaluate_constraints(Iterate& iterate, std::vector<double>& constraints) const override;
       void evaluate_constraint_jacobian(Iterate& iterate, RectangularMatrix<double>& constraint_jacobian) const override;
       void evaluate_lagrangian_hessian(const Vector<double>& x, const Vector<double>& multipliers, SymmetricMatrix<size_t, double>& hessian) const override;
+      void compute_hessian_vector_product(const Vector<double>& x, const Vector<double>& multipliers, Vector<double>& result) const override;
 
       [[nodiscard]] double variable_lower_bound(size_t variable_index) const override { return this->model.variable_lower_bound(variable_index); }
       [[nodiscard]] double variable_upper_bound(size_t variable_index) const override { return this->model.variable_upper_bound(variable_index); }
diff --git a/uno/reformulation/OptimizationProblem.hpp b/uno/reformulation/OptimizationProblem.hpp
index cb01c884..f5f8dee7 100644
--- a/uno/reformulation/OptimizationProblem.hpp
+++ b/uno/reformulation/OptimizationProblem.hpp
@@ -42,6 +42,7 @@ namespace uno {
       virtual void evaluate_constraints(Iterate& iterate, std::vector<double>& constraints) const = 0;
       virtual void evaluate_constraint_jacobian(Iterate& iterate, RectangularMatrix<double>& constraint_jacobian) const = 0;
       virtual void evaluate_lagrangian_hessian(const Vector<double>& x, const Vector<double>& multipliers, SymmetricMatrix<size_t, double>& hessian) const = 0;
+      virtual void compute_hessian_vector_product(const Vector<double>& x, const Vector<double>& multipliers, Vector<double>& result) const = 0;
 
       [[nodiscard]] size_t get_number_original_variables() const;
       [[nodiscard]] virtual double variable_lower_bound(size_t variable_index) const = 0;
diff --git a/uno/reformulation/l1RelaxedProblem.cpp b/uno/reformulation/l1RelaxedProblem.cpp
index d8085d0b..da9ae3cc 100644
--- a/uno/reformulation/l1RelaxedProblem.cpp
+++ b/uno/reformulation/l1RelaxedProblem.cpp
@@ -113,6 +113,19 @@ namespace uno {
       }
    }
 
+   void l1RelaxedProblem::compute_hessian_vector_product(const Vector<double>& x, const Vector<double>& multipliers, Vector<double>& result) const {
+      this->model.compute_hessian_vector_product(x, this->objective_multiplier, multipliers, result);
+
+      // proximal contribution
+      if (this->proximal_center != nullptr && this->proximal_coefficient != 0.) {
+         for (size_t variable_index: Range(this->model.number_variables)) {
+            const double scaling = std::min(1., 1./std::abs(this->proximal_center[variable_index]));
+            const double proximal_term = this->proximal_coefficient * scaling * scaling;
+            result[variable_index] += proximal_term * x[variable_index] * x[variable_index];
+         }
+      }
+   }
+
    // Lagrangian gradient split in two parts: objective contribution and constraints' contribution
    void l1RelaxedProblem::evaluate_lagrangian_gradient(LagrangianGradient<double>& lagrangian_gradient, Iterate& iterate,
          const Multipliers& multipliers) const {
diff --git a/uno/reformulation/l1RelaxedProblem.hpp b/uno/reformulation/l1RelaxedProblem.hpp
index e5ab02ed..d9d927c3 100644
--- a/uno/reformulation/l1RelaxedProblem.hpp
+++ b/uno/reformulation/l1RelaxedProblem.hpp
@@ -19,6 +19,7 @@ namespace uno {
       void evaluate_constraints(Iterate& iterate, std::vector<double>& constraints) const override;
       void evaluate_constraint_jacobian(Iterate& iterate, RectangularMatrix<double>& constraint_jacobian) const override;
       void evaluate_lagrangian_hessian(const Vector<double>& x, const Vector<double>& multipliers, SymmetricMatrix<size_t, double>& hessian) const override;
+      void compute_hessian_vector_product(const Vector<double>& x, const Vector<double>& multipliers, Vector<double>& result) const override;
 
       [[nodiscard]] double variable_lower_bound(size_t variable_index) const override;
       [[nodiscard]] double variable_upper_bound(size_t variable_index) const override;
diff --git a/uno/solvers/BQPD/BQPDSolver.cpp b/uno/solvers/BQPD/BQPDSolver.cpp
index ee3aceac..91c21c52 100644
--- a/uno/solvers/BQPD/BQPDSolver.cpp
+++ b/uno/solvers/BQPD/BQPDSolver.cpp
@@ -317,7 +317,9 @@ namespace uno {
 } // namespace
 
 void hessian_vector_product(int *n, const double x[], const double ws[], const int lws[], double v[]) {
-   for (size_t i = 0; i < *n; i++) {
+   assert(n != nullptr && "BQPDSolver::hessian_vector_product: the dimension n passed by pointer is NULL");
+
+   for (size_t i = 0; i < static_cast<size_t>(*n); i++) {
       v[i] = 0.;
    }