What I did was to preprocess the data. I created a @bools2 that maps each boolean to the sets it is in. @sets2 contains the number of booleans in each set. And instead of splicing a set from @sets, I delete it:

use strict;
use warnings;

use List::Util qw[shuffle];

our $N //= 1e3;
our $S //= 1e4;
our $Z //= 5;

my @bools = (0) x $N;
    
my @sets = map [map int(rand $N), 1 .. 1 + int(rand $Z)], 1 .. $S;
my @someOrder = shuffle 0 .. $N - 1; 

#
# Preprocessing
#
my @bool2;
my @sets2;
for (my $i = 0; $i < @sets; $i ++) {
    push @{$bool2[$_]}, $i for @{$sets[$i]};
    $sets2[$i] = @{$sets[$i]};
}

#
# Main loop
#
foreach my $next (@someOrder) {
    foreach my $set_nr (@{$bool2[$next]}) {
        delete $sets[$set_nr] unless --$sets2[$set_nr];
    }
}
[download]

I think this algorithm is O(k), where k is the sum of the number of elements in each set.