diff --git a/mne/cov.py b/mne/cov.py
index df5918cf727..6d440322a2c 100644
--- a/mne/cov.py
+++ b/mne/cov.py
@@ -553,16 +553,17 @@ def make_ad_hoc_cov(info, std=None, *, verbose=None):
     return Covariance(data, ch_names, info["bads"], info["projs"], nfree=0)
 
 
-def _check_n_samples(n_samples, n_chan):
+def _check_n_samples(n_samples, n_chan, on_few_samples="warn"):
     """Check to see if there are enough samples for reliable cov calc."""
     n_samples_min = 10 * (n_chan + 1) // 2
     if n_samples <= 0:
         raise ValueError("No samples found to compute the covariance matrix")
     if n_samples < n_samples_min:
-        warn(
+        msg = (
             f"Too few samples (required : {n_samples_min} got : {n_samples}), "
             "covariance estimate may be unreliable"
         )
+        _on_missing(on_few_samples, msg, "on_few_samples")
 
 
 @verbose
@@ -582,6 +583,7 @@ def compute_raw_covariance(
     return_estimators=False,
     reject_by_annotation=True,
     rank=None,
+    on_few_samples="warn",
     verbose=None,
 ):
     """Estimate noise covariance matrix from a continuous segment of raw data.
@@ -662,6 +664,10 @@ def compute_raw_covariance(
 
         .. versionadded:: 0.18
            Support for 'info' mode.
+    on_few_samples : str
+        Can be 'warn' (default), 'ignore', or 'raise' to control behavior when
+        there are fewer samples than channels, which can lead to inaccurate
+        covariance or rank estimates.
     %(verbose)s
 
     Returns
@@ -736,7 +742,7 @@ def compute_raw_covariance(
             mu += raw_segment.sum(axis=1)
             data += np.dot(raw_segment, raw_segment.T)
             n_samples += raw_segment.shape[1]
-        _check_n_samples(n_samples, len(picks))
+        _check_n_samples(n_samples, len(picks), on_few_samples)
         data -= mu[:, None] * (mu[None, :] / n_samples)
         data /= n_samples - 1.0
         logger.info("Number of samples used : %d", n_samples)
@@ -872,6 +878,7 @@ def compute_covariance(
     return_estimators=False,
     on_mismatch="raise",
     rank=None,
+    on_few_samples="warn",
     verbose=None,
 ):
     """Estimate noise covariance matrix from epochs.
@@ -966,6 +973,10 @@ def compute_covariance(
 
         .. versionadded:: 0.18
            Support for 'info' mode.
+    on_few_samples : str
+        Can be 'warn' (default), 'ignore', or 'raise' to control behavior when
+        there are fewer samples than channels, which can lead to inaccurate
+        covariance or rank estimates.
     %(verbose)s
 
     Returns
@@ -1144,7 +1155,7 @@ def _unpack_epochs(epochs):
 
     epochs = np.hstack(epochs)
     n_samples_tot = epochs.shape[-1]
-    _check_n_samples(n_samples_tot, len(picks_meeg))
+    _check_n_samples(n_samples_tot, len(picks_meeg), on_few_samples)
 
     epochs = epochs.T  # sklearn | C-order
     cov_data = _compute_covariance_auto(
@@ -1158,6 +1169,7 @@ def _unpack_epochs(epochs):
         picks_list=picks_list,
         scalings=scalings,
         rank=rank,
+        on_few_samples=on_few_samples,
     )
 
     if keep_sample_mean is False:
@@ -1221,7 +1233,7 @@ def _eigvec_subspace(eig, eigvec, mask):
 
 @verbose
 def _compute_rank_raw_array(
-    data, info, rank, scalings, *, log_ch_type=None, verbose=None
+    data, info, rank, scalings, *, log_ch_type=None, on_few_samples="warn", verbose=None
 ):
     from .io import RawArray
 
@@ -1231,6 +1243,7 @@ def _compute_rank_raw_array(
         scalings,
         info,
         log_ch_type=log_ch_type,
+        on_few_samples=on_few_samples,
     )
 
 
@@ -1249,6 +1262,7 @@ def _compute_covariance_auto(
     cov_kind="",
     log_ch_type=None,
     log_rank=True,
+    on_few_samples="warn",
 ):
     """Compute covariance auto mode."""
     # rescale to improve numerical stability
@@ -1258,6 +1272,7 @@ def _compute_covariance_auto(
         info,
         rank=rank,
         scalings=scalings,
+        on_few_samples=on_few_samples,
         verbose=_verbose_safe_false(),
     )
     with _scaled_array(data.T, picks_list, scalings):
@@ -1268,6 +1283,7 @@ def _compute_covariance_auto(
             rank,
             proj_subspace=True,
             do_compute_rank=False,
+            on_few_samples=on_few_samples,
             log_ch_type=log_ch_type,
             verbose=None if log_rank else _verbose_safe_false(),
         )
@@ -1729,6 +1745,7 @@ def prepare_noise_cov(
     rank=None,
     scalings=None,
     on_rank_mismatch="ignore",
+    on_few_samples="warn",
     verbose=None,
 ):
     """Prepare noise covariance matrix.
@@ -1751,6 +1768,10 @@ def prepare_noise_cov(
 
             dict(mag=1e12, grad=1e11, eeg=1e5)
     %(on_rank_mismatch)s
+    on_few_samples : str
+        Can be 'warn' (default), 'ignore', or 'raise' to control behavior when
+        there are fewer samples than channels, which can lead to inaccurate
+        covariance or rank estimates.
     %(verbose)s
 
     Returns
@@ -1792,6 +1813,7 @@ def prepare_noise_cov(
         projs,
         ch_names,
         on_rank_mismatch=on_rank_mismatch,
+        on_few_samples=on_few_samples,
     )
     noise_cov.update(eig=eig, eigvec=eigvec)
     return noise_cov
@@ -1808,6 +1830,7 @@ def _smart_eigh(
     proj_subspace=False,
     do_compute_rank=True,
     on_rank_mismatch="ignore",
+    on_few_samples="warn",
     *,
     log_ch_type=None,
     verbose=None,
@@ -1838,6 +1861,7 @@ def _smart_eigh(
             scalings,
             info,
             on_rank_mismatch=on_rank_mismatch,
+            on_few_samples=on_few_samples,
             log_ch_type=log_ch_type,
         )
     assert C.ndim == 2 and C.shape[0] == C.shape[1]
@@ -1916,6 +1940,7 @@ def regularize(
     dbs=0.1,
     rank=None,
     scalings=None,
+    on_few_samples="warn",
     verbose=None,
 ):
     """Regularize noise covariance matrix.
@@ -1978,6 +2003,10 @@ def regularize(
         See :func:`mne.compute_covariance`.
 
         .. versionadded:: 0.17
+    on_few_samples : str
+        Can be 'warn' (default), 'ignore', or 'raise' to control behavior when
+        there are fewer samples than channels, which can lead to inaccurate
+        covariance or rank estimates.
     %(verbose)s
 
     Returns
@@ -2032,7 +2061,7 @@ def regularize(
     else:
         regs.update(mag=mag, grad=grad)
     if rank != "full":
-        rank = _compute_rank(cov, rank, scalings, info)
+        rank = _compute_rank(cov, rank, scalings, info, on_few_samples=on_few_samples)
 
     info_ch_names = info["ch_names"]
     ch_names_by_type = dict()
@@ -2092,7 +2121,9 @@ def regularize(
             this_info = pick_info(info, this_picks)
             # Here we could use proj_subspace=True, but this should not matter
             # since this is already in a loop over channel types
-            _, eigvec, mask = _smart_eigh(this_C, this_info, rank)
+            _, eigvec, mask = _smart_eigh(
+                this_C, this_info, rank, on_few_samples=on_few_samples
+            )
             U = eigvec[mask].T
         this_C = np.dot(U.T, np.dot(this_C, U))
 
@@ -2119,6 +2150,7 @@ def _regularized_covariance(
     log_ch_type=None,
     log_rank=None,
     cov_kind="",
+    on_few_samples="warn",
     verbose=None,
 ):
     """Compute a regularized covariance from data using sklearn.
@@ -2166,6 +2198,7 @@ def _regularized_covariance(
         cov_kind=cov_kind,
         log_ch_type=log_ch_type,
         log_rank=log_rank,
+        on_few_samples=on_few_samples,
     )[reg]["data"]
     return cov
 
diff --git a/mne/decoding/_covs_ged.py b/mne/decoding/_covs_ged.py
index 074b4c15625..e6421d94757 100644
--- a/mne/decoding/_covs_ged.py
+++ b/mne/decoding/_covs_ged.py
@@ -29,6 +29,7 @@ def _concat_cov(x_class, *, cov_kind, log_rank, reg, cov_method_params, info, ra
         cov_kind=cov_kind,
         log_rank=log_rank,
         log_ch_type="data",
+        on_few_samples="ignore",
     )
 
     return cov, n_channels  # the weight here is just the number of channels
@@ -88,6 +89,7 @@ def _csp_estimate(X, y, reg, cov_method_params, cov_est, info, rank, norm_trace)
             rank=rank,
             scalings=None,
             log_ch_type="data",
+            on_few_samples="ignore",
         )
 
     covs = []
@@ -130,7 +132,12 @@ def _xdawn_estimate(
     # Retrieve or compute whitening covariance
     if R is None:
         R = _regularized_covariance(
-            np.hstack(X), reg, cov_method_params, info, rank=rank
+            np.hstack(X),
+            reg,
+            cov_method_params,
+            info,
+            rank=rank,
+            on_few_samples="ignore",
         )
     elif isinstance(R, Covariance):
         R = R.data
@@ -146,7 +153,9 @@ def _xdawn_estimate(
     for evo, toeplitz in zip(evokeds, toeplitzs):
         # Estimate covariance matrix of the prototype response
         evo = np.dot(evo, toeplitz)
-        evo_cov = _regularized_covariance(evo, reg, cov_method_params, info, rank=rank)
+        evo_cov = _regularized_covariance(
+            evo, reg, cov_method_params, info, rank=rank, on_few_samples="ignore"
+        )
         covs.append(evo_cov)
 
     covs.append(R)
@@ -158,6 +167,7 @@ def _xdawn_estimate(
             rank=rank,
             scalings=None,
             log_ch_type="data",
+            on_few_samples="ignore",
         )
     return covs, C_ref, info, rank, dict()
 
@@ -197,6 +207,7 @@ def _ssd_estimate(
         method_params=cov_method_params,
         rank="full",
         info=picked_info,
+        on_few_samples="ignore",
     )
     R = _regularized_covariance(
         X_noise,
@@ -204,6 +215,7 @@ def _ssd_estimate(
         method_params=cov_method_params,
         rank="full",
         info=picked_info,
+        on_few_samples="ignore",
     )
     covs = [S, R]
     C_ref = S
@@ -223,6 +235,7 @@ def _ssd_estimate(
                 rank,
                 _handle_default("scalings_cov_rank", None),
                 info,
+                on_few_samples="ignore",
             ).values()
         )[0]
     all_ranks.append(r)
@@ -266,6 +279,7 @@ def _spoc_estimate(X, y, reg, cov_method_params, info, rank):
             rank=rank,
             log_ch_type="data",
             log_rank=ii == 0,
+            on_few_samples="ignore",
         )
 
     S = np.mean(covs * target[:, np.newaxis, np.newaxis], axis=0)
@@ -280,5 +294,6 @@ def _spoc_estimate(X, y, reg, cov_method_params, info, rank):
             rank=rank,
             scalings=None,
             log_ch_type="data",
+            on_few_samples="ignore",
         )
     return covs, C_ref, info, rank, dict()
diff --git a/mne/decoding/tests/test_csp.py b/mne/decoding/tests/test_csp.py
index 0f26b7c30c0..4411267b407 100644
--- a/mne/decoding/tests/test_csp.py
+++ b/mne/decoding/tests/test_csp.py
@@ -431,7 +431,7 @@ def test_spoc():
     # check y
     pytest.raises(ValueError, spoc.fit, X, y * 0)
 
-    # Check that doesn't take CSP-spcific input
+    # Check that doesn't take CSP-specific input
     pytest.raises(TypeError, SPoC, cov_est="epoch")
 
     # Check mixing matrix on simulated data
diff --git a/mne/decoding/tests/test_ged.py b/mne/decoding/tests/test_ged.py
index 4781553fd01..c6935dba362 100644
--- a/mne/decoding/tests/test_ged.py
+++ b/mne/decoding/tests/test_ged.py
@@ -89,13 +89,17 @@ def _mock_cov_callable(X, y, cov_method_params=None, compute_C_ref=True):
     for ci, this_class in enumerate(classes):
         class_data = X[y == this_class]
         class_data = class_data.transpose(1, 0, 2).reshape(n_channels, -1)
-        cov = _regularized_covariance(class_data, **cov_method_params)
+        cov = _regularized_covariance(
+            class_data, on_few_samples="ignore", **cov_method_params
+        )
         covs.append(cov)
         sample_weights.append(class_data.shape[0])
 
     ref_data = X.transpose(1, 0, 2).reshape(n_channels, -1)
     if compute_C_ref:
-        C_ref = _regularized_covariance(ref_data, **cov_method_params)
+        C_ref = _regularized_covariance(
+            ref_data, on_few_samples="ignore", **cov_method_params
+        )
     else:
         C_ref = None
     info = _mock_info(n_channels)
diff --git a/mne/rank.py b/mne/rank.py
index 22ae5068732..e25ad30f395 100644
--- a/mne/rank.py
+++ b/mne/rank.py
@@ -130,14 +130,27 @@ def _estimate_rank_from_s(s, tol="auto", tol_kind="absolute"):
 
 
 def _estimate_rank_raw(
-    raw, picks=None, tol=1e-4, scalings="norm", with_ref_meg=False, tol_kind="absolute"
+    raw,
+    picks=None,
+    tol=1e-4,
+    scalings="norm",
+    with_ref_meg=False,
+    tol_kind="absolute",
+    on_few_samples="warn",
 ):
     """Aid the transition away from raw.estimate_rank."""
     if picks is None:
         picks = _picks_to_idx(raw.info, picks, with_ref_meg=with_ref_meg)
     # conveniency wrapper to expose the expert "tol" option + scalings options
     return _estimate_rank_meeg_signals(
-        raw[picks][0], pick_info(raw.info, picks), scalings, tol, False, tol_kind
+        raw[picks][0],
+        pick_info(raw.info, picks),
+        scalings,
+        tol,
+        False,
+        tol_kind,
+        log_ch_type=None,
+        on_few_samples=on_few_samples,
     )
 
 
@@ -150,6 +163,7 @@ def _estimate_rank_meeg_signals(
     return_singular=False,
     tol_kind="absolute",
     log_ch_type=None,
+    on_few_samples="warn",
 ):
     """Estimate rank for M/EEG data.
 
@@ -173,6 +187,10 @@ def _estimate_rank_meeg_signals(
         to determine the rank.
     tol_kind : str
         Tolerance kind. See ``estimate_rank``.
+    on_few_samples : str
+        Can be 'warn' (default), 'ignore', or 'raise' to control behavior when
+        there are fewer samples than channels, which can lead to inaccurate rank
+        estimates.
 
     Returns
     -------
@@ -184,10 +202,11 @@ def _estimate_rank_meeg_signals(
     """
     picks_list = _picks_by_type(info)
     if data.shape[1] < data.shape[0]:
-        ValueError(
+        msg = (
             "You've got fewer samples than channels, your "
             "rank estimate might be inaccurate."
         )
+        _on_missing(on_few_samples, msg, "on_few_samples")
     with _scaled_array(data, picks_list, scalings):
         out = estimate_rank(
             data,
@@ -214,6 +233,7 @@ def _estimate_rank_meeg_cov(
     return_singular=False,
     *,
     log_ch_type=None,
+    on_few_samples="warn",
     verbose=None,
 ):
     """Estimate rank of M/EEG covariance data, given the covariance.
@@ -236,6 +256,10 @@ def _estimate_rank_meeg_cov(
     return_singular : bool
         If True, also return the singular values that were used
         to determine the rank.
+    on_few_samples : str
+        Can be 'warn' (default), 'ignore', or 'raise' to control behavior when
+        there are fewer samples than channels, which can lead to inaccurate rank
+        estimates.
 
     Returns
     -------
@@ -249,10 +273,11 @@ def _estimate_rank_meeg_cov(
     scalings = _handle_default("scalings_cov_rank", scalings)
     _apply_scaling_cov(data, picks_list, scalings)
     if data.shape[1] < data.shape[0]:
-        ValueError(
+        msg = (
             "You've got fewer samples than channels, your "
             "rank estimate might be inaccurate."
         )
+        _on_missing(on_few_samples, msg, "on_few_samples")
     out = estimate_rank(data, tol=tol, norm=False, return_singular=return_singular)
     rank = out[0] if isinstance(out, tuple) else out
     if log_ch_type is None:
@@ -338,6 +363,7 @@ def compute_rank(
     proj=True,
     tol_kind="absolute",
     on_rank_mismatch="ignore",
+    on_few_samples="warn",
     verbose=None,
 ):
     """Compute the rank of data or noise covariance.
@@ -363,6 +389,10 @@ def compute_rank(
         considered when ``rank=None`` or ``rank='info'``.
     %(tol_kind_rank)s
     %(on_rank_mismatch)s
+    on_few_samples : str
+        Can be 'warn' (default), 'ignore', or 'raise' to control behavior when
+        there are fewer samples than channels, which can lead to inaccurate rank
+        estimates.
     %(verbose)s
 
     Returns
@@ -384,6 +414,7 @@ def compute_rank(
         proj=proj,
         tol_kind=tol_kind,
         on_rank_mismatch=on_rank_mismatch,
+        on_few_samples=on_few_samples,
     )
 
 
@@ -398,6 +429,7 @@ def _compute_rank(
     proj=True,
     tol_kind="absolute",
     on_rank_mismatch="ignore",
+    on_few_samples="warn",
     log_ch_type=None,
     verbose=None,
 ):
@@ -503,6 +535,7 @@ def _compute_rank(
                     False,
                     tol_kind,
                     log_ch_type=log_ch_type,
+                    on_few_samples=on_few_samples,
                 )
             else:
                 assert isinstance(inst, Covariance)
@@ -520,6 +553,7 @@ def _compute_rank(
                         tol,
                         return_singular=True,
                         log_ch_type=log_ch_type,
+                        on_few_samples=on_few_samples,
                         verbose=est_verbose,
                     )
                     if ch_type in rank:
diff --git a/mne/simulation/tests/test_evoked.py b/mne/simulation/tests/test_evoked.py
index 003b16ef47d..167d01c1843 100644
--- a/mne/simulation/tests/test_evoked.py
+++ b/mne/simulation/tests/test_evoked.py
@@ -100,6 +100,7 @@ def test_simulate_evoked():
 @pytest.mark.filterwarnings("ignore:No average EEG reference present")
 @pytest.mark.filterwarnings("ignore:Too few samples")
 @pytest.mark.filterwarnings("ignore:Epochs are not baseline corrected")
+@pytest.mark.filterwarnings("ignore:.*You've got fewer samples than channels.*")
 def test_add_noise():
     """Test noise addition."""
     rng = np.random.default_rng(0)
diff --git a/mne/tests/test_cov.py b/mne/tests/test_cov.py
index ad05f5d8b1f..a29c5a41952 100644
--- a/mne/tests/test_cov.py
+++ b/mne/tests/test_cov.py
@@ -72,8 +72,7 @@ def test_compute_whitener(proj, pca):
         raw.apply_proj()
     else:
         raw.del_proj()
-    with pytest.warns(RuntimeWarning, match="Too few samples"):
-        cov = compute_raw_covariance(raw)
+    cov = compute_raw_covariance(raw, on_few_samples="ignore")
     assert cov["names"] == raw.ch_names
     W, _, C = compute_whitener(
         cov, raw.info, pca=pca, return_colorer=True, verbose="error"
@@ -96,9 +95,8 @@ def test_compute_whitener(proj, pca):
 
     raw.info["bads"] = [raw.ch_names[0]]
     picks = pick_types(raw.info, meg=True, eeg=True, exclude=[])
-    with pytest.warns(RuntimeWarning, match="Too few samples"):
-        cov2 = compute_raw_covariance(raw, picks=picks)
-        cov3 = compute_raw_covariance(raw, picks=None)
+    cov2 = compute_raw_covariance(raw, picks=picks, on_few_samples="ignore")
+    cov3 = compute_raw_covariance(raw, picks=None, on_few_samples="ignore")
     assert_allclose(cov2["data"][1:, 1:], cov3["data"])
     W2, _, C2 = compute_whitener(
         cov2, raw.info, pca=pca, return_colorer=True, picks=picks, verbose="error"
@@ -434,14 +432,16 @@ def test_cov_estimation_with_triggers(rank, tmp_path):
         preload=True,
     )
 
-    cov = compute_covariance(epochs, keep_sample_mean=True)
+    cov = compute_covariance(epochs, keep_sample_mean=True, on_few_samples="ignore")
     cov_km = read_cov(cov_km_fname)
     # adjust for nfree bug
     cov_km["nfree"] -= 1
     _assert_cov(cov, cov_km)
 
     # Test with tmin and tmax (different but not too much)
-    cov_tmin_tmax = compute_covariance(epochs, tmin=-0.19, tmax=-0.01)
+    cov_tmin_tmax = compute_covariance(
+        epochs, tmin=-0.19, tmax=-0.01, on_few_samples="ignore"
+    )
     assert np.all(cov.data != cov_tmin_tmax.data)
     err = np.linalg.norm(cov.data - cov_tmin_tmax.data) / np.linalg.norm(
         cov_tmin_tmax.data
@@ -462,12 +462,12 @@ def test_cov_estimation_with_triggers(rank, tmp_path):
         )
         for ev_id in event_ids
     ]
-    cov2 = compute_covariance(epochs, keep_sample_mean=True)
+    cov2 = compute_covariance(epochs, keep_sample_mean=True, on_few_samples="ignore")
     assert_array_almost_equal(cov.data, cov2.data)
     assert cov.ch_names == cov2.ch_names
 
     # cov with keep_sample_mean=False using a list of epochs
-    cov = compute_covariance(epochs, keep_sample_mean=False)
+    cov = compute_covariance(epochs, keep_sample_mean=False, on_few_samples="ignore")
     assert cov_km.nfree == cov.nfree
     _assert_cov(cov, read_cov(cov_fname), nfree=False)
 
@@ -506,10 +506,10 @@ def test_cov_estimation_with_triggers(rank, tmp_path):
     pytest.raises(ValueError, compute_covariance, epochs)
     pytest.raises(ValueError, compute_covariance, epochs, projs=None)
     # these should work, but won't be equal to above
-    with pytest.warns(RuntimeWarning, match="Too few samples"):
-        cov = compute_covariance(epochs, projs=epochs[0].info["projs"])
-    with pytest.warns(RuntimeWarning, match="Too few samples"):
-        cov = compute_covariance(epochs, projs=[])
+    cov = compute_covariance(
+        epochs, projs=epochs[0].info["projs"], on_few_samples="ignore"
+    )
+    cov = compute_covariance(epochs, projs=[], on_few_samples="ignore")
 
     # test new dict support
     epochs = Epochs(
@@ -522,10 +522,8 @@ def test_cov_estimation_with_triggers(rank, tmp_path):
         reject=reject,
         preload=True,
     )
-    with pytest.warns(RuntimeWarning, match="Too few samples"):
-        compute_covariance(epochs)
-    with pytest.warns(RuntimeWarning, match="Too few samples"):
-        compute_covariance(epochs, projs=[])
+    compute_covariance(epochs, on_few_samples="ignore")
+    compute_covariance(epochs, projs=[], on_few_samples="ignore")
     pytest.raises(TypeError, compute_covariance, epochs, projs="foo")
     pytest.raises(TypeError, compute_covariance, epochs, projs=["foo"])
 
@@ -799,10 +797,14 @@ def test_low_rank_methods(rank, raw_epochs_events):
             empirical=(-15000, -5000), diagonal_fixed=(-700, -600), oas=(-700, -600)
         ),
     }
-    with pytest.warns(RuntimeWarning, match="Too few samples"):
-        covs = compute_covariance(
-            epochs, method=methods, return_estimators=True, rank=rank, verbose=True
-        )
+    covs = compute_covariance(
+        epochs,
+        method=methods,
+        return_estimators=True,
+        rank=rank,
+        verbose=True,
+        on_few_samples="ignore",
+    )
     for cov in covs:
         method = cov["method"]
         these_bounds = bounds[str(rank)][method]
@@ -821,23 +823,33 @@ def test_low_rank_cov(raw_epochs_events):
     sss_proj_rank = 139  # 80 MEG + 60 EEG - 1 proj
     n_ch = 366
     proj_rank = 365  # one EEG proj
-    with pytest.warns(RuntimeWarning, match="Too few samples"):
-        emp_cov = compute_covariance(epochs)
+    emp_cov = compute_covariance(epochs, on_few_samples="ignore")
     # Test equivalence with mne.cov.regularize subspace
     with pytest.raises(ValueError, match="are dependent.*must equal"):
         regularize(emp_cov, epochs.info, rank=None, mag=0.1, grad=0.2)
     assert _cov_rank(emp_cov, epochs.info) == sss_proj_rank
-    reg_cov = regularize(emp_cov, epochs.info, proj=True, rank="full")
+    reg_cov = regularize(
+        emp_cov, epochs.info, proj=True, rank="full", on_few_samples="ignore"
+    )
     assert _cov_rank(reg_cov, epochs.info) == proj_rank
     with pytest.warns(RuntimeWarning, match="exceeds the theoretical"):
         _compute_rank_int(reg_cov, info=epochs.info)
     del reg_cov
     with catch_logging() as log:
-        reg_r_cov = regularize(emp_cov, epochs.info, proj=True, rank=None, verbose=True)
+        reg_r_cov = regularize(
+            emp_cov,
+            epochs.info,
+            proj=True,
+            rank=None,
+            verbose=True,
+            on_few_samples="ignore",
+        )
     log = log.getvalue()
     assert "jointly" in log
     assert _cov_rank(reg_r_cov, epochs.info) == sss_proj_rank
-    reg_r_only_cov = regularize(emp_cov, epochs.info, proj=False, rank=None)
+    reg_r_only_cov = regularize(
+        emp_cov, epochs.info, proj=False, rank=None, on_few_samples="ignore"
+    )
     assert _cov_rank(reg_r_only_cov, epochs.info) == sss_proj_rank
     assert_allclose(reg_r_only_cov["data"], reg_r_cov["data"])
     del reg_r_only_cov, reg_r_cov
@@ -848,15 +860,20 @@ def test_low_rank_cov(raw_epochs_events):
     with epochs_meg.info._unlock():
         epochs_meg.info.update(bads=[], projs=[])
     cov_full = compute_covariance(
-        epochs_meg, method="oas", rank="full", verbose="error"
+        epochs_meg, method="oas", rank="full", verbose="error", on_few_samples="ignore"
     )
     assert _cov_rank(cov_full, epochs_meg.info) == 306
-    with pytest.warns(RuntimeWarning, match="few samples"):
-        cov_dict = compute_covariance(epochs_meg, method="oas", rank=dict(meg=306))
+    cov_dict = compute_covariance(
+        epochs_meg, method="oas", rank=dict(meg=306), on_few_samples="ignore"
+    )
     assert _cov_rank(cov_dict, epochs_meg.info) == 306
     assert_allclose(cov_full["data"], cov_dict["data"])
     cov_dict = compute_covariance(
-        epochs_meg, method="oas", rank=dict(meg=306), verbose="error"
+        epochs_meg,
+        method="oas",
+        rank=dict(meg=306),
+        verbose="error",
+        on_few_samples="ignore",
     )
     assert _cov_rank(cov_dict, epochs_meg.info) == 306
     assert_allclose(cov_full["data"], cov_dict["data"])
@@ -870,14 +887,24 @@ def test_low_rank_cov(raw_epochs_events):
     assert len(raw.info["projs"]) == 3
     epochs = Epochs(raw, events, tmin=-0.2, tmax=0, preload=True)
     assert len(raw.ch_names) == n_ch
-    emp_cov = compute_covariance(epochs, rank="full", verbose="error")
+    emp_cov = compute_covariance(
+        epochs, rank="full", verbose="error", on_few_samples="ignore"
+    )
     assert _cov_rank(emp_cov, epochs.info) == rank
-    reg_cov = regularize(emp_cov, epochs.info, proj=True, rank="full")
+    reg_cov = regularize(
+        emp_cov, epochs.info, proj=True, rank="full", on_few_samples="ignore"
+    )
     assert _cov_rank(reg_cov, epochs.info) == rank
-    reg_r_cov = regularize(emp_cov, epochs.info, proj=False, rank=None)
+    reg_r_cov = regularize(
+        emp_cov, epochs.info, proj=False, rank=None, on_few_samples="ignore"
+    )
     assert _cov_rank(reg_r_cov, epochs.info) == rank
     dia_cov = compute_covariance(
-        epochs, rank=None, method="diagonal_fixed", verbose="error"
+        epochs,
+        rank=None,
+        method="diagonal_fixed",
+        verbose="error",
+        on_few_samples="ignore",
     )
     assert _cov_rank(dia_cov, epochs.info) == rank
     assert_allclose(dia_cov["data"], reg_cov["data"])
@@ -904,20 +931,24 @@ def test_cov_ctf():
     for comp in [0, 1]:
         raw.apply_gradient_compensation(comp)
         epochs = Epochs(raw, events, None, -0.2, 0.2, preload=True)
-        with _record_warnings(), pytest.warns(RuntimeWarning, match="Too few samples"):
-            noise_cov = compute_covariance(epochs, tmax=0.0, method=["empirical"])
+        with _record_warnings():
+            noise_cov = compute_covariance(
+                epochs, tmax=0.0, method=["empirical"], on_few_samples="ignore"
+            )
         with pytest.warns(RuntimeWarning, match="orders of magnitude"):
-            prepare_noise_cov(noise_cov, raw.info, ch_names)
+            prepare_noise_cov(noise_cov, raw.info, ch_names, on_few_samples="ignore")
 
     raw.apply_gradient_compensation(0)
     epochs = Epochs(raw, events, None, -0.2, 0.2, preload=True)
-    with _record_warnings(), pytest.warns(RuntimeWarning, match="Too few samples"):
-        noise_cov = compute_covariance(epochs, tmax=0.0, method=["empirical"])
+    with _record_warnings():
+        noise_cov = compute_covariance(
+            epochs, tmax=0.0, method=["empirical"], on_few_samples="ignore"
+        )
     raw.apply_gradient_compensation(1)
 
     # TODO This next call in principle should fail.
     with pytest.warns(RuntimeWarning, match="orders of magnitude"):
-        prepare_noise_cov(noise_cov, raw.info, ch_names)
+        prepare_noise_cov(noise_cov, raw.info, ch_names, on_few_samples="ignore")
 
     # make sure comps matrices was not removed from raw
     assert raw.info["comps"], "Comps matrices removed"