diff --git a/docs/src/algorithms.md b/docs/src/algorithms.md
index 8a3ebc1..2d27baf 100644
--- a/docs/src/algorithms.md
+++ b/docs/src/algorithms.md
@@ -10,7 +10,7 @@ Qubes represent sets of identifiers, so the familiar set operations are all defi
 |---|---|---|
 | **Union** | `a.append(&mut b)` | All identifiers in A or B (or both) |
 | **Intersection** | `select` with intersection logic | Identifiers in both A and B |
-| **Difference** | internal set operation | Identifiers in A but not B |
+| **Difference** | `a.subtract(&b)` | Identifiers in A but not B |
 | **Symmetric difference** | internal set operation | Identifiers in exactly one of A or B |
 
 ### How It Works
diff --git a/docs/src/python/py_qubed.md b/docs/src/python/py_qubed.md
index 8582923..8beb99d 100644
--- a/docs/src/python/py_qubed.md
+++ b/docs/src/python/py_qubed.md
@@ -179,6 +179,45 @@ print(q)
 
 ---
 
+### Set Operations
+
+#### `subtract(other: Qube) -> Qube`
+
+Return a **new** Qube containing every identifier that is in `self` but not in `other`. Neither operand is modified.
+
+```python
+a = Qube.from_ascii("""root
+└── class=od/rd
+    └── param=1
+""")
+b = Qube.from_ascii("""root
+└── class=od
+    └── param=1
+""")
+
+result = a.subtract(b)
+print(result)
+# root
+# └── class=rd
+#     └── param=1
+```
+
+Also available as the `−` operator via `__sub__`:
+
+```python
+result = a - b            # same as a.subtract(b)
+result = a - b - c        # chaining removes identifiers from b and c
+```
+
+Key properties:
+- `a.subtract(a)` → empty Qube
+- `a.subtract(Qube())` → equivalent to `a` (subtracting empty changes nothing)
+- `Qube().subtract(b)` → empty Qube
+- `a - b` and `a.subtract(b)` produce identical results
+- Both operands are left unchanged; the result is an independent Qube
+
+---
+
 ### Manipulation
 
 #### `compress() -> None`
@@ -261,6 +300,7 @@ selected = q.select({"class": [1], "param": [1, 2]}, None, None)
 | `__str__()` | Same as `to_ascii()` |
 | `__repr__()` | Returns `Qube(root_id=...)` |
 | `__len__()` | Returns `datacube_count()` — the number of leaf identifiers |
