Skip to content

RSL algorithms implementation details

Common functionality for Recursive Skeleton Learning algorithms.

Parameters

ci_test : Callable[[int, int, list[int], np.ndarray], bool] Conditional independence test. The callable receives the indices of two variables, a conditioning set, and the data matrix, and returns True when the variables are conditionally independent given the set. find_markov_boundary_matrix_fun : Callable[[np.ndarray], np.ndarray], optional Custom routine that estimates the Markov boundary matrix. If omitted, a Gaussian partial-correlation based estimator is used.

Notes

Concrete subclasses must implement :meth:find_neighborhood and :meth:is_removable.

Source code in rcd/rsl/rsl_base.py 
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 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
272
273
274
275
276
277
278
279
280
281
282
283
class _RSLBase:
    """Common functionality for Recursive Skeleton Learning algorithms.

    Parameters
    ----------
    ci_test : Callable[[int, int, list[int], np.ndarray], bool]
        Conditional independence test. The callable receives the indices of two
        variables, a conditioning set, and the data matrix, and returns ``True``
        when the variables are conditionally independent given the set.
    find_markov_boundary_matrix_fun : Callable[[np.ndarray], np.ndarray], optional
        Custom routine that estimates the Markov boundary matrix. If omitted, a
        Gaussian partial-correlation based estimator is used.

    Notes
    -----
    Concrete subclasses must implement :meth:`find_neighborhood` and
    :meth:`is_removable`.
    """

    def __init__(
        self,
        ci_test: Callable[[int, int, list[int], np.ndarray], bool],
        find_markov_boundary_matrix_fun: Callable[[np.ndarray], np.ndarray] | None = None,
    ) -> None:
        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: int | None = None
        self.data: np.ndarray | None = 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 a variable only when the corresponding flag is False)
        self.skip_rem_check_vec: np.ndarray | None = None  # SkipCheck_VEC in the paper
        self.markov_boundary_matrix: np.ndarray | None = None
        self.learned_skeleton: nx.Graph | None = None
        self.is_rsl_d: bool = False
        self.clique_num: int | None = None

    def learn_and_get_skeleton(
        self,
        data: np.ndarray,
        clique_num: int | None = None,
    ) -> nx.Graph:
        """Learn and return the skeleton implied by the data.

        Parameters
        ----------
        data : np.ndarray
            Matrix of shape ``(n_samples, n_vars)`` whose columns correspond to
            variables.
        clique_num : int, optional
            Upper bound on the clique number. Required for algorithms that are
            not diamond-free variants.

        Returns
        -------
        nx.Graph
            Undirected skeleton learned by the recursive elimination procedure.

        Raises
        ------
        ValueError
            If ``clique_num`` is not provided for algorithms that require it.
        """

        self._prepare_learning_state(data, clique_num)

        if self.num_vars is None:
            raise RuntimeError("Learning state failed to initialize number of variables.")

        skeleton = nx.Graph()
        skeleton.add_nodes_from(range(self.num_vars))

        self._run_recursive_elimination(build_skeleton=True, skeleton=skeleton)
        self.learned_skeleton = skeleton
        return skeleton

    def compute_removal_order(
        self,
        data: np.ndarray,
        clique_num: int | None = None,
    ) -> np.ndarray:
        """Compute the removal (r-) order without constructing the skeleton.

        Parameters
        ----------
        data : np.ndarray
            Matrix of shape ``(n_samples, n_vars)`` whose columns correspond to
            variables.
        clique_num : int, optional
            Upper bound on the clique number. Required for non diamond-free
            variants.

        Returns
        -------
        np.ndarray
            Integer array containing the order in which variables are removed.

        Raises
        ------
        ValueError
            If ``clique_num`` is not provided for algorithms that require it.
        """

        self._prepare_learning_state(data, clique_num)
        return self._run_recursive_elimination(build_skeleton=False)

    def find_neighborhood(self, var: int) -> np.ndarray:
        """Return the neighborhood of ``var`` in the current skeleton estimate.

        Parameters
        ----------
        var : int
            Index of the variable whose neighborhood is requested.

        Returns
        -------
        np.ndarray
            Indices that are neighbors of ``var``.
        """

        raise NotImplementedError

    def is_removable(self, var: int) -> bool:
        """Determine whether ``var`` can be removed without violating constraints.

        Parameters
        ----------
        var : int
            Index of the candidate variable.

        Returns
        -------
        bool
            ``True`` if the variable can be removed, ``False`` otherwise.
        """

        raise NotImplementedError

    def find_removable(self, var_arr: np.ndarray) -> int:
        """Locate the first removable variable in ``var_arr``.

        Parameters
        ----------
        var_arr : np.ndarray
            Candidate variable indices ordered by preference.

        Returns
        -------
        int
            Index of the first removable variable, or ``REMOVABLE_NOT_FOUND`` if
            none qualify.
        """

        if self.skip_rem_check_vec is None:
            raise RuntimeError("Learning state has not been initialized.")

        for var in var_arr:
            if self.is_removable(var):
                return var
            self.skip_rem_check_vec[var] = True
        return REMOVABLE_NOT_FOUND

    # ---------------------------------------------------------------------
    # Internal helpers
    # ---------------------------------------------------------------------

    def _prepare_learning_state(
        self,
        data: np.ndarray,
        clique_num: int | None,
    ) -> None:
        """Validate inputs and initialize shared state for a learning run."""

        if not self.is_rsl_d and clique_num is None:
            raise ValueError("Clique number not given!")

        self.num_vars = data.shape[1]
        self.data = data
        self.clique_num = clique_num

        self.skip_rem_check_vec = np.zeros(self.num_vars, dtype=bool)
        self.markov_boundary_matrix = self.find_markov_boundary_matrix(self.data)
        self.learned_skeleton = None

    def _run_recursive_elimination(
        self,
        build_skeleton: bool,
        skeleton: nx.Graph | None = None,
    ) -> np.ndarray:
        """Run the recursive elimination routine and optionally build the skeleton."""

        if self.data is None or self.markov_boundary_matrix is None:
            raise RuntimeError("Learning state has not been initialized.")

        if build_skeleton and skeleton is None:
            raise ValueError("Skeleton must be provided when build_skeleton is True.")

        if self.skip_rem_check_vec is None:
            raise RuntimeError("Skip-check vector is not initialized.")

        if self.num_vars is None:
            raise RuntimeError("Number of variables is unknown.")

        num_vars = self.num_vars
        r_order = np.zeros(num_vars, dtype=int)
        var_arr = np.arange(num_vars)
        var_left_mask = np.ones(num_vars, dtype=bool)

        def data_included_ci_test(x: int, y: int, z: list[int]) -> bool:
            """Closure that injects ``self.data`` into the CI test signature."""

            return self.ci_test(x, y, z, self.data)

        for iter_idx in range(num_vars - 1):
            # only consider variables that are still left and whose removal needs checking
            candidate_vars = var_arr[var_left_mask & ~self.skip_rem_check_vec]

            # sort the remaining candidates by the size of their Markov boundary
            mb_size = np.sum(self.markov_boundary_matrix[candidate_vars], axis=1)
            sorted_candidates = candidate_vars[np.argsort(mb_size, kind="stable")]

            removable_var = self.find_removable(sorted_candidates)

            if removable_var == REMOVABLE_NOT_FOUND:
                # if no removable variable was found, fall back to the smallest Markov boundary
                remaining_vars = np.flatnonzero(var_left_mask)
                mb_size_all = np.sum(self.markov_boundary_matrix[remaining_vars], axis=1)
                removable_var = remaining_vars[np.argmin(mb_size_all)]
                self.skip_rem_check_vec[:] = False

            # find the neighbors of the removable variable using the subclass rule
            neighbors = self.find_neighborhood(removable_var)

            # update the Markov boundary matrix to reflect the removal
            update_markov_boundary_matrix(
                self.markov_boundary_matrix,
                data_included_ci_test,
                removable_var,
                neighbors,
                self.is_rsl_d,
                skip_check=self.skip_rem_check_vec,
            )

            if build_skeleton and skeleton is not None:
                # add edges between the removable variable and its neighbors
                for neighbor_idx in neighbors:
                    skeleton.add_edge(removable_var, neighbor_idx)

            # mark the variable as removed and record it in the r-order
            var_left_mask[removable_var] = False
            r_order[iter_idx] = removable_var

        remaining_var = np.flatnonzero(var_left_mask)
        if remaining_var.size != 1:
            raise RuntimeError("Recursive elimination did not terminate correctly.")
        r_order[-1] = remaining_var[0]
        return r_order

