Skip to content

Exposed public API for L-Marvel

learn_and_get_skeleton(ci_test, data, find_markov_boundary_matrix_fun=None)

Learn the skeleton of a causal graph with latent variables using the L-MARVEL algorithm.

Parameters:

Name Type Description Default
ci_test Callable[[int, int, List[int], ndarray], bool]

A conditional independence test function that takes in the indices of two variables and a list of variable indices as the conditioning set, and returns True if the two variables are independent given the conditioning set, and False otherwise.

 required
data_matrix ndarray

The data matrix with shape (num_samples, num_vars), where each column corresponds to a variable and each row corresponds to a sample.

 required

Returns:

Type Description
Graph

nx.Graph: A networkx graph representing the learned skeleton.
Source code in rcd/l_marvel/l_marvel.py 
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
def learn_and_get_skeleton(ci_test: Callable[[int, int, List[int], np.ndarray], bool], data,
                           find_markov_boundary_matrix_fun=None) -> nx.Graph:
    """
    Learn the skeleton of a causal graph with latent variables using the L-MARVEL algorithm.

    Args:
        ci_test (Callable[[int, int, List[int], np.ndarray], bool]):
            A conditional independence test function that takes in the indices of two variables
            and a list of variable indices as the conditioning set, and returns True if the two
            variables are independent given the conditioning set, and False otherwise.
        data_matrix (np.ndarray): The data matrix with shape (num_samples, num_vars), where each column corresponds
                                  to a variable and each row corresponds to a sample.

    Returns:
        nx.Graph: A networkx graph representing the learned skeleton.
    """

    data_mat = sanitize_data(data)
    l_marvel = _LMarvel(ci_test, find_markov_boundary_matrix_fun)
    learned_skeleton = l_marvel.learn_and_get_skeleton(data_mat)
    return learned_skeleton

L-Marvel implementation details

Implementation for the L-MARVEL algorithm for learning causal graphs with latent variables.

This class is initialized with a conditional independence test function, which determines whether two variables are independent given another set of variables, using the data provided.

The class has a learn_and_get_skeleton function that takes in a data matrix (numpy array), where each column corresponds to a variable and each row corresponds to a sample, and returns a networkx graph representing the learned skeleton.