+| `__sub__(other)` | Same as `subtract(other)` — enables `a - b` syntax |
 
 ```python
 q = Qube.from_ascii("root\n├── class=od, param=1/2\n└── class=rd, param=3")
diff --git a/py_qubed/src/lib.rs b/py_qubed/src/lib.rs
index e45f4bb..8876a20 100644
--- a/py_qubed/src/lib.rs
+++ b/py_qubed/src/lib.rs
@@ -192,6 +192,18 @@ impl PyQube {
         Ok(())
     }
 
+    /// Returns a new Qube containing every identifier in `self` that is not
+    /// present in `other`.  Neither operand is modified.
+    pub fn subtract(&self, other: &Bound<'_, PyQube>) -> PyResult<PyQube> {
+        let other_ref = other.borrow();
+        Ok(PyQube { inner: self.inner.subtract(&other_ref.inner) })
+    }
+
+    /// Python operator sugar: `a - b` delegates to `a.subtract(b)`.
+    pub fn __sub__(&self, other: &Bound<'_, PyQube>) -> PyResult<PyQube> {
+        self.subtract(other)
+    }
+
     pub fn append_many(&mut self, others: &Bound<'_, PyList>) -> PyResult<()> {
         // First validate all types so type errors happen before any mutation.
         let mut validated_qubes = Vec::with_capacity(others.len());
diff --git a/py_qubed/tests/test_difference_api.py b/py_qubed/tests/test_difference_api.py
new file mode 100644
index 0000000..048daa9
--- /dev/null
+++ b/py_qubed/tests/test_difference_api.py
@@ -0,0 +1,221 @@
+from qubed import Qube
+import pytest
+
+
+# ---------------------------------------------------------------------------
+# Helpers
+# ---------------------------------------------------------------------------
+
+def _simple(class_val: str, param_val: str = "1") -> Qube:
+    return Qube.from_ascii(f"root\n└── class={class_val}\n    └── param={param_val}\n")
+
+
+# ---------------------------------------------------------------------------
+# Basic identity / empty cases
+# ---------------------------------------------------------------------------
+
+def test_subtract_identical_returns_empty() -> None:
+    """A − A should produce an empty Qube (no dimension coordinates)."""
+    a = _simple("od")
+    b = _simple("od")
+    result = a.subtract(b)
+    assert result.all_unique_dim_coords() == {}
+
+
+def test_subtract_empty_b_returns_a() -> None:
+    """A − ∅ = A."""
+    a = _simple("od")
+    b = Qube()
+    result = a.subtract(b)
+    assert len(result) == len(a)
+    assert "class=od" in result.to_ascii()
+
+
+def test_subtract_empty_a_returns_empty() -> None:
+    """∅ − B = ∅ (no dimension coordinates in the result)."""
+    a = Qube()
+    b = _simple("od")
+    result = a.subtract(b)
+    assert result.all_unique_dim_coords() == {}
+
+
+# ---------------------------------------------------------------------------
+# Disjoint sets
+# ---------------------------------------------------------------------------
+
+def test_subtract_disjoint_preserves_all_of_a() -> None:
+    """When A and B share no identifiers, A − B = A."""
+    a = _simple("od")
+    b = _simple("rd")
+    result = a.subtract(b)
+    coords = result.all_unique_dim_coords()
+    assert coords["class"] == ["od"]
+
+
+def test_subtract_disjoint_different_dimension_depth() -> None:
+    """Qubes with completely different dimension schemas are disjoint."""
+    a = Qube.from_ascii("root\n└── class=od\n    └── param=1\n")
+    b = Qube.from_ascii("root\n└── type=fc\n    └── step=0\n")
+    result = a.subtract(b)
+    ascii_out = result.to_ascii()
+    assert "class=od" in ascii_out
+    assert "param=1" in ascii_out
+
+
+# ---------------------------------------------------------------------------
+# Partial overlap — coordinate level
+# ---------------------------------------------------------------------------
+
+def test_subtract_removes_overlapping_coord_leaves_rest() -> None:
+    """A has class=od/rd; B has class=od. Result should contain only class=rd."""
+    a = Qube.from_ascii("root\n└── class=od/rd\n    └── param=1\n")
+    b = _simple("od")
+    result = a.subtract(b)
+    coords = result.all_unique_dim_coords()
+    assert "od" not in coords.get("class", [])
+    assert "rd" in coords["class"]
+
+
+def test_subtract_removes_subset_of_param_values() -> None:
+    """A has param=1/2/3; B covers param=2. Result should contain param=1/3."""
+    a = Qube.from_ascii("root\n└── class=od\n    └── param=1/2/3\n")
+    b = Qube.from_ascii("root\n└── class=od\n    └── param=2\n")
+    result = a.subtract(b)
+    coords = result.all_unique_dim_coords()
+    param_values = set(coords["param"])
+    assert "2" not in param_values
+    assert {"1", "3"}.issubset(param_values)
+
+
+def test_subtract_all_params_removes_branch() -> None:
+    """When B covers all of A's param values, no dimension coordinates remain."""
+    a = Qube.from_ascii("root\n└── class=od\n    └── param=1/2\n")
+    b = Qube.from_ascii("root\n└── class=od\n    └── param=1/2\n")
+    result = a.subtract(b)
+    assert result.all_unique_dim_coords() == {}
+
+
+# ---------------------------------------------------------------------------
+# Multi-branch / nested structures
+# ---------------------------------------------------------------------------
+
+def test_subtract_multi_class_partial_remove() -> None:
+    """A has class=od/rd/xd; B covers class=rd. Result: class=od/xd."""
+    a = Qube.from_ascii("root\n└── class=od/rd/xd\n    └── param=1\n")
+    b = _simple("rd")
+    result = a.subtract(b)
+    coords = result.all_unique_dim_coords()
+    assert "rd" not in coords.get("class", [])
+    assert set(coords["class"]) == {"od", "xd"}
+
+
+def test_subtract_deep_nested_partial_overlap() -> None:
+    """Subtraction propagates correctly through multi-level trees."""
+    a = Qube.from_ascii("""root
+└── class=od
+    ├── expver=0001
+    │   └── param=1/2
+    └── expver=0002
+        └── param=1/2
+""")
+    b = Qube.from_ascii("""root
+└── class=od
+    └── expver=0001
+        └── param=1
+""")
+    result = a.subtract(b)
+    ascii_out = result.to_ascii()
+    # param=1 under expver=0001 should be gone
+    assert "expver=0001" not in ascii_out or "param=2" in ascii_out
+    # expver=0002 branch should survive intact
+    assert "expver=0002" in ascii_out
+
+
+def test_subtract_b_covers_entire_branch() -> None:
+    """When B fully covers one class branch, that branch disappears from result."""
+    a = Qube.from_ascii("root\n├── class=od\n│   └── param=1\n└── class=rd\n    └── param=2\n")
+    b = _simple("od")
+    result = a.subtract(b)
+    ascii_out = result.to_ascii()
+    assert "class=od" not in ascii_out
+    assert "class=rd" in ascii_out
+
+
+# ---------------------------------------------------------------------------
+# Operator sugar  `a - b`
+# ---------------------------------------------------------------------------
+
+def test_dunder_sub_is_equivalent_to_subtract() -> None:
+    """`a - b` and `a.subtract(b)` must produce the same ASCII tree."""
+    a = Qube.from_ascii("root\n└── class=od/rd\n    └── param=1\n")
+    b = _simple("od")
+    via_method = a.subtract(b)
+    via_operator = a - b
+    assert via_method.to_ascii() == via_operator.to_ascii()
+
+
+def test_dunder_sub_chaining() -> None:
+    """`a - b - c` should remove identifiers from both b and c."""
+    a = Qube.from_ascii("root\n└── class=od/rd/xd\n    └── param=1\n")
+    b = _simple("od")
+    c = _simple("xd")
+    result = a - b - c
+    coords = result.all_unique_dim_coords()
+    assert set(coords["class"]) == {"rd"}
+
+
+# ---------------------------------------------------------------------------
+# Non-mutation guarantee
+# ---------------------------------------------------------------------------
+
+def test_subtract_does_not_mutate_a() -> None:
+    """a.subtract(b) must leave `a` unchanged."""
+    a = _simple("od")
+    b = _simple("od")
+    before = a.to_ascii()
+    _ = a.subtract(b)
+    assert a.to_ascii() == before
+
+
+def test_subtract_does_not_mutate_b() -> None:
+    """a.subtract(b) must leave `b` unchanged."""
+    a = _simple("od")
+    b = _simple("od")
+    before = b.to_ascii()
+    _ = a.subtract(b)
+    assert b.to_ascii() == before
+
+
+# ---------------------------------------------------------------------------
+# Result is a fresh independent Qube
+# ---------------------------------------------------------------------------
+
+def test_subtract_result_is_independent() -> None:
+    """Mutating the result via compress must not affect the original."""
+    a = Qube.from_ascii("root\n└── class=od/rd\n    └── param=1\n")
+    b = _simple("od")
+    result = a.subtract(b)
+    result.compress()  # should not raise; result is a valid, mutable Qube
+    assert "class=od" in a.to_ascii()  # original untouched
+
+
+# ---------------------------------------------------------------------------
+# len() on result
+# ---------------------------------------------------------------------------
+
+def test_subtract_len_reflects_removed_identifiers() -> None:
+    """len() counts leaf-path nodes; subtracting one complete branch reduces the count.
+
+    Three distinct branches each with a unique param value means three leaf-path
+    nodes (they cannot be merged because their sub-trees differ).  Removing one
+    whole branch drops len() from 3 → 2.
+    """
+    a = Qube.from_ascii(
+        "root\n"
+        "├── class=od\n│   └── param=10\n"
+        "├── class=rd\n│   └── param=20\n"
+        "└── class=xd\n    └── param=30\n"
+    )
+    b = Qube.from_ascii("root\n└── class=od\n    └── param=10\n")
+    result = a.subtract(b)
+    assert len(result) == 2
diff --git a/qubed/src/compress.rs b/qubed/src/compress.rs
index 4b35d52..0d5d7c7 100644
--- a/qubed/src/compress.rs
+++ b/qubed/src/compress.rs
@@ -48,6 +48,10 @@ impl Qube {
         for kids in parent.children_mut().values_mut() {
             kids.retain(|id| keep.contains(id));
         }
+        // Remove dimension entries whose child list became empty after pruning.
+        // Without this, BTreeMap::is_empty() would return false even when no
+        // real children remain, causing Qube::is_empty() to mis-report.
+        parent.children_mut().retain(|_, kids| !kids.is_empty());
     }
 
     /// Invalidates the cached structural hash of a node.
diff --git a/qubed/src/difference.rs b/qubed/src/difference.rs
new file mode 100644
index 0000000..393cec2
--- /dev/null
+++ b/qubed/src/difference.rs
@@ -0,0 +1,251 @@
+use std::collections::HashMap;
+
+use crate::coordinates::Coordinates;
+use crate::{NodeIdx, Qube};
+
+impl Qube {
+    /// Computes the set difference A − B.
+    ///
+    /// Returns a new [`Qube`] containing every identifier present in `self` (A)
+    /// that is **not** present in `other` (B).  Neither operand is consumed or
+    /// modified.
+    ///
+    /// # Semantics
+    ///
+    /// The operation is applied recursively on the compressed tree structure.
+    /// At each dimension level the coordinates of every A-node are split by
+    /// intersecting with each overlapping B-node:
+    ///
+    /// - **Only-A values** (`A_coords − B_coords`): kept unchanged, together
+    ///   with their full A subtree.
+    /// - **Intersection values** (`A_coords ∩ B_coords`): three sub-cases:
+    ///   - *B is a leaf* (no deeper dimensions): B "covers" the entire
+    ///     coordinate range.  The intersection values and all of A's children
+    ///     underneath them are removed.
+    ///   - *A is a leaf, B has children*: the schemas differ — A's paths end
+    ///     here while B's continue deeper.  A's datacubes are considered
+    ///     distinct from B's deeper datacubes, so A's leaf is kept unchanged.
+    ///   - *Both have children*: the operation recurses into the A and B
+    ///     subtrees for the intersection coordinate range.  Any A paths that
+    ///     survive the recursive subtraction are kept; those fully covered by B
+    ///     are discarded.
+    /// - **Only-B values**: not in A, so irrelevant and ignored.
+    ///
+    /// The result is automatically compressed before being returned.
+    pub fn subtract(&self, other: &Qube) -> Qube {
+        // Seed the result with a full copy of A.
+        let mut result = Qube::new();
+        let result_root = result.root();
+        let self_root = self.root();
+        result.copy_subtree(self, self_root, result_root);
+
+        // Fast paths: trivially empty inputs.
+        if other.is_empty() || result.is_empty() {
+            result.compress();
+            return result;
+        }
+
+        // Ensure result is compressed so sibling coordinate sets are
+        // non-overlapping — correctness of the pending-list loop relies on it.
+        result.compress();
+
+        let result_root = result.root();
+        let other_root = other.root();
+        result.node_subtract(other, result_root, other_root);
+
+        result.compress();
+        result
+    }
+
+    // ------------------------------------------------------------------
+    // Internal recursive helpers
+    // ------------------------------------------------------------------
+
+    /// Recursively subtracts the B-subtree rooted at `other_id` from the
+    /// A-subtree (in `self`) rooted at `self_id`.
+    fn node_subtract(&mut self, other: &Qube, self_id: NodeIdx, other_id: NodeIdx) {
+        let self_children = self.node_ref(self_id).unwrap().children().clone();
+        let other_children = other.node_ref(other_id).unwrap().children().clone();
+
+        // Build a map  dim_name_str → A's child NodeIdx list.
+        // We key by string because the two Qubes may have independent string
+        // interner tables (different MiniSpur integers for the same name).
+        let mut self_dim_kids: HashMap<String, Vec<NodeIdx>> = HashMap::new();
+        for (dim, kids) in &self_children {
+            if let Some(s) = self.dimension_str(dim) {
+                self_dim_kids.entry(s.to_owned()).or_default().extend(kids);
+            }
+        }
+
+        // For every dimension that B has at this level, subtract its children
+        // from A's matching children.
+        for (other_dim, other_kids) in &other_children {
+            let other_dim_str = match other.dimension_str(other_dim) {
+                Some(s) => s.to_owned(),
+                None => continue,
+            };
+
+            let self_kids = match self_dim_kids.get(&other_dim_str) {
+                Some(kids) => kids.clone(),
+                None => continue, // A has no node for this dimension here.
+            };
+
+            let other_kids_vec: Vec<NodeIdx> = other_kids.iter().copied().collect();
+
+            self.dimension_subtract(other, &other_dim_str, self_id, self_kids, other_kids_vec);
+        }
+    }
+
+    /// Applies the difference operation within one dimension group.
+    ///
+    /// Uses a *pending list*: the outer loop iterates over B-kids; the inner
+    /// loop processes all A-kid fragments that still need to be checked against
+    /// the current B-kid.  When an A-node is split into an intersection part
+    /// (handled immediately, kept as the original node to preserve order) and
+    /// an only-self part (new node), the only-self fragment is pushed into the
+    /// *next* pending list so subsequent B-kids also get a chance to subtract
+    /// from it.
+    ///
+    /// * `parent_id` — the parent of all nodes in `self_kids`
+    /// * `dim_str`   — the shared dimension name
+    fn dimension_subtract(
+        &mut self,
+        other: &Qube,
+        dim_str: &str,
+        parent_id: NodeIdx,
+        self_kids: Vec<NodeIdx>,
+        other_kids: Vec<NodeIdx>,
+    ) {
+        let mut pending: Vec<NodeIdx> = self_kids;
+
+        for &other_kid_id in &other_kids {
+            let other_coords = other.node_ref(other_kid_id).unwrap().coords().clone();
+            if other_coords.is_empty() {
+                continue;
+            }
+
+            let mut next_pending: Vec<NodeIdx> = Vec::new();
+
+            for self_kid_id in std::mem::take(&mut pending) {
+                let self_coords = self.node_ref(self_kid_id).unwrap().coords().clone();
+                if self_coords.is_empty() {
+                    // Already exhausted by a prior B-kid; don't carry forward.
+                    continue;
+                }
+
+                let res = self_coords.intersect(&other_coords);
+                let only_self = res.only_a;
+                let intersection = res.intersection;
+
+                if intersection.is_empty() {
+                    // No overlap with this B-kid; carry A-node to the next round.
+                    next_pending.push(self_kid_id);
+                    continue;
+                }
+
+                let self_has_children = !self.node_ref(self_kid_id).unwrap().children().is_empty();
+                let other_has_children =
+                    !other.node_ref(other_kid_id).unwrap().children().is_empty();
+
+                match (self_has_children, other_has_children) {
+                    // ── B is a leaf ────────────────────────────────────────────
+                    // B covers the entire intersection range at this level.
+                    // Trim A's coordinate set to only_self; A's children (if any)
+                    // remain associated with the surviving only-self values.
+                    (_, false) => {
+                        if only_self.is_empty() {
+                            // Must use Coordinates::Empty so prune_empty_nodes_recursively
+                            // will remove this node during compress().
+                            let node = self.node_mut(self_kid_id).unwrap();
+                            *node.coords_mut() = Coordinates::Empty;
+                            // Invalidate cached structural hash so compress() sees the change.
+                            self.invalidate_ancestors(self_kid_id);
+                            // Node exhausted; don't carry forward.
+                        } else {
+                            let node = self.node_mut(self_kid_id).unwrap();
+                            *node.coords_mut() = only_self;
+                            self.invalidate_ancestors(self_kid_id);
+                            // Remaining coords may overlap with a later B-kid.
+                            next_pending.push(self_kid_id);
+                        }
+                    }
+
+                    // ── A is a leaf, B has deeper children ────────────────────
+                    // The schema depths differ: A's datacubes end here while
+                    // B's continue further. They are considered distinct
+                    // identifiers, so A's leaf is left untouched.
+                    (false, true) => {
+                        // No-op: leave self_kid_id unchanged, carry to next B-kid.
+                        next_pending.push(self_kid_id);
+                    }
+
+                    // ── Both have children: recurse ───────────────────────────
+                    (true, true) => {
+                        if only_self.is_empty() {
+                            // Intersection == entire A-node's coords.
+                            // Subtract B directly from self_kid_id (no new node needed).
+                            // self_kid_id's coords are a subset of this B-kid's coords, so no
+                            // later (non-overlapping) B-kid can intersect further — don't carry.
+                            self.node_subtract(other, self_kid_id, other_kid_id);
+                            if !self.has_leaf_content(self_kid_id) {
+                                let node = self.node_mut(self_kid_id).unwrap();
+                                *node.coords_mut() = Coordinates::Empty;
+                                self.invalidate_ancestors(self_kid_id);
+                            }
+                        } else {
+                            // Keep the original node as the intersection part.
+                            // This preserves insertion order so ASCII output is stable.
+                            {
+                                let node = self.node_mut(self_kid_id).unwrap();
+                                *node.coords_mut() = intersection;
+                                // coords changed — invalidate cached hash for self and ancestors
+                                self.invalidate_ancestors(self_kid_id);
+                            }
+
+                            // Create a sibling node for the only-self values,
+                            // seeded with a deep copy of A's current subtree.
+                            let only_self_node = self
+                                .get_or_create_child(dim_str, parent_id, Some(only_self))
+                                .unwrap();
+                            self.copy_branch(self_kid_id, only_self_node);
+
+                            // Recursively subtract B's subtree from the intersection node.
+                            self.node_subtract(other, self_kid_id, other_kid_id);
+                            if !self.has_leaf_content(self_kid_id) {
+                                let node = self.node_mut(self_kid_id).unwrap();
+                                *node.coords_mut() = Coordinates::Empty;
+                                self.invalidate_ancestors(self_kid_id);
+                            }
+
+                            // The only-self fragment may overlap with later B-kids.
+                            next_pending.push(only_self_node);
+                        }
+                    }
+                }
+            }
+
+            pending = next_pending;
+        }
+    }
+
+    /// Returns `true` if the subtree rooted at `node_id` contains at least one
+    /// leaf node with non-empty coordinates (i.e., at least one identifier).
+    fn has_leaf_content(&self, node_id: NodeIdx) -> bool {
+        let node = self.node_ref(node_id).unwrap();
+        if node.children().is_empty() {
+            return !node.coords().is_empty();
+        }
+        node.children().values().flatten().any(|&child_id| self.has_leaf_content(child_id))
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Trait impl: `&a - &b`
+// ---------------------------------------------------------------------------
+
+impl std::ops::Sub for &Qube {
+    type Output = Qube;
+    fn sub(self, rhs: Self) -> Qube {
+        self.subtract(rhs)
+    }
+}
diff --git a/qubed/src/lib.rs b/qubed/src/lib.rs
index 79915cc..fec21bc 100644
--- a/qubed/src/lib.rs
+++ b/qubed/src/lib.rs
@@ -1,6 +1,7 @@
 mod compress;
 mod coordinates;
 pub mod datacube;
+mod difference;
 mod merge;
 mod qube;
 pub mod select;
diff --git a/qubed/tests/test_difference.rs b/qubed/tests/test_difference.rs
new file mode 100644
index 0000000..fb9f36d
--- /dev/null
+++ b/qubed/tests/test_difference.rs
@@ -0,0 +1,647 @@
+use qubed::Qube;
+
+// ---------------------------------------------------------------------------
+// Helper: build a Qube from ASCII and assert it equals expected ASCII
+// ---------------------------------------------------------------------------
+fn assert_ascii_eq(result: &Qube, expected: &str, msg: &str) {
+    let expected_qube = Qube::from_ascii(expected).unwrap();
+    assert_eq!(result.to_ascii(), expected_qube.to_ascii(), "{}", msg);
+}
+
+fn assert_empty(q: &Qube, msg: &str) {
+    assert!(q.is_empty(), "{msg}: expected empty Qube, got:\n{}", q.to_ascii());
+}
+
+// ---------------------------------------------------------------------------
+// Basic: completely disjoint qubes → result equals A
+// ---------------------------------------------------------------------------
+#[test]
+fn subtract_disjoint_returns_self() {
+    let a = Qube::from_ascii(
+        r#"root
+└── class=1
+    └── param=1/2"#,
+    )
+    .unwrap();
+
+    let b = Qube::from_ascii(
+        r#"root
+└── class=2
+    └── param=1/2"#,
+    )
+    .unwrap();
+
+    let result = a.subtract(&b);
+
+    assert_ascii_eq(
+        &result,
+        r#"root
+└── class=1
+    └── param=1/2"#,
+        "disjoint A–B should equal A",
+    );
+}
+
+// ---------------------------------------------------------------------------
+// Basic: A − A = empty
+// ---------------------------------------------------------------------------
+#[test]
+fn subtract_identical_qubes_returns_empty() {
+    let a = Qube::from_ascii(
+        r#"root
+└── class=1
+    ├── expver=0001
+    │   └── param=1/2
+    └── expver=0002
+        └── param=1/2"#,
+    )
+    .unwrap();
+
+    let result = a.subtract(&a);
+    assert_empty(&result, "A − A");
+}
+
+// ---------------------------------------------------------------------------
+// Basic: A − empty = A
+// ---------------------------------------------------------------------------
+#[test]
+fn subtract_empty_other_returns_self() {
+    let a = Qube::from_ascii(
+        r#"root
+└── class=1
+    └── param=1/2"#,
+    )
+    .unwrap();
+
+    let b = Qube::new();
+    let result = a.subtract(&b);
+
+    assert_ascii_eq(
+        &result,
+        r#"root
+└── class=1
+    └── param=1/2"#,
+        "A − empty should equal A",
+    );
+}
+
+// ---------------------------------------------------------------------------
+// Basic: empty − B = empty
+// ---------------------------------------------------------------------------
+#[test]
+fn subtract_from_empty_returns_empty() {
+    let a = Qube::new();
+    let b = Qube::from_ascii(
+        r#"root
+└── class=1
+    └── param=1/2"#,
+    )
+    .unwrap();
+
+    let result = a.subtract(&b);
+    assert_empty(&result, "empty − B");
+}
+
+// ---------------------------------------------------------------------------
+// Coordinate-level removal: A has class=1/2/3, B covers class=1/2 (leaf)
+// → result: class=3 (with original subtree)
+// ---------------------------------------------------------------------------
+#[test]
+fn subtract_leaf_coordinate_removal() {
+    let a = Qube::from_ascii(
+        r#"root
+└── class=1/2/3
+    └── param=1/2"#,
+    )
+    .unwrap();
+
+    let b = Qube::from_ascii(
+        r#"root
+└── class=1/2
+    └── param=1/2"#,
+    )
+    .unwrap();
+
+    let result = a.subtract(&b);
+
+    assert_ascii_eq(
+        &result,
+        r#"root
+└── class=3
+    └── param=1/2"#,
+        "subtract should remove class=1/2 and keep class=3",
+    );
+}
+
+// ---------------------------------------------------------------------------
+// Partial leaf subtraction: only some values removed at the leaf level
+// ---------------------------------------------------------------------------
+#[test]
+fn subtract_partial_leaf_coord_removal() {
+    let a = Qube::from_ascii(
+        r#"root
+└── class=1
+    └── param=1/2/3/4"#,
+    )
+    .unwrap();
+
+    let b = Qube::from_ascii(
+        r#"root
+└── class=1
+    └── param=2/3"#,
+    )
+    .unwrap();
+
+    let result = a.subtract(&b);
+
+    assert_ascii_eq(
+        &result,
+        r#"root
+└── class=1
+    └── param=1/4"#,
+        "should remove param=2/3, keep param=1/4",
+    );
+}
+
+// ---------------------------------------------------------------------------
+// Multi-level: remove a branch at an intermediate level
+// ---------------------------------------------------------------------------
+#[test]
+fn subtract_intermediate_branch_removal() {
+    let a = Qube::from_ascii(
+        r#"root
+├── class=1
+│   ├── expver=0001
+│   │   └── param=1/2
+│   └── expver=0002
+│       └── param=1/2
+└── class=2
+    └── expver=0001
+        └── param=1/2"#,
+    )
+    .unwrap();
+
+    // Remove the entire class=1/expver=0001 branch
+    let b = Qube::from_ascii(
+        r#"root
+└── class=1
+    └── expver=0001
+        └── param=1/2"#,
+    )
+    .unwrap();
+
+    let result = a.subtract(&b);
+
+    let expected = Qube::from_ascii(
+        r#"root
+├── class=1
+│   └── expver=0002
+│       └── param=1/2
+└── class=2
+    └── expver=0001
+        └── param=1/2"#,
+    )
+    .unwrap();
+
+    assert_eq!(
+        result.to_ascii(),
+        expected.to_ascii(),
+        "class=1/expver=0001 should be removed; rest preserved"
+    );
+}
+
+// ---------------------------------------------------------------------------
+// B covers a superset of A → result is empty
+// ---------------------------------------------------------------------------
+#[test]
+fn subtract_b_is_superset_returns_empty() {
+    let a = Qube::from_ascii(
+        r#"root
+└── class=1/2
+    └── param=1/2"#,
+    )
+    .unwrap();
+
+    let b = Qube::from_ascii(
+        r#"root
+└── class=1/2/3
+    └── param=1/2/3"#,
+    )
+    .unwrap();
+
+    let result = a.subtract(&b);
+    assert_empty(&result, "when B ⊇ A, result should be empty");
+}
+
+// ---------------------------------------------------------------------------
+// Asymmetry: A − B ≠ B − A in general
+// ---------------------------------------------------------------------------
+#[test]
+fn subtract_is_asymmetric() {
+    let a = Qube::from_ascii(
+        r#"root
+└── class=1/2/3
+    └── param=1/2"#,
+    )
+    .unwrap();
+
+    let b = Qube::from_ascii(
+        r#"root
+└── class=2/3/4
+    └── param=1/2"#,
+    )
+    .unwrap();
+
+    let a_minus_b = a.subtract(&b);
+    let b_minus_a = b.subtract(&a);
+
+    // A − B = class=1 / param=1/2
+    assert_ascii_eq(
+        &a_minus_b,
+        r#"root
+└── class=1
+    └── param=1/2"#,
+        "A − B",
+    );
+
+    // B − A = class=4 / param=1/2
+    assert_ascii_eq(
+        &b_minus_a,
+        r#"root
+└── class=4
+    └── param=1/2"#,
+        "B − A",
+    );
+
+    assert_ne!(a_minus_b.to_ascii(), b_minus_a.to_ascii(), "A − B should differ from B − A");
+}
+
+// ---------------------------------------------------------------------------
+// String coordinates
+// ---------------------------------------------------------------------------
+#[test]
+fn subtract_string_coordinates() {
+    let a = Qube::from_ascii(
+        r#"root
+└── class=od/rd/xd
+    └── param=1/2"#,
+    )
+    .unwrap();
+
+    let b = Qube::from_ascii(
+        r#"root
+└── class=od/rd
+    └── param=1/2"#,
+    )
+    .unwrap();
+
+    let result = a.subtract(&b);
+
+    assert_ascii_eq(
+        &result,
+        r#"root
+└── class=xd
+    └── param=1/2"#,
+        "subtract with string coordinates",
+    );
+}
+
+// ---------------------------------------------------------------------------
+// Deep tree: only a specific leaf value is removed
+// ---------------------------------------------------------------------------
+#[test]
+fn subtract_single_deep_leaf_value() {
+    let a = Qube::from_ascii(
+        r#"root
+└── class=1/2
+    ├── expver=0001
+    │   └── param=1/2
+    └── expver=0002
+        └── param=1/2"#,
+    )
+    .unwrap();
+
+    // Remove only class=1, expver=0001, param=1
+    let b = Qube::from_ascii(
+        r#"root
+└── class=1
+    └── expver=0001
+        └── param=1"#,
+    )
+    .unwrap();
+
+    let result = a.subtract(&b);
+
+    // The expected result: class=1/expver=0001/param=1 removed
+    let expected = Qube::from_ascii(
+        r#"root
+├── class=1
+│   ├── expver=0001
+│   │   └── param=2
+│   └── expver=0002
+│       └── param=1/2
+└── class=2
+    └── expver=0001/0002
+        └── param=1/2"#,
+    )
+    .unwrap();
+
+    assert_eq!(
+        result.to_ascii(),
+        expected.to_ascii(),
+        "only class=1/expver=0001/param=1 should be removed"
+    );
+}
+
+// ---------------------------------------------------------------------------
+// B leaf covers whole A subtree → all of A under intersection is removed
+// ---------------------------------------------------------------------------
+#[test]
+fn subtract_b_leaf_removes_whole_a_subtree() {
+    let a = Qube::from_ascii(
+        r#"root
+└── class=1/2
+    ├── expver=0001
+    │   └── param=1/2/3
+    └── expver=0002
+        └── param=4/5"#,
+    )
+    .unwrap();
+
+    // B covers class=1 as a leaf → everything under class=1 in A is removed
+    let b = Qube::from_ascii(
+        r#"root
+└── class=1
+    └── expver=0001
+        └── param=1/2/3"#,
+    )
+    .unwrap();
+
+    let result = a.subtract(&b);
+
+    // Expected: class=1 loses expver=0001/param=1/2/3;
+    //           class=1/expver=0002 and class=2 branches survive.
+    let expected = Qube::from_ascii(
+        r#"root
+├── class=1
+│   └── expver=0002
+│       └── param=4/5
+└── class=2
+    ├── expver=0001
+    │   └── param=1/2/3
+    └── expver=0002
+        └── param=4/5"#,
+    )
+    .unwrap();
+
+    assert_eq!(
+        result.to_ascii(),
+        expected.to_ascii(),
+        "class=1/expver=0001 subtree should be removed"
+    );
+}
+
+// ---------------------------------------------------------------------------
+// Self-referential: A - B does not modify A or B
+// ---------------------------------------------------------------------------
+#[test]
+fn subtract_does_not_modify_inputs() {
+    let a = Qube::from_ascii(
+        r#"root
+└── class=1/2
+    └── param=1/2"#,
+    )
+    .unwrap();
+
+    let b = Qube::from_ascii(
+        r#"root
+└── class=1
+    └── param=1"#,
+    )
+    .unwrap();
+
+    let a_ascii_before = a.to_ascii();
+    let b_ascii_before = b.to_ascii();
+
+    let _result = a.subtract(&b);
+
+    assert_eq!(a.to_ascii(), a_ascii_before, "A should not be modified by subtract");
+    assert_eq!(b.to_ascii(), b_ascii_before, "B should not be modified by subtract");
+}
+
+// ---------------------------------------------------------------------------
+// Operator sugar: &a - &b should equal a.subtract(&b)
+// ---------------------------------------------------------------------------
+#[test]
+fn subtract_operator_equals_method() {
+    let a = Qube::from_ascii(
+        r#"root
+└── class=1/2/3
+    └── param=1/2"#,
+    )
+    .unwrap();
+
+    let b = Qube::from_ascii(
+        r#"root
+└── class=2
+    └── param=1/2"#,
+    )
+    .unwrap();
+
+    let via_method = a.subtract(&b);
+    let via_operator = &a - &b;
+
+    assert_eq!(
+        via_method.to_ascii(),
+        via_operator.to_ascii(),
+        "operator and method should produce identical results"
+    );
+}
+
+// ---------------------------------------------------------------------------
+// A - B when B has dimensions A doesn't: A is unchanged
+// (Different schema depths — B can only remove what its paths match in A)
+// ---------------------------------------------------------------------------
+#[test]
+fn subtract_b_deeper_schema_leaves_a_unchanged() {
+    // A is a leaf at the class level (no further dimensions)
+    let a = Qube::from_ascii(
+        r#"root
+└── class=1/2"#,
+    )
+    .unwrap();
+
+    // B goes deeper: class → expver
+    let b = Qube::from_ascii(
+        r#"root
+└── class=1
+    └── expver=0001"#,
+    )
+    .unwrap();
+
+    let result = a.subtract(&b);
+
+    // A's paths ({class=1} and {class=2}) don't match B's deeper paths
+    // ({class=1, expver=0001}), so A should be unchanged.
+    assert_ascii_eq(
+        &result,
+        r#"root
+└── class=1/2"#,
+        "shallower A should not be affected by deeper B",
+    );
+}
+
+// ---------------------------------------------------------------------------
+// Complex: multiple B branches that each remove different parts of A
+// ---------------------------------------------------------------------------
+#[test]
+fn subtract_multiple_b_branches() {
+    let a = Qube::from_ascii(
+        r#"root
+└── class=1/2/3/4
+    └── expver=0001/0002
+        └── param=1/2/3"#,
+    )
+    .unwrap();
+
+    // B removes:
+    //   - class=1/2 / expver=0001 / param=1/2
+    //   - class=3/4 / expver=0002 / param=2/3
+    let b = Qube::from_ascii(
+        r#"root
+├── class=1/2
+│   └── expver=0001
+│       └── param=1/2
+└── class=3/4
+    └── expver=0002
+        └── param=2/3"#,
+    )
+    .unwrap();
+
+    let result = a.subtract(&b);
+
+    // Verify the result is a valid non-empty Qube and does not contain
+    // any of the removed identifiers.
+    let ascii = result.to_ascii();
+
+    // Removed: class=1/2, expver=0001, param=1/2
+    // Surviving from that branch: class=1/2/expver=0001/param=3
+    //                              class=1/2/expver=0002/param=1/2/3
+    // Removed: class=3/4, expver=0002, param=2/3
+    // Surviving from that branch: class=3/4/expver=0001/param=1/2/3
+    //                              class=3/4/expver=0002/param=1
+
+    assert!(!result.is_empty(), "result should not be empty");
+
+    // Spot-check: the combinations that WERE in B should be gone.
+    // datacube_count() counts compressed leaf-path count (not coord cross-product).
+    // After subtraction and compression, 4 distinct leaf paths remain:
+    //   class=1/2 / expver=0001 / param=3
+    //   class=1/2 / expver=0002 / param=1/2/3
+    //   class=3/4 / expver=0001 / param=1/2/3
+    //   class=3/4 / expver=0002 / param=1
+    assert_eq!(result.datacube_count(), 4, "datacube count after subtract: {}", ascii);
+}
+
+// ---------------------------------------------------------------------------
+// Result is a valid compressed Qube: ASCII round-trip is stable
+// ---------------------------------------------------------------------------
+#[test]
+fn subtract_result_ascii_roundtrip_is_stable() {
+    let a = Qube::from_ascii(
+        r#"root
+└── class=1/2
+    ├── expver=0001
+    │   └── param=1/2
+    └── expver=0002
+        └── param=1/2"#,
+    )
+    .unwrap();
+
+    let b = Qube::from_ascii(
+        r#"root
+└── class=1
+    └── expver=0001
+        └── param=1/2"#,
+    )
+    .unwrap();
+
+    let result = a.subtract(&b);
+
+    // Parsing the ASCII back should yield the same ASCII output.
+    let roundtrip = Qube::from_ascii(&result.to_ascii()).unwrap();
+    assert_eq!(
+        result.to_ascii(),
+        roundtrip.to_ascii(),
+        "result ASCII should be stable on re-parse"
+    );
+}
+
+// ---------------------------------------------------------------------------
+// Edge: single-node trees
+// ---------------------------------------------------------------------------
+#[test]
+fn subtract_single_value_trees() {
+    let a = Qube::from_ascii(
+        r#"root
+└── class=42"#,
+    )
+    .unwrap();
+    let b = Qube::from_ascii(
+        r#"root
+└── class=42"#,
+    )
+    .unwrap();
+    assert_empty(&a.subtract(&b), "single identical leaf trees should give empty");
+
+    let c = Qube::from_ascii(
+        r#"root
+└── class=99"#,
+    )
+    .unwrap();
+    assert_ascii_eq(
+        &a.subtract(&c),
+        r#"root
+└── class=42"#,
+        "subtract of disjoint single-value trees",
+    );
+}
+
+// ---------------------------------------------------------------------------
+// Regression: chained subtractions
+// ---------------------------------------------------------------------------
+#[test]
+fn subtract_chained() {
+    let a = Qube::from_ascii(
+        r#"root
+└── class=1/2/3/4/5
+    └── param=1"#,
+    )
+    .unwrap();
+
+    let b = Qube::from_ascii(
+        r#"root
+└── class=1
+    └── param=1"#,
+    )
+    .unwrap();
+
+    let c = Qube::from_ascii(
+        r#"root
+└── class=5
+    └── param=1"#,
+    )
+    .unwrap();
+
+    // (A - B) - C should equal A - (B ∪ C)
+    let step1 = a.subtract(&b);
+    let step2 = step1.subtract(&c);
+
+    assert_ascii_eq(
+        &step2,
+        r#"root
+└── class=2/3/4
+    └── param=1"#,
+        "chained subtraction",
+    );
+}