compute_removal_order(data, clique_num=None)

Compute the removal (r-) order without constructing the skeleton.

Parameters

data : np.ndarray Matrix of shape (n_samples, n_vars) whose columns correspond to variables. clique_num : int, optional Upper bound on the clique number. Required for non diamond-free variants.

Returns

np.ndarray Integer array containing the order in which variables are removed.

Raises

ValueError If clique_num is not provided for algorithms that require it.

Source code in rcd/rsl/rsl_base.py 
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
def compute_removal_order(
    self,
    data: np.ndarray,
    clique_num: int | None = None,
) -> np.ndarray:
    """Compute the removal (r-) order without constructing the skeleton.

    Parameters
    ----------
    data : np.ndarray
        Matrix of shape ``(n_samples, n_vars)`` whose columns correspond to
        variables.
    clique_num : int, optional
        Upper bound on the clique number. Required for non diamond-free
        variants.

    Returns
    -------
    np.ndarray
        Integer array containing the order in which variables are removed.

    Raises
    ------
    ValueError
        If ``clique_num`` is not provided for algorithms that require it.
    """

    self._prepare_learning_state(data, clique_num)
    return self._run_recursive_elimination(build_skeleton=False)

find_neighborhood(var)

Return the neighborhood of var in the current skeleton estimate.

Parameters

var : int Index of the variable whose neighborhood is requested.

Returns

np.ndarray Indices that are neighbors of var.

Source code in rcd/rsl/rsl_base.py 
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
def find_neighborhood(self, var: int) -> np.ndarray:
    """Return the neighborhood of ``var`` in the current skeleton estimate.

    Parameters
    ----------
    var : int
        Index of the variable whose neighborhood is requested.

    Returns
    -------
    np.ndarray
        Indices that are neighbors of ``var``.
    """

    raise NotImplementedError

find_removable(var_arr)

Locate the first removable variable in var_arr.

Parameters

var_arr : np.ndarray Candidate variable indices ordered by preference.

Returns

int Index of the first removable variable, or REMOVABLE_NOT_FOUND if none qualify.

Source code in rcd/rsl/rsl_base.py 
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
def find_removable(self, var_arr: np.ndarray) -> int:
    """Locate the first removable variable in ``var_arr``.

    Parameters
    ----------
    var_arr : np.ndarray
        Candidate variable indices ordered by preference.

    Returns
    -------
    int
        Index of the first removable variable, or ``REMOVABLE_NOT_FOUND`` if
        none qualify.
    """

    if self.skip_rem_check_vec is None:
        raise RuntimeError("Learning state has not been initialized.")

    for var in var_arr:
        if self.is_removable(var):
            return var
        self.skip_rem_check_vec[var] = True
    return REMOVABLE_NOT_FOUND

is_removable(var)

Determine whether var can be removed without violating constraints.

Parameters

var : int Index of the candidate variable.

Returns

bool True if the variable can be removed, False otherwise.

Source code in rcd/rsl/rsl_base.py 
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
def is_removable(self, var: int) -> bool:
    """Determine whether ``var`` can be removed without violating constraints.

    Parameters
    ----------
    var : int
        Index of the candidate variable.

    Returns
    -------
    bool
        ``True`` if the variable can be removed, ``False`` otherwise.
    """

    raise NotImplementedError

learn_and_get_skeleton(data, clique_num=None)

Learn and return the skeleton implied by the data.

Parameters

data : np.ndarray Matrix of shape (n_samples, n_vars) whose columns correspond to variables. clique_num : int, optional Upper bound on the clique number. Required for algorithms that are not diamond-free variants.

Returns

nx.Graph Undirected skeleton learned by the recursive elimination procedure.

Raises

ValueError If clique_num is not provided for algorithms that require it.

Source code in rcd/rsl/rsl_base.py 
 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
def learn_and_get_skeleton(
    self,
    data: np.ndarray,
    clique_num: int | None = None,
) -> nx.Graph:
    """Learn and return the skeleton implied by the data.

    Parameters
    ----------
    data : np.ndarray
        Matrix of shape ``(n_samples, n_vars)`` whose columns correspond to
        variables.
    clique_num : int, optional
        Upper bound on the clique number. Required for algorithms that are
        not diamond-free variants.

    Returns
    -------
    nx.Graph
        Undirected skeleton learned by the recursive elimination procedure.

    Raises
    ------
    ValueError
        If ``clique_num`` is not provided for algorithms that require it.
    """

    self._prepare_learning_state(data, clique_num)

    if self.num_vars is None:
        raise RuntimeError("Learning state failed to initialize number of variables.")

    skeleton = nx.Graph()
    skeleton.add_nodes_from(range(self.num_vars))

    self._run_recursive_elimination(build_skeleton=True, skeleton=skeleton)
    self.learned_skeleton = skeleton
    return skeleton