diff --git a/linear_algebra/src/conjugate_gradient.py b/linear_algebra/src/conjugate_gradient.py
index 4c0b58deb978..45da35813978 100644
--- a/linear_algebra/src/conjugate_gradient.py
+++ b/linear_algebra/src/conjugate_gradient.py
@@ -61,7 +61,8 @@ def _create_spd_matrix(dimension: int) -> Any:
     >>> _is_matrix_spd(spd_matrix)
     True
     """
-    random_matrix = np.random.randn(dimension, dimension)
+    rng = np.random.default_rng()
+    random_matrix = rng.normal(size=(dimension, dimension))
     spd_matrix = np.dot(random_matrix, random_matrix.T)
     assert _is_matrix_spd(spd_matrix)
     return spd_matrix
@@ -157,7 +158,8 @@ def test_conjugate_gradient() -> None:
     # Create linear system with SPD matrix and known solution x_true.
     dimension = 3
     spd_matrix = _create_spd_matrix(dimension)
-    x_true = np.random.randn(dimension, 1)
+    rng = np.random.default_rng()
+    x_true = rng.normal(size=(dimension, 1))
     b = np.dot(spd_matrix, x_true)
 
     # Numpy solution.
diff --git a/machine_learning/decision_tree.py b/machine_learning/decision_tree.py
index 7f129919a3ce..e48905eeac6a 100644
--- a/machine_learning/decision_tree.py
+++ b/machine_learning/decision_tree.py
@@ -187,7 +187,8 @@ def main():
     tree = DecisionTree(depth=10, min_leaf_size=10)
     tree.train(x, y)
 
-    test_cases = (np.random.rand(10) * 2) - 1
+    rng = np.random.default_rng()
+    test_cases = (rng.random(10) * 2) - 1
     predictions = np.array([tree.predict(x) for x in test_cases])
     avg_error = np.mean((predictions - test_cases) ** 2)
 
diff --git a/machine_learning/k_means_clust.py b/machine_learning/k_means_clust.py
index 9f6646944458..a926362fc18b 100644
--- a/machine_learning/k_means_clust.py
+++ b/machine_learning/k_means_clust.py
@@ -55,12 +55,12 @@
 
 def get_initial_centroids(data, k, seed=None):
     """Randomly choose k data points as initial centroids"""
-    if seed is not None:  # useful for obtaining consistent results
-        np.random.seed(seed)
+    # useful for obtaining consistent results
+    rng = np.random.default_rng(seed)
     n = data.shape[0]  # number of data points
 
     # Pick K indices from range [0, N).
-    rand_indices = np.random.randint(0, n, k)
+    rand_indices = rng.integers(0, n, k)
 
     # Keep centroids as dense format, as many entries will be nonzero due to averaging.
     # As long as at least one document in a cluster contains a word,
diff --git a/machine_learning/sequential_minimum_optimization.py b/machine_learning/sequential_minimum_optimization.py
index be16baca1a4c..408d59ab5d29 100644
--- a/machine_learning/sequential_minimum_optimization.py
+++ b/machine_learning/sequential_minimum_optimization.py
@@ -289,12 +289,13 @@ def _choose_a2(self, i1):
             if cmd is None:
                 return
 
-        for i2 in np.roll(self.unbound, np.random.choice(self.length)):
+        rng = np.random.default_rng()
+        for i2 in np.roll(self.unbound, rng.choice(self.length)):
             cmd = yield i1, i2
             if cmd is None:
                 return
 
-        for i2 in np.roll(self._all_samples, np.random.choice(self.length)):
+        for i2 in np.roll(self._all_samples, rng.choice(self.length)):
             cmd = yield i1, i2
             if cmd is None:
                 return
diff --git a/neural_network/back_propagation_neural_network.py b/neural_network/back_propagation_neural_network.py
index 7e0bdbbe2857..6131a13e945e 100644
--- a/neural_network/back_propagation_neural_network.py
+++ b/neural_network/back_propagation_neural_network.py
@@ -51,8 +51,9 @@ def __init__(
         self.is_input_layer = is_input_layer
 
     def initializer(self, back_units):
-        self.weight = np.asmatrix(np.random.normal(0, 0.5, (self.units, back_units)))
-        self.bias = np.asmatrix(np.random.normal(0, 0.5, self.units)).T
+        rng = np.random.default_rng()
+        self.weight = np.asmatrix(rng.normal(0, 0.5, (self.units, back_units)))
+        self.bias = np.asmatrix(rng.normal(0, 0.5, self.units)).T
         if self.activation is None:
             self.activation = sigmoid
 
@@ -174,7 +175,8 @@ def plot_loss(self):
 
 
 def example():
-    x = np.random.randn(10, 10)
+    rng = np.random.default_rng()
+    x = rng.normal(size=(10, 10))
     y = np.asarray(
         [
             [0.8, 0.4],
diff --git a/neural_network/convolution_neural_network.py b/neural_network/convolution_neural_network.py
index 07cc456b7466..3c551924442d 100644
--- a/neural_network/convolution_neural_network.py
+++ b/neural_network/convolution_neural_network.py
@@ -41,15 +41,16 @@ def __init__(
         self.size_pooling1 = size_p1
         self.rate_weight = rate_w
         self.rate_thre = rate_t
+        rng = np.random.default_rng()
         self.w_conv1 = [
-            np.asmatrix(-1 * np.random.rand(self.conv1[0], self.conv1[0]) + 0.5)
+            np.asmatrix(-1 * rng.random((self.conv1[0], self.conv1[0])) + 0.5)
             for i in range(self.conv1[1])
         ]
-        self.wkj = np.asmatrix(-1 * np.random.rand(self.num_bp3, self.num_bp2) + 0.5)
-        self.vji = np.asmatrix(-1 * np.random.rand(self.num_bp2, self.num_bp1) + 0.5)
-        self.thre_conv1 = -2 * np.random.rand(self.conv1[1]) + 1
-        self.thre_bp2 = -2 * np.random.rand(self.num_bp2) + 1
-        self.thre_bp3 = -2 * np.random.rand(self.num_bp3) + 1
+        self.wkj = np.asmatrix(-1 * rng.random((self.num_bp3, self.num_bp2)) + 0.5)
+        self.vji = np.asmatrix(-1 * rng.random((self.num_bp2, self.num_bp1)) + 0.5)
+        self.thre_conv1 = -2 * rng.random(self.conv1[1]) + 1
+        self.thre_bp2 = -2 * rng.random(self.num_bp2) + 1
+        self.thre_bp3 = -2 * rng.random(self.num_bp3) + 1
 
     def save_model(self, save_path):
         # save model dict with pickle
diff --git a/neural_network/input_data.py b/neural_network/input_data.py
index 9d4195487dbb..d189e3f9e0d9 100644
--- a/neural_network/input_data.py
+++ b/neural_network/input_data.py
@@ -153,7 +153,7 @@ def __init__(
         """
         seed1, seed2 = random_seed.get_seed(seed)
         # If op level seed is not set, use whatever graph level seed is returned
-        np.random.seed(seed1 if seed is None else seed2)
+        self._rng = np.random.default_rng(seed1 if seed is None else seed2)
         dtype = dtypes.as_dtype(dtype).base_dtype
         if dtype not in (dtypes.uint8, dtypes.float32):
             raise TypeError("Invalid image dtype %r, expected uint8 or float32" % dtype)
@@ -211,7 +211,7 @@ def next_batch(self, batch_size, fake_data=False, shuffle=True):
         # Shuffle for the first epoch
         if self._epochs_completed == 0 and start == 0 and shuffle:
             perm0 = np.arange(self._num_examples)
-            np.random.shuffle(perm0)
+            self._rng.shuffle(perm0)
             self._images = self.images[perm0]
             self._labels = self.labels[perm0]
         # Go to the next epoch
@@ -225,7 +225,7 @@ def next_batch(self, batch_size, fake_data=False, shuffle=True):
             # Shuffle the data
             if shuffle:
                 perm = np.arange(self._num_examples)
-                np.random.shuffle(perm)
+                self._rng.shuffle(perm)
                 self._images = self.images[perm]
                 self._labels = self.labels[perm]
             # Start next epoch
diff --git a/neural_network/two_hidden_layers_neural_network.py b/neural_network/two_hidden_layers_neural_network.py
index dea7e2342d9f..d488de590cc2 100644
--- a/neural_network/two_hidden_layers_neural_network.py
+++ b/neural_network/two_hidden_layers_neural_network.py
@@ -28,19 +28,20 @@ def __init__(self, input_array: np.ndarray, output_array: np.ndarray) -> None:
         # Random initial weights are assigned.
         # self.input_array.shape[1] is used to represent number of nodes in input layer.
         # First hidden layer consists of 4 nodes.
-        self.input_layer_and_first_hidden_layer_weights = np.random.rand(
-            self.input_array.shape[1], 4
+        rng = np.random.default_rng()
+        self.input_layer_and_first_hidden_layer_weights = rng.random(
+            (self.input_array.shape[1], 4)
         )
 
         # Random initial values for the first hidden layer.
         # First hidden layer has 4 nodes.
         # Second hidden layer has 3 nodes.
-        self.first_hidden_layer_and_second_hidden_layer_weights = np.random.rand(4, 3)
+        self.first_hidden_layer_and_second_hidden_layer_weights = rng.random((4, 3))
 
         # Random initial values for the second hidden layer.
         # Second hidden layer has 3 nodes.
         # Output layer has 1 node.
-        self.second_hidden_layer_and_output_layer_weights = np.random.rand(3, 1)
+        self.second_hidden_layer_and_output_layer_weights = rng.random((3, 1))
 
         # Real output values provided.
         self.output_array = output_array
diff --git a/pyproject.toml b/pyproject.toml
index 22da7cb777b5..827cea1fcf37 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -7,7 +7,6 @@ lint.ignore = [    # `ruff rule S101` for a description of that rule
   "EXE001",   # Shebang is present but file is not executable" -- FIX ME
   "G004",     # Logging statement uses f-string
   "INP001",   # File `x/y/z.py` is part of an implicit namespace package. Add an `__init__.py`. -- FIX ME
-  "NPY002",   # Replace legacy `np.random.choice` call with `np.random.Generator` -- FIX ME
   "PGH003",   # Use specific rule codes when ignoring type issues -- FIX ME
   "PLC1901",  # `{}` can be simplified to `{}` as an empty string is falsey
   "PLW060",   # Using global for `{name}` but no assignment is done -- DO NOT FIX