diff --git a/inflation/lp/InflationLP.py b/inflation/lp/InflationLP.py
index 5461d8de..f5bbbdd1 100644
--- a/inflation/lp/InflationLP.py
+++ b/inflation/lp/InflationLP.py
@@ -27,7 +27,7 @@
 from ..sdp.fast_npa import nb_is_knowable as is_knowable
 from .monomial_classes import InternalAtomicMonomial, CompoundMoment
 from ..sdp.quantum_tools import flatten_symbolic_powers
-from .lp_utils import solveLP
+from .lp_utils import solveLP, solve_Gurobi
 from functools import reduce
 from ..utils import clean_coefficients, eprint, partsextractor, \
     expand_sparse_vec
@@ -490,6 +490,7 @@ def set_extra_inequalities(self,
                                    for ineq in extra_inequalities]
 
     def solve(self,
+              interpreter: str = "solveLP_sparse",
               relax_known_vars: bool = False,
               relax_inequalities: bool = False,
               solve_dual: bool = True,
@@ -503,6 +504,8 @@ def solve(self,
 
         Parameters
         ----------
+        interpreter : str, optional
+            The solver to be called. By default, ``"solveLP_sparse"``.
         relax_known_vars : bool, optional
             Do feasibility as optimization where each known value equality
             becomes two relaxed inequality constraints. E.g., P(A) = 0.7
@@ -531,6 +534,10 @@ def solve(self,
                  "relax_known_vars=False, relax_inequalities=False, and "
                  "optimizing the objective...")
             relax_known_vars = relax_inequalities = False
+        if interpreter == "solve_Gurobi" and solve_dual is True:
+            warn("Gurobi does not support solving the dual problem. Use "
+                 "solveLP_sparse if you want to solve the dual problem. "
+                 "Setting solve_dual=False...")
         if verbose is None:
             real_verbose = self.verbose
         else:
@@ -539,8 +546,14 @@ def solve(self,
             real_default_non_negative = self.default_non_negative
         else:
             real_default_non_negative = default_non_negative
-
-        args = self._prepare_solver_matrices()
+            
+        qcp_interpreters = {"solve_Gurobi"}
+        sparse_interpreters = {"solveLP", "solveLP_sparse", "solve_Gurobi"}
+        is_qcp = (interpreter in qcp_interpreters)
+        if interpreter in sparse_interpreters:
+            args = self._prepare_solver_matrices(include_nonlinear=is_qcp)
+        else:
+            args = self._prepare_solver_arguments(include_nonlinear=is_qcp)
         
         # Still allow for 'feas_as_optim' as an argument
         if 'feas_as_optim' in solver_arguments:
@@ -559,7 +572,10 @@ def solve(self,
                      "solverparameters": solverparameters,
                      "solve_dual": solve_dual})
 
-        self.solution_object = solveLP(**args)
+        if is_qcp:
+            self.solution_object = solve_Gurobi(**args)
+        else:
+            self.solution_object = solveLP(**args)
         self.success = self.solution_object["success"]
         self.status = self.solution_object["status"]
         if self.success:
@@ -1672,7 +1688,8 @@ def semiknown_by_name(self) -> Dict:
                                          in self.semiknown_moments.items()}
 
     def _prepare_solver_matrices(self,
-                                 separate_bounds: bool = True) -> dict:
+                                 separate_bounds: bool = True,
+                                 include_nonlinear: bool = False) -> dict:
         """Convert arguments from dictionaries to sparse coo_matrix form to
         pass to the solver.
 
@@ -1682,6 +1699,9 @@ def _prepare_solver_matrices(self,
             Whether to have variable bounds as a separate item in
             ``solverargs`` (True) or absorb them into the inequalities (False).
             By default, ``True``.
+        include_nonlinear : bool, optional
+            Whether to include quadratic factorization conditions in the
+            solver arguments, by default, ``False``.
 
         Returns
         -------
@@ -1710,6 +1730,9 @@ def _prepare_solver_matrices(self,
                       "equalities": internal_equalities,
                       "inequalities": self.sparse_inequalities,
                       "variables": self.monomial_names}
+        if include_nonlinear:
+            solverargs["factorization_conditions"] = \
+                self.sparse_quadratic_factorization_conditions
         if separate_bounds:
             solverargs["lower_bounds"] = self.sparse_lowerbounds
             solverargs["upper_bounds"] = self.sparse_upperbounds
@@ -1723,7 +1746,8 @@ def _prepare_solver_matrices(self,
         return solverargs
 
     def _prepare_solver_arguments(self, 
-                                  separate_bounds: bool = True) -> dict:
+                                  separate_bounds: bool = True,
+                                  include_nonlinear: bool = False) -> dict:
         """Prepare arguments to pass to the solver.
 
         The solver takes as input the following arguments, which are all
@@ -1744,6 +1768,9 @@ def _prepare_solver_arguments(self,
             Whether to have variable bounds as a separate item in
             ``solverargs`` (True) or absorb them into the inequalities (False).
             By default, ``True``.
+        include_nonlinear : bool, optional
+            Whether to include quadratic factorization conditions in the
+            solver arguments, by default, ``False``.
 
         Returns
         -------
@@ -1770,6 +1797,9 @@ def _prepare_solver_arguments(self,
                       "equalities": self.moment_equalities_by_name,
                       "inequalities": self.moment_inequalities_by_name
                       }
+        if include_nonlinear:
+            solverargs["factorization_conditions"] = \
+                self.quadratic_factorization_conditions_by_name
         # Add the constant 1 in case of unnormalized problems removed it
         solverargs["known_vars"][self.constant_term_name] = 1.
         if separate_bounds:
diff --git a/inflation/lp/lp_utils.py b/inflation/lp/lp_utils.py
index 59e867fe..0f5ca279 100644
--- a/inflation/lp/lp_utils.py
+++ b/inflation/lp/lp_utils.py
@@ -13,6 +13,10 @@
 from time import perf_counter
 from gc import collect
 from inflation.utils import partsextractor, expand_sparse_vec, vstack
+try:
+    import gurobipy as gp
+except ModuleNotFoundError:
+    pass
 
 
 def drop_zero_rows(coo_mat: coo_matrix):
@@ -532,6 +536,142 @@ def solveLP_sparse(objective: coo_matrix = blank_coo_matrix,
             }
 
 
+def solve_Gurobi(objective: coo_matrix = blank_coo_matrix,
+                 known_vars: coo_matrix = blank_coo_matrix,
+                 inequalities: coo_matrix = blank_coo_matrix,
+                 equalities: coo_matrix = blank_coo_matrix,
+                 lower_bounds: coo_matrix = blank_coo_matrix,
+                 upper_bounds: coo_matrix = blank_coo_matrix,
+                 factorization_conditions: dict = None,
+                 default_non_negative: bool = True,
+                 verbose: int = 0,
+                 **kwargs
+                 ) -> Dict:
+    """Internal function to solve an LP with the Gurobi Optimizer API using
+    sparse matrices. Columns of each matrix correspond to a fixed order of
+    variables in the LP.
+    
+    Parameters
+    ----------
+    objective : coo_matrix, optional
+        Objective function with coefficients as matrix entries.
+    known_vars : coo_matrix, optional
+        Known values of the monomials with values as matrix entries.
+    inequalities : coo_matrix, optional
+        Inequality constraints in matrix form.
+    equalities : coo_matrix, optional
+        Equality constraints in matrix form.
+    lower_bounds : coo_matrix, optional
+        Lower bounds of variables with bounds as matrix entries.
+    upper_bounds : coo_matrix, optional
+        Upper bounds of variables with bounds as matrix entries.
+    factorization_conditions : dict, optional
+        Allows one to specify that certain variables are equal to the product of other variables
+    default_non_negative : bool, optional
+        Whether to set default primal variables as non-negative. By default,
+        ``True``.
+    verbose : int, optional
+        Verbosity. Higher means more messages. By default, 0.
+        
+    Returns
+    -------
+    dict
+        Primal objective value, problem status, success status, x values
+    """
+    if verbose > 1:
+        t0 = perf_counter()
+        t_total = perf_counter()
+        print("Starting pre-processing for the LP solver...")
+
+    env = gp.Env(empty=True)
+    if verbose == 0:
+        env.setParam('OutputFlag', 0)  # Supress output from solver
+    env.start()
+
+    # Create new model
+    m = gp.Model(env=env)
+
+    (nof_primal_inequalities, nof_primal_variables) = inequalities.shape
+    nof_primal_equalities = equalities.shape[0]
+    kv_matrix = expand_sparse_vec(known_vars)
+
+    # Create variables and set variable bounds
+    x = m.addMVar(shape=nof_primal_variables)
+    m.update()
+    if not default_non_negative:
+        x.lb = -gp.GRB.INFINITY
+    if lower_bounds.nnz > 0:
+        new_lb = x.lb
+        new_lb[lower_bounds.col] = lower_bounds.data
+        x.lb = new_lb
+    if upper_bounds.nnz > 0:
+        new_ub = x.ub
+        new_ub[upper_bounds.col] = upper_bounds.data
+        x.ub = new_ub
+
+    # Set objective
+    objective_vector = objective.toarray().ravel()
+    m.setObjective(objective_vector @ x, gp.GRB.MAXIMIZE)
+
+    # Add inequality constraint matrix
+    if inequalities.shape[0]:
+        m.addConstr(inequalities @ x >= 0, name="ineq")
+
+    # Add equality constraint matrix
+    if equalities.shape[0]:
+        m.addConstr(equalities @ x == 0, name="eq")
+
+    # Add known value constraint matrix
+    rhs_kv = known_vars.data
+    if len(rhs_kv):
+        m.addConstr(kv_matrix @ x == rhs_kv, name="kv")
+
+    # Add quadratic constraints
+    if factorization_conditions:
+        m.setParam("NonConvex", 2)
+        m.addConstrs((x[mon] == x[factors[0]] * x[factors[1]]
+                     for mon, factors in factorization_conditions.items()),
+                     name="fac")
+
+    collect()
+    if verbose > 1:
+        print("Pre-processing took",
+              format(perf_counter() - t0, ".4f"), "seconds.\n")
+        t0 = perf_counter()
+
+    # Solve the problem
+    if verbose > 0:
+        print("\nSolving the problem...\n")
+    m.optimize()
+    if verbose > 1:
+        print("Solving took", format(perf_counter() - t0, ".4f"),
+              "seconds.")
+
+    sc = gp.StatusConstClass
+    solution_statuses = {sc.__dict__[k]: k for k in sc.__dict__.keys()
+                         if 'A' <= k[0] <= 'Z'}
+    status_str = solution_statuses[m.Status]
+    if status_str == "OPTIMAL":
+        success = True
+        primal = m.ObjVal
+        x_values = x.X
+    else:
+        success = False
+        primal = None
+        x_values = None
+
+    if verbose > 1:
+        print("\nTotal execution time:",
+              format(perf_counter() - t_total, ".4f"), "seconds.")
+
+    return {
+        "primal_value": primal,
+        "status": status_str,
+        "success": success,
+        "x": x_values,
+    }
+
+
 def streamprinter(text: str) -> None:
     """A stream printer to get output from Mosek.