diff --git a/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs
index 825b2877..a2cf4afa 100644
--- a/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs
@@ -1152,43 +1152,52 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// <exception cref="ArgumentOutOfRangeException">If the two matrices don't have the same dimensions.</exception>
         protected override void DoAdd(Matrix<Complex> other, Matrix<Complex> result)
         {
-            result.Clear();
-
+            var sparseOther = other as SparseMatrix;
             var sparseResult = result as SparseMatrix;
-            if (sparseResult == null)
+            if (sparseOther == null || sparseResult == null)
             {
-                for (var i = 0; i < other.RowCount; i++)
-                {
-                    // Get the begin / end index for the current row
-                    var startIndex = _rowIndex[i];
-                    var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount;
+                base.DoAdd(other, result);
+                return;
+            }
 
-                    for (var j = startIndex; j < endIndex; j++)
-                    {
-                        var resVal = _nonZeroValues[j] + other.At(i, _columnIndices[j]);
-                        if (resVal != 0.0)
-                        {
-                            result.At(i, _columnIndices[j], resVal);
-                        }
-                    }
+            if (ReferenceEquals(this, other))
+            {
+                if (!ReferenceEquals(this, result))
+                {
+                    CopyTo(result);
                 }
+
+                Control.LinearAlgebraProvider.ScaleArray(2.0, _nonZeroValues, _nonZeroValues);
+                return;
+            }
+
+            SparseMatrix left;
+
+            if (ReferenceEquals(sparseOther, sparseResult))
+            {
+                left = this;
+            }
+            else if (ReferenceEquals(this, sparseResult))
+            {
+                left = sparseOther;
             }
             else
             {
-                for (var i = 0; i < other.RowCount; i++)
-                {
-                    // Get the begin / end index for the current row
-                    var startIndex = _rowIndex[i];
-                    var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount;
+                CopyTo(sparseResult);
+                left = sparseOther;
+            }
 
-                    for (var j = startIndex; j < endIndex; j++)
-                    {
-                        var resVal = _nonZeroValues[j] + other.At(i, _columnIndices[j]);
-                        if (resVal != 0.0)
-                        {
-                            sparseResult.SetValueAt(i, _columnIndices[j], resVal);
-                        }
-                    }
+            for (var i = 0; i < left.RowCount; i++)
+            {
+                // Get the begin / end index for the current row
+                var startIndex = left._rowIndex[i];
+                var endIndex = i < left._rowIndex.Length - 1 ? left._rowIndex[i + 1] : left.NonZerosCount;
+
+                for (var j = startIndex; j < endIndex; j++)
+                {
+                    var columnIndex = left._columnIndices[j];
+                    var resVal = left._nonZeroValues[j] + result.At(i, columnIndex);
+                    result.At(i, columnIndex, resVal);
                 }
             }
         }
@@ -1202,30 +1211,46 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// <exception cref="ArgumentOutOfRangeException">If the two matrices don't have the same dimensions.</exception>
         protected override void DoSubtract(Matrix<Complex> other, Matrix<Complex> result)
         {
-            result.Clear();
-
+            var sparseOther = other as SparseMatrix;
             var sparseResult = result as SparseMatrix;
-            if (sparseResult == null)
+            if (sparseOther == null || sparseResult == null)
+            {
+                base.DoSubtract(other, result);
+                return;
+            }
+
+            if (ReferenceEquals(this, other))
             {
-                for (var i = 0; i < other.RowCount; i++)
+                result.Clear();
+                return;
+            }
+
+            if (ReferenceEquals(this, sparseResult))
+            {
+                for (var i = 0; i < sparseOther.RowCount; i++)
                 {
                     // Get the begin / end index for the current row
-                    var startIndex = _rowIndex[i];
-                    var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount;
+                    var startIndex = sparseOther._rowIndex[i];
+                    var endIndex = i < sparseOther._rowIndex.Length - 1 ? sparseOther._rowIndex[i + 1] : sparseOther.NonZerosCount;
 
                     for (var j = startIndex; j < endIndex; j++)
                     {
-                        var resVal = _nonZeroValues[j] - other.At(i, _columnIndices[j]);
-                        if (resVal != 0.0)
-                        {
-                            result.At(i, _columnIndices[j], resVal);
-                        }
+                        var columnIndex = sparseOther._columnIndices[j];
+                        var resVal = sparseResult.At(i, columnIndex) - sparseOther._nonZeroValues[j];
+                        result.At(i, columnIndex, resVal);
                     }
                 }
             }
             else
             {
-                for (var i = 0; i < other.RowCount; i++)
+                if (!ReferenceEquals(sparseOther, sparseResult))
+                {
+                    sparseOther.CopyTo(sparseResult);
+                }
+
+                sparseResult.Negate(sparseResult);
+
+                for (var i = 0; i < RowCount; i++)
                 {
                     // Get the begin / end index for the current row
                     var startIndex = _rowIndex[i];
@@ -1233,11 +1258,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
 
                     for (var j = startIndex; j < endIndex; j++)
                     {
-                        var resVal = _nonZeroValues[j] - other.At(i, _columnIndices[j]);
-                        if (resVal != 0.0)
-                        {
-                            sparseResult.SetValueAt(i, _columnIndices[j], resVal);
-                        }
+                        var columnIndex = _columnIndices[j];
+                        var resVal = sparseResult.At(i, columnIndex) + _nonZeroValues[j];
+                        result.At(i, columnIndex, resVal);
                     }
                 }
             }
diff --git a/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs
index 1e6d03c6..f5c3c4c0 100644
--- a/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs
@@ -1152,43 +1152,52 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// <exception cref="ArgumentOutOfRangeException">If the two matrices don't have the same dimensions.</exception>
         protected override void DoAdd(Matrix<Complex32> other, Matrix<Complex32> result)
         {
-            result.Clear();
-
+            var sparseOther = other as SparseMatrix;
             var sparseResult = result as SparseMatrix;
-            if (sparseResult == null)
+            if (sparseOther == null || sparseResult == null)
             {
-                for (var i = 0; i < other.RowCount; i++)
-                {
-                    // Get the begin / end index for the current row
-                    var startIndex = _rowIndex[i];
-                    var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount;
+                base.DoAdd(other, result);
+                return;
+            }
 
-                    for (var j = startIndex; j < endIndex; j++)
-                    {
-                        var resVal = _nonZeroValues[j] + other.At(i, _columnIndices[j]);
-                        if (resVal != 0.0f)
-                        {
-                            result.At(i, _columnIndices[j], resVal);
-                        }
-                    }
+            if (ReferenceEquals(this, other))
+            {
+                if (!ReferenceEquals(this, result))
+                {
+                    CopyTo(result);
                 }
+
+                Control.LinearAlgebraProvider.ScaleArray(2.0f, _nonZeroValues, _nonZeroValues);
+                return;
+            }
+
+            SparseMatrix left;
+
+            if (ReferenceEquals(sparseOther, sparseResult))
+            {
+                left = this;
+            }
+            else if (ReferenceEquals(this, sparseResult))
+            {
+                left = sparseOther;
             }
             else
             {
-                for (var i = 0; i < other.RowCount; i++)
-                {
-                    // Get the begin / end index for the current row
-                    var startIndex = _rowIndex[i];
-                    var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount;
+                CopyTo(sparseResult);
+                left = sparseOther;
+            }
 
-                    for (var j = startIndex; j < endIndex; j++)
-                    {
-                        var resVal = _nonZeroValues[j] + other.At(i, _columnIndices[j]);
-                        if (resVal != 0.0f)
-                        {
-                            sparseResult.SetValueAt(i, _columnIndices[j], resVal);
-                        }
-                    }
+            for (var i = 0; i < left.RowCount; i++)
+            {
+                // Get the begin / end index for the current row
+                var startIndex = left._rowIndex[i];
+                var endIndex = i < left._rowIndex.Length - 1 ? left._rowIndex[i + 1] : left.NonZerosCount;
+
+                for (var j = startIndex; j < endIndex; j++)
+                {
+                    var columnIndex = left._columnIndices[j];
+                    var resVal = left._nonZeroValues[j] + result.At(i, columnIndex);
+                    result.At(i, columnIndex, resVal);
                 }
             }
         }
@@ -1202,30 +1211,46 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// <exception cref="ArgumentOutOfRangeException">If the two matrices don't have the same dimensions.</exception>
         protected override void DoSubtract(Matrix<Complex32> other, Matrix<Complex32> result)
         {
-            result.Clear();
-
+            var sparseOther = other as SparseMatrix;
             var sparseResult = result as SparseMatrix;
-            if (sparseResult == null)
+            if (sparseOther == null || sparseResult == null)
+            {
+                base.DoSubtract(other, result);
+                return;
+            }
+
+            if (ReferenceEquals(this, other))
             {
-                for (var i = 0; i < other.RowCount; i++)
+                result.Clear();
+                return;
+            }
+
+            if (ReferenceEquals(this, sparseResult))
+            {
+                for (var i = 0; i < sparseOther.RowCount; i++)
                 {
                     // Get the begin / end index for the current row
-                    var startIndex = _rowIndex[i];
-                    var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount;
+                    var startIndex = sparseOther._rowIndex[i];
+                    var endIndex = i < sparseOther._rowIndex.Length - 1 ? sparseOther._rowIndex[i + 1] : sparseOther.NonZerosCount;
 
                     for (var j = startIndex; j < endIndex; j++)
                     {
-                        var resVal = _nonZeroValues[j] - other.At(i, _columnIndices[j]);
-                        if (resVal != 0.0f)
-                        {
-                            result.At(i, _columnIndices[j], resVal);
-                        }
+                        var columnIndex = sparseOther._columnIndices[j];
+                        var resVal = sparseResult.At(i, columnIndex) - sparseOther._nonZeroValues[j];
+                        result.At(i, columnIndex, resVal);
                     }
                 }
             }
             else
             {
-                for (var i = 0; i < other.RowCount; i++)
+                if (!ReferenceEquals(sparseOther, sparseResult))
+                {
+                    sparseOther.CopyTo(sparseResult);
+                }
+
+                sparseResult.Negate(sparseResult);
+
+                for (var i = 0; i < RowCount; i++)
                 {
                     // Get the begin / end index for the current row
                     var startIndex = _rowIndex[i];
@@ -1233,11 +1258,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
 
                     for (var j = startIndex; j < endIndex; j++)
                     {
-                        var resVal = _nonZeroValues[j] - other.At(i, _columnIndices[j]);
-                        if (resVal != 0.0f)
-                        {
-                            sparseResult.SetValueAt(i, _columnIndices[j], resVal);
-                        }
+                        var columnIndex = _columnIndices[j];
+                        var resVal = sparseResult.At(i, columnIndex) + _nonZeroValues[j];
+                        result.At(i, columnIndex, resVal);
                     }
                 }
             }
diff --git a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs
index 43a10b57..1a1def63 100644
--- a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs
@@ -1163,43 +1163,52 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <exception cref="ArgumentOutOfRangeException">If the two matrices don't have the same dimensions.</exception>
         protected override void DoAdd(Matrix<double> other, Matrix<double> result)
         {
-            result.Clear();
-
+            var sparseOther = other as SparseMatrix;
             var sparseResult = result as SparseMatrix;
-            if (sparseResult == null)
+            if (sparseOther == null || sparseResult == null)
             {
-                for (var i = 0; i < other.RowCount; i++)
-                {
-                    // Get the begin / end index for the current row
-                    var startIndex = _rowIndex[i];
-                    var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount;
+                base.DoAdd(other, result);
+                return;
+            }
 
-                    for (var j = startIndex; j < endIndex; j++)
-                    {
-                        var resVal = _nonZeroValues[j] + other.At(i, _columnIndices[j]);
-                        if (resVal != 0.0)
-                        {
-                            result.At(i, _columnIndices[j], resVal);
-                        }
-                    }
+            if (ReferenceEquals(this, other))
+            {
+                if (!ReferenceEquals(this, result))
+                {
+                    CopyTo(result);
                 }
+
+                Control.LinearAlgebraProvider.ScaleArray(2.0, _nonZeroValues, _nonZeroValues);
+                return;
+            }
+
+            SparseMatrix left;
+
+            if (ReferenceEquals(sparseOther, sparseResult))
+            {
+                left = this;
+            }
+            else if (ReferenceEquals(this, sparseResult))
+            {
+                left = sparseOther;
             }
             else
             {
-                for (var i = 0; i < other.RowCount; i++)
-                {
-                    // Get the begin / end index for the current row
-                    var startIndex = _rowIndex[i];
-                    var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount;
+                CopyTo(sparseResult);
+                left = sparseOther;
+            }
 
-                    for (var j = startIndex; j < endIndex; j++)
-                    {
-                        var resVal = _nonZeroValues[j] + other.At(i, _columnIndices[j]);
-                        if (resVal != 0.0)
-                        {
-                            sparseResult.SetValueAt(i, _columnIndices[j], resVal);
-                        }
-                    }
+            for (var i = 0; i < left.RowCount; i++)
+            {
+                // Get the begin / end index for the current row
+                var startIndex = left._rowIndex[i];
+                var endIndex = i < left._rowIndex.Length - 1 ? left._rowIndex[i + 1] : left.NonZerosCount;
+
+                for (var j = startIndex; j < endIndex; j++)
+                {
+                    var columnIndex = left._columnIndices[j];
+                    var resVal = left._nonZeroValues[j] + result.At(i, columnIndex);
+                    result.At(i, columnIndex, resVal);
                 }
             }
         }
@@ -1213,30 +1222,46 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <exception cref="ArgumentOutOfRangeException">If the two matrices don't have the same dimensions.</exception>
         protected override void DoSubtract(Matrix<double> other, Matrix<double> result)
         {
-            result.Clear();
-
+            var sparseOther = other as SparseMatrix;
             var sparseResult = result as SparseMatrix;
-            if (sparseResult == null)
+            if (sparseOther == null || sparseResult == null)
+            {
+                base.DoSubtract(other, result);
+                return;
+            }
+
+            if (ReferenceEquals(this, other))
             {
-                for (var i = 0; i < other.RowCount; i++)
+                result.Clear();
+                return;
+            }
+
+            if (ReferenceEquals(this, sparseResult))
+            {
+                for (var i = 0; i < sparseOther.RowCount; i++)
                 {
                     // Get the begin / end index for the current row
-                    var startIndex = _rowIndex[i];
-                    var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount;
+                    var startIndex = sparseOther._rowIndex[i];
+                    var endIndex = i < sparseOther._rowIndex.Length - 1 ? sparseOther._rowIndex[i + 1] : sparseOther.NonZerosCount;
 
                     for (var j = startIndex; j < endIndex; j++)
                     {
-                        var resVal = _nonZeroValues[j] - other.At(i, _columnIndices[j]);
-                        if (resVal != 0.0)
-                        {
-                            result.At(i, _columnIndices[j], resVal);
-                        }
+                        var columnIndex = sparseOther._columnIndices[j];
+                        var resVal = sparseResult.At(i, columnIndex) - sparseOther._nonZeroValues[j];
+                        result.At(i, columnIndex, resVal);
                     }
                 }
             }
             else
             {
-                for (var i = 0; i < other.RowCount; i++)
+                if (!ReferenceEquals(sparseOther, sparseResult))
+                {
+                    sparseOther.CopyTo(sparseResult);
+                }
+
+                sparseResult.Negate(sparseResult);
+
+                for (var i = 0; i < RowCount; i++)
                 {
                     // Get the begin / end index for the current row
                     var startIndex = _rowIndex[i];
@@ -1244,11 +1269,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double
 
                     for (var j = startIndex; j < endIndex; j++)
                     {
-                        var resVal = _nonZeroValues[j] - other.At(i, _columnIndices[j]);
-                        if (resVal != 0.0)
-                        {
-                            sparseResult.SetValueAt(i, _columnIndices[j], resVal);
-                        }
+                        var columnIndex = _columnIndices[j];
+                        var resVal = sparseResult.At(i, columnIndex) + _nonZeroValues[j];
+                        result.At(i, columnIndex, resVal);
                     }
                 }
             }
diff --git a/src/Numerics/LinearAlgebra/Generic/Matrix.Arithmetic.cs b/src/Numerics/LinearAlgebra/Generic/Matrix.Arithmetic.cs
index c0ab7de3..9525ced1 100644
--- a/src/Numerics/LinearAlgebra/Generic/Matrix.Arithmetic.cs
+++ b/src/Numerics/LinearAlgebra/Generic/Matrix.Arithmetic.cs
@@ -66,7 +66,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
             }
 
             var result = CreateMatrix(RowCount, ColumnCount);
-            Add(other, result);
+            DoAdd(other, result);
             return result;
         }
 
@@ -177,7 +177,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
             }
 
             var result = CreateMatrix(RowCount, ColumnCount);
-            Multiply(scalar, result);
+            DoMultiply(scalar, result);
             return result;
         }
 
@@ -238,7 +238,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
             }
 
             var result = CreateMatrix(RowCount, ColumnCount);
-            Divide(scalar, result);
+            DoDivide(scalar, result);
             return result;
         }
 
@@ -577,7 +577,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         public virtual Matrix<T> Negate()
         {
             var result = CreateMatrix(RowCount, ColumnCount);
-            Negate(result);
+            DoNegate(result);
             return result;
         }
 
diff --git a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs
index f2efa00c..945bf895 100644
--- a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs
@@ -1146,43 +1146,52 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// <exception cref="ArgumentOutOfRangeException">If the two matrices don't have the same dimensions.</exception>
         protected override void DoAdd(Matrix<float> other, Matrix<float> result)
         {
-            result.Clear();
-
+            var sparseOther = other as SparseMatrix;
             var sparseResult = result as SparseMatrix;
-            if (sparseResult == null)
+            if (sparseOther == null || sparseResult == null)
             {
-                for (var i = 0; i < other.RowCount; i++)
-                {
-                    // Get the begin / end index for the current row
-                    var startIndex = _rowIndex[i];
-                    var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount;
+                base.DoAdd(other, result);
+                return;
+            }
 
-                    for (var j = startIndex; j < endIndex; j++)
-                    {
-                        var resVal = _nonZeroValues[j] + other.At(i, _columnIndices[j]);
-                        if (resVal != 0.0f)
-                        {
-                            result.At(i, _columnIndices[j], resVal);
-                        }
-                    }
+            if (ReferenceEquals(this, other))
+            {
+                if (!ReferenceEquals(this, result))
+                {
+                    CopyTo(result);
                 }
+
+                Control.LinearAlgebraProvider.ScaleArray(2.0f, _nonZeroValues, _nonZeroValues);
+                return;
+            }
+
+            SparseMatrix left;
+
+            if (ReferenceEquals(sparseOther, sparseResult))
+            {
+                left = this;
+            }
+            else if (ReferenceEquals(this, sparseResult))
+            {
+                left = sparseOther;
             }
             else
             {
-                for (var i = 0; i < other.RowCount; i++)
-                {
-                    // Get the begin / end index for the current row
-                    var startIndex = _rowIndex[i];
-                    var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount;
+                CopyTo(sparseResult);
+                left = sparseOther;
+            }
 
-                    for (var j = startIndex; j < endIndex; j++)
-                    {
-                        var resVal = _nonZeroValues[j] + other.At(i, _columnIndices[j]);
-                        if (resVal != 0.0f)
-                        {
-                            sparseResult.SetValueAt(i, _columnIndices[j], resVal);
-                        }
-                    }
+            for (var i = 0; i < left.RowCount; i++)
+            {
+                // Get the begin / end index for the current row
+                var startIndex = left._rowIndex[i];
+                var endIndex = i < left._rowIndex.Length - 1 ? left._rowIndex[i + 1] : left.NonZerosCount;
+
+                for (var j = startIndex; j < endIndex; j++)
+                {
+                    var columnIndex = left._columnIndices[j];
+                    var resVal = left._nonZeroValues[j] + result.At(i, columnIndex);
+                    result.At(i, columnIndex, resVal);
                 }
             }
         }
@@ -1196,30 +1205,46 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// <exception cref="ArgumentOutOfRangeException">If the two matrices don't have the same dimensions.</exception>
         protected override void DoSubtract(Matrix<float> other, Matrix<float> result)
         {
-            result.Clear();
-
+            var sparseOther = other as SparseMatrix;
             var sparseResult = result as SparseMatrix;
-            if (sparseResult == null)
+            if (sparseOther == null || sparseResult == null)
+            {
+                base.DoSubtract(other, result);
+                return;
+            }
+
+            if (ReferenceEquals(this, other))
+            {
+                result.Clear();
+                return;
+            }
+
+            if (ReferenceEquals(this, sparseResult))
             {
-                for (var i = 0; i < other.RowCount; i++)
+                for (var i = 0; i < sparseOther.RowCount; i++)
                 {
                     // Get the begin / end index for the current row
-                    var startIndex = _rowIndex[i];
-                    var endIndex = i < _rowIndex.Length - 1 ? _rowIndex[i + 1] : NonZerosCount;
+                    var startIndex = sparseOther._rowIndex[i];
+                    var endIndex = i < sparseOther._rowIndex.Length - 1 ? sparseOther._rowIndex[i + 1] : sparseOther.NonZerosCount;
 
                     for (var j = startIndex; j < endIndex; j++)
                     {
-                        var resVal = _nonZeroValues[j] - other.At(i, _columnIndices[j]);
-                        if (resVal != 0.0f)
-                        {
-                            result.At(i, _columnIndices[j], resVal);
-                        }
+                        var columnIndex = sparseOther._columnIndices[j];
+                        var resVal = sparseResult.At(i, columnIndex) - sparseOther._nonZeroValues[j];
+                        result.At(i, columnIndex, resVal);
                     }
                 }
             }
             else
             {
-                for (var i = 0; i < other.RowCount; i++)
+                if (!ReferenceEquals(sparseOther, sparseResult))
+                {
+                    sparseOther.CopyTo(sparseResult);
+                }
+
+                sparseResult.Negate(sparseResult);
+
+                for (var i = 0; i < RowCount; i++)
                 {
                     // Get the begin / end index for the current row
                     var startIndex = _rowIndex[i];
@@ -1227,16 +1252,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single
 
                     for (var j = startIndex; j < endIndex; j++)
                     {
-                        var resVal = _nonZeroValues[j] - other.At(i, _columnIndices[j]);
-                        if (resVal != 0.0f)
-                        {
-                            sparseResult.SetValueAt(i, _columnIndices[j], resVal);
-                        }
+                        var columnIndex = _columnIndices[j];
+                        var resVal = sparseResult.At(i, columnIndex) + _nonZeroValues[j];
+                        result.At(i, columnIndex, resVal);
                     }
                 }
             }
         }
-
+        
         /// <summary>
         /// Multiplies each element of the matrix by a scalar and places results into the result matrix.
         /// </summary>
diff --git a/src/UnitTests/LinearAlgebraTests/Complex/SparseMatrixTests.cs b/src/UnitTests/LinearAlgebraTests/Complex/SparseMatrixTests.cs
index 610026a2..886e7022 100644
--- a/src/UnitTests/LinearAlgebraTests/Complex/SparseMatrixTests.cs
+++ b/src/UnitTests/LinearAlgebraTests/Complex/SparseMatrixTests.cs
@@ -209,5 +209,93 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex
 
             Assert.AreEqual(matrix.NonZerosCount, nonzero);
         }
+
+        /// <summary>
+        /// Test whether order matters when adding sparse matrices.
+        /// </summary>
+        [Test]
+        public void CanAddSparseMatricesBothWays()
+        {
+            var m1 = new SparseMatrix(1, 3);
+            var m2 = new SparseMatrix(new Complex[,] { { 0, 1, 1 } });
+            var sum1 = m1 + m2;
+            var sum2 = m2 + m1;
+            Assert.IsTrue(sum1.Equals(m2));
+            Assert.IsTrue(sum1.Equals(sum2));
+
+            var sparseResult = new SparseMatrix(1, 3);
+            sparseResult.Add(m2, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(sum1));
+
+            sparseResult = new SparseMatrix(new Complex[,] { { 0, 1, 1 } });
+            sparseResult.Add(m1, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(sum1));
+
+            sparseResult = new SparseMatrix(new Complex[,] { { 0, 1, 1 } });
+            m1.Add(sparseResult, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(sum1));
+
+            sparseResult = new SparseMatrix(new Complex[,] { { 0, 1, 1 } });
+            sparseResult.Add(sparseResult, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(2 * sum1));
+
+            var denseResult = new DenseMatrix(1, 3);
+            denseResult.Add(m2, denseResult);
+            Assert.IsTrue(denseResult.Equals(sum1));
+
+            denseResult = new DenseMatrix(new Complex[,] { { 0, 1, 1 } });
+            denseResult.Add(m1, denseResult);
+            Assert.IsTrue(denseResult.Equals(sum1));
+
+            var m3 = new DenseMatrix(new Complex[,] { { 0, 1, 1 } });
+            var sum3 = m1 + m3;
+            var sum4 = m3 + m1;
+            Assert.IsTrue(sum3.Equals(m3));
+            Assert.IsTrue(sum3.Equals(sum4));
+        }
+
+        /// <summary>
+        /// Test whether order matters when subtracting sparse matrices.
+        /// </summary>
+        [Test]
+        public void CanSubtractSparseMatricesBothWays()
+        {
+            var m1 = new SparseMatrix(1, 3);
+            var m2 = new SparseMatrix(new Complex[,] { { 0, 1, 1 } });
+            var diff1 = m1 - m2;
+            var diff2 = m2 - m1;
+            Assert.IsTrue(diff1.Equals(m2.Negate()));
+            Assert.IsTrue(diff1.Equals(diff2.Negate()));
+
+            var sparseResult = new SparseMatrix(1, 3);
+            sparseResult.Subtract(m2, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(diff1));
+
+            sparseResult = new SparseMatrix(new Complex[,] { { 0, 1, 1 } });
+            sparseResult.Subtract(m1, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(diff2));
+
+            sparseResult = new SparseMatrix(new Complex[,] { { 0, 1, 1 } });
+            m1.Subtract(sparseResult, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(diff1));
+
+            sparseResult = new SparseMatrix(new Complex[,] { { 0, 1, 1 } });
+            sparseResult.Subtract(sparseResult, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(0 * diff1));
+
+            var denseResult = new DenseMatrix(1, 3);
+            denseResult.Subtract(m2, denseResult);
+            Assert.IsTrue(denseResult.Equals(diff1));
+
+            denseResult = new DenseMatrix(new Complex[,] { { 0, 1, 1 } });
+            denseResult.Subtract(m1, denseResult);
+            Assert.IsTrue(denseResult.Equals(diff2));
+
+            var m3 = new DenseMatrix(new Complex[,] { { 0, 1, 1 } });
+            var diff3 = m1 - m3;
+            var diff4 = m3 - m1;
+            Assert.IsTrue(diff3.Equals(m3.Negate()));
+            Assert.IsTrue(diff3.Equals(diff4.Negate()));
+        }
     }
 }
diff --git a/src/UnitTests/LinearAlgebraTests/Complex32/SparseMatrixTests.cs b/src/UnitTests/LinearAlgebraTests/Complex32/SparseMatrixTests.cs
index d204b878..62c63a73 100644
--- a/src/UnitTests/LinearAlgebraTests/Complex32/SparseMatrixTests.cs
+++ b/src/UnitTests/LinearAlgebraTests/Complex32/SparseMatrixTests.cs
@@ -209,5 +209,93 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32
 
             Assert.AreEqual(matrix.NonZerosCount, nonzero);
         }
+
+        /// <summary>
+        /// Test whether order matters when adding sparse matrices.
+        /// </summary>
+        [Test]
+        public void CanAddSparseMatricesBothWays()
+        {
+            var m1 = new SparseMatrix(1, 3);
+            var m2 = new SparseMatrix(new Complex32[,] { { 0, 1, 1 } });
+            var sum1 = m1 + m2;
+            var sum2 = m2 + m1;
+            Assert.IsTrue(sum1.Equals(m2));
+            Assert.IsTrue(sum1.Equals(sum2));
+
+            var sparseResult = new SparseMatrix(1, 3);
+            sparseResult.Add(m2, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(sum1));
+
+            sparseResult = new SparseMatrix(new Complex32[,] { { 0, 1, 1 } });
+            sparseResult.Add(m1, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(sum1));
+
+            sparseResult = new SparseMatrix(new Complex32[,] { { 0, 1, 1 } });
+            m1.Add(sparseResult, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(sum1));
+
+            sparseResult = new SparseMatrix(new Complex32[,] { { 0, 1, 1 } });
+            sparseResult.Add(sparseResult, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(2 * sum1));
+
+            var denseResult = new DenseMatrix(1, 3);
+            denseResult.Add(m2, denseResult);
+            Assert.IsTrue(denseResult.Equals(sum1));
+
+            denseResult = new DenseMatrix(new Complex32[,] { { 0, 1, 1 } });
+            denseResult.Add(m1, denseResult);
+            Assert.IsTrue(denseResult.Equals(sum1));
+
+            var m3 = new DenseMatrix(new Complex32[,] { { 0, 1, 1 } });
+            var sum3 = m1 + m3;
+            var sum4 = m3 + m1;
+            Assert.IsTrue(sum3.Equals(m3));
+            Assert.IsTrue(sum3.Equals(sum4));
+        }
+
+        /// <summary>
+        /// Test whether order matters when subtracting sparse matrices.
+        /// </summary>
+        [Test]
+        public void CanSubtractSparseMatricesBothWays()
+        {
+            var m1 = new SparseMatrix(1, 3);
+            var m2 = new SparseMatrix(new Complex32[,] { { 0, 1, 1 } });
+            var diff1 = m1 - m2;
+            var diff2 = m2 - m1;
+            Assert.IsTrue(diff1.Equals(m2.Negate()));
+            Assert.IsTrue(diff1.Equals(diff2.Negate()));
+
+            var sparseResult = new SparseMatrix(1, 3);
+            sparseResult.Subtract(m2, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(diff1));
+
+            sparseResult = new SparseMatrix(new Complex32[,] { { 0, 1, 1 } });
+            sparseResult.Subtract(m1, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(diff2));
+
+            sparseResult = new SparseMatrix(new Complex32[,] { { 0, 1, 1 } });
+            m1.Subtract(sparseResult, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(diff1));
+
+            sparseResult = new SparseMatrix(new Complex32[,] { { 0, 1, 1 } });
+            sparseResult.Subtract(sparseResult, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(0 * diff1));
+
+            var denseResult = new DenseMatrix(1, 3);
+            denseResult.Subtract(m2, denseResult);
+            Assert.IsTrue(denseResult.Equals(diff1));
+
+            denseResult = new DenseMatrix(new Complex32[,] { { 0, 1, 1 } });
+            denseResult.Subtract(m1, denseResult);
+            Assert.IsTrue(denseResult.Equals(diff2));
+
+            var m3 = new DenseMatrix(new Complex32[,] { { 0, 1, 1 } });
+            var diff3 = m1 - m3;
+            var diff4 = m3 - m1;
+            Assert.IsTrue(diff3.Equals(m3.Negate()));
+            Assert.IsTrue(diff3.Equals(diff4.Negate()));
+        }
     }
 }
