diff --git a/docs/core_concepts/outputs_and_inspection.md b/docs/core_concepts/outputs_and_inspection.md
index b8b4d80..87c2889 100644
--- a/docs/core_concepts/outputs_and_inspection.md
+++ b/docs/core_concepts/outputs_and_inspection.md
@@ -94,7 +94,9 @@ Two parameters control the physical layout of the stacked output.
 
 `output_format`
 : One of `"dense"`, `"sparse"`, or `"auto"`. `"auto"` picks sparse when it saves memory (for
-example wide one-hot blocks) and dense otherwise. Default `"dense"`.
+example wide one-hot blocks) and dense otherwise. The choice is made once at `fit`, from the
+training output, and stored in `output_format_`, so every `transform` returns the same container,
+even for a single row. Default `"dense"`.
 
 `dtype`
 : The floating-point precision of the output, for example `numpy.float32` to halve memory.
diff --git a/pretab/compose/output.py b/pretab/compose/output.py
index 5a0223c..b670f48 100644
--- a/pretab/compose/output.py
+++ b/pretab/compose/output.py
@@ -22,6 +22,7 @@
     "compute_output_report",
     "format_output",
     "resolve_embedding_dimensions",
+    "resolve_output_format",
     "to_dataframe_output",
     "validate_embedding_request",
 ]
@@ -62,6 +63,26 @@ def validate_embedding_request(embeddings, *, expected: bool, output_kind: str =
         raise IncompatibleParamsError(_EMBEDDINGS_DICT_ONLY)
 
 
+def _density(array) -> float:
+    """Fraction of non-zero entries in a dense or sparse array (0.0 when empty)."""
+    size = int(np.prod(array.shape, dtype=np.int64))
+    nonzero = int(array.count_nonzero()) if sp.issparse(array) else int(np.count_nonzero(array))
+    return float(nonzero) / size if size else 0.0
+
+
+def resolve_output_format(array, output_format, *, threshold=_SPARSE_AUTO_THRESHOLD) -> str:
+    """Resolve ``output_format`` to ``"dense"`` or ``"sparse"`` for ``array``.
+
+    ``"auto"`` picks ``"sparse"`` when the density of ``array`` is below
+    ``threshold``. The Preprocessor calls this once at ``fit`` on the training
+    output, so the container is fixed per fit rather than chosen per batch (a
+    sparse single row would break a dense-only downstream estimator).
+    """
+    if output_format == "auto":
+        return "sparse" if _density(array) < threshold else "dense"
+    return output_format
+
+
 def compute_output_report(array, output_format, *, threshold=_SPARSE_AUTO_THRESHOLD):
     """Resolve the concrete output format and build the memory report.
 
@@ -71,8 +92,9 @@ def compute_output_report(array, output_format, *, threshold=_SPARSE_AUTO_THRESH
         The stacked output. Sparse inputs are inspected through their shape,
         dtype, and stored values without converting them to a dense array.
     output_format : {"auto", "dense", "sparse"}
-        Requested format. ``"auto"`` picks ``"sparse"`` when the density is below
-        ``threshold``.
+        Requested format. ``"auto"`` picks ``"sparse"`` when the density of this
+        ``array`` is below ``threshold``; pass the format resolved at fit (see
+        :func:`resolve_output_format`) to keep the container stable across calls.
     threshold : float, default=0.3
         Density cut-off for the ``"auto"`` decision.
 
@@ -84,16 +106,9 @@ def compute_output_report(array, output_format, *, threshold=_SPARSE_AUTO_THRESH
         ``memory_saved_bytes``.
     """
     size = int(np.prod(array.shape, dtype=np.int64))
-    nonzero = int(array.count_nonzero()) if sp.issparse(array) else int(np.count_nonzero(array))
-    density = float(nonzero) / size if size else 0.0
+    density = _density(array)
     dense_bytes = size * int(array.dtype.itemsize)
-
-    if output_format == "sparse":
-        use_sparse = True
-    elif output_format == "auto":
-        use_sparse = density < threshold
-    else:  # "dense"
-        use_sparse = False
+    use_sparse = resolve_output_format(array, output_format, threshold=threshold) == "sparse"
 
     if use_sparse:
         csr = array.tocsr(copy=False) if sp.issparse(array) else sp.csr_matrix(array)
diff --git a/pretab/preprocessor.py b/pretab/preprocessor.py
index a8193b8..3f0d9f0 100644
--- a/pretab/preprocessor.py
+++ b/pretab/preprocessor.py
@@ -25,6 +25,7 @@
     compute_output_report,
     format_output,
     resolve_embedding_dimensions,
+    resolve_output_format,
     to_dataframe_output,
     validate_embedding_request,
 )
