diff --git a/docs/src/reference/standard_form.md b/docs/src/reference/standard_form.md
index ffa3efac0c..0c2485540c 100644
--- a/docs/src/reference/standard_form.md
+++ b/docs/src/reference/standard_form.md
@@ -12,28 +12,30 @@ DocTestFilters = [r"MathOptInterface|MOI"]
 
 ```@docs
 AbstractFunction
+output_dimension
+constant
+```
+
+## Scalar functions
+
+```@docs
 AbstractScalarFunction
-AbstractVectorFunction
 VariableIndex
-VectorOfVariables
 ScalarAffineTerm
 ScalarAffineFunction
-VectorAffineTerm
-VectorAffineFunction
 ScalarQuadraticTerm
 ScalarQuadraticFunction
-VectorQuadraticTerm
-VectorQuadraticFunction
 ```
 
-### Utilities
+## Vector functions
 
 ```@docs
-output_dimension
-constant(f::Union{ScalarAffineFunction, ScalarQuadraticFunction})
-constant(f::Union{VectorAffineFunction, VectorQuadraticFunction})
-constant(f::VariableIndex, ::Type)
-constant(f::VectorOfVariables, T::Type)
+AbstractVectorFunction
+VectorOfVariables
+VectorAffineTerm
+VectorAffineFunction
+VectorQuadraticTerm
+VectorQuadraticFunction
 ```
 
 ## Sets
diff --git a/src/functions.jl b/src/functions.jl
index b320c2e552..2807574e38 100644
--- a/src/functions.jl
+++ b/src/functions.jl
@@ -7,72 +7,353 @@
 """
     output_dimension(f::AbstractFunction)
 
-Return 1 if `f` has a scalar output and the number of output components if `f`
-has a vector output.
+Return 1 if `f` is an [`AbstractScalarFunction`](@ref), or the number of output
+components if `f` is an [`AbstractVectorFunction`](@ref).
 """
 function output_dimension end
 
 output_dimension(::AbstractScalarFunction) = 1
 
-Base.broadcastable(f::AbstractScalarFunction) = Ref(f)
-Base.ndims(::Type{<:AbstractScalarFunction}) = 0
-Base.ndims(::AbstractScalarFunction) = 0
+"""
+    constant(f::AbstractFunction[, ::Type{T}]) where {T}
+
+Returns the constant term of a scalar-valued function, or the constant vector of
+a vector-valued function.
 
+If `f` is untyped and `T` is provided, returns `zero(T)`.
 """
-    VectorOfVariables(variables)
+constant(f::AbstractFunction, ::Type{T}) where {T} = constant(f)
 
-The function that extracts the vector of variables referenced by `variables`, a `Vector{VariableIndex}`.
-This function is naturally be used for constraints that apply to groups of variables, such as an "all different" constraint, an indicator constraint, or a complementarity constraint.
 """
-struct VectorOfVariables <: AbstractVectorFunction
-    variables::Vector{VariableIndex}
-end
-output_dimension(f::VectorOfVariables) = length(f.variables)
+    coefficient(t::ScalarAffineTerm)
+    coefficient(t::ScalarQuadraticTerm)
+    coefficient(t::VectorAffineTerm)
+    coefficient(t::VectorQuadraticTerm)
 
+Finds the coefficient stored in the term `t`.
 """
-    struct ScalarAffineTerm{T}
-        coefficient::T
-        variable::VariableIndex
-    end
+function coefficient end
+
+"""
+    term_indices(t::ScalarAffineTerm)
+    term_indices(t::ScalarQuadraticTerm)
+    term_indices(t::VectorAffineTerm)
+    term_indices(t::VectorQuadraticTerm)
+
+Returns the indices of the input term `t` as a tuple of `Int`s.
+
+* For `t::ScalarAffineTerm`, this is a 1-tuple of the variable index.
+* For `t::ScalarQuadraticTerm`, this is a 2-tuple of the variable indices
+  in non-decreasing order.
+* For `t::VectorAffineTerm`, this is a 2-tuple of the row/output and
+  variable indices.
+* For `t::VectorQuadraticTerm`, this is a 3-tuple of the row/output and
+  variable indices in non-decreasing order.
+"""
+function term_indices end
+
+"""
+    term_pair(t::ScalarAffineTerm)
+    term_pair(t::ScalarQuadraticTerm)
+    term_pair(t::VectorAffineTerm)
+    term_pair(t::VectorQuadraticTerm)
 
-Represents ``c x_i`` where ``c`` is `coefficient` and ``x_i`` is the variable
-identified by `variable`.
+Returns the pair [`term_indices`](@ref) `=>` [`coefficient`](@ref) of the term.
+"""
+function term_pair end
+
+# VariableIndex is defined in indextypes.jl
+
+constant(f::VariableIndex, ::Type{T}) where {T} = zero(T)
+
+Base.copy(x::VariableIndex) = x
+
+Base.isapprox(x::VariableIndex, y::VariableIndex; kwargs...) = x == y
+
+"""
+    ScalarAffineTerm{T}(coefficient::T, variable::VariableIndex) where {T}
+
+Represents the scalar-valued term `coefficient * variable`.
+
+## Example
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> x = MOI.VariableIndex(1)
+MathOptInterface.VariableIndex(1)
+
+julia> MOI.ScalarAffineTerm(2.0, x)
+MathOptInterface.ScalarAffineTerm{Float64}(2.0, MathOptInterface.VariableIndex(1))
+```
 """
 struct ScalarAffineTerm{T}
     coefficient::T
     variable::VariableIndex
 end
 
-# Note: ScalarAffineFunction is mutable because its `constant` field is likely of an immutable
-# type, while its `terms` field is of a mutable type, meaning that creating a `ScalarAffineFunction`
-# allocates, and it is desirable to provide a zero-allocation option for working with
-# ScalarAffineFunctions. See https://github.com/jump-dev/MathOptInterface.jl/pull/343.
+coefficient(t::ScalarAffineTerm) = t.coefficient
+
+term_indices(t::ScalarAffineTerm) = (t.variable.value,)
+
+term_pair(t::ScalarAffineTerm) = term_indices(t) => coefficient(t)
+
+# !!! developer note
+#
+#     ScalarAffineFunction is mutable because its `constant` field is likely of
+#     an immutable type, while its `terms` field is of a mutable type, meaning
+#     that creating a `ScalarAffineFunction` allocates, and it is desirable to
+#     provide a zero-allocation option for working with ScalarAffineFunctions.
+#
+#     See https://github.com/jump-dev/MathOptInterface.jl/pull/343.
+
 """
-    ScalarAffineFunction{T}(terms, constant)
+    ScalarAffineFunction{T}(terms::ScalarAffineTerm{T}, constant::T) where {T}
 
-The scalar-valued affine function ``a^T x + b``, where:
-* ``a`` is a sparse vector specified by a list of
-  [`ScalarAffineTerm`](@ref) structs.
-* ``b`` is a scalar specified by `constant::T`
+Represents the scalar-valued affine function ``a^\\top x + b``, where:
+
+ * ``a^\\top x`` is represented by the vector of [`ScalarAffineTerm`](@ref)s
+ * ``b`` is a scalar `constant::T`
+
+## Duplicates
 
 Duplicate variable indices in `terms` are accepted, and the corresponding
 coefficients are summed together.
+
+## Example
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> x = MOI.VariableIndex(1)
+MathOptInterface.VariableIndex(1)
+
+julia> terms = [MOI.ScalarAffineTerm(2.0, x), MOI.ScalarAffineTerm(3.0, x)]
+2-element Vector{MathOptInterface.ScalarAffineTerm{Float64}}:
+ MathOptInterface.ScalarAffineTerm{Float64}(2.0, MathOptInterface.VariableIndex(1))
+ MathOptInterface.ScalarAffineTerm{Float64}(3.0, MathOptInterface.VariableIndex(1))
+
+julia> f = MOI.ScalarAffineFunction(terms, 4.0)
+MathOptInterface.ScalarAffineFunction{Float64}(MathOptInterface.ScalarAffineTerm{Float64}[MathOptInterface.ScalarAffineTerm{Float64}(2.0, MathOptInterface.VariableIndex(1)), MathOptInterface.ScalarAffineTerm{Float64}(3.0, MathOptInterface.VariableIndex(1))], 4.0)
+```
 """
 mutable struct ScalarAffineFunction{T} <: AbstractScalarFunction
     terms::Vector{ScalarAffineTerm{T}}
     constant::T
 end
 
+constant(f::ScalarAffineFunction) = f.constant
+
+function Base.copy(f::ScalarAffineFunction)
+    return ScalarAffineFunction(copy(f.terms), copy(f.constant))
+end
+
+function ScalarAffineFunction{T}(x::VariableIndex) where {T}
+    return ScalarAffineFunction([ScalarAffineTerm(one(T), x)], zero(T))
+end
+
 """
-    struct VectorAffineTerm{T}
-        output_index::Int64
-        scalar_term::ScalarAffineTerm{T}
-    end
+    ScalarQuadraticTerm{T}(
+        coefficient::T,
+        variable_1::VariableIndex,
+        variable_2::VariableIndex,
+    ) where {T}
+
+Represents the scalar-valued term ``c x_i x_j`` where ``c`` is `coefficient`,
+``x_i`` is `variable_1` and ``x_j`` is `variable_2`.
+
+## Example
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> x = MOI.VariableIndex(1)
+MathOptInterface.VariableIndex(1)
+
+julia> MOI.ScalarQuadraticTerm(2.0, x, x)
+MathOptInterface.ScalarQuadraticTerm{Float64}(2.0, MathOptInterface.VariableIndex(1), MathOptInterface.VariableIndex(1))
+```
+"""
+struct ScalarQuadraticTerm{T}
+    coefficient::T
+    variable_1::VariableIndex
+    variable_2::VariableIndex
+end
+
+coefficient(t::ScalarQuadraticTerm) = t.coefficient
+
+function term_indices(t::ScalarQuadraticTerm)
+    return minmax(t.variable_1.value, t.variable_2.value)
+end
+
+term_pair(t::ScalarQuadraticTerm) = term_indices(t) => coefficient(t)
+
+# !!! developer note
+#
+#     ScalarQuadraticFunction is mutable because its `constant` field is likely
+#     of an immutable type, while its `terms` field is of a mutable type,
+#     meaning that creating a `ScalarQuadraticFunction` allocates, and it is
+#     desirable to provide a zero-allocation option for working with
+#     ScalarQuadraticFunctions.
+#
+#     See https://github.com/jump-dev/MathOptInterface.jl/pull/343.
+
+"""
+    ScalarQuadraticFunction{T}(
+        quadratic_terms::Vector{ScalarQuadraticTerm{T}},
+        affine_terms::Vector{ScalarAffineTerm{T}},
+        constant::T,
+    ) wher {T}
+
+The scalar-valued quadratic function ``\\frac{1}{2}x^\\top Q x + a^\\top x + b``,
+where:
+
+ * ``Q`` is the symmetric matrix given by the vector of [`ScalarQuadraticTerm`](@ref)s
+ * ``a^\\top x`` is a sparse vector given by the vector of [`ScalarAffineTerm`](@ref)s
+ * ``b`` is the scalar `constant::T`.
+
+## Duplicates
+
+Duplicate indices in `quadratic_terms` or `affine_terms` are accepted, and the
+corresponding coefficients are summed together.
+
+In `quadratic_terms`, "mirrored" indices, `(q, r)` and `(r, q)` where `r` and
+`q` are [`VariableIndex`](@ref)es, are considered duplicates; only one needs to
+be specified.
+
+## The 0.5 factor
+
+Coupled with the interpretation of mirrored indices, the `0.5` factor in front
+of the ``Q`` matrix is a common source of bugs.
+
+As a rule, to represent ``a * x^2 + b * x * y``:
+
+ * The coefficient ``a`` in front of squared variables (diagonal elements in
+   ``Q``) must be doubled when creating a [`ScalarQuadraticTerm`](@ref)
+ * The coefficient ``b`` in front of off-diagonal elements in ``Q`` should be
+   left as ``b``, be cause the mirrored index ``b * y * x`` will be implicitly
+   added.
+
+## Example
+
+To represent the function ``f(x, y) = 2 * x^2 + 3 * x * y + 4 * x + 5``, do:
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> x = MOI.VariableIndex(1);
+
+julia> y = MOI.VariableIndex(2);
+
+julia> constant = 5.0;
+
+julia> affine_terms = [MOI.ScalarAffineTerm(4.0, x)];
+
+julia> quadratic_terms = [
+           MOI.ScalarQuadraticTerm(4.0, x, x),  # Note the changed coefficient
+           MOI.ScalarQuadraticTerm(3.0, x, y),
+       ]
+2-element Vector{MathOptInterface.ScalarQuadraticTerm{Float64}}:
+ MathOptInterface.ScalarQuadraticTerm{Float64}(4.0, MathOptInterface.VariableIndex(1), MathOptInterface.VariableIndex(1))
+ MathOptInterface.ScalarQuadraticTerm{Float64}(3.0, MathOptInterface.VariableIndex(1), MathOptInterface.VariableIndex(2))
+
+julia> f = MOI.ScalarQuadraticFunction(quadratic_terms, affine_terms, constant)
+MathOptInterface.ScalarQuadraticFunction{Float64}(MathOptInterface.ScalarQuadraticTerm{Float64}[MathOptInterface.ScalarQuadraticTerm{Float64}(4.0, MathOptInterface.VariableIndex(1), MathOptInterface.VariableIndex(1)), MathOptInterface.ScalarQuadraticTerm{Float64}(3.0, MathOptInterface.VariableIndex(1), MathOptInterface.VariableIndex(2))], MathOptInterface.ScalarAffineTerm{Float64}[MathOptInterface.ScalarAffineTerm{Float64}(4.0, MathOptInterface.VariableIndex(1))], 5.0)
+```
+
+"""
+mutable struct ScalarQuadraticFunction{T} <: AbstractScalarFunction
+    quadratic_terms::Vector{ScalarQuadraticTerm{T}}
+    affine_terms::Vector{ScalarAffineTerm{T}}
+    constant::T
+end
+
+constant(f::ScalarQuadraticFunction) = f.constant
+
+function Base.copy(f::ScalarQuadraticFunction)
+    return ScalarQuadraticFunction(
+        copy(f.quadratic_terms),
+        copy(f.affine_terms),
+        copy(f.constant),
+    )
+end
+
+"""
+    abstract type AbstractVectorFunction <: AbstractFunction
+
+Abstract supertype for vector-valued [`AbstractFunction`](@ref)s.
+
+## Required methods
+
+All subtypes of `AbstractVectorFunction` must implement:
+
+ * [`output_dimension`](@ref)
+"""
+abstract type AbstractVectorFunction <: AbstractFunction end
+
+"""
+    VectorOfVariables(variables::Vector{VariableIndex}) <: AbstractVectorFunction
+
+The vector-valued function `f(x) = variables`, where `variables` is a subset of
+[`VariableIndex`](@ref)es in the model.
+
+The list of `variables` may contain duplicates.
+
+## Example
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> x = MOI.VariableIndex.(1:2)
+2-element Vector{MathOptInterface.VariableIndex}:
+ MathOptInterface.VariableIndex(1)
+ MathOptInterface.VariableIndex(2)
+
+julia> f = MOI.VectorOfVariables([x[1], x[2], x[1]])
+MathOptInterface.VectorOfVariables(MathOptInterface.VariableIndex[MathOptInterface.VariableIndex(1), MathOptInterface.VariableIndex(2), MathOptInterface.VariableIndex(1)])
+
+julia> MOI.output_dimension(f)
+3
+```
+"""
+struct VectorOfVariables <: AbstractVectorFunction
+    variables::Vector{VariableIndex}
+end
+
+output_dimension(f::VectorOfVariables) = length(f.variables)
+
+function constant(f::VectorOfVariables, ::Type{T}) where {T}
+    return zeros(T, output_dimension(f))
+end
+
+Base.copy(f::VectorOfVariables) = VectorOfVariables(copy(f.variables))
+
+function Base.:(==)(f::VectorOfVariables, g::VectorOfVariables)
+    return f.variables == g.variables
+end
+
+Base.isapprox(x::VectorOfVariables, y::VectorOfVariables; kwargs...) = x == y
 
-A `ScalarAffineTerm` plus its index of the output component of a
-`VectorAffineFunction` or `VectorQuadraticFunction`.
-`output_index` can also be interpreted as a row index into a sparse matrix,
-where the `scalar_term` contains the column index and coefficient.
+"""
+    VectorAffineTerm{T}(
+        output_index::Int64,
+        scalar_term::ScalarAffineTerm{T},
+    ) where {T}
+
+A `VectorAffineTerm` is a `scalar_term` that appears in the `output_index` row
+of the vector-valued [`VectorAffineFunction`](@ref) or
+[`VectorQuadraticFunction`](@ref).
+
+## Example
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> x = MOI.VariableIndex(1);
+
+julia> MOI.VectorAffineTerm(Int64(2), MOI.ScalarAffineTerm(3.0, x))
+MathOptInterface.VectorAffineTerm{Float64}(2, MathOptInterface.ScalarAffineTerm{Float64}(3.0, MathOptInterface.VariableIndex(1)))
+```
 """
 struct VectorAffineTerm{T}
     output_index::Int64
@@ -86,80 +367,89 @@ function VectorAffineTerm(
     return VectorAffineTerm(convert(Int64, output_index), scalar_term)
 end
 
+coefficient(t::VectorAffineTerm) = t.scalar_term.coefficient
+
+function term_indices(t::VectorAffineTerm)
+    return (t.output_index, term_indices(t.scalar_term)...)
+end
+
+term_pair(t::VectorAffineTerm) = term_indices(t) => coefficient(t)
+
 """
-    VectorAffineFunction{T}(terms, constants)
+    VectorAffineFunction{T}(
+        terms::Vector{VectorAffineTerm{T}},
+        constants::Vector{T},
+    ) where {T}
 
 The vector-valued affine function ``A x + b``, where:
-* ``A`` is a sparse matrix specified by a list of `VectorAffineTerm` objects.
-* ``b`` is a vector specified by `constants`
+
+ * ``A x`` is the sparse matrix given by the vector of [`VectorAffineTerm`](@ref)s
+ * ``b`` is the vector `constants`
+
+## Duplicates
 
 Duplicate indices in the ``A`` are accepted, and the corresponding coefficients
 are summed together.
+
+## Example
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> x = MOI.VariableIndex(1);
+
+julia> terms = [
+           MOI.VectorAffineTerm(Int64(1), MOI.ScalarAffineTerm(2.0, x)),
+           MOI.VectorAffineTerm(Int64(2), MOI.ScalarAffineTerm(3.0, x)),
+       ];
+
+julia> f = MOI.VectorAffineFunction(terms, [4.0, 5.0])
+MathOptInterface.VectorAffineFunction{Float64}(MathOptInterface.VectorAffineTerm{Float64}[MathOptInterface.VectorAffineTerm{Float64}(1, MathOptInterface.ScalarAffineTerm{Float64}(2.0, MathOptInterface.VariableIndex(1))), MathOptInterface.VectorAffineTerm{Float64}(2, MathOptInterface.ScalarAffineTerm{Float64}(3.0, MathOptInterface.VariableIndex(1)))], [4.0, 5.0])
+
+julia> MOI.output_dimension(f)
+2
+```
 """
 struct VectorAffineFunction{T} <: AbstractVectorFunction
     terms::Vector{VectorAffineTerm{T}}
     constants::Vector{T}
 end
+
 output_dimension(f::VectorAffineFunction) = length(f.constants)
 
-"""
-    struct ScalarQuadraticTerm{T}
-        coefficient::T
-        variable_1::VariableIndex
-        variable_2::VariableIndex
-    end
+constant(f::VectorAffineFunction) = f.constants
 
-Represents ``c x_i x_j`` where ``c`` is `coefficient`, ``x_i`` is the variable
-identified by `variable_1` and ``x_j`` is the variable identified by
-`variable_2`.
-"""
-struct ScalarQuadraticTerm{T}
-    coefficient::T
-    variable_1::VariableIndex
-    variable_2::VariableIndex
+function Base.copy(f::VectorAffineFunction)
+    return VectorAffineFunction(copy(f.terms), copy(f.constants))
 end
 
-# Note: ScalarQuadraticFunction is mutable because its `constant` field is
-# likely of an immutable type, while its other fields are of mutable types,
-# meaning that creating a `ScalarQuadraticFunction` allocates, and it is
-# desirable to provide a zero-allocation option for working with
-# ScalarQuadraticFunctions.
-#
-# See https://github.com/jump-dev/MathOptInterface.jl/pull/343.
+function VectorAffineFunction{T}(f::VectorOfVariables) where {T}
+    terms = map(1:output_dimension(f)) do i
+        return VectorAffineTerm(i, ScalarAffineTerm(one(T), f.variables[i]))
+    end
+    return VectorAffineFunction(terms, zeros(T, output_dimension(f)))
+end
 
 """
-    ScalarQuadraticFunction{T}(quadratic_terms, affine_terms, constant)
+    VectorQuadraticTerm{T}(
+        output_index::Int64,
+        scalar_term::ScalarQuadraticTerm{T},
+    ) where {T}
 
-The scalar-valued quadratic function ``\\frac{1}{2}x^TQx + a^T x + b``, where:
-* ``a`` is a sparse vector specified by a list of `ScalarAffineTerm` structs.
-* ``b`` is a scalar specified by `constant`.
-* ``Q`` is a symmetric matrix specified by a list of `ScalarQuadraticTerm`
-  structs.
+A `VectorQuadraticTerm` is a [`ScalarQuadraticTerm`](@ref) `scalar_term` that
+appears in the `output_index` row of the vector-valued
+[`VectorQuadraticFunction`](@ref).
 
-Duplicate indices in ``a`` or ``Q`` are accepted, and the corresponding
-coefficients are summed together. "Mirrored" indices `(q,r)` and `(r,q)` (where
-`r` and `q` are `VariableIndex`es) are considered duplicates; only one need be
-specified.
+## Example
 
-For example, for two scalar variables ``y, z``, the quadratic expression
-``yz + y^2`` is represented by the terms
-`ScalarQuadraticTerm.([1.0, 2.0], [y, y], [z, y])`.
-"""
-mutable struct ScalarQuadraticFunction{T} <: AbstractScalarFunction
-    quadratic_terms::Vector{ScalarQuadraticTerm{T}}
-    affine_terms::Vector{ScalarAffineTerm{T}}
-    constant::T
-end
+```jldoctest
+julia> import MathOptInterface as MOI
 
-"""
-    struct VectorQuadraticTerm{T}
-        output_index::Int64
-        scalar_term::ScalarQuadraticTerm{T}
-    end
+julia> x = MOI.VariableIndex(1);
 
-A [`ScalarQuadraticTerm`](@ref) plus its index of the output component of a
-`VectorQuadraticFunction`. Each output component corresponds to a
-distinct sparse matrix ``Q_i``.
+julia> MOI.VectorQuadraticTerm(Int64(2), MOI.ScalarQuadraticTerm(3.0, x, x))
+MathOptInterface.VectorQuadraticTerm{Float64}(2, MathOptInterface.ScalarQuadraticTerm{Float64}(3.0, MathOptInterface.VariableIndex(1), MathOptInterface.VariableIndex(1)))
+```
 """
 struct VectorQuadraticTerm{T}
     output_index::Int64
@@ -173,21 +463,64 @@ function VectorQuadraticTerm(
     return VectorQuadraticTerm(convert(Int64, output_index), scalar_term)
 end
 
+coefficient(t::VectorQuadraticTerm) = t.scalar_term.coefficient
+
+function term_indices(t::VectorQuadraticTerm)
+    return (t.output_index, term_indices(t.scalar_term)...)
+end
+
+term_pair(t::VectorQuadraticTerm) = term_indices(t) => coefficient(t)
+
 """
-    VectorQuadraticFunction{T}(quadratic_terms, affine_terms, constants)
+    VectorQuadraticFunction{T}(
+        quadratic_terms::Vector{VectorQuadraticTerm{T}},
+        affine_terms::Vector{VectorAffineTerm{T}},
+        constants::Vector{T},
+    ) where {T}
 
 The vector-valued quadratic function with i`th` component ("output index")
-defined as ``\\frac{1}{2}x^TQ_ix + a_i^T x + b_i``, where:
-* ``a_i`` is a sparse vector specified by the `VectorAffineTerm`s with
-  `output_index == i`.
-* ``b_i`` is a scalar specified by `constants[i]`
-* ``Q_i`` is a symmetric matrix specified by the `VectorQuadraticTerm` with
-  `output_index == i`.
-
-Duplicate indices in ``a_i`` or ``Q_i`` are accepted, and the corresponding
-coefficients are summed together. "Mirrored" indices `(q,r)` and `(r,q)` (where
-`r` and `q` are `VariableIndex`es) are considered duplicates; only one need be
-specified.
+defined as ``\\frac{1}{2}x^\\top Q_i x + a_i^\\top x + b_i``, where:
+
+ * ``\\frac{1}{2}x^\\top Q_i x`` is the symmetric matrix given by the
+   [`VectorQuadraticTerm`](@ref) elements in `quadratic_terms` with
+   `output_index == i`
+ * ``a_i^\\top x`` is the sparse vector given by the [`VectorAffineTerm`](@ref)
+   elements in `affine_terms` with `output_index == i`
+ * ``b_i`` is a scalar given by `constants[i]`
+
+## Duplicates
+
+Duplicate indices in `quadratic_terms` and `affine_terms` with the same
+`output_index` are handled in the same manner as duplicates in
+[`ScalarQuadraticFunction`](@ref).
+
+## Example
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> x = MOI.VariableIndex(1);
+
+julia> y = MOI.VariableIndex(2);
+
+julia> constants = [4.0, 5.0];
+
+julia> affine_terms = [
+           MOI.VectorAffineTerm(Int64(1), MOI.ScalarAffineTerm(2.0, x)),
+           MOI.VectorAffineTerm(Int64(2), MOI.ScalarAffineTerm(3.0, x)),
+       ];
+
+julia> quad_terms = [
+        MOI.VectorQuadraticTerm(Int64(1), MOI.ScalarQuadraticTerm(2.0, x, x)),
+        MOI.VectorQuadraticTerm(Int64(2), MOI.ScalarQuadraticTerm(3.0, x, y)),
+           ];
+
+julia> f = MOI.VectorQuadraticFunction(quad_terms, affine_terms, constants)
+MathOptInterface.VectorQuadraticFunction{Float64}(MathOptInterface.VectorQuadraticTerm{Float64}[MathOptInterface.VectorQuadraticTerm{Float64}(1, MathOptInterface.ScalarQuadraticTerm{Float64}(2.0, MathOptInterface.VariableIndex(1), MathOptInterface.VariableIndex(1))), MathOptInterface.VectorQuadraticTerm{Float64}(2, MathOptInterface.ScalarQuadraticTerm{Float64}(3.0, MathOptInterface.VariableIndex(1), MathOptInterface.VariableIndex(2)))], MathOptInterface.VectorAffineTerm{Float64}[MathOptInterface.VectorAffineTerm{Float64}(1, MathOptInterface.ScalarAffineTerm{Float64}(2.0, MathOptInterface.VariableIndex(1))), MathOptInterface.VectorAffineTerm{Float64}(2, MathOptInterface.ScalarAffineTerm{Float64}(3.0, MathOptInterface.VariableIndex(1)))], [4.0, 5.0])
+
+julia> MOI.output_dimension(f)
+2
+```
 """
 struct VectorQuadraticFunction{T} <: AbstractVectorFunction
     quadratic_terms::Vector{VectorQuadraticTerm{T}}
@@ -197,20 +530,34 @@ end
 
 output_dimension(f::VectorQuadraticFunction) = length(f.constants)
 
+constant(f::VectorQuadraticFunction) = f.constants
+
+function Base.copy(f::VectorQuadraticFunction)
+    return VectorQuadraticFunction(
+        copy(f.quadratic_terms),
+        copy(f.affine_terms),
+        copy(f.constants),
+    )
+end
+
 # Function modifications
 
 """
     AbstractFunctionModification
 
-An abstract supertype for structs which specify partial modifications to functions, to be used for making small modifications instead of replacing the functions entirely.
+An abstract supertype for structs which specify partial modifications to
+functions, to be used for making small modifications instead of replacing the
+functions entirely.
 """
 abstract type AbstractFunctionModification end
 
 """
     ScalarConstantChange{T}(new_constant::T)
 
-A struct used to request a change in the constant term of a scalar-valued function.
-Applicable to `ScalarAffineFunction` and `ScalarQuadraticFunction`.
+A struct used to request a change in the constant term of a scalar-valued
+function.
+
+Applicable to [`ScalarAffineFunction`](@ref) and [`ScalarQuadraticFunction`](@ref).
 """
 struct ScalarConstantChange{T} <: AbstractFunctionModification
     new_constant::T
@@ -219,8 +566,10 @@ end
 """
     VectorConstantChange{T}(new_constant::Vector{T})
 
-A struct used to request a change in the constant vector of a vector-valued function.
-Applicable to `VectorAffineFunction` and `VectorQuadraticFunction`.
+A struct used to request a change in the constant vector of a vector-valued
+function.
+
+Applicable to [`VectorAffineFunction`](@ref) and [`VectorQuadraticFunction`](@ref).
 """
 struct VectorConstantChange{T} <: AbstractFunctionModification
     new_constant::Vector{T}
@@ -231,24 +580,34 @@ end
 
 A struct used to request a change in the linear coefficient of a single variable
 in a scalar-valued function.
-Applicable to `ScalarAffineFunction` and `ScalarQuadraticFunction`.
+
+Applicable to [`ScalarAffineFunction`](@ref) and [`ScalarQuadraticFunction`](@ref).
 """
 struct ScalarCoefficientChange{T} <: AbstractFunctionModification
     variable::VariableIndex
     new_coefficient::T
 end
 
-# Note: MultiRowChange is mutable because its `variable` field of an immutable
-# type, while `new_coefficients` is of a mutable type, meaning that creating a `MultiRowChange`
-# allocates, and it is desirable to provide a zero-allocation option for working with
-# MultiRowChanges. See https://github.com/jump-dev/MathOptInterface.jl/pull/343.
+# !!! developer note
+#     MultiRowChange is mutable because its `variable` field of an immutable
+#     type, while `new_coefficients` is of a mutable type, meaning that creating
+#     a `MultiRowChange` allocates, and it is desirable to provide a
+#     zero-allocation option for working with MultiRowChanges.
+#
+#     See https://github.com/jump-dev/MathOptInterface.jl/pull/343.
+
 """
-    MultirowChange{T}(variable::VariableIndex, new_coefficients::Vector{Tuple{Int64, T}})
+    MultirowChange{T}(
+        variable::VariableIndex,
+        new_coefficients::Vector{Tuple{Int64,T}},
+    ) where {T}
 
 A struct used to request a change in the linear coefficients of a single
-variable in a vector-valued function. New coefficients are specified by
-`(output_index, coefficient)` tuples. Applicable to `VectorAffineFunction` and
-`VectorQuadraticFunction`.
+variable in a vector-valued function.
+
+New coefficients are specified by `(output_index, coefficient)` tuples.
+
+Applicable to [`VectorAffineFunction`](@ref) and [`VectorQuadraticFunction`](@ref).
 """
 mutable struct MultirowChange{T} <: AbstractFunctionModification
     variable::VariableIndex
@@ -265,18 +624,7 @@ function MultirowChange(
     )
 end
 
-# Implementation of comparison for functions
-function Base.:(==)(f::VectorOfVariables, g::VectorOfVariables)
-    return f.variables == g.variables
-end
-
-function Base.isapprox(
-    f::Union{VariableIndex,VectorOfVariables},
-    g::Union{VariableIndex,VectorOfVariables};
-    kwargs...,
-)
-    return f == g
-end
+# isapprox
 
 # For affine and quadratic functions, terms are compressed in a dictionary using
 # `_dicts` and then the dictionaries are compared with `dict_compare`
@@ -292,67 +640,12 @@ end
 # Build a dictionary where the duplicate keys are summed
 function sum_dict(kvs::Vector{Pair{K,V}}) where {K,V}
     d = Dict{K,V}()
-    for kv in kvs
-        key = kv.first
-        d[key] = kv.second + Base.get(d, key, zero(V))
+    for (key, value) in kvs
+        d[key] = value + Base.get(d, key, zero(V))
     end
     return d
 end
 
-"""
-    coefficient(t::Union{ScalarAffineTerm, ScalarQuadraticTerm
-                         VectorAffineTerm, VectorQuadraticTerm})
-
-Finds the coefficient stored in the term `t`.
-"""
-function coefficient end
-
-function coefficient(t::Union{ScalarAffineTerm,ScalarQuadraticTerm})
-    return t.coefficient
-end
-function coefficient(t::Union{VectorAffineTerm,VectorQuadraticTerm})
-    return t.scalar_term.coefficient
-end
-
-"""
-    term_indices(t::Union{ScalarAffineTerm, ScalarQuadraticTerm,
-                         VectorAffineTerm, VectorQuadraticTerm})
-
-Returns the indices of the input term `t` as a tuple of `Int`s.
-
-* For `t::ScalarAffineTerm`, this is a 1-tuple of the variable index.
-* For `t::ScalarQuadraticTerm`, this is a 2-tuple of the variable indices
-  in non-decreasing order.
-* For `t::VectorAffineTerm`, this is a 2-tuple of the row/output and
-  variable indices.
-* For `t::VectorQuadraticTerm`, this is a 3-tuple of the row/output and
-  variable indices in non-decreasing order.
-"""
-term_indices(t::ScalarAffineTerm) = (t.variable.value,)
-function term_indices(t::ScalarQuadraticTerm)
-    return minmax(t.variable_1.value, t.variable_2.value)
-end
-function term_indices(t::Union{VectorAffineTerm,VectorQuadraticTerm})
-    return (t.output_index, term_indices(t.scalar_term)...)
-end
-
-"""
-    term_pair(t::Union{ScalarAffineTerm, ScalarQuadraticTerm,
-                       VectorAffineTerm, VectorQuadraticTerm})
-
-Returns the pair [`term_indices`](@ref) `=>` [`coefficient`](@ref) of the term.
-"""
-function term_pair(
-    t::Union{
-        ScalarAffineTerm,
-        ScalarQuadraticTerm,
-        VectorAffineTerm,
-        VectorQuadraticTerm,
-    },
-)
-    return term_indices(t) => coefficient(t)
-end
-
 function _dicts(f::Union{ScalarAffineFunction,VectorAffineFunction})
     return (sum_dict(term_pair.(f.terms)),)
 end
@@ -364,20 +657,6 @@ function _dicts(f::Union{ScalarQuadraticFunction,VectorQuadraticFunction})
     )
 end
 
-"""
-    constant(f::Union{ScalarAffineFunction, ScalarQuadraticFunction})
-
-Returns the constant term of the scalar function
-"""
-constant(f::Union{ScalarAffineFunction,ScalarQuadraticFunction}) = f.constant
-
-"""
-    constant(f::Union{VectorAffineFunction, VectorQuadraticFunction})
-
-Returns the vector of constant terms of the vector function
-"""
-constant(f::Union{VectorAffineFunction,VectorQuadraticFunction}) = f.constants
-
 function Base.isapprox(
     f::F,
     g::G;
@@ -405,84 +684,12 @@ function Base.isapprox(
     )
 end
 
-function constant(f::Union{ScalarAffineFunction,ScalarQuadraticFunction}, ::Any)
-    return constant(f)
-end
-
-function constant(f::Union{VectorAffineFunction,VectorQuadraticFunction}, ::Any)
-    return constant(f)
-end
-
-"""
-    constant(f::VariableIndex, ::Type{T}) where {T}
-
-The constant term of a `VariableIndex` function is
-the zero value of the specified type `T`.
-"""
-constant(f::VariableIndex, ::Type{T}) where {T} = zero(T)
-
-"""
-    constant(f::VectorOfVariables, ::Type{T}) where {T}
-
-The constant term of a `VectorOfVariables` function is a
-vector of zero values of the specified type `T`.
-"""
-function constant(f::VectorOfVariables, ::Type{T}) where {T}
-    return zeros(T, length(f.variables))
-end
-
-# isbits type, nothing to copy
-Base.copy(func::VariableIndex) = func
-
-Base.copy(func::VectorOfVariables) = VectorOfVariables(copy(func.variables))
-
-"""
-    copy(func::Union{ScalarAffineFunction, VectorAffineFunction})
-
-Return a new affine function with a shallow copy of the terms and constant(s)
-from `func`.
-"""
-function Base.copy(
-    func::F,
-) where {F<:Union{ScalarAffineFunction,VectorAffineFunction}}
-    return F(copy(func.terms), copy(constant(func)))
-end
-
-"""
-    copy(func::Union{ScalarQuadraticFunction, VectorQuadraticFunction})
-
-Return a new quadratic function with a shallow copy of the terms and constant(s)
-from `func`.
-"""
-function Base.copy(
-    func::F,
-) where {F<:Union{ScalarQuadraticFunction,VectorQuadraticFunction}}
-    return F(
-        copy(func.quadratic_terms),
-        copy(func.affine_terms),
-        copy(constant(func)),
-    )
-end
+###
+### Base.convert
+###
 
-# Define shortcuts for
-# VariableIndex -> ScalarAffineFunction
-function ScalarAffineFunction{T}(f::VariableIndex) where {T}
-    return ScalarAffineFunction([ScalarAffineTerm(one(T), f)], zero(T))
-end
-# VectorOfVariables -> VectorAffineFunction
-function VectorAffineFunction{T}(f::VectorOfVariables) where {T}
-    n = length(f.variables)
-    return VectorAffineFunction(
-        map(
-            i -> VectorAffineTerm(i, ScalarAffineTerm(one(T), f.variables[i])),
-            1:n,
-        ),
-        zeros(T, n),
-    )
-end
+# VariableIndex
 
-# Conversion between scalar functions
-# Conversion to VariableIndex
 function Base.convert(::Type{VariableIndex}, f::ScalarAffineFunction)
     if (
         !iszero(f.constant) ||
@@ -501,7 +708,8 @@ function Base.convert(
     return convert(VariableIndex, convert(ScalarAffineFunction{T}, f))
 end
 
-# Conversion to ScalarAffineFunction
+# ScalarAffineFunction
+
 function Base.convert(::Type{ScalarAffineFunction{T}}, α::T) where {T}
     return ScalarAffineFunction{T}(ScalarAffineTerm{T}[], α)
 end
@@ -537,7 +745,8 @@ function Base.convert(
     return ScalarAffineFunction{T}(f.affine_terms, f.constant)
 end
 
-# Conversion to ScalarQuadraticFunction
+# ScalarQuadraticFunction
+
 function Base.convert(::Type{ScalarQuadraticFunction{T}}, α::T) where {T}
     return ScalarQuadraticFunction{T}(
         ScalarQuadraticTerm{T}[],
@@ -567,26 +776,37 @@ function Base.convert(
     )
 end
 
-function Base.convert(::Type{VectorOfVariables}, g::VariableIndex)
-    return VectorOfVariables([g])
+function Base.convert(
+    ::Type{ScalarQuadraticTerm{T}},
+    f::ScalarQuadraticTerm,
+) where {T}
+    return ScalarQuadraticTerm{T}(f.coefficient, f.variable_1, f.variable_2)
 end
 
 function Base.convert(
-    ::Type{VectorAffineFunction{T}},
-    g::VariableIndex,
+    ::Type{ScalarQuadraticFunction{T}},
+    f::ScalarQuadraticFunction,
 ) where {T}
-    return VectorAffineFunction{T}(
-        [VectorAffineTerm(1, ScalarAffineTerm(one(T), g))],
-        [zero(T)],
+    return ScalarQuadraticFunction{T}(
+        f.quadratic_terms,
+        f.affine_terms,
+        f.constant,
     )
 end
 
+# VectorOfVariables
+
+function Base.convert(::Type{VectorOfVariables}, g::VariableIndex)
+    return VectorOfVariables([g])
+end
+
+# VectorAffineFunction
+
 function Base.convert(
-    ::Type{VectorQuadraticFunction{T}},
+    ::Type{VectorAffineFunction{T}},
     g::VariableIndex,
 ) where {T}
-    return VectorQuadraticFunction{T}(
-        VectorQuadraticTerm{T}[],
+    return VectorAffineFunction{T}(
         [VectorAffineTerm(1, ScalarAffineTerm(one(T), g))],
         [zero(T)],
     )
@@ -602,6 +822,33 @@ function Base.convert(
     )
 end
 
+function Base.convert(
+    ::Type{VectorAffineTerm{T}},
+    f::VectorAffineTerm,
+) where {T}
+    return VectorAffineTerm{T}(f.output_index, f.scalar_term)
+end
+
+function Base.convert(
+    ::Type{VectorAffineFunction{T}},
+    f::VectorAffineFunction,
+) where {T}
+    return VectorAffineFunction{T}(f.terms, f.constants)
+end
+
+# VectorQuadraticFunction
+
+function Base.convert(
+    ::Type{VectorQuadraticFunction{T}},
+    g::VariableIndex,
+) where {T}
+    return VectorQuadraticFunction{T}(
+        VectorQuadraticTerm{T}[],
+        [VectorAffineTerm(1, ScalarAffineTerm(one(T), g))],
+        [zero(T)],
+    )
+end
+
 function Base.convert(
     ::Type{VectorQuadraticFunction{T}},
     g::ScalarAffineFunction,
@@ -628,38 +875,6 @@ function Base.convert(
     )
 end
 
-function Base.convert(
-    ::Type{ScalarQuadraticTerm{T}},
-    f::ScalarQuadraticTerm,
-) where {T}
-    return ScalarQuadraticTerm{T}(f.coefficient, f.variable_1, f.variable_2)
-end
-
-function Base.convert(
-    ::Type{ScalarQuadraticFunction{T}},
-    f::ScalarQuadraticFunction,
-) where {T}
-    return ScalarQuadraticFunction{T}(
-        f.quadratic_terms,
-        f.affine_terms,
-        f.constant,
-    )
-end
-
-function Base.convert(
-    ::Type{VectorAffineTerm{T}},
-    f::VectorAffineTerm,
-) where {T}
-    return VectorAffineTerm{T}(f.output_index, f.scalar_term)
-end
-
-function Base.convert(
-    ::Type{VectorAffineFunction{T}},
-    f::VectorAffineFunction,
-) where {T}
-    return VectorAffineFunction{T}(f.terms, f.constants)
-end
-
 function Base.convert(
     ::Type{VectorQuadraticTerm{T}},
     f::VectorQuadraticTerm,
diff --git a/src/indextypes.jl b/src/indextypes.jl
index b3fd526e11..7720169af5 100644
--- a/src/indextypes.jl
+++ b/src/indextypes.jl
@@ -8,22 +8,32 @@
     AbstractFunction
 
 Abstract supertype for function objects.
+
+## Required methods
+
+All functions must implement:
+
+ * `Base.copy`
+ * `Base.isapprox`
+ * [`constant`](@ref)
+
+Abstract subtypes of `AbstractFunction` may require additional methods to be
+implemented.
 """
 abstract type AbstractFunction <: MA.AbstractMutable end
 
 """
-    AbstractScalarFunction
+    abstract type AbstractScalarFunction <: AbstractFunction
 
-Abstract supertype for scalar-valued function objects.
+Abstract supertype for scalar-valued [`AbstractFunction`](@ref)s.
 """
 abstract type AbstractScalarFunction <: AbstractFunction end
 
-"""
-    AbstractVectorFunction
+Base.broadcastable(f::AbstractScalarFunction) = Ref(f)
 
-Abstract supertype for vector-valued function objects.
-"""
-abstract type AbstractVectorFunction <: AbstractFunction end
+Base.ndims(::Type{<:AbstractScalarFunction}) = 0
+
+Base.ndims(::AbstractScalarFunction) = 0
 
 """
     VariableIndex