diff --git a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
index 7deaf8d1..63ea0f19 100644
--- a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
@@ -506,11 +506,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
                 return;
             }
 
-            // dense + diagonal = matrix
+            // dense + diagonal = any
             var diagonalOther = other.Storage as DiagonalMatrixStorage<Complex>;
             if (diagonalOther != null)
             {
-                CopyTo(result);
+                Storage.CopyToUnchecked(result.Storage);
                 var diagonal = diagonalOther.Data;
                 for (int i = 0; i < diagonal.Length; i++)
                 {
diff --git a/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs
index d18bdff1..094d8986 100644
--- a/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs
@@ -202,16 +202,19 @@ 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)
         {
+            // diagonal + diagonal = diagonal
             var diagOther = other as DiagonalMatrix;
             var diagResult = result as DiagonalMatrix;
-
-            if (diagOther == null || diagResult == null)
+            if (diagOther != null && diagResult != null)
             {
-                base.DoAdd(other, result);
+                Control.LinearAlgebraProvider.AddArrays(_data, diagOther._data, diagResult._data);
+                return;
             }
-            else
+
+            other.CopyTo(result);
+            for (int i = 0; i < _data.Length; i++)
             {
-                Control.LinearAlgebraProvider.AddArrays(_data, diagOther._data, diagResult._data);
+                result.At(i, i, result.At(i, i) + _data[i]);
             }
         }
 
@@ -223,16 +226,19 @@ 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)
         {
+            // diagonal - diagonal = diagonal
             var diagOther = other as DiagonalMatrix;
             var diagResult = result as DiagonalMatrix;
-
-            if (diagOther == null || diagResult == null)
+            if (diagOther != null && diagResult != null)
             {
-                base.DoSubtract(other, result);
+                Control.LinearAlgebraProvider.SubtractArrays(_data, diagOther._data, diagResult._data);
+                return;
             }
-            else
+
+            other.Negate(result);
+            for (int i = 0; i < _data.Length; i++)
             {
-                Control.LinearAlgebraProvider.SubtractArrays(_data, diagOther._data, diagResult._data);
+                result.At(i, i, result.At(i, i) + _data[i]);
             }
         }
 
diff --git a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
index 1337b9aa..cfe659f0 100644
--- a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
@@ -501,11 +501,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
                 return;
             }
 
-            // dense + diagonal = matrix
+            // dense + diagonal = any
             var diagonalOther = other.Storage as DiagonalMatrixStorage<Complex32>;
             if (diagonalOther != null)
             {
-                CopyTo(result);
+                Storage.CopyToUnchecked(result.Storage);
                 var diagonal = diagonalOther.Data;
                 for (int i = 0; i < diagonal.Length; i++)
                 {
diff --git a/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs
index 233b51d0..4c55048f 100644
--- a/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs
@@ -197,16 +197,19 @@ 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)
         {
+            // diagonal + diagonal = diagonal
             var diagOther = other as DiagonalMatrix;
             var diagResult = result as DiagonalMatrix;
-
-            if (diagOther == null || diagResult == null)
+            if (diagOther != null && diagResult != null)
             {
-                base.DoAdd(other, result);
+                Control.LinearAlgebraProvider.AddArrays(_data, diagOther._data, diagResult._data);
+                return;
             }
-            else
+
+            other.CopyTo(result);
+            for (int i = 0; i < _data.Length; i++)
             {
-                Control.LinearAlgebraProvider.AddArrays(_data, diagOther._data, diagResult._data);
+                result.At(i, i, result.At(i, i) + _data[i]);
             }
         }
 
@@ -218,16 +221,19 @@ 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)
         {
+            // diagonal - diagonal = diagonal
             var diagOther = other as DiagonalMatrix;
             var diagResult = result as DiagonalMatrix;
-
-            if (diagOther == null || diagResult == null)
+            if (diagOther != null && diagResult != null)
             {
-                base.DoSubtract(other, result);
+                Control.LinearAlgebraProvider.SubtractArrays(_data, diagOther._data, diagResult._data);
+                return;
             }
-            else
+
+            other.Negate(result);
+            for (int i = 0; i < _data.Length; i++)
             {
-                Control.LinearAlgebraProvider.SubtractArrays(_data, diagOther._data, diagResult._data);
+                result.At(i, i, result.At(i, i) + _data[i]);
             }
         }
 
diff --git a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
index a45cd15f..abb831bc 100644
--- a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
@@ -481,11 +481,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 return;
             }
 
-            // dense + diagonal = matrix
+            // dense + diagonal = any
             var diagonalOther = other.Storage as DiagonalMatrixStorage<double>;
             if (diagonalOther != null)
             {
-                CopyTo(result);
+                Storage.CopyToUnchecked(result.Storage);
                 var diagonal = diagonalOther.Data;
                 for (int i = 0; i < diagonal.Length; i++)
                 {
diff --git a/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
index 9d735cb5..2daca108 100644
--- a/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
@@ -196,16 +196,19 @@ 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)
         {
+            // diagonal + diagonal = diagonal
             var diagOther = other as DiagonalMatrix;
             var diagResult = result as DiagonalMatrix;
-
-            if (diagOther == null || diagResult == null)
+            if (diagOther != null && diagResult != null)
             {
-                base.DoAdd(other, result);
+                Control.LinearAlgebraProvider.AddArrays(_data, diagOther._data, diagResult._data);
+                return;
             }
-            else
+
+            other.CopyTo(result);
+            for (int i = 0; i < _data.Length; i++)
             {
-                Control.LinearAlgebraProvider.AddArrays(_data, diagOther._data, diagResult._data);
+                result.At(i, i, result.At(i, i) + _data[i]);
             }
         }
 
@@ -217,16 +220,19 @@ 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)
         {
+            // diagonal - diagonal = diagonal
             var diagOther = other as DiagonalMatrix;
             var diagResult = result as DiagonalMatrix;
-
-            if (diagOther == null || diagResult == null)
+            if (diagOther != null && diagResult != null)
             {
-                base.DoSubtract(other, result);
+                Control.LinearAlgebraProvider.SubtractArrays(_data, diagOther._data, diagResult._data);
+                return;
             }
-            else
+
+            other.Negate(result);
+            for (int i = 0; i < _data.Length; i++)
             {
-                Control.LinearAlgebraProvider.SubtractArrays(_data, diagOther._data, diagResult._data);
+                result.At(i, i, result.At(i, i) + _data[i]);
             }
         }
 
diff --git a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
index 61dc0c34..4691979b 100644
--- a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
@@ -481,11 +481,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single
                 return;
             }
 
-            // dense + diagonal = matrix
+            // dense + diagonal = any
             var diagonalOther = other.Storage as DiagonalMatrixStorage<float>;
             if (diagonalOther != null)
             {
-                CopyTo(result);
+                Storage.CopyToUnchecked(result.Storage);
                 var diagonal = diagonalOther.Data;
                 for (int i = 0; i < diagonal.Length; i++)
                 {
diff --git a/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
index fa9aa7c7..5b84ee75 100644
--- a/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
@@ -196,16 +196,19 @@ 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)
         {
+            // diagonal + diagonal = diagonal
             var diagOther = other as DiagonalMatrix;
             var diagResult = result as DiagonalMatrix;
-
-            if (diagOther == null || diagResult == null)
+            if (diagOther != null && diagResult != null)
             {
-                base.DoAdd(other, result);
+                Control.LinearAlgebraProvider.AddArrays(_data, diagOther._data, diagResult._data);
+                return;
             }
-            else
+
+            other.CopyTo(result);
+            for (int i = 0; i < _data.Length; i++)
             {
-                Control.LinearAlgebraProvider.AddArrays(_data, diagOther._data, diagResult._data);
+                result.At(i, i, result.At(i, i) + _data[i]);
             }
         }
 
@@ -217,16 +220,19 @@ 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)
         {
+            // diagonal - diagonal = diagonal
             var diagOther = other as DiagonalMatrix;
             var diagResult = result as DiagonalMatrix;
-
-            if (diagOther == null || diagResult == null)
+            if (diagOther != null && diagResult != null)
             {
-                base.DoSubtract(other, result);
+                Control.LinearAlgebraProvider.SubtractArrays(_data, diagOther._data, diagResult._data);
+                return;
             }
-            else
+
+            other.Negate(result);
+            for (int i = 0; i < _data.Length; i++)
             {
-                Control.LinearAlgebraProvider.SubtractArrays(_data, diagOther._data, diagResult._data);
+                result.At(i, i, result.At(i, i) + _data[i]);
             }
         }