diff --git a/docs/pages.jl b/docs/pages.jl
index 2a7c9501..6fd3cb79 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -1,5 +1,6 @@
 pages = [
     "Home" => "index.md",
-    "Functions" => "functions.md",
     "Tutorials" => Any["tutorials/deim_FitzHugh-Nagumo.md"],
+    "Functions" => "functions.md",
+    "Developer Interface" => "interfaces.md",
 ]
diff --git a/docs/src/interfaces.md b/docs/src/interfaces.md
new file mode 100644
index 00000000..a355add8
--- /dev/null
+++ b/docs/src/interfaces.md
@@ -0,0 +1,54 @@
+# Developer Interface
+
+The supported user workflow is deliberately small: construct a [`POD`](@ref) from
+state snapshots, then call [`reduce!`](@ref) with one of the built-in backends,
+[`SVD`](@ref), [`TSVD`](@ref), or [`RSVD`](@ref). The public API does not require users
+to interact with the abstract types described below.
+
+## POD reduction contract
+
+`POD` accepts either an `AbstractMatrix{T}` whose columns are state snapshots or an
+`AbstractVector{<:AbstractVector{T}}` containing state snapshots, with `T <: AbstractFloat`.
+All snapshots must describe the same state dimension. The positional `nmodes` form fixes
+the retained rank. The keyword form selects the rank from `min_renergy`, `min_nmodes`, and
+`max_nmodes`.
+
+The reduction contract is:
+
+1. Construct `pod = POD(snapshots, nmodes)` or `POD(snapshots; kwargs...)`.
+2. Call `reduce!(pod, SVD())`, `reduce!(pod, TSVD())`, or `reduce!(pod, RSVD())`.
+3. Read `pod.rbasis`, `pod.spectrum`, `pod.nmodes`, and `pod.renergy` after the call.
+
+`reduce!` mutates `pod` and returns `nothing`. The three concrete backends are the
+supported choices; their keyword arguments are forwarded to their respective SVD
+implementations as described in the API reference.
+
+```jldoctest
+julia> using ModelOrderReduction
+
+julia> snapshots = [1.0 0.0; 0.0 1.0];
+
+julia> pod = POD(snapshots, 1);
+
+julia> reduce!(pod, SVD()); pod.nmodes
+1
+```
+
+## Internal type hierarchy
+
+The following abstract types are implementation details. They are documented so that
+contributors can understand the source hierarchy, but they are not exported and are not
+stable extension points:
+
+```@docs
+ModelOrderReduction.AbstractReductionProblem
+ModelOrderReduction.AbstractMORProblem
+ModelOrderReduction.AbstractDRProblem
+ModelOrderReduction.AbstractSVD
+```
+
+Do not subtype these types, dispatch on them from downstream packages, or add new
+`reduce!` methods for custom algorithm types. Such methods would depend on internal
+dispatch and field invariants that are not covered by the public API. If a supported
+extension point is needed, it should first be designed and documented as a separate
+interface with generic tests.
diff --git a/src/DataReduction/POD.jl b/src/DataReduction/POD.jl
index 41936977..ed7cdf04 100644
--- a/src/DataReduction/POD.jl
+++ b/src/DataReduction/POD.jl
@@ -27,30 +27,44 @@ end
 _rsvd(data, n::Int, p::Int) = rsvd(data, n, p)
 
 """
-    POD(snapshots; min_renergy = 1.0, min_nmodes = 1, max_nmodes = length(snapshots[1]))
-    POD(snapshots, nmodes)
+    POD(snapshots; min_renergy = 1.0, min_nmodes = 1,
+        max_nmodes = length(snapshots[1])) -> POD
+    POD(snapshots, nmodes::Int) -> POD
 
 Proper orthogonal decomposition reduction problem built from state snapshots. Call
 [`reduce!`](@ref) with an SVD backend to compute its basis and spectrum.
 
 # Arguments
-- `snapshots`: a state-by-snapshot matrix or a vector of state vectors.
+- `snapshots::AbstractMatrix{T}`: a state-by-snapshot matrix, where each column is a
+  state snapshot.
+- `snapshots::AbstractVector{<:AbstractVector{T}}`: a vector of state snapshots.
 - `nmodes::Int`: the fixed number of retained modes. This positional form disables
   energy-based truncation.
 
-# Keyword Arguments
-- `min_renergy = 1.0`: minimum captured relative spectral energy when selecting modes.
+# Keywords
+- `min_renergy::T = 1.0`: minimum captured relative spectral energy when selecting modes.
 - `min_nmodes::Int = 1`: lower bound on the number of retained modes.
 - `max_nmodes::Int = length(snapshots[1])`: upper bound on the number of retained modes.
 
 # Fields
-- `snapshots`: the input snapshot data.
-- `min_renergy`, `min_nmodes`, `max_nmodes`: the truncation policy.
-- `nmodes`: the selected number of retained modes.
-- `rbasis`: the reduced basis after [`reduce!`](@ref), or `missing` before reduction.
-- `renergy`: the captured relative spectral energy.
-- `spectrum`: the singular-value spectrum after [`reduce!`](@ref), or `missing` before
-  reduction.
+- `snapshots::S`: the input snapshot data.
+- `min_renergy::T`: the minimum relative spectral energy used for truncation.
+- `min_nmodes::Int`: the lower bound on the number of retained modes.
+- `max_nmodes::Int`: the upper bound on the number of retained modes.
+- `nmodes::Int`: the selected number of retained modes.
+- `rbasis::Union{Missing, Matrix{T}}`: the reduced basis after [`reduce!`](@ref), or
+  `missing` before reduction.
+- `renergy::T`: the captured relative spectral energy.
+- `spectrum::Union{Missing, Vector{T}}`: the singular-value spectrum after
+  [`reduce!`](@ref), or `missing` before reduction.
+
+# Throws
+- `AssertionError`: if the snapshots are not vector-valued or the truncation bounds are
+  invalid.
+
+# Returns
+- `POD`: an initialized reduction problem. The basis and spectrum are populated by
+  [`reduce!`](@ref).
 
 # Examples
 ```jldoctest
@@ -126,6 +140,9 @@ Compute the reduced basis and full singular-value spectrum for `pod` using dense
 - `pod::POD`: reduction problem to update in place.
 - `alg::SVD`: dense singular value decomposition backend.
 
+# Returns
+- `nothing`: `pod` is updated in place.
+
 # Examples
 ```jldoctest
 julia> using ModelOrderReduction
@@ -157,6 +174,23 @@ truncated SVD.
 # Arguments
 - `pod::POD`: reduction problem to update in place.
 - `alg::TSVD`: truncated singular value decomposition backend.
+
+# Returns
+- `nothing`: `pod` is updated in place.
+
+# Throws
+- An exception propagated by `TSVD.tsvd` if the snapshot matrix is not accepted by the
+  truncated SVD backend.
+
+# Examples
+```jldoctest
+julia> using ModelOrderReduction
+
+julia> pod = POD([3.0 0.0; 0.0 1.0], 1);
+
+julia> reduce!(pod, TSVD()); pod.nmodes
+1
+```
 """
 function reduce!(pod::POD{S, T}, alg::TSVD)::Nothing where {S, T}
     u, s, v = _tsvd(pod.snapshots, pod.nmodes; alg.kwargs...)
@@ -176,6 +210,23 @@ randomized SVD.
 # Arguments
 - `pod::POD`: reduction problem to update in place.
 - `alg::RSVD`: randomized singular value decomposition backend.
+
+# Returns
+- `nothing`: `pod` is updated in place.
+
+# Throws
+- An exception propagated by `RandomizedLinAlg.rsvd` if the snapshot matrix is not
+  accepted by the randomized SVD backend.
+
+# Examples
+```jldoctest
+julia> using ModelOrderReduction
+
+julia> pod = POD([3.0 0.0; 0.0 1.0], 1);
+
+julia> reduce!(pod, RSVD()); pod.nmodes
+1
+```
 """
 function reduce!(pod::POD{S, T}, alg::RSVD)::Nothing where {S, T}
     u, s, v = _rsvd(pod.snapshots, pod.nmodes, alg.p)
diff --git a/src/Types.jl b/src/Types.jl
index 05ce4f9f..6249a0c5 100644
--- a/src/Types.jl
+++ b/src/Types.jl
@@ -1,15 +1,50 @@
+"""
+    AbstractReductionProblem
+
+Internal common supertype for model-reduction problem types.
+
+This type is not exported and is not a supported extension point. User code should
+construct one of the documented concrete problem types instead.
+"""
 abstract type AbstractReductionProblem end
+
+"""
+    AbstractMORProblem <: AbstractReductionProblem
+
+Internal marker type for model-order-reduction problems.
+
+This type is not exported and is not a supported extension point. Its hierarchy may
+change without a public API guarantee.
+"""
 abstract type AbstractMORProblem <: AbstractReductionProblem end
+
+"""
+    AbstractDRProblem <: AbstractReductionProblem
+
+Internal marker type for data-reduction problems such as [`POD`](@ref).
+
+This type is not exported and is not a supported extension point. Use the documented
+[`POD`](@ref) constructors and [`reduce!`](@ref) methods instead.
+"""
 abstract type AbstractDRProblem <: AbstractReductionProblem end
 
+"""
+    AbstractSVD
+
+Internal supertype for the built-in singular-value-decomposition backends.
+
+This type is not exported and is not a supported extension point. The supported
+backends are [`SVD`](@ref), [`TSVD`](@ref), and [`RSVD`](@ref); custom subtypes and
+additional `reduce!` dispatches are not part of the public API.
+"""
 abstract type AbstractSVD end
 
 """
-    SVD(; kwargs...)
+    SVD(; kwargs...) -> SVD
 
 Dense singular value decomposition backend for [`reduce!`](@ref).
 
-# Keyword Arguments
+# Keywords
 - `kwargs...`: keyword arguments forwarded to `LinearAlgebra.svd`.
 
 # Fields
@@ -32,12 +67,12 @@ struct SVD{K <: NamedTuple} <: AbstractSVD
 end
 
 """
-    TSVD(; kwargs...)
+    TSVD(; kwargs...) -> TSVD
 
 Truncated singular value decomposition backend for [`reduce!`](@ref). Use this backend
 when only the requested reduced modes should be computed.
 
-# Keyword Arguments
+# Keywords
 - `kwargs...`: keyword arguments forwarded to `TSVD.tsvd`.
 
 # Fields
@@ -60,16 +95,19 @@ struct TSVD{K <: NamedTuple} <: AbstractSVD
 end
 
 """
-    RSVD([p = 0])
+    RSVD([p::Int = 0]) -> RSVD
 
 Randomized singular value decomposition backend for [`reduce!`](@ref).
 
 # Arguments
-- `p::Integer = 0`: number of oversampling vectors used by `RandomizedLinAlg.rsvd`.
+- `p::Int = 0`: number of oversampling vectors used by `RandomizedLinAlg.rsvd`.
 
 # Fields
 - `p::Int`: the oversampling parameter.
 
+# Throws
+- `MethodError`: if `p` is not an `Int`.
+
 # Examples
 ```jldoctest
 julia> using ModelOrderReduction
diff --git a/src/deim.jl b/src/deim.jl
index 1b97f2df..bbce3651 100644
--- a/src/deim.jl
+++ b/src/deim.jl
@@ -3,7 +3,8 @@ $(TYPEDSIGNATURES)
 
 Compute the DEIM interpolation indices for the given projection basis.
 
-The orthonormal `basis` should not be a sparse matrix.
+The orthonormal `basis` should not be a sparse matrix. This helper is internal and is
+not part of the exported API.
 """
 function deim_interpolation_indices(basis::AbstractMatrix)::Vector{Int}
     dim = size(basis, 2)
@@ -88,22 +89,34 @@ algorithm, and ``\\mathbf e_{\\rho_i}=[0,\\ldots,0,1,0,\\ldots,0]^T\\in\\mathbb
 the ``\\rho_i``-th column of the identity matrix ``I_n\\in\\mathbb R^{n\\times n}``.
 
 # Arguments
-- `full_vars::AbstractVector`: the dependent variables ``\\underset{n\\times 1}{\\mathbf y}`` in FOM.
-- `linear_coeffs::AbstractMatrix`: the coefficient matrix ``\\underset{n\\times n}A`` of linear terms in FOM.
-- `constant_part::AbstractVector`: the constant terms ``\\underset{n\\times 1}{\\mathbf g}`` in FOM.
-- `nonlinear_part::AbstractVector`: the nonlinear functions ``\\underset{n\\times 1}{\\mathbf F}`` in FOM.
-- `reduced_vars::AbstractVector`: the dependent variables ``\\underset{k\\times 1}{\\hat{\\mathbf y}}`` in the reduced-order model.
-- `linear_projection_matrix::AbstractMatrix`: the projection matrix ``\\underset{n\\times k}V`` for the dependent variables ``\\mathbf y``.
-- `nonlinear_projection_matrix::AbstractMatrix`: the projection matrix ``\\underset{n\\times m}U`` for the nonlinear functions ``\\mathbf F``.
-
-# Keyword Arguments
+- `full_vars::AbstractVector`: the dependent variables
+  ``\\underset{n\\times 1}{\\mathbf y}`` in FOM.
+- `linear_coeffs::AbstractMatrix`: the coefficient matrix
+  ``\\underset{n\\times n}A`` of linear terms in FOM.
+- `constant_part::AbstractVector`: the constant terms
+  ``\\underset{n\\times 1}{\\mathbf g}`` in FOM.
+- `nonlinear_part::AbstractVector`: the nonlinear functions
+  ``\\underset{n\\times 1}{\\mathbf F}`` in FOM.
+- `reduced_vars::AbstractVector`: the dependent variables
+  ``\\underset{k\\times 1}{\\hat{\\mathbf y}}`` in the reduced-order model.
+- `linear_projection_matrix::AbstractMatrix`: the projection matrix
+  ``\\underset{n\\times k}V`` for the dependent variables ``\\mathbf y``.
+- `nonlinear_projection_matrix::AbstractMatrix`: the projection matrix
+  ``\\underset{n\\times m}U`` for the nonlinear functions ``\\mathbf F``.
+
+# Keywords
 - `kwargs...`: keyword arguments forwarded to `Symbolics.substitute` when constructing
   the nonlinear reduced model.
 
-# Return
-- `reduced_rhss`: the right-hand side of ROM.
+# Returns
+- `Tuple{AbstractVector, Vector{Equation}}`: the reduced right-hand side and
+  linear projection equations.
 - `linear_projection_eqs`: the linear projection mapping ``\\mathbf y=V\\hat{\\mathbf y}``.
 
+# Throws
+- An exception from the matrix operations if the projection matrices have incompatible
+  dimensions.
+
 # Examples
 ```julia
 reduced_rhs, projection_equations = deim(
@@ -174,7 +187,7 @@ nonlinear expressions.
 - `snapshot::AbstractMatrix`: state-by-time snapshot matrix for `sys`.
 - `pod_dim::Integer`: number of POD state modes to retain.
 
-# Keyword Arguments
+# Keywords
 - `deim_dim::Integer = pod_dim`: number of DEIM modes for nonlinear terms.
 - `name::Symbol = Symbol(nameof(sys), :_deim)`: name assigned to the reduced system.
 - `kwargs...`: keyword arguments forwarded to ModelingToolkit transformations and
@@ -184,6 +197,12 @@ nonlinear expressions.
 ```julia
 reduced_system = deim(compiled_system, snapshots, 4; deim_dim = 6)
 ```
+
+# Returns
+- `ModelingToolkit.ODESystem`: the reduced and completed system.
+
+# Throws
+- `ArgumentError`: if the observed equations contain a dependency cycle.
 """
 function deim(
         sys::ODESystem, snapshot::AbstractMatrix, pod_dim::Integer;
diff --git a/test/interface.jl b/test/interface.jl
new file mode 100644
index 00000000..9654044c
--- /dev/null
+++ b/test/interface.jl
@@ -0,0 +1,17 @@
+using Test, ModelOrderReduction
+
+@testset "Public POD interface" begin
+    matrix_snapshots = [3.0 0.0; 0.0 1.0]
+    vector_snapshots = [[3.0, 0.0], [0.0, 1.0]]
+
+    for snapshots in (matrix_snapshots, vector_snapshots)
+        for alg in (SVD(), TSVD(), RSVD())
+            pod = POD(snapshots, 1)
+            @test reduce!(pod, alg) === nothing
+            @test size(pod.rbasis) == (2, 1)
+            @test length(pod.spectrum) >= pod.nmodes
+            @test pod.nmodes == 1
+            @test 0.0 < pod.renergy <= 1.0
+        end
+    end
+end