Source code in rcd/l_marvel/l_marvel.py 
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
class _LMarvel:
    """
    Implementation for the L-MARVEL algorithm for learning causal graphs with latent variables.

    This class is initialized with a conditional independence test function, which determines whether two variables are
    independent given another set of variables, using the data provided.

    The class has a learn_and_get_skeleton function that takes in a data matrix (numpy array), where each column
    corresponds to a variable and each row corresponds to a sample, and returns a networkx graph representing the
    learned skeleton.
    """

    def __init__(self, ci_test: Callable[[int, int, List[int], np.ndarray], bool], find_markov_boundary_matrix_fun: Callable[[np.ndarray], np.ndarray] = None):
        """
        Initialize the L-MARVEL algorithm with the conditional independence test to use.

        Args:
            ci_test (Callable[[int, int, List[int], np.ndarray], bool]):
                A conditional independence test function. It takes the indices of two variables
                and a list of variable indices as the conditioning set. It returns True if the two
                variables are independent given the conditioning set, and False otherwise.
                Signature:
                    ci_test(var1: int, var2: int, cond_set: List[int], data: np.ndarray) -> bool.
            find_markov_boundary_matrix_fun (Callable[[np.ndarray], np.ndarray], optional):
                A function to find the Markov boundary matrix. It takes a numpy array of data
                and returns a 2D numpy array. The (i, j)th entry is True if the jth variable is in the
                Markov boundary of the ith variable, and False otherwise.
                Signature:
                    find_markov_boundary_matrix_fun(data: np.ndarray) -> np.ndarray.
        """
        if find_markov_boundary_matrix_fun is None:
            self.find_markov_boundary_matrix = compute_mb_gaussian
        else:
            self.find_markov_boundary_matrix = find_markov_boundary_matrix_fun

        self.num_vars = None
        self.data = None
        self.ci_test = ci_test

        # we use a flag array to keep track of which variables need to be checked for removal (i.e., we check if False)
        self.skip_rem_check_vec = None  # SkipCheck_VEC in the paper

        # we use a set to keep track of which Y and Z pairs have been checked for a given X (see IsRemovable in the paper)
        self.skip_rem_check_set = None  # SkipCheck_MAT in the paper

        # we use a flag array to keep track of which variables' neighbors need to be learned (i.e., we learn if False)
        self.neighbor_learned_arr = None
        self.var_idx_set = None
        self.markov_boundary_matrix = None
        self.learned_skeleton = None

    def learn_and_get_skeleton(self, data: np.ndarray) -> nx.Graph:
        """
        Run the L-MARVEL algorithm on the data to learn and return the learned skeleton graph.

        Args:
            data (np.ndarray): The data matrix with shape (num_samples, num_vars).

        Returns:
            nx.Graph: A networkx graph representing the learned skeleton.
        """
        # Initialize algorithm state
        self.num_vars = data.shape[1]
        self.data = data

        self.skip_rem_check_vec = np.zeros(self.num_vars, dtype=bool)
        self.skip_rem_check_set = set()
        self.neighbor_learned_arr = np.zeros(self.num_vars, dtype=bool)
        self.var_idx_set = set(range(self.num_vars))
        self.markov_boundary_matrix = self.find_markov_boundary_matrix(self.data)
        self.learned_skeleton = nx.Graph()
        self.learned_skeleton.add_nodes_from(range(self.num_vars))

        data_included_ci_test = lambda x, y, z: self.ci_test(x, y, z, self.data)

        var_arr = np.arange(self.num_vars)
        var_left_bool_arr = np.ones(self.num_vars, dtype=bool)  # Indicates if variable is left

        for _ in range(self.num_vars - 1):
            # sort variables by decreasing Markov boundary size
            # only sort variables that are still left and whose removability has NOT been checked
            var_to_sort_bool_arr = var_left_bool_arr & ~self.skip_rem_check_vec
            var_to_sort_arr = var_arr[var_to_sort_bool_arr]
            sorted_var_arr = sort_vars_by_mkb_size(self.markov_boundary_matrix[var_to_sort_bool_arr], var_to_sort_arr)

            removable_var = REMOVABLE_NOT_FOUND
            for var in sorted_var_arr:
                # Check whether we need to learn the neighbors of var
                if not self.neighbor_learned_arr[var]:
                    neighbors = self.find_neighborhood(var)
                    self.neighbor_learned_arr[var] = True

                    # Add edges between the variable and its neighbors
                    for neighbor in neighbors:
                        self.learned_skeleton.add_edge(var, neighbor)
                else:
                    # If neighbors already learned, get them from the graph
                    neighbors = list(self.learned_skeleton.neighbors(var))

                    # Ensure only to include neighbors that are still left
                    neighbors = [neighbor for neighbor in neighbors if var_left_bool_arr[neighbor]]

                # Check if variable is removable
                if self.is_removable(var, neighbors):
                    removable_var = var
                    break
                else:
                    self.skip_rem_check_vec[var] = True

            if removable_var == REMOVABLE_NOT_FOUND:
                # If no removable found, pick the variable with the smallest Markov boundary from var_left_bool_arr
                var_left_arr = np.flatnonzero(var_left_bool_arr)
                mb_size_all = np.sum(self.markov_boundary_matrix[var_left_arr], axis=1)
                removable_var = var_left_arr[np.argmin(mb_size_all)]

                self.skip_rem_check_vec[:] = False
            else:
                # Remove the removable variable from the set of variables left
                var_left_bool_arr[removable_var] = False

            # Make sure to only include neighbors that are still left
            neighbors = [neighbor for neighbor in self.learned_skeleton.neighbors(removable_var) if var_left_bool_arr[neighbor]]

            # Update the Markov boundary matrix
            update_markov_boundary_matrix(
                self.markov_boundary_matrix,
                data_included_ci_test,
                removable_var,
                neighbors,
                skip_check=self.skip_rem_check_vec,
            )

        return self.learned_skeleton

    def find_neighborhood(self, var: int) -> np.ndarray:
        """
        Find the neighborhood of a variable using Lemma 27.

        Args:
            var (int): The variable whose neighborhood we want to find.

        Returns:
            np.ndarray: 1D numpy array containing the variables in the neighborhood.
        """
        var_mk_bool_arr = self.markov_boundary_matrix[var]
        var_mk_arr = np.flatnonzero(var_mk_bool_arr)
        var_mk_set = set(var_mk_arr)

        neighbor_bool_arr = np.copy(var_mk_bool_arr)

        for var_y in var_mk_arr:
            # Check if Y is already a neighbor of X
            if not self.learned_skeleton.has_edge(var, var_y):
                if not self.is_neighbor(var, var_y, var_mk_set):
                    neighbor_bool_arr[var_y] = False

        # Remove all variables that are not neighbors
        neighbors = np.flatnonzero(neighbor_bool_arr)
        return neighbors

    def is_neighbor(self, var: int, var_y: int, var_mk_set: Set[int]) -> bool:
        """
        Check if var_y is a neighbor of variable var using Lemma 27.

        Args:
            var (int): Index of the variable.
            var_y (int): The variable to check.
            var_mk_set (Set[int]): Set of the variables in the Markov boundary of var.

        Returns:
            bool: True if var_y is a neighbor, False otherwise.
        """
        var_mk_left_list = list(var_mk_set - {var_y})
        # Use Lemma 27 and check all proper subsets of Mb(X) - {Y}
        for cond_set_size in range(len(var_mk_left_list) + 1):
            for var_s in itertools.combinations(var_mk_left_list, cond_set_size):
                cond_set = list(var_s)
                if self.ci_test(var, var_y, cond_set, self.data):
                    # Y is a co-parent and thus NOT a neighbor
                    return False
        return True

    def is_removable(self, var: int, neighbors: np.ndarray) -> bool:
        """
        Check whether a variable is removable using Theorem 32.

        Args:
            var (int): Index of the variable.
            neighbors (np.ndarray): Neighbors of the variable.

        Returns:
            bool: True if the variable is removable, False otherwise.
        """
        var_mk_bool_arr = self.markov_boundary_matrix[var]
        var_mk_arr = np.flatnonzero(var_mk_bool_arr)
        var_mk_set = set(var_mk_arr)

        def cond_1(var_y, var_z):
            # there exists subset W in Mb(X) - {Y, Z}, s.t. Y ind. Z | W
            var_mk_left_list = list(var_mk_set - {var_y, var_z})
            for cond_set_size in range(len(var_mk_left_list) + 1):
                for var_s in itertools.combinations(var_mk_left_list, cond_set_size):
                    cond_set = list(var_s)
                    if self.ci_test(var_y, var_z, cond_set, self.data):
                        return True
            return False

        def cond_2(var_y, var_z):
            # for all subset W in Mb(X) - {Y, Z}, s.t. Y NOT ind. Z | W + {X}
            var_mk_left_left = list(var_mk_set - {var_y, var_z})
            for cond_set_size in range(len(var_mk_left_left) + 1):
                for var_s in itertools.combinations(var_mk_left_left, cond_set_size):
                    cond_set = list(var_s) + [var]
                    if self.ci_test(var_y, var_z, cond_set, self.data):
                        return False
            return True

        # Use Theorem 32 to check if X is removable. Loop over Y in Mb(X) and Z in Ne(X)
        for var_y in var_mk_arr:
            for var_z in neighbors:
                if var_y == var_z:
                    continue
                xyz_tuple = (var, min(var_y, var_z), max(var_y, var_z))
                if xyz_tuple in self.skip_rem_check_set:
                    continue
                if not (cond_1(var_y, var_z) or cond_2(var_y, var_z)):
                    return False
                self.skip_rem_check_set.add(xyz_tuple)
        return True