@@ -269,10 +270,12 @@ class Preprocessor(TransformerMixin, BaseEstimator):
         Container used for the transformed output. ``"dense"`` (the default, for
         backward compatibility) returns NumPy arrays; ``"sparse"`` returns SciPy
         CSR matrices (a single stacked CSR when ``return_array=True``, otherwise CSR
-        blocks in the output dict); ``"auto"`` selects ``"sparse"`` when the output
-        density falls below ``0.3`` and ``"dense"`` otherwise. Ignored when
+        blocks in the output dict); ``"auto"`` selects ``"sparse"`` when the density of
+        the training output falls below ``0.3`` and ``"dense"`` otherwise. The choice is
+        made once at ``fit`` (stored in ``output_format_``), so every ``transform``
+        returns the same container regardless of the batch. Ignored when
         :meth:`set_output` requests a pandas or polars DataFrame. Every ``transform``
-        records the resolved choice and its memory footprint in ``output_report_``.
+        records the format and its memory footprint in ``output_report_``.
     dtype : numpy dtype or None, default=None
         Optional dtype to cast the transformed output to (e.g. ``numpy.float32`` to
         halve memory). ``None`` keeps the native ``float64`` output.
@@ -314,6 +317,9 @@ class Preprocessor(TransformerMixin, BaseEstimator):
     output_dims\_ : dict
         Per-feature expanded output-column counts, keyed by input feature name.
         The values sum to ``total_output_dim_``.
+    output_format\_ : {"dense", "sparse"}
+        The output container resolved at ``fit`` from ``output_format`` (``"auto"``
+        is decided from the training output's density).
     output_report\_ : dict
         Memory report for the most recent ``transform``, with keys ``format``
         (``"dense"`` or ``"sparse"``), ``shape``, ``density``, ``dense_bytes``,
@@ -580,7 +586,12 @@ def fit(self, X, y=None, embeddings=None):
             categorical_features,
             sparse_threshold=sparse_threshold,
         )
-        self.column_transformer_.fit(X, y)
+        # ColumnTransformer.fit runs fit_transform internally; keep the training
+        # output to resolve output_format="auto" once, from the training density.
+        training_output = self.column_transformer_.fit_transform(X, y)
+        if not sp.issparse(training_output):
+            training_output = np.asarray(training_output)
+        self.output_format_ = resolve_output_format(training_output, self.output_format)
         self.n_features_in_ = X.shape[1]
 
         self._enforce_output_budget(X.shape[0])
@@ -653,7 +664,9 @@ def transform(self, X, embeddings=None, return_array=None):
         if self.dtype is not None:
             transformed_X = transformed_X.astype(self.dtype, copy=False)
 
-        fmt, self.output_report_ = compute_output_report(transformed_X, self.output_format)
+        # Spec files written before output_format_ existed keep the per-call choice.
+        output_format = getattr(self, "output_format_", self.output_format)
+        fmt, self.output_report_ = compute_output_report(transformed_X, output_format)
 
         if container in ("pandas", "polars"):
             return to_dataframe_output(transformed_X, self.get_feature_names_out(), container)
diff --git a/tests/integration/test_output_format.py b/tests/integration/test_output_format.py
index abe2613..b16ccd3 100644
--- a/tests/integration/test_output_format.py
+++ b/tests/integration/test_output_format.py
@@ -5,6 +5,8 @@
 wrapping.
 """
 
+from typing import cast
+
 import numpy as np
 import pandas as pd
 import pytest
@@ -241,3 +243,36 @@ def test_invalid_output_format_raises(frame, y):
 
     with pytest.raises(InvalidParamError):
         _bspline(output_format="nope").fit(frame, y)
+
+
+# --- "auto" is resolved once per fit (issue #73) ----------------------------------
+
+
+def test_auto_format_is_fixed_at_fit_for_dense_training_output():
+    """Regression guard for issue #73: low-density single rows came back sparse."""
+    from sklearn.datasets import load_diabetes
+
+    X, y = cast("tuple[pd.DataFrame, pd.Series]", load_diabetes(return_X_y=True, as_frame=True))
+    pre = Preprocessor(output_format="auto", random_state=0).fit(X, y)
+    assert pre.output_format_ == "dense"
+    for i in range(len(X)):
+        assert isinstance(pre.transform(X.iloc[[i]]), np.ndarray)
+    # The report still describes each call's own density.
+    pre.transform(X.iloc[[0]])
+    assert pre.output_report_["format"] == "dense"
+    assert pre.output_report_["shape"] == (1, pre.total_output_dim_)
+
+
+def test_auto_format_is_fixed_at_fit_for_sparse_training_output():
+    rng = np.random.default_rng(0)
+    frame = pd.DataFrame({"city": rng.choice([f"c{i}" for i in range(20)], 200)})
+    pre = Preprocessor(categorical_method="one-hot", output_format="auto").fit(frame)
+    assert pre.output_format_ == "sparse"
+    assert sp.issparse(pre.transform(frame.iloc[:1]))
+    assert sp.issparse(pre.transform(frame))
+
+
+@pytest.mark.parametrize("output_format", ["dense", "sparse"])
+def test_explicit_formats_are_stored_as_resolved(frame, y, output_format):
+    pre = Preprocessor(numerical_method="minmax", output_format=output_format).fit(frame, y)
+    assert pre.output_format_ == output_format