diff --git a/src/aff_expr.jl b/src/aff_expr.jl
index 079a3e7aa32..0580654e79c 100644
--- a/src/aff_expr.jl
+++ b/src/aff_expr.jl
@@ -185,8 +185,11 @@ function add_to_expression! end
 
 # TODO: add deprecations for Base.push! and Base.append!
 
-function add_to_expression!(aff::GenericAffExpr{C,V}, new_coef::C, new_var::V) where {C,V}
-    _add_or_set!(aff.terms, new_var, new_coef)
+# With one factor.
+
+function add_to_expression!(aff::GenericAffExpr{C,V},
+                            other::Real) where {C,V}
+    aff.constant += other
     return aff
 end
 
@@ -195,23 +198,48 @@ function add_to_expression!(aff::GenericAffExpr{C,V}, new_var::V) where {C,V}
     return aff
 end
 
-function add_to_expression!(aff::GenericAffExpr{C,V}, other::GenericAffExpr{C,V}) where {C,V}
+function add_to_expression!(aff::GenericAffExpr{C,V},
+                            other::GenericAffExpr{C,V}) where {C,V}
+    # Note: merge!() doesn't appear to call sizehint!(). Is this important?
     merge!(+, aff.terms, other.terms)
     aff.constant += other.constant
     return aff
 end
 
-function add_to_expression!(aff::GenericAffExpr{C,V}, other::C) where {C,V}
-    aff.constant += other
+# With two factors.
+
+function add_to_expression!(aff::GenericAffExpr{C,V}, new_coef::Real,
+                            new_var::V) where {C,V}
+    _add_or_set!(aff.terms, new_var, convert(C, new_coef))
     return aff
 end
-function add_to_expression!(aff::GenericAffExpr{C,V}, other::Real) where {C,V}
-    aff.constant += other
+
+function add_to_expression!(aff::GenericAffExpr{C,V}, new_var::V,
+                            new_coef::Real) where {C,V}
+    return add_to_expression!(aff, new_coef, new_var)
+end
+
+function add_to_expression!(aff::GenericAffExpr{C,V}, coef::Real,
+                            other::GenericAffExpr{C,V}) where {C,V}
+    sizehint!(aff, length(linear_terms(aff)) + length(linear_terms(other)))
+    for (term_coef, var) in linear_terms(other)
+        _add_or_set!(aff.terms, var, coef * term_coef)
+    end
+    aff.constant += coef * other.constant
     return aff
 end
 
-function Base.isequal(aff::GenericAffExpr{C,V},other::GenericAffExpr{C,V}) where {C,V}
-    return isequal(aff.constant, other.constant) && isequal(aff.terms, other.terms)
+function add_to_expression!(aff::GenericAffExpr{C,V},
+                            other::GenericAffExpr{C,V},
+                            coef::Real) where {C,V}
+    return add_to_expression!(aff, coef, other)
+end
+
+
+function Base.isequal(aff::GenericAffExpr{C,V},
+                      other::GenericAffExpr{C,V}) where {C,V}
+    return isequal(aff.constant, other.constant) &&
+        isequal(aff.terms, other.terms)
 end
 
 Base.hash(aff::GenericAffExpr, h::UInt) = hash(aff.constant, hash(aff.terms, h))
diff --git a/src/operators.jl b/src/operators.jl
index b6a6a02ae5c..aed50405c50 100644
--- a/src/operators.jl
+++ b/src/operators.jl
@@ -156,11 +156,7 @@ function Base.:+(lhs::GenericAffExpr{C, V}, rhs::GenericAffExpr{C, V}) where {C,
             operator_warn(owner_model(first(linear_terms(lhs))[2]))
         end
     end
-    result_terms = copy(lhs.terms)
-    # merge() returns a Dict(), so we need to call merge!() instead.
-    # Note: merge!() doesn't appear to call sizehint!(). Is this important?
-    merge!(+, result_terms, rhs.terms)
-    return GenericAffExpr(lhs.constant + rhs.constant, result_terms)
+    return add_to_expression!(copy(lhs), rhs)
 end
 
 function Base.:-(lhs::GenericAffExpr{C, V}, rhs::GenericAffExpr{C, V}) where {C, V<:_JuMPTypes}
@@ -175,44 +171,7 @@ end
 
 function Base.:*(lhs::GenericAffExpr{C, V}, rhs::GenericAffExpr{C, V}) where {C, V<:_JuMPTypes}
     result = zero(GenericQuadExpr{C, V})
-
-    lhs_length = length(linear_terms(lhs))
-    rhs_length = length(linear_terms(rhs))
-
-    # Quadratic terms
-    for (lhscoef, lhsvar) in linear_terms(lhs)
-        for (rhscoef, rhsvar) in linear_terms(rhs)
-            add_to_expression!(result, lhscoef*rhscoef, lhsvar, rhsvar)
-        end
-    end
-
-    # Try to preallocate space for aff
-    if !iszero(lhs.constant) && !iszero(rhs.constant)
-        sizehint!(result.aff, lhs_length + rhs_length)
-    elseif !iszero(lhs.constant)
-        sizehint!(result.aff, rhs_length)
-    elseif !iszero(rhs.constant)
-        sizehint!(result.aff, lhs_length)
-    end
-
-    # [LHS constant] * [RHS linear terms]
-    if !iszero(lhs.constant)
-        c = lhs.constant
-        for (rhscoef, rhsvar) in linear_terms(rhs)
-            add_to_expression!(result.aff, c*rhscoef, rhsvar)
-        end
-    end
-
-    # [RHS constant] * [LHS linear terms]
-    if !iszero(rhs.constant)
-        c = rhs.constant
-        for (lhscoef, lhsvar) in linear_terms(lhs)
-            add_to_expression!(result.aff, c*lhscoef, lhsvar)
-        end
-    end
-
-    result.aff.constant = lhs.constant * rhs.constant
-
+    add_to_expression!(result, lhs, rhs)
     return result
 end
 # GenericAffExpr--GenericQuadExpr
@@ -413,11 +372,10 @@ function _A_mul_B!(ret::AbstractArray{T}, A, B) where {T <: _JuMPTypes}
         q = ret[i, j]
         _sizehint_expr!(q, size(A, 2))
         for k ∈ 1:size(A, 2)
-            tmp = convert(T, A[i, k] * B[k, j])
-            add_to_expression!(q, tmp)
+            add_to_expression!(q, A[i, k], B[k, j])
         end
     end
-    ret
+    return ret
 end
 
 function _A_mul_B!(ret::AbstractArray{<:GenericAffOrQuadExpr}, A::SparseMatrixCSC, B)
@@ -426,11 +384,11 @@ function _A_mul_B!(ret::AbstractArray{<:GenericAffOrQuadExpr}, A::SparseMatrixCS
     for col ∈ 1:size(A, 2)
         for k ∈ 1:size(ret, 2)
             for j ∈ nzrange(A, col)
-                add_to_expression!(ret[rv[j], k], nzv[j] * B[col, k])
+                add_to_expression!(ret[rv[j], k], nzv[j], B[col, k])
             end
         end
     end
-    ret
+    return ret
 end
 
 function _A_mul_B!(ret::AbstractArray{<:GenericAffOrQuadExpr}, A::AbstractMatrix, B::SparseMatrixCSC)
@@ -442,11 +400,11 @@ function _A_mul_B!(ret::AbstractArray{<:GenericAffOrQuadExpr}, A::AbstractMatrix
             q = ret[multivec_row, col]
             _sizehint_expr!(q, length(idxset))
             for k ∈ idxset
-                add_to_expression!(q, A[multivec_row, rowval[k]] * nzval[k])
+                add_to_expression!(q, A[multivec_row, rowval[k]], nzval[k])
             end
         end
     end
-    ret
+    return ret
 end
 
 # TODO: Implement sparse * sparse code as in base/sparse/linalg.jl (spmatmul).
@@ -471,8 +429,7 @@ function _At_mul_B!(ret::AbstractArray{T}, A, B) where {T <: _JuMPTypes}
         q = ret[i, j]
         _sizehint_expr!(q, size(A, 1))
         for k ∈ 1:size(A, 1)
-            tmp = convert(T, A[k, i] * B[k, j]) # transpose
-            add_to_expression!(q, tmp)
+            add_to_expression!(q, A[k, i], B[k, j]) # transpose
         end
     end
     ret
@@ -517,23 +474,23 @@ function _At_mul_B!(ret::StridedVecOrMat{<:GenericAffOrQuadExpr}, A::SparseMatri
     size(B, 2) == size(ret, 2) || throw(DimensionMismatch())
     nzv = A.nzval
     rv = A.rowval
-    # ret is already filled with zeros by _return_arrayt.
+    # `ret` is already filled with zeros by `_return_arrayt`.
     for k = 1:size(ret, 2)
         @inbounds for col = 1:A.n
             tmp = zero(eltype(ret))
             for j = A.colptr[col]:(A.colptr[col + 1] - 1)
-                tmp += adjoint(nzv[j]) * B[rv[j],k]
+                add_to_expression!(tmp, adjoint(nzv[j]), B[rv[j],k])
             end
-            ret[col,k] += tmp
+            ret[col, k] += tmp
         end
     end
-    ret
+    return ret
 end
 function _At_mul_B(A, B)
     size(A, 1) == size(B, 1) || error("Incompatible sizes")
     ret = _A_mul_B_ret(transpose(A), B)
     _At_mul_B!(ret, A, B)
-    ret
+    return ret
 end
 
 # TODO: Implement sparse * sparse code as in base/sparse/linalg.jl (spmatmul).
@@ -617,7 +574,7 @@ function Base.:-(x::AbstractArray{T}) where {T <: _JuMPTypes}
     for I in eachindex(ret)
         ret[I] = -x[I]
     end
-    ret
+    return ret
 end
 
 ###############################################################################
diff --git a/src/quad_expr.jl b/src/quad_expr.jl
index 9b79aee976f..20b6b9277e2 100644
--- a/src/quad_expr.jl
+++ b/src/quad_expr.jl
@@ -126,33 +126,137 @@ function Base.eltype(qti::QuadTermIterator{GenericQuadExpr{C, V}}
     return Tuple{C, V, V}
 end
 
-function add_to_expression!(quad::GenericQuadExpr{C,V}, new_coef::C, new_var1::V, new_var2::V) where {C,V}
-    # Node: OrderedDict updates the *key* as well. That is, if there was a
-    # previous value for UnorderedPair(new_var2, new_var1), it's key will now be
-    # UnorderedPair(new_var1, new_var2) (because these are defined as equal).
-    key = UnorderedPair(new_var1, new_var2)
-    _add_or_set!(quad.terms, key, new_coef)
-    return quad
-end
+# With one factor.
 
-function add_to_expression!(quad::GenericQuadExpr{C, V}, new_coef::C, new_var::V) where {C,V}
-    add_to_expression!(quad.aff, new_coef, new_var)
-    return quad
+function add_to_expression!(quad::GenericQuadExpr{C}, other::C) where C
+    return add_to_expression!(quad.aff, other)
 end
 
-function add_to_expression!(q::GenericQuadExpr{T,S}, other::GenericAffExpr{T,S}) where {T,S}
+function add_to_expression!(q::GenericQuadExpr{T,S},
+                            other::GenericAffExpr{T,S}) where {T,S}
     add_to_expression!(q.aff, other)
     return q
 end
 
-function add_to_expression!(q::GenericQuadExpr{T,S}, other::GenericQuadExpr{T,S}) where {T,S}
+function add_to_expression!(q::GenericQuadExpr{T,S},
+                            other::GenericQuadExpr{T,S}) where {T,S}
     merge!(+, q.terms, other.terms)
     add_to_expression!(q.aff, other.aff)
     return q
 end
 
-function add_to_expression!(quad::GenericQuadExpr{C}, other::C) where C
-    return add_to_expression!(quad.aff, other)
+# With two factors.
+
+function add_to_expression!(quad::GenericQuadExpr{C, V},
+                            new_coef::Real,
+                            new_var::V) where {C,V}
+    add_to_expression!(quad.aff, new_coef, new_var)
+    return quad
+end
+
+function add_to_expression!(quad::GenericQuadExpr{C, V},
+                            new_var::Union{V, GenericAffExpr{C, V}},
+                            new_coef::Real) where {C,V}
+    return add_to_expression!(quad, new_coef, new_var)
+end
+
+function add_to_expression!(quad::GenericQuadExpr{C},
+                            new_coef::Real,
+                            new_aff::GenericAffExpr{C}) where {C}
+    add_to_expression!(quad.aff, new_coef, new_aff)
+    return quad
+end
+
+function add_to_expression!(quad::GenericQuadExpr{C, V}, coef::Real,
+                            other::GenericQuadExpr{C, V}) where {C, V}
+    for (key, term_coef) in other.terms
+        _add_or_set!(quad.terms, key, coef * term_coef)
+    end
+    return add_to_expression!(quad, coef, other.aff)
+end
+
+function add_to_expression!(quad::GenericQuadExpr{C, V},
+                            other::GenericQuadExpr{C, V},
+                            coef::Real) where {C, V}
+    return add_to_expression!(quad, coef, other)
+end
+
+function add_to_expression!(quad::GenericQuadExpr{C},
+                            var_1::AbstractVariableRef,
+                            var_2::AbstractVariableRef) where {C}
+    return add_to_expression!(quad, one(C), var_1, var_2)
+end
+
+function add_to_expression!(quad::GenericQuadExpr{C,V},
+                            var::V,
+                            aff::GenericAffExpr{C,V}) where {C,V}
+    for (coef, term_var) in linear_terms(aff)
+        key = UnorderedPair(var, term_var)
+        _add_or_set!(quad.terms, key, coef)
+    end
+    return add_to_expression!(quad, var, aff.constant)
+end
+
+function add_to_expression!(quad::GenericQuadExpr{C,V},
+                            aff::GenericAffExpr{C,V},
+                            var::V) where {C,V}
+    return add_to_expression!(quad, var, aff)
+end
+
+function add_to_expression!(quad::GenericQuadExpr{C,V},
+                            lhs::GenericAffExpr{C,V},
+                            rhs::GenericAffExpr{C,V}) where {C,V}
+    lhs_length = length(linear_terms(lhs))
+    rhs_length = length(linear_terms(rhs))
+
+    # Quadratic terms
+    for (lhscoef, lhsvar) in linear_terms(lhs)
+        for (rhscoef, rhsvar) in linear_terms(rhs)
+            add_to_expression!(quad, lhscoef*rhscoef, lhsvar, rhsvar)
+        end
+    end
+
+    # Try to preallocate space for aff
+    cur = length(linear_terms(quad))
+    if !iszero(lhs.constant) && !iszero(rhs.constant)
+        sizehint!(quad.aff, cur + lhs_length + rhs_length)
+    elseif !iszero(lhs.constant)
+        sizehint!(quad.aff, cur + rhs_length)
+    elseif !iszero(rhs.constant)
+        sizehint!(quad.aff, cur + lhs_length)
+    end
+
+    # [LHS constant] * [RHS linear terms]
+    if !iszero(lhs.constant)
+        c = lhs.constant
+        for (rhscoef, rhsvar) in linear_terms(rhs)
+            add_to_expression!(quad.aff, c*rhscoef, rhsvar)
+        end
+    end
+
+    # [RHS constant] * [LHS linear terms]
+    if !iszero(rhs.constant)
+        c = rhs.constant
+        for (lhscoef, lhsvar) in linear_terms(lhs)
+            add_to_expression!(quad.aff, c*lhscoef, lhsvar)
+        end
+    end
+
+    quad.aff.constant += lhs.constant * rhs.constant
+
+    return quad
+end
+
+# With three factors.
+
+function add_to_expression!(quad::GenericQuadExpr{C,V}, new_coef::C,
+                            new_var1::V, new_var2::V) where {C,V}
+    # Node: OrderedDict updates the *key* as well. That is, if there was a
+    # previous value for UnorderedPair(new_var2, new_var1), it's key will now be
+    # UnorderedPair(new_var1, new_var2) (because these are defined as equal).
+    key = UnorderedPair(new_var1, new_var2)
+    _add_or_set!(quad.terms, key, new_coef)
+    return quad
 end
 
 
diff --git a/test/perf/matrix_product.jl b/test/perf/matrix_product.jl
new file mode 100644
index 00000000000..103b5a5b113
--- /dev/null
+++ b/test/perf/matrix_product.jl
@@ -0,0 +1,32 @@
+using LinearAlgebra
+using JuMP
+
+function test(n::Int, verbose::Bool)
+    model = Model()
+    a = rand(n, n)
+    @variable(model, x[1:n, 1:n])
+
+    # Compile it
+    x * a
+    a * x * a
+
+    elapsed = @elapsed x * a
+    if verbose
+        println("2D Matrix product `x * a`: $(elapsed)")
+    end
+
+    elapsed = @elapsed a * x * a
+    if verbose
+        println("2D Matrix product `a * x * a`: $(elapsed)")
+    end
+
+    return
+end
+
+# Compile the methods needed.
+test(2, false)
+
+for n in [10, 20, 50, 100]
+    println("n = $n")
+    test(n, true)
+end
diff --git a/test/perf/vector_speedtest.jl b/test/perf/vector_speedtest.jl
index cb630d62e94..bf669884303 100644
--- a/test/perf/vector_speedtest.jl
+++ b/test/perf/vector_speedtest.jl
@@ -1,36 +1,37 @@
+using LinearAlgebra
 using JuMP
 
 function test(n::Int)
-    m = Model()
-    a = rand(n,n)
-    @variable(m,x[1:n,1:n])
-    b = rand(n,n,n)
-    @variable(m,y[1:n,1:n,1:n])
+    model = Model()
+    a = rand(n, n)
+    @variable(model, x[1:n, 1:n])
+    b = rand(n, n, n)
+    @variable(model, y[1:n, 1:n, 1:n])
     c = rand(n)
-    @variable(m,z[1:n])
+    @variable(model, z[1:n])
 
     #initialize
-    @constraint(m,sum(c[i]*z[i] for i=1:n)<=0)
+    @constraint(model, sum(c[i] * z[i] for i=1:n) <= 0)
 
     #Vector
-    elapsed = @elapsed @constraint(m,sum(c[i]*z[i] for i=1:n)<=1)
+    elapsed = @elapsed @constraint(model, sum(c[i] * z[i] for i=1:n) <= 1)
     println("Vector with sum(): $(elapsed)")
 
-    elapsed = @elapsed @constraint(m,dot(c,z) <= 1)
+    elapsed = @elapsed @constraint(model, dot(c, z) <= 1)
     println("Vector with vecdot() : $(elapsed)")
 
     #2D Matrix
-    elapsed = @elapsed @constraint(m,sum(a[i,j]*x[i,j] for i=1:n,j=1:n)<=1)
+    elapsed = @elapsed @constraint(model, sum(a[i, j] * x[i, j] for i=1:n, j=1:n) <= 1)
     println("2D Matrix with sum(): $(elapsed)")
 
-    elapsed = @elapsed @constraint(m,dot(a,x)<=1)
+    elapsed = @elapsed @constraint(model, dot(a, x) <= 1)
     println("2D Matrix with bigvecdot(): $(elapsed)")
 
     #3D Matrix
-    elapsed = @elapsed @constraint(m,sum(b[i,j,k]*y[i,j,k] for i=1:n,j=1:n,k=1:n)<=1)
+    elapsed = @elapsed @constraint(model, sum(b[i, j, k] * y[i, j, k] for i=1:n, j=1:n, k=1:n) <= 1)
     println("3D Matrix with sum(): $(elapsed)")
 
-    elapsed = @elapsed @constraint(m,dot(b,y)<=1)
+    elapsed = @elapsed @constraint(model, dot(b, y) <= 1)
     println("3D Matrix with vecdot(): $(elapsed)")
     return 0
 end