diff --git a/src/UnitTests/LinearAlgebraTests/Double/SparseMatrixTests.cs b/src/UnitTests/LinearAlgebraTests/Double/SparseMatrixTests.cs
index 027446e4..bb7bd52a 100644
--- a/src/UnitTests/LinearAlgebraTests/Double/SparseMatrixTests.cs
+++ b/src/UnitTests/LinearAlgebraTests/Double/SparseMatrixTests.cs
@@ -208,5 +208,93 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
 
             Assert.AreEqual(matrix.NonZerosCount, nonzero);
         }
+
+        /// <summary>
+        /// Test whether order matters when adding sparse matrices.
+        /// </summary>
+        [Test]
+        public void CanAddSparseMatricesBothWays()
+        {
+            var m1 = new SparseMatrix(1, 3);
+            var m2 = new SparseMatrix(new double[,] { { 0, 1, 1 } });
+            var sum1 = m1 + m2;
+            var sum2 = m2 + m1;
+            Assert.IsTrue(sum1.Equals(m2));
+            Assert.IsTrue(sum1.Equals(sum2));
+
+            var sparseResult = new SparseMatrix(1, 3);
+            sparseResult.Add(m2, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(sum1));
+
+            sparseResult = new SparseMatrix(new double[,] { { 0, 1, 1 } });
+            sparseResult.Add(m1, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(sum1));
+
+            sparseResult = new SparseMatrix(new double[,] { { 0, 1, 1 } });
+            m1.Add(sparseResult, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(sum1));
+
+            sparseResult = new SparseMatrix(new double[,] { { 0, 1, 1 } });
+            sparseResult.Add(sparseResult, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(2 * sum1));
+            
+            var denseResult = new DenseMatrix(1, 3);
+            denseResult.Add(m2, denseResult);
+            Assert.IsTrue(denseResult.Equals(sum1));
+
+            denseResult = new DenseMatrix(new double[,] { { 0, 1, 1 } });
+            denseResult.Add(m1, denseResult);
+            Assert.IsTrue(denseResult.Equals(sum1));
+
+            var m3 = new DenseMatrix(new double[,] { { 0, 1, 1 } });
+            var sum3 = m1 + m3;
+            var sum4 = m3 + m1;
+            Assert.IsTrue(sum3.Equals(m3));
+            Assert.IsTrue(sum3.Equals(sum4));
+        }
+
+        /// <summary>
+        /// Test whether order matters when subtracting sparse matrices.
+        /// </summary>
+        [Test]
+        public void CanSubtractSparseMatricesBothWays()
+        {
+            var m1 = new SparseMatrix(1, 3);
+            var m2 = new SparseMatrix(new double[,] { { 0, 1, 1 } });
+            var diff1 = m1 - m2;
+            var diff2 = m2 - m1;
+            Assert.IsTrue(diff1.Equals(m2.Negate()));
+            Assert.IsTrue(diff1.Equals(diff2.Negate()));
+
+            var sparseResult = new SparseMatrix(1, 3);
+            sparseResult.Subtract(m2, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(diff1));
+
+            sparseResult = new SparseMatrix(new double[,] { { 0, 1, 1 } });
+            sparseResult.Subtract(m1, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(diff2));
+
+            sparseResult = new SparseMatrix(new double[,] { { 0, 1, 1 } });
+            m1.Subtract(sparseResult, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(diff1));
+
+            sparseResult = new SparseMatrix(new double[,] { { 0, 1, 1 } });
+            sparseResult.Subtract(sparseResult, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(0 * diff1));
+
+            var denseResult = new DenseMatrix(1, 3);
+            denseResult.Subtract(m2, denseResult);
+            Assert.IsTrue(denseResult.Equals(diff1));
+
+            denseResult = new DenseMatrix(new double[,] { { 0, 1, 1 } });
+            denseResult.Subtract(m1, denseResult);
+            Assert.IsTrue(denseResult.Equals(diff2));
+
+            var m3 = new DenseMatrix(new double[,] { { 0, 1, 1 } });
+            var diff3 = m1 - m3;
+            var diff4 = m3 - m1;
+            Assert.IsTrue(diff3.Equals(m3.Negate()));
+            Assert.IsTrue(diff3.Equals(diff4.Negate()));
+        }
     }
 }