__init__(ci_test, find_markov_boundary_matrix_fun=None)

Initialize the L-MARVEL algorithm with the conditional independence test to use.

Parameters:

Name Type Description Default
ci_test Callable[[int, int, List[int], ndarray], bool]

A conditional independence test function. It takes the indices of two variables and a list of variable indices as the conditioning set. It returns True if the two variables are independent given the conditioning set, and False otherwise. Signature: ci_test(var1: int, var2: int, cond_set: List[int], data: np.ndarray) -> bool.

 required
find_markov_boundary_matrix_fun Callable[[ndarray], ndarray]

A function to find the Markov boundary matrix. It takes a numpy array of data and returns a 2D numpy array. The (i, j)th entry is True if the jth variable is in the Markov boundary of the ith variable, and False otherwise. Signature: find_markov_boundary_matrix_fun(data: np.ndarray) -> np.ndarray.

 None
Source code in rcd/l_marvel/l_marvel.py 
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
def __init__(self, ci_test: Callable[[int, int, List[int], np.ndarray], bool], find_markov_boundary_matrix_fun: Callable[[np.ndarray], np.ndarray] = None):
    """
    Initialize the L-MARVEL algorithm with the conditional independence test to use.

    Args:
        ci_test (Callable[[int, int, List[int], np.ndarray], bool]):
            A conditional independence test function. It takes the indices of two variables
            and a list of variable indices as the conditioning set. It returns True if the two
            variables are independent given the conditioning set, and False otherwise.
            Signature:
                ci_test(var1: int, var2: int, cond_set: List[int], data: np.ndarray) -> bool.
        find_markov_boundary_matrix_fun (Callable[[np.ndarray], np.ndarray], optional):
            A function to find the Markov boundary matrix. It takes a numpy array of data
            and returns a 2D numpy array. The (i, j)th entry is True if the jth variable is in the
            Markov boundary of the ith variable, and False otherwise.
            Signature:
                find_markov_boundary_matrix_fun(data: np.ndarray) -> np.ndarray.
    """
    if find_markov_boundary_matrix_fun is None:
        self.find_markov_boundary_matrix = compute_mb_gaussian
    else:
        self.find_markov_boundary_matrix = find_markov_boundary_matrix_fun

    self.num_vars = None
    self.data = None
    self.ci_test = ci_test

    # we use a flag array to keep track of which variables need to be checked for removal (i.e., we check if False)
    self.skip_rem_check_vec = None  # SkipCheck_VEC in the paper

    # we use a set to keep track of which Y and Z pairs have been checked for a given X (see IsRemovable in the paper)
    self.skip_rem_check_set = None  # SkipCheck_MAT in the paper

    # we use a flag array to keep track of which variables' neighbors need to be learned (i.e., we learn if False)
    self.neighbor_learned_arr = None
    self.var_idx_set = None
    self.markov_boundary_matrix = None
    self.learned_skeleton = None

find_neighborhood(var)

Find the neighborhood of a variable using Lemma 27.

Parameters:

Name Type Description Default
var int

The variable whose neighborhood we want to find.

 required

Returns:

Type Description
ndarray

np.ndarray: 1D numpy array containing the variables in the neighborhood.
Source code in rcd/l_marvel/l_marvel.py 
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
def find_neighborhood(self, var: int) -> np.ndarray:
    """
    Find the neighborhood of a variable using Lemma 27.

    Args:
        var (int): The variable whose neighborhood we want to find.

    Returns:
        np.ndarray: 1D numpy array containing the variables in the neighborhood.
    """
    var_mk_bool_arr = self.markov_boundary_matrix[var]
    var_mk_arr = np.flatnonzero(var_mk_bool_arr)
    var_mk_set = set(var_mk_arr)

    neighbor_bool_arr = np.copy(var_mk_bool_arr)

    for var_y in var_mk_arr:
        # Check if Y is already a neighbor of X
        if not self.learned_skeleton.has_edge(var, var_y):
            if not self.is_neighbor(var, var_y, var_mk_set):
                neighbor_bool_arr[var_y] = False

    # Remove all variables that are not neighbors
    neighbors = np.flatnonzero(neighbor_bool_arr)
    return neighbors

is_neighbor(var, var_y, var_mk_set)

Check if var_y is a neighbor of variable var using Lemma 27.

Parameters:

Name Type Description Default
var int

Index of the variable.

 required
var_y int

The variable to check.

 required
var_mk_set Set[int]

Set of the variables in the Markov boundary of var.

 required

Returns:

Name Type Description
bool bool

True if var_y is a neighbor, False otherwise.
Source code in rcd/l_marvel/l_marvel.py 
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
def is_neighbor(self, var: int, var_y: int, var_mk_set: Set[int]) -> bool:
    """
    Check if var_y is a neighbor of variable var using Lemma 27.

    Args:
        var (int): Index of the variable.
        var_y (int): The variable to check.
        var_mk_set (Set[int]): Set of the variables in the Markov boundary of var.

    Returns:
        bool: True if var_y is a neighbor, False otherwise.
    """
    var_mk_left_list = list(var_mk_set - {var_y})
    # Use Lemma 27 and check all proper subsets of Mb(X) - {Y}
    for cond_set_size in range(len(var_mk_left_list) + 1):
        for var_s in itertools.combinations(var_mk_left_list, cond_set_size):
            cond_set = list(var_s)
            if self.ci_test(var, var_y, cond_set, self.data):
                # Y is a co-parent and thus NOT a neighbor
                return False
    return True

is_removable(var, neighbors)

Check whether a variable is removable using Theorem 32.

Parameters:

Name Type Description Default
var int

Index of the variable.

 required
neighbors ndarray

Neighbors of the variable.

 required

Returns:

Name Type Description
bool bool

True if the variable is removable, False otherwise.
Source code in rcd/l_marvel/l_marvel.py 
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
def is_removable(self, var: int, neighbors: np.ndarray) -> bool:
    """
    Check whether a variable is removable using Theorem 32.

    Args:
        var (int): Index of the variable.
        neighbors (np.ndarray): Neighbors of the variable.

    Returns:
        bool: True if the variable is removable, False otherwise.
    """
    var_mk_bool_arr = self.markov_boundary_matrix[var]
    var_mk_arr = np.flatnonzero(var_mk_bool_arr)
    var_mk_set = set(var_mk_arr)

    def cond_1(var_y, var_z):
        # there exists subset W in Mb(X) - {Y, Z}, s.t. Y ind. Z | W
        var_mk_left_list = list(var_mk_set - {var_y, var_z})
        for cond_set_size in range(len(var_mk_left_list) + 1):
            for var_s in itertools.combinations(var_mk_left_list, cond_set_size):
                cond_set = list(var_s)
                if self.ci_test(var_y, var_z, cond_set, self.data):
                    return True
        return False

    def cond_2(var_y, var_z):
        # for all subset W in Mb(X) - {Y, Z}, s.t. Y NOT ind. Z | W + {X}
        var_mk_left_left = list(var_mk_set - {var_y, var_z})
        for cond_set_size in range(len(var_mk_left_left) + 1):
            for var_s in itertools.combinations(var_mk_left_left, cond_set_size):
                cond_set = list(var_s) + [var]
                if self.ci_test(var_y, var_z, cond_set, self.data):
                    return False
        return True

    # Use Theorem 32 to check if X is removable. Loop over Y in Mb(X) and Z in Ne(X)
    for var_y in var_mk_arr:
        for var_z in neighbors:
            if var_y == var_z:
                continue
            xyz_tuple = (var, min(var_y, var_z), max(var_y, var_z))
            if xyz_tuple in self.skip_rem_check_set:
                continue
            if not (cond_1(var_y, var_z) or cond_2(var_y, var_z)):
                return False
            self.skip_rem_check_set.add(xyz_tuple)
    return True

learn_and_get_skeleton(data)

Run the L-MARVEL algorithm on the data to learn and return the learned skeleton graph.

Parameters:

Name Type Description Default
data ndarray

The data matrix with shape (num_samples, num_vars).

 required

Returns:

Type Description
Graph

nx.Graph: A networkx graph representing the learned skeleton.
Source code in rcd/l_marvel/l_marvel.py 
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
def learn_and_get_skeleton(self, data: np.ndarray) -> nx.Graph:
    """
    Run the L-MARVEL algorithm on the data to learn and return the learned skeleton graph.

    Args:
        data (np.ndarray): The data matrix with shape (num_samples, num_vars).

    Returns:
        nx.Graph: A networkx graph representing the learned skeleton.
    """
    # Initialize algorithm state
    self.num_vars = data.shape[1]
    self.data = data

    self.skip_rem_check_vec = np.zeros(self.num_vars, dtype=bool)
    self.skip_rem_check_set = set()
    self.neighbor_learned_arr = np.zeros(self.num_vars, dtype=bool)
    self.var_idx_set = set(range(self.num_vars))
    self.markov_boundary_matrix = self.find_markov_boundary_matrix(self.data)
    self.learned_skeleton = nx.Graph()
    self.learned_skeleton.add_nodes_from(range(self.num_vars))

    data_included_ci_test = lambda x, y, z: self.ci_test(x, y, z, self.data)

    var_arr = np.arange(self.num_vars)
    var_left_bool_arr = np.ones(self.num_vars, dtype=bool)  # Indicates if variable is left

    for _ in range(self.num_vars - 1):
        # sort variables by decreasing Markov boundary size
        # only sort variables that are still left and whose removability has NOT been checked
        var_to_sort_bool_arr = var_left_bool_arr & ~self.skip_rem_check_vec
        var_to_sort_arr = var_arr[var_to_sort_bool_arr]
        sorted_var_arr = sort_vars_by_mkb_size(self.markov_boundary_matrix[var_to_sort_bool_arr], var_to_sort_arr)

        removable_var = REMOVABLE_NOT_FOUND
        for var in sorted_var_arr:
            # Check whether we need to learn the neighbors of var
            if not self.neighbor_learned_arr[var]:
                neighbors = self.find_neighborhood(var)
                self.neighbor_learned_arr[var] = True

                # Add edges between the variable and its neighbors
                for neighbor in neighbors:
                    self.learned_skeleton.add_edge(var, neighbor)
            else:
                # If neighbors already learned, get them from the graph
                neighbors = list(self.learned_skeleton.neighbors(var))

                # Ensure only to include neighbors that are still left
                neighbors = [neighbor for neighbor in neighbors if var_left_bool_arr[neighbor]]

            # Check if variable is removable
            if self.is_removable(var, neighbors):
                removable_var = var
                break
            else:
                self.skip_rem_check_vec[var] = True

        if removable_var == REMOVABLE_NOT_FOUND:
            # If no removable found, pick the variable with the smallest Markov boundary from var_left_bool_arr
            var_left_arr = np.flatnonzero(var_left_bool_arr)
            mb_size_all = np.sum(self.markov_boundary_matrix[var_left_arr], axis=1)
            removable_var = var_left_arr[np.argmin(mb_size_all)]

            self.skip_rem_check_vec[:] = False
        else:
            # Remove the removable variable from the set of variables left
            var_left_bool_arr[removable_var] = False

        # Make sure to only include neighbors that are still left
        neighbors = [neighbor for neighbor in self.learned_skeleton.neighbors(removable_var) if var_left_bool_arr[neighbor]]

        # Update the Markov boundary matrix
        update_markov_boundary_matrix(
            self.markov_boundary_matrix,
            data_included_ci_test,
            removable_var,
            neighbors,
            skip_check=self.skip_rem_check_vec,
        )

    return self.learned_skeleton