diff --git a/CHANGELOG.md b/CHANGELOG.md
index cc7dd4b59..f771fb4c8 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -19,7 +19,6 @@
 - changed default value of enablepricing flag to True
 - Speed up MatrixExpr.add.reduce via quicksum
 - Speed up np.ndarray(..., dtype=np.float64) @ MatrixExpr
-- Minimum numpy version increased from 1.16.0 to 1.19.0
 - MatrixExpr and MatrixExprCons use `__array_ufunc__` protocol to control all numpy.ufunc inputs and outputs
 ### Removed
 
diff --git a/pyproject.toml b/pyproject.toml
index 9ac13da7f..d990f8d50 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -8,7 +8,7 @@ description = "Python interface and modeling environment for SCIP"
 authors = [
     {name = "Zuse Institute Berlin", email = "scip@zib.de"},
 ]
-dependencies = ['numpy >=1.19.0']
+dependencies = ['numpy >=1.16.0']
 requires-python = ">=3.8"
 readme = "README.md"
 license = {text = "MIT License"}
diff --git a/src/pyscipopt/matrix.pxi b/src/pyscipopt/matrix.pxi
index 84fb09b77..0043b7638 100644
--- a/src/pyscipopt/matrix.pxi
+++ b/src/pyscipopt/matrix.pxi
@@ -1,7 +1,6 @@
 import operator
 from typing import Literal, Optional, Tuple, Union
 import numpy as np
-from numpy.typing import NDArray
 try:
     # NumPy 2.x location
     from numpy.lib.array_utils import normalize_axis_tuple
@@ -73,7 +72,6 @@ class MatrixExpr(np.ndarray):
             res = super().__array_ufunc__(ufunc, method, *args, **kwargs)
         return res.view(MatrixExpr) if isinstance(res, np.ndarray) else res
 
-
     def _evaluate(self, Solution sol) -> NDArray[np.float64]:
         return _vec_evaluate(self, sol).view(np.ndarray)
 
@@ -100,7 +98,7 @@ _vec_eq = np.frompyfunc(operator.eq, 2, 1)
 _vec_evaluate = np.frompyfunc(lambda expr, sol: expr._evaluate(sol), 2, 1)
 
 
-cdef inline _ensure_array(arg, bool convert_scalar = True):
+cdef inline _ensure_array(arg, bint convert_scalar = True):
     if isinstance(arg, np.ndarray):
         return arg.view(np.ndarray)
     elif isinstance(arg, (list, tuple)):
@@ -127,14 +125,14 @@ def _core_dot(cnp.ndarray a, cnp.ndarray b) -> Union[Expr, np.ndarray]:
         If both `a` and `b` are 1-D arrays, return an `Expr`, otherwise return a
         `np.ndarray` of type `object` and containing `Expr` objects.
     """
-    cdef bool a_is_1d = a.ndim == 1
-    cdef bool b_is_1d = b.ndim == 1
+    cdef bint a_is_1d = a.ndim == 1
+    cdef bint b_is_1d = b.ndim == 1
     cdef cnp.ndarray a_nd = a[..., np.newaxis, :] if a_is_1d else a
     cdef cnp.ndarray b_nd = b[..., :, np.newaxis] if b_is_1d else b
-    cdef bool a_is_num = a_nd.dtype.kind in "fiub"
+    cdef bint a_is_num = a_nd.dtype.kind in "fiub"
 
     if a_is_num ^ (b_nd.dtype.kind in "fiub"):
-        res = _core_dot_2d(a_nd, b_nd) if a_is_num else _core_dot_2d(b_nd.T, a_nd.T).T
+        res = _core_dot_nd(a_nd, b_nd) if a_is_num else _core_dot_nd(b_nd.T, a_nd.T).T
         if a_is_1d and b_is_1d:
             return res.item()
         if a_is_1d:
@@ -145,7 +143,6 @@ def _core_dot(cnp.ndarray a, cnp.ndarray b) -> Union[Expr, np.ndarray]:
     return NotImplemented
 
 
-@np.vectorize(otypes=[object], signature="(m,n),(n,p)->(m,p)")
 def _core_dot_2d(cnp.ndarray a, cnp.ndarray x) -> np.ndarray:
     """
     Perform matrix multiplication between a 2-Demension constant array and a 2-Demension
@@ -183,6 +180,13 @@ def _core_dot_2d(cnp.ndarray a, cnp.ndarray x) -> np.ndarray:
     return res
 
 
+_core_dot_nd = np.vectorize(
+    _core_dot_2d,
+    otypes=[object],
+    signature="(m,n),(n,p)->(m,p)",
+)
+
+
 def _core_sum(
     cnp.ndarray a,
     axis: Optional[Union[int, Tuple[int, ...]]] = None,
diff --git a/src/pyscipopt/scip.pxi b/src/pyscipopt/scip.pxi
index 87757fb3c..f4f0ecd35 100644
--- a/src/pyscipopt/scip.pxi
+++ b/src/pyscipopt/scip.pxi
@@ -21,7 +21,6 @@ from dataclasses import dataclass
 from typing import Union
 
 import numpy as np
-from numpy.typing import NDArray
 
 include "expr.pxi"
 include "lp.pxi"
@@ -10735,7 +10734,7 @@ cdef class Model:
         self,
         Solution sol,
         expr: Union[Expr, GenExpr],
-    ) -> Union[float, NDArray[np.float64]]:
+    ) -> Union[float, np.ndarray]:
         """
         Retrieve value of given variable or expression in the given solution or in
         the LP/pseudo solution if sol == None