diff --git a/src/UnitTests/LinearAlgebraTests/Single/SparseMatrixTests.cs b/src/UnitTests/LinearAlgebraTests/Single/SparseMatrixTests.cs
index 09b10117..54b14d6d 100644
--- a/src/UnitTests/LinearAlgebraTests/Single/SparseMatrixTests.cs
+++ b/src/UnitTests/LinearAlgebraTests/Single/SparseMatrixTests.cs
@@ -208,5 +208,93 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single
 
             Assert.AreEqual(matrix.NonZerosCount, nonzero);
         }
+
+        /// <summary>
+        /// Test whether order matters when adding sparse matrices.
+        /// </summary>
+        [Test]
+        public void CanAddSparseMatricesBothWays()
+        {
+            var m1 = new SparseMatrix(1, 3);
+            var m2 = new SparseMatrix(new float[,] { { 0, 1, 1 } });
+            var sum1 = m1 + m2;
+            var sum2 = m2 + m1;
+            Assert.IsTrue(sum1.Equals(m2));
+            Assert.IsTrue(sum1.Equals(sum2));
+
+            var sparseResult = new SparseMatrix(1, 3);
+            sparseResult.Add(m2, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(sum1));
+
+            sparseResult = new SparseMatrix(new float[,] { { 0, 1, 1 } });
+            sparseResult.Add(m1, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(sum1));
+
+            sparseResult = new SparseMatrix(new float[,] { { 0, 1, 1 } });
+            m1.Add(sparseResult, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(sum1));
+
+            sparseResult = new SparseMatrix(new float[,] { { 0, 1, 1 } });
+            sparseResult.Add(sparseResult, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(2 * sum1));
+
+            var denseResult = new DenseMatrix(1, 3);
+            denseResult.Add(m2, denseResult);
+            Assert.IsTrue(denseResult.Equals(sum1));
+
+            denseResult = new DenseMatrix(new float[,] { { 0, 1, 1 } });
+            denseResult.Add(m1, denseResult);
+            Assert.IsTrue(denseResult.Equals(sum1));
+
+            var m3 = new DenseMatrix(new float[,] { { 0, 1, 1 } });
+            var sum3 = m1 + m3;
+            var sum4 = m3 + m1;
+            Assert.IsTrue(sum3.Equals(m3));
+            Assert.IsTrue(sum3.Equals(sum4));
+        }
+
+        /// <summary>
+        /// Test whether order matters when subtracting sparse matrices.
+        /// </summary>
+        [Test]
+        public void CanSubtractSparseMatricesBothWays()
+        {
+            var m1 = new SparseMatrix(1, 3);
+            var m2 = new SparseMatrix(new float[,] { { 0, 1, 1 } });
+            var diff1 = m1 - m2;
+            var diff2 = m2 - m1;
+            Assert.IsTrue(diff1.Equals(m2.Negate()));
+            Assert.IsTrue(diff1.Equals(diff2.Negate()));
+
+            var sparseResult = new SparseMatrix(1, 3);
+            sparseResult.Subtract(m2, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(diff1));
+
+            sparseResult = new SparseMatrix(new float[,] { { 0, 1, 1 } });
+            sparseResult.Subtract(m1, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(diff2));
+
+            sparseResult = new SparseMatrix(new float[,] { { 0, 1, 1 } });
+            m1.Subtract(sparseResult, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(diff1));
+
+            sparseResult = new SparseMatrix(new float[,] { { 0, 1, 1 } });
+            sparseResult.Subtract(sparseResult, sparseResult);
+            Assert.IsTrue(sparseResult.Equals(0 * diff1));
+
+            var denseResult = new DenseMatrix(1, 3);
+            denseResult.Subtract(m2, denseResult);
+            Assert.IsTrue(denseResult.Equals(diff1));
+
+            denseResult = new DenseMatrix(new float[,] { { 0, 1, 1 } });
+            denseResult.Subtract(m1, denseResult);
+            Assert.IsTrue(denseResult.Equals(diff2));
+
+            var m3 = new DenseMatrix(new float[,] { { 0, 1, 1 } });
+            var diff3 = m1 - m3;
+            var diff4 = m3 - m1;
+            Assert.IsTrue(diff3.Equals(m3.Negate()));
+            Assert.IsTrue(diff3.Equals(diff4.Negate()));
+        }
     }
 }