diff --git a/src/nlp_expr.jl b/src/nlp_expr.jl
index 13f572c2c7f..fbe7e631753 100644
--- a/src/nlp_expr.jl
+++ b/src/nlp_expr.jl
@@ -405,10 +405,10 @@ end
 
 function _MA.operate!!(
     ::typeof(_MA.add_mul),
-    x::GenericNonlinearExpr,
-    ::AbstractArray,
+    ::GenericNonlinearExpr,
+    x::AbstractArray,
 )
-    return _throw_operator_error(_MA.add_mul)
+    return _throw_operator_error(_MA.add_mul, x)
 end
 
 """
diff --git a/src/operators.jl b/src/operators.jl
index 76a06742e2e..8150658fdb2 100644
--- a/src/operators.jl
+++ b/src/operators.jl
@@ -452,37 +452,53 @@ function LinearAlgebra.issymmetric(x::Matrix{T}) where {T<:_JuMPTypes}
     return true
 end
 
-function _throw_operator_error(::Union{typeof(+),typeof(_MA.add_mul)})
-    return error(
+function _throw_operator_error(
+    ::Union{typeof(+),typeof(_MA.add_mul)},
+    x::AbstractArray,
+)
+    msg =
         "Addition between an array and a JuMP scalar is not supported: " *
-        "instead of `x + y`, do `x .+ y` for element-wise addition.",
-    )
+        "instead of `x + y`, do `x .+ y` for element-wise addition."
+    if ndims(x) == 2 && size(x, 1) == size(x, 2)
+        msg *=
+            " If you are modifying the diagonal entries of a square matrix, " *
+            "do `x + y * LinearAlgebra.I(n)`, where `n` is the side length."
+    end
+    return error(msg)
 end
 
-function _throw_operator_error(::Union{typeof(-),typeof(_MA.sub_mul)})
-    return error(
-        "Subtraction between a Matrix and a JuMP scalar is not supported: " *
-        "instead of `x - y`, do `x .- y` for element-wise subtraction.",
-    )
+function _throw_operator_error(
+    ::Union{typeof(-),typeof(_MA.sub_mul)},
+    x::AbstractArray,
+)
+    msg =
+        "Subtraction between an array and a JuMP scalar is not supported: " *
+        "instead of `x - y`, do `x .- y` for element-wise subtraction."
+    if ndims(x) == 2 && size(x, 1) == size(x, 2)
+        msg *=
+            " If you are modifying the diagonal entries of a square matrix, " *
+            "do `x - y * LinearAlgebra.I(n)`, where `n` is the side length."
+    end
+    return error(msg)
 end
 
-Base.:+(::AbstractJuMPScalar, ::AbstractMatrix) = _throw_operator_error(+)
-Base.:+(::AbstractMatrix, ::AbstractJuMPScalar) = _throw_operator_error(+)
-Base.:-(::AbstractJuMPScalar, ::AbstractMatrix) = _throw_operator_error(-)
-Base.:-(::AbstractMatrix, ::AbstractJuMPScalar) = _throw_operator_error(-)
+Base.:+(::AbstractJuMPScalar, x::AbstractArray) = _throw_operator_error(+, x)
+Base.:+(x::AbstractArray, ::AbstractJuMPScalar) = _throw_operator_error(+, x)
+Base.:-(::AbstractJuMPScalar, x::AbstractArray) = _throw_operator_error(-, x)
+Base.:-(x::AbstractArray, ::AbstractJuMPScalar) = _throw_operator_error(-, x)
 
 function _MA.operate!!(
     op::Union{typeof(_MA.add_mul),typeof(_MA.sub_mul)},
-    ::AbstractArray,
+    x::AbstractArray,
     ::AbstractJuMPScalar,
 )
-    return _throw_operator_error(op)
+    return _throw_operator_error(op, x)
 end
 
 function _MA.operate!!(
     op::Union{typeof(_MA.add_mul),typeof(_MA.sub_mul)},
     ::AbstractJuMPScalar,
-    ::AbstractArray,
+    x::AbstractArray,
 )
-    return _throw_operator_error(op)
+    return _throw_operator_error(op, x)
 end
diff --git a/test/test_operator.jl b/test/test_operator.jl
index c79862aa733..58e8b03e3fb 100644
--- a/test/test_operator.jl
+++ b/test/test_operator.jl
@@ -628,16 +628,41 @@ end
 function test_matrix_abstractscalar_add()
     model = Model()
     @variable(model, x)
-    A = rand(Float64, 2, 2)
+    A = rand(Float64, 3, 2)
+    B = rand(Float64, 3)
     err_add = ErrorException(
         "Addition between an array and a JuMP scalar is not supported: " *
         "instead of `x + y`, do `x .+ y` for element-wise addition.",
     )
     err_sub = ErrorException(
-        "Subtraction between a Matrix and a JuMP scalar is not supported: " *
+        "Subtraction between an array and a JuMP scalar is not supported: " *
         "instead of `x - y`, do `x .- y` for element-wise subtraction.",
     )
-    for lhs in (A, A'), rhs in (x, 1.0 * x, x^2, sin(x))
+    for lhs in (A, A', B, B'), rhs in (x, 1.0 * x, x^2, sin(x))
+        @test_throws(err_add, lhs + rhs)
+        @test_throws(err_add, rhs + lhs)
+        @test_throws(err_add, @expression(model, lhs + rhs))
+        @test_throws(err_add, @expression(model, rhs + lhs))
+        @test_throws(err_sub, lhs - rhs)
+        @test_throws(err_sub, rhs - lhs)
+        @test_throws(err_sub, @expression(model, lhs - rhs))
+        @test_throws(err_sub, @expression(model, rhs - lhs))
+    end
+    C = rand(Float64, 2, 2)
+    err_add = ErrorException(
+        "Addition between an array and a JuMP scalar is not supported: " *
+        "instead of `x + y`, do `x .+ y` for element-wise addition." *
+        " If you are modifying the diagonal entries of a square matrix, " *
+        "do `x + y * LinearAlgebra.I(n)`, where `n` is the side length.",
+    )
+    err_sub = ErrorException(
+        "Subtraction between an array and a JuMP scalar is not supported: " *
+        "instead of `x - y`, do `x .- y` for element-wise subtraction." *
+        " If you are modifying the diagonal entries of a square matrix, " *
+        "do `x - y * LinearAlgebra.I(n)`, where `n` is the side length.",
+
+    )
+    for lhs in (C, C'), rhs in (x, 1.0 * x, x^2, sin(x))
         @test_throws(err_add, lhs + rhs)
         @test_throws(err_add, rhs + lhs)
         @test_throws(err_add, @expression(model, lhs + rhs))