diff --git a/src/Numerics/Generate.cs b/src/Numerics/Generate.cs
index 4c4b4ab6..ddb3d20c 100644
--- a/src/Numerics/Generate.cs
+++ b/src/Numerics/Generate.cs
@@ -793,6 +793,24 @@ namespace MathNet.Numerics
             return distribution.Samples();
         }
 
+        /// <summary>
+        /// Create random samples.
+        /// </summary>
+        public static float[] RandomSingle(int length, IContinuousDistribution distribution)
+        {
+            var samples = new double[length];
+            distribution.Samples(samples);
+            return Map(samples, v => (float)v);
+        }
+
+        /// <summary>
+        /// Create an infinite random sample sequence.
+        /// </summary>
+        public static IEnumerable<float> RandomSingle(IContinuousDistribution distribution)
+        {
+            return distribution.Samples().Select(v => (float)v);
+        }
+
         /// <summary>
         /// Create random samples.
         /// </summary>
@@ -809,6 +827,22 @@ namespace MathNet.Numerics
             return RandomMap2Sequence(distribution, (r, i) => new Complex(r, i));
         }
 
+        /// <summary>
+        /// Create random samples.
+        /// </summary>
+        public static Complex32[] RandomComplex32(int length, IContinuousDistribution distribution)
+        {
+            return RandomMap2(length, distribution, (r, i) => new Complex32((float)r, (float)i));
+        }
+
+        /// <summary>
+        /// Create an infinite random sample sequence.
+        /// </summary>
+        public static IEnumerable<Complex32> RandomComplex32(IContinuousDistribution distribution)
+        {
+            return RandomMap2Sequence(distribution, (r, i) => new Complex32((float)r, (float)i));
+        }
+
         /// <summary>
         /// Generate samples by sampling a function at samples from a probability distribution.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Builder.cs b/src/Numerics/LinearAlgebra/Builder.cs
index 05ef6c43..88bc4fa9 100644
--- a/src/Numerics/LinearAlgebra/Builder.cs
+++ b/src/Numerics/LinearAlgebra/Builder.cs
@@ -67,7 +67,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
 
         public override Matrix<double> Random(int rows, int columns, IContinuousDistribution distribution)
         {
-            return Dense(rows, columns, (i, j) => distribution.Sample());
+            return Dense(rows, columns, Generate.Random(rows*columns, distribution));
         }
 
         public override IIterationStopCriterion<double>[] IterativeSolverStopCriteria(int maxIterations = 1000)
@@ -106,7 +106,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
 
         public override Vector<double> Random(int length, IContinuousDistribution distribution)
         {
-            return Dense(length, i => distribution.Sample());
+            return Dense(Generate.Random(length, distribution));
         }
     }
 }
@@ -142,7 +142,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
 
         public override Matrix<float> Random(int rows, int columns, IContinuousDistribution distribution)
         {
-            return Dense(rows, columns, (i, j) => (float) distribution.Sample());
+            return Dense(rows, columns, Generate.RandomSingle(rows*columns, distribution));
         }
 
         public override IIterationStopCriterion<float>[] IterativeSolverStopCriteria(int maxIterations = 1000)
@@ -181,7 +181,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
 
         public override Vector<float> Random(int length, IContinuousDistribution distribution)
         {
-            return Dense(length, i => (float) distribution.Sample());
+            return Dense(Generate.RandomSingle(length, distribution));
         }
     }
 }
@@ -223,7 +223,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
 
         public override Matrix<Complex> Random(int rows, int columns, IContinuousDistribution distribution)
         {
-            return Dense(rows, columns, (i, j) => new Complex(distribution.Sample(), distribution.Sample()));
+            return Dense(rows, columns, Generate.RandomComplex(rows*columns, distribution));
         }
 
         public override IIterationStopCriterion<Complex>[] IterativeSolverStopCriteria(int maxIterations = 1000)
@@ -262,7 +262,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
 
         public override Vector<Complex> Random(int length, IContinuousDistribution distribution)
         {
-            return Dense(length, i => new Complex(distribution.Sample(), distribution.Sample()));
+            return Dense(Generate.RandomComplex(length, distribution));
         }
     }
 }
@@ -298,7 +298,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
 
         public override Matrix<Numerics.Complex32> Random(int rows, int columns, IContinuousDistribution distribution)
         {
-            return Dense(rows, columns, (i, j) => new Numerics.Complex32((float) distribution.Sample(), (float) distribution.Sample()));
+            return Dense(rows, columns, Generate.RandomComplex32(rows*columns, distribution));
         }
 
         public override IIterationStopCriterion<Numerics.Complex32>[] IterativeSolverStopCriteria(int maxIterations = 1000)
@@ -337,7 +337,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
 
         public override Vector<Numerics.Complex32> Random(int length, IContinuousDistribution distribution)
         {
-            return Dense(length, i => new Numerics.Complex32((float) distribution.Sample(), (float) distribution.Sample()));
+            return Dense(Generate.RandomComplex32(length, distribution));
         }
     }
 }
@@ -526,7 +526,7 @@ namespace MathNet.Numerics.LinearAlgebra
         public Matrix<T> Dense(int rows, int columns, T value)
         {
             if (Zero.Equals(value)) return Dense(rows, columns);
-            return Dense(DenseColumnMajorMatrixStorage<T>.OfInit(rows, columns, (i, j) => value));
+            return Dense(DenseColumnMajorMatrixStorage<T>.OfValue(rows, columns, value));
         }
 
         /// <summary>
@@ -860,7 +860,7 @@ namespace MathNet.Numerics.LinearAlgebra
         public Matrix<T> Sparse(int rows, int columns, T value)
         {
             if (Zero.Equals(value)) return Sparse(rows, columns);
-            return Sparse(SparseCompressedRowMatrixStorage<T>.OfInit(rows, columns, (i, j) => value));
+            return Sparse(SparseCompressedRowMatrixStorage<T>.OfValue(rows, columns, value));
         }
 
         /// <summary>
@@ -1229,7 +1229,7 @@ namespace MathNet.Numerics.LinearAlgebra
         public Matrix<T> Diagonal(int rows, int columns, T value)
         {
             if (Zero.Equals(value)) return Diagonal(rows, columns);
-            return Diagonal(DiagonalMatrixStorage<T>.OfInit(rows, columns, i => value));
+            return Diagonal(DiagonalMatrixStorage<T>.OfValue(rows, columns, value));
         }
 
         /// <summary>
@@ -1245,7 +1245,7 @@ namespace MathNet.Numerics.LinearAlgebra
         /// </summary>
         public Matrix<T> DiagonalIdentity(int rows, int columns)
         {
-            return Diagonal(DiagonalMatrixStorage<T>.OfInit(rows, columns, i => One));
+            return Diagonal(DiagonalMatrixStorage<T>.OfValue(rows, columns, One));
         }
 
         /// <summary>
@@ -1253,7 +1253,7 @@ namespace MathNet.Numerics.LinearAlgebra
         /// </summary>
         public Matrix<T> DiagonalIdentity(int order)
         {
-            return Diagonal(DiagonalMatrixStorage<T>.OfInit(order, order, i => One));
+            return Diagonal(DiagonalMatrixStorage<T>.OfValue(order, order, One));
         }
 
 
@@ -1445,7 +1445,7 @@ namespace MathNet.Numerics.LinearAlgebra
         public Vector<T> Dense(int length, T value)
         {
             if (Zero.Equals(value)) return Dense(length);
-            return Dense(DenseVectorStorage<T>.OfInit(length, i => value));
+            return Dense(DenseVectorStorage<T>.OfValue(length, value));
         }
 
         /// <summary>
@@ -1520,7 +1520,7 @@ namespace MathNet.Numerics.LinearAlgebra
         public Vector<T> Sparse(int length, T value)
         {
             if (Zero.Equals(value)) return Sparse(length);
-            return Sparse(SparseVectorStorage<T>.OfInit(length, i => value));
+            return Sparse(SparseVectorStorage<T>.OfValue(length, value));
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
index 8e20573a..4f9a26ab 100644
--- a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
@@ -353,7 +353,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         public static DenseMatrix Create(int rows, int columns, Complex value)
         {
             if (value == Complex.Zero) return new DenseMatrix(rows, columns);
-            return new DenseMatrix(DenseColumnMajorMatrixStorage<Complex>.OfInit(rows, columns, (i, j) => value));
+            return new DenseMatrix(DenseColumnMajorMatrixStorage<Complex>.OfValue(rows, columns, value));
         }
 
         /// <summary>
@@ -394,8 +394,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// </summary>
         public static DenseMatrix CreateRandom(int rows, int columns, IContinuousDistribution distribution)
         {
-            return new DenseMatrix(DenseColumnMajorMatrixStorage<Complex>.OfInit(rows, columns,
-                (i, j) => new Complex(distribution.Sample(), distribution.Sample())));
+            return new DenseMatrix(new DenseColumnMajorMatrixStorage<Complex>(rows, columns, Generate.RandomComplex(rows*columns, distribution)));
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Complex/DenseVector.cs b/src/Numerics/LinearAlgebra/Complex/DenseVector.cs
index a2a6eb1a..c5bdc6bd 100644
--- a/src/Numerics/LinearAlgebra/Complex/DenseVector.cs
+++ b/src/Numerics/LinearAlgebra/Complex/DenseVector.cs
@@ -143,7 +143,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         public static DenseVector Create(int length, Complex value)
         {
             if (value == Complex.Zero) return new DenseVector(length);
-            return new DenseVector(DenseVectorStorage<Complex>.OfInit(length, i => value));
+            return new DenseVector(DenseVectorStorage<Complex>.OfValue(length, value));
         }
 
         /// <summary>
@@ -159,8 +159,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// </summary>
         public static DenseVector CreateRandom(int length, IContinuousDistribution distribution)
         {
-            return new DenseVector(DenseVectorStorage<Complex>.OfInit(length,
-                i => new Complex(distribution.Sample(), distribution.Sample())));
+            var samples = Generate.RandomComplex(length, distribution);
+            return new DenseVector(new DenseVectorStorage<Complex>(length, samples));
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs
index fce362c0..0688a802 100644
--- a/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs
@@ -180,7 +180,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// </summary>
         public static DiagonalMatrix CreateIdentity(int order)
         {
-            return new DiagonalMatrix(DiagonalMatrixStorage<Complex>.OfInit(order, order, i => One));
+            return new DiagonalMatrix(DiagonalMatrixStorage<Complex>.OfValue(order, order, One));
         }
 
         /// <summary>
@@ -188,8 +188,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// </summary>
         public static DiagonalMatrix CreateRandom(int rows, int columns, IContinuousDistribution distribution)
         {
-            return new DiagonalMatrix(DiagonalMatrixStorage<Complex>.OfInit(rows, columns,
-                i => new Complex(distribution.Sample(), distribution.Sample())));
+            return new DiagonalMatrix(new DiagonalMatrixStorage<Complex>(rows, columns, Generate.RandomComplex(Math.Min(rows, columns), distribution)));
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs
index 2ae00785..d76c7b34 100644
--- a/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs
@@ -340,7 +340,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         public static SparseMatrix Create(int rows, int columns, Complex value)
         {
             if (value == Complex.Zero) return new SparseMatrix(rows, columns);
-            return new SparseMatrix(SparseCompressedRowMatrixStorage<Complex>.OfInit(rows, columns, (i, j) => value));
+            return new SparseMatrix(SparseCompressedRowMatrixStorage<Complex>.OfValue(rows, columns, value));
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Complex/SparseVector.cs b/src/Numerics/LinearAlgebra/Complex/SparseVector.cs
index d81325c3..d5384960 100644
--- a/src/Numerics/LinearAlgebra/Complex/SparseVector.cs
+++ b/src/Numerics/LinearAlgebra/Complex/SparseVector.cs
@@ -121,8 +121,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// </summary>
         public static SparseVector Create(int length, Complex value)
         {
-            if (value == Complex.Zero) return new SparseVector(length);
-            return new SparseVector(SparseVectorStorage<Complex>.OfInit(length, i => value));
+            return new SparseVector(SparseVectorStorage<Complex>.OfValue(length, value));
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
index 83a48790..dd812358 100644
--- a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
@@ -348,7 +348,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         public static DenseMatrix Create(int rows, int columns, Complex32 value)
         {
             if (value == Complex32.Zero) return new DenseMatrix(rows, columns);
-            return new DenseMatrix(DenseColumnMajorMatrixStorage<Complex32>.OfInit(rows, columns, (i, j) => value));
+            return new DenseMatrix(DenseColumnMajorMatrixStorage<Complex32>.OfValue(rows, columns, value));
         }
 
         /// <summary>
@@ -389,8 +389,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// </summary>
         public static DenseMatrix CreateRandom(int rows, int columns, IContinuousDistribution distribution)
         {
-            return new DenseMatrix(DenseColumnMajorMatrixStorage<Complex32>.OfInit(rows, columns,
-                (i, j) => new Complex32((float) distribution.Sample(), (float) distribution.Sample())));
+            return new DenseMatrix(new DenseColumnMajorMatrixStorage<Complex32>(rows, columns, Generate.RandomComplex32(rows*columns, distribution)));
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Complex32/DenseVector.cs b/src/Numerics/LinearAlgebra/Complex32/DenseVector.cs
index 5bba0a10..d41c99bb 100644
--- a/src/Numerics/LinearAlgebra/Complex32/DenseVector.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/DenseVector.cs
@@ -138,7 +138,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         public static DenseVector Create(int length, Complex32 value)
         {
             if (value == Complex32.Zero) return new DenseVector(length);
-            return new DenseVector(DenseVectorStorage<Complex32>.OfInit(length, i => value));
+            return new DenseVector(DenseVectorStorage<Complex32>.OfValue(length, value));
         }
 
         /// <summary>
@@ -154,8 +154,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// </summary>
         public static DenseVector CreateRandom(int length, IContinuousDistribution distribution)
         {
-            return new DenseVector(DenseVectorStorage<Complex32>.OfInit(length,
-                i => new Complex32((float)distribution.Sample(), (float)distribution.Sample())));
+            var samples = Generate.RandomComplex32(length, distribution);
+            return new DenseVector(new DenseVectorStorage<Complex32>(length, samples));
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs
index ca47a63e..b06bad4d 100644
--- a/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs
@@ -175,7 +175,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// </summary>
         public static DiagonalMatrix CreateIdentity(int order)
         {
-            return new DiagonalMatrix(DiagonalMatrixStorage<Complex32>.OfInit(order, order, i => One));
+            return new DiagonalMatrix(DiagonalMatrixStorage<Complex32>.OfValue(order, order, One));
         }
 
         /// <summary>
@@ -183,8 +183,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// </summary>
         public static DiagonalMatrix CreateRandom(int rows, int columns, IContinuousDistribution distribution)
         {
-            return new DiagonalMatrix(DiagonalMatrixStorage<Complex32>.OfInit(rows, columns,
-                i => new Complex32((float) distribution.Sample(), (float) distribution.Sample())));
+            return new DiagonalMatrix(new DiagonalMatrixStorage<Complex32>(rows, columns, Generate.RandomComplex32(Math.Min(rows, columns), distribution)));
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs
index 3e779217..ff9fefe2 100644
--- a/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs
@@ -335,7 +335,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         public static SparseMatrix Create(int rows, int columns, Complex32 value)
         {
             if (value == Complex32.Zero) return new SparseMatrix(rows, columns);
-            return new SparseMatrix(SparseCompressedRowMatrixStorage<Complex32>.OfInit(rows, columns, (i, j) => value));
+            return new SparseMatrix(SparseCompressedRowMatrixStorage<Complex32>.OfValue(rows, columns, value));
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Complex32/SparseVector.cs b/src/Numerics/LinearAlgebra/Complex32/SparseVector.cs
index 48c2f13c..b82aa7a4 100644
--- a/src/Numerics/LinearAlgebra/Complex32/SparseVector.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/SparseVector.cs
@@ -116,8 +116,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// </summary>
         public static SparseVector Create(int length, Complex32 value)
         {
-            if (value == Complex32.Zero) return new SparseVector(length);
-            return new SparseVector(SparseVectorStorage<Complex32>.OfInit(length, i => value));
+            return new SparseVector(SparseVectorStorage<Complex32>.OfValue(length, value));
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
index 496984cb..3d576fc3 100644
--- a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
@@ -346,7 +346,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         public static DenseMatrix Create(int rows, int columns, double value)
         {
             if (value == 0d) return new DenseMatrix(rows, columns);
-            return new DenseMatrix(DenseColumnMajorMatrixStorage<double>.OfInit(rows, columns, (i, j) => value));
+            return new DenseMatrix(DenseColumnMajorMatrixStorage<double>.OfValue(rows, columns, value));
         }
 
         /// <summary>
@@ -387,7 +387,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// </summary>
         public static DenseMatrix CreateRandom(int rows, int columns, IContinuousDistribution distribution)
         {
-            return new DenseMatrix(DenseColumnMajorMatrixStorage<double>.OfInit(rows, columns, (i, j) => distribution.Sample()));
+            return new DenseMatrix(new DenseColumnMajorMatrixStorage<double>(rows, columns, Generate.Random(rows*columns, distribution)));
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Double/DenseVector.cs b/src/Numerics/LinearAlgebra/Double/DenseVector.cs
index 12b96828..7d7dc876 100644
--- a/src/Numerics/LinearAlgebra/Double/DenseVector.cs
+++ b/src/Numerics/LinearAlgebra/Double/DenseVector.cs
@@ -138,7 +138,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         public static DenseVector Create(int length, double value)
         {
             if (value == 0d) return new DenseVector(length);
-            return new DenseVector(DenseVectorStorage<double>.OfInit(length, i => value));
+            return new DenseVector(DenseVectorStorage<double>.OfValue(length, value));
         }
 
         /// <summary>
@@ -154,8 +154,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// </summary>
         public static DenseVector CreateRandom(int length, IContinuousDistribution distribution)
         {
-            return new DenseVector(DenseVectorStorage<double>.OfInit(length,
-                i => distribution.Sample()));
+            var samples = Generate.Random(length, distribution);
+            return new DenseVector(new DenseVectorStorage<double>(length, samples));
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
index 1d809f74..adf278e7 100644
--- a/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
@@ -173,7 +173,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// </summary>
         public static DiagonalMatrix CreateIdentity(int order)
         {
-            return new DiagonalMatrix(DiagonalMatrixStorage<double>.OfInit(order, order, i => One));
+            return new DiagonalMatrix(DiagonalMatrixStorage<double>.OfValue(order, order, One));
         }
 
         /// <summary>
@@ -181,8 +181,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// </summary>
         public static DiagonalMatrix CreateRandom(int rows, int columns, IContinuousDistribution distribution)
         {
-            return new DiagonalMatrix(DiagonalMatrixStorage<double>.OfInit(rows, columns,
-                i => distribution.Sample()));
+            return new DiagonalMatrix(new DiagonalMatrixStorage<double>(rows, columns, Generate.Random(Math.Min(rows, columns), distribution)));
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs
index eb2cbdb1..1de758ad 100644
--- a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs
@@ -333,7 +333,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         public static SparseMatrix Create(int rows, int columns, double value)
         {
             if (value == 0d) return new SparseMatrix(rows, columns);
-            return new SparseMatrix(SparseCompressedRowMatrixStorage<double>.OfInit(rows, columns, (i, j) => value));
+            return new SparseMatrix(SparseCompressedRowMatrixStorage<double>.OfValue(rows, columns, value));
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Double/SparseVector.cs b/src/Numerics/LinearAlgebra/Double/SparseVector.cs
index 9ec18d65..7bb2fc9a 100644
--- a/src/Numerics/LinearAlgebra/Double/SparseVector.cs
+++ b/src/Numerics/LinearAlgebra/Double/SparseVector.cs
@@ -116,8 +116,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// </summary>
         public static SparseVector Create(int length, double value)
         {
-            if (value == 0d) return new SparseVector(new SparseVectorStorage<double>(length));
-            return new SparseVector(SparseVectorStorage<double>.OfInit(length, i => value));
+            return new SparseVector(SparseVectorStorage<double>.OfValue(length, value));
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
index 7ec7905a..bc488a59 100644
--- a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
@@ -346,7 +346,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         public static DenseMatrix Create(int rows, int columns, float value)
         {
             if (value == 0f) return new DenseMatrix(rows, columns);
-            return new DenseMatrix(DenseColumnMajorMatrixStorage<float>.OfInit(rows, columns, (i, j) => value));
+            return new DenseMatrix(DenseColumnMajorMatrixStorage<float>.OfValue(rows, columns, value));
         }
 
         /// <summary>
@@ -387,7 +387,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// </summary>
         public static DenseMatrix CreateRandom(int rows, int columns, IContinuousDistribution distribution)
         {
-            return new DenseMatrix(DenseColumnMajorMatrixStorage<float>.OfInit(rows, columns, (i, j) => (float) distribution.Sample()));
+            return new DenseMatrix(new DenseColumnMajorMatrixStorage<float>(rows, columns, Generate.RandomSingle(rows*columns, distribution)));
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Single/DenseVector.cs b/src/Numerics/LinearAlgebra/Single/DenseVector.cs
index 4b69bd9f..be93a48a 100644
--- a/src/Numerics/LinearAlgebra/Single/DenseVector.cs
+++ b/src/Numerics/LinearAlgebra/Single/DenseVector.cs
@@ -137,7 +137,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         public static DenseVector Create(int length, float value)
         {
             if (value == 0f) return new DenseVector(length);
-            return new DenseVector(DenseVectorStorage<float>.OfInit(length, i => value));
+            return new DenseVector(DenseVectorStorage<float>.OfValue(length, value));
         }
 
         /// <summary>
@@ -153,8 +153,8 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// </summary>
         public static DenseVector CreateRandom(int length, IContinuousDistribution distribution)
         {
-            return new DenseVector(DenseVectorStorage<float>.OfInit(length,
-                i => (float)distribution.Sample()));
+            var samples = Generate.RandomSingle(length, distribution);
+            return new DenseVector(new DenseVectorStorage<float>(length, samples));
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
index dcd5056a..7c058b58 100644
--- a/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
@@ -173,7 +173,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// </summary>
         public static DiagonalMatrix CreateIdentity(int order)
         {
-            return new DiagonalMatrix(DiagonalMatrixStorage<float>.OfInit(order, order, i => One));
+            return new DiagonalMatrix(DiagonalMatrixStorage<float>.OfValue(order, order, One));
         }
 
         /// <summary>
@@ -181,8 +181,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// </summary>
         public static DiagonalMatrix CreateRandom(int rows, int columns, IContinuousDistribution distribution)
         {
-            return new DiagonalMatrix(DiagonalMatrixStorage<float>.OfInit(rows, columns,
-                i => (float) distribution.Sample()));
+            return new DiagonalMatrix(new DiagonalMatrixStorage<float>(rows, columns, Generate.RandomSingle(Math.Min(rows, columns), distribution)));
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs
index 615f63f5..5e01087f 100644
--- a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs
@@ -333,7 +333,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         public static SparseMatrix Create(int rows, int columns, float value)
         {
             if (value == 0f) return new SparseMatrix(rows, columns);
-            return new SparseMatrix(SparseCompressedRowMatrixStorage<float>.OfInit(rows, columns, (i, j) => value));
+            return new SparseMatrix(SparseCompressedRowMatrixStorage<float>.OfValue(rows, columns, value));
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Single/SparseVector.cs b/src/Numerics/LinearAlgebra/Single/SparseVector.cs
index 2bc0fc85..1f78d586 100644
--- a/src/Numerics/LinearAlgebra/Single/SparseVector.cs
+++ b/src/Numerics/LinearAlgebra/Single/SparseVector.cs
@@ -116,8 +116,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// </summary>
         public static SparseVector Create(int length, float value)
         {
-            if (value == 0f) return new SparseVector(length);
-            return new SparseVector(SparseVectorStorage<float>.OfInit(length, i => value));
+            return new SparseVector(SparseVectorStorage<float>.OfValue(length, value));
         }
 
         /// <summary>
diff --git a/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs
index 29a3e660..4fedf2fb 100644
--- a/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs
+++ b/src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs
@@ -156,6 +156,20 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             return storage;
         }
 
+        public static DenseColumnMajorMatrixStorage<T> OfValue(int rows, int columns, T value)
+        {
+            var storage = new DenseColumnMajorMatrixStorage<T>(rows, columns);
+            var data = storage.Data;
+            CommonParallel.For(0, data.Length, 4096, (a, b) =>
+            {
+                for (int i = a; i < b; i++)
+                {
+                    data[i] = value;
+                }
+            });
+            return storage;
+        }
+
         public static DenseColumnMajorMatrixStorage<T> OfInit(int rows, int columns, Func<int, int, T> init)
         {
             var storage = new DenseColumnMajorMatrixStorage<T>(rows, columns);
diff --git a/src/Numerics/LinearAlgebra/Storage/DenseVectorStorage.cs b/src/Numerics/LinearAlgebra/Storage/DenseVectorStorage.cs
index f1c6293f..d472ca6b 100644
--- a/src/Numerics/LinearAlgebra/Storage/DenseVectorStorage.cs
+++ b/src/Numerics/LinearAlgebra/Storage/DenseVectorStorage.cs
@@ -109,6 +109,24 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             return storage;
         }
 
+        public static DenseVectorStorage<T> OfValue(int length, T value)
+        {
+            if (length < 1)
+            {
+                throw new ArgumentOutOfRangeException("length", string.Format(Resources.ArgumentLessThanOne, length));
+            }
+
+            var data = new T[length];
+            CommonParallel.For(0, data.Length, 4096, (a, b) =>
+            {
+                for (int i = a; i < b; i++)
+                {
+                    data[i] = value;
+                }
+            });
+            return new DenseVectorStorage<T>(length, data);
+        }
+
         public static DenseVectorStorage<T> OfInit(int length, Func<int, T> init)
         {
             if (length < 1)
diff --git a/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs
index 3bd8c9ec..1101d474 100644
--- a/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs
+++ b/src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs
@@ -233,6 +233,16 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             return storage;
         }
 
+        public static DiagonalMatrixStorage<T> OfValue(int rows, int columns, T diagonalValue)
+        {
+            var storage = new DiagonalMatrixStorage<T>(rows, columns);
+            for (var i = 0; i < storage.Data.Length; i++)
+            {
+                storage.Data[i] = diagonalValue;
+            }
+            return storage;
+        }
+
         public static DiagonalMatrixStorage<T> OfInit(int rows, int columns, Func<int, T> init)
         {
             var storage = new DiagonalMatrixStorage<T>(rows, columns);
diff --git a/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs b/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs
index 25417f61..152d9520 100644
--- a/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs
+++ b/src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs
@@ -444,6 +444,44 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             return storage;
         }
 
+        public static SparseCompressedRowMatrixStorage<T> OfValue(int rows, int columns, T value)
+        {
+            if (Zero.Equals(value))
+            {
+                return new SparseCompressedRowMatrixStorage<T>(rows, columns);
+            }
+
+            var storage = new SparseCompressedRowMatrixStorage<T>(rows, columns);
+
+            var values = new T[rows * columns];
+            for (int i = 0; i < values.Length; i++)
+            {
+                values[i] = value;
+            }
+
+            var rowPointers = storage.RowPointers;
+            for (int i = 0; i <= rows; i++)
+            {
+                rowPointers[i] = i*columns;
+            }
+
+            var columnIndices = new int[values.Length];
+            for (int row = 0; row < rows; row++)
+            {
+                int offset = row*columns;
+                for (int col = 0; col < columns; col++)
+                {
+                    columnIndices[offset + col] = col;
+                }
+            }
+
+            rowPointers[rows] = values.Length;
+            storage.ColumnIndices = columnIndices;
+            storage.Values = values;
+            return storage;
+        }
+
+
         public static SparseCompressedRowMatrixStorage<T> OfInit(int rows, int columns, Func<int, int, T> init)
         {
             var storage = new SparseCompressedRowMatrixStorage<T>(rows, columns);
diff --git a/src/Numerics/LinearAlgebra/Storage/SparseVectorStorage.cs b/src/Numerics/LinearAlgebra/Storage/SparseVectorStorage.cs
index 4145eb5e..0f919562 100644
--- a/src/Numerics/LinearAlgebra/Storage/SparseVectorStorage.cs
+++ b/src/Numerics/LinearAlgebra/Storage/SparseVectorStorage.cs
@@ -302,6 +302,34 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             return storage;
         }
 
+        public static SparseVectorStorage<T> OfValue(int length, T value)
+        {
+            if (Zero.Equals(value))
+            {
+                return new SparseVectorStorage<T>(length);
+            }
+
+            if (length < 1)
+            {
+                throw new ArgumentOutOfRangeException("length", string.Format(Resources.ArgumentLessThanOne, length));
+            }
+
+            var indices = new int[length];
+            var values = new T[length];
+            for (int i = 0; i < indices.Length; i++)
+            {
+                indices[i] = i;
+                values[i] = value;
+            }
+
+            return new SparseVectorStorage<T>(length)
+            {
+                Indices = indices,
+                Values = values,
+                ValueCount = length
+            };
+        }
+
         public static SparseVectorStorage<T> OfInit(int length, Func<int, T> init)
         {
             if (length < 1)
diff --git a/src/UnitTests/LinearAlgebraTests/Complex/DenseMatrixTests.cs b/src/UnitTests/LinearAlgebraTests/Complex/DenseMatrixTests.cs
index 13b40876..e7b90861 100644
--- a/src/UnitTests/LinearAlgebraTests/Complex/DenseMatrixTests.cs
+++ b/src/UnitTests/LinearAlgebraTests/Complex/DenseMatrixTests.cs
@@ -140,7 +140,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex
         [Test]
         public void CanCreateMatrixWithUniformValues()
         {
-            var matrix = DenseMatrix.Create(10, 10, (i, j) => new Complex(10.0, 1));
+            var matrix = DenseMatrix.Create(10, 10, new Complex(10.0, 1));
             for (var i = 0; i < matrix.RowCount; i++)
             {
                 for (var j = 0; j < matrix.ColumnCount; j++)
diff --git a/src/UnitTests/LinearAlgebraTests/Complex/DenseVectorTests.cs b/src/UnitTests/LinearAlgebraTests/Complex/DenseVectorTests.cs
index dd0015b1..2b53d5de 100644
--- a/src/UnitTests/LinearAlgebraTests/Complex/DenseVectorTests.cs
+++ b/src/UnitTests/LinearAlgebraTests/Complex/DenseVectorTests.cs
@@ -145,7 +145,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex
         [Test]
         public void CanCreateDenseVectorWithConstantValues()
         {
-            var vector = DenseVector.Create(5, i => 5);
+            var vector = DenseVector.Create(5, 5);
             foreach (var t in vector)
             {
                 Assert.AreEqual(t, new Complex(5.0, 0));
diff --git a/src/UnitTests/LinearAlgebraTests/Complex/Solvers/Iterative/BiCgStabTest.cs b/src/UnitTests/LinearAlgebraTests/Complex/Solvers/Iterative/BiCgStabTest.cs
index 820f65ae..3d4b0b95 100644
--- a/src/UnitTests/LinearAlgebraTests/Complex/Solvers/Iterative/BiCgStabTest.cs
+++ b/src/UnitTests/LinearAlgebraTests/Complex/Solvers/Iterative/BiCgStabTest.cs
@@ -96,7 +96,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Solvers.Iterativ
             var matrix = SparseMatrix.CreateIdentity(100);
 
             // Create the y vector
-            var y = DenseVector.Create(matrix.RowCount, i => Complex.One);
+            var y = DenseVector.Create(matrix.RowCount, Complex.One);
 
             // Create an iteration monitor which will keep track of iterative convergence
             var monitor = new Iterator<Complex>(
@@ -140,7 +140,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Solvers.Iterativ
             matrix.Multiply(new Complex(Math.PI, Math.PI), matrix);
 
             // Create the y vector
-            var y = DenseVector.Create(matrix.RowCount, i => Complex.One);
+            var y = DenseVector.Create(matrix.RowCount, Complex.One);
 
             // Create an iteration monitor which will keep track of iterative convergence
             var monitor = new Iterator<Complex>(new IterationCountStopCriterion<Complex>(MaximumIterations),
@@ -216,7 +216,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Solvers.Iterativ
             }
 
             // Create the y vector
-            var y = DenseVector.Create(matrix.RowCount, i => Complex.One);
+            var y = DenseVector.Create(matrix.RowCount, Complex.One);
 
             // Create an iteration monitor which will keep track of iterative convergence
             var monitor = new Iterator<Complex>(new IterationCountStopCriterion<Complex>(MaximumIterations),
diff --git a/src/UnitTests/LinearAlgebraTests/Complex32/DenseVectorTests.cs b/src/UnitTests/LinearAlgebraTests/Complex32/DenseVectorTests.cs
index e14f74d3..f5c7869d 100644
--- a/src/UnitTests/LinearAlgebraTests/Complex32/DenseVectorTests.cs
+++ b/src/UnitTests/LinearAlgebraTests/Complex32/DenseVectorTests.cs
@@ -141,7 +141,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32
         [Test]
         public void CanCreateDenseVectorWithConstantValues()
         {
-            var vector = DenseVector.Create(5, i => 5);
+            var vector = DenseVector.Create(5, 5);
             foreach (var t in vector)
             {
                 Assert.AreEqual(t, new Complex32(5.0f, 0));
diff --git a/src/UnitTests/LinearAlgebraTests/Complex32/Solvers/Iterative/BiCgStabTest.cs b/src/UnitTests/LinearAlgebraTests/Complex32/Solvers/Iterative/BiCgStabTest.cs
index 56845219..ad8ed33b 100644
--- a/src/UnitTests/LinearAlgebraTests/Complex32/Solvers/Iterative/BiCgStabTest.cs
+++ b/src/UnitTests/LinearAlgebraTests/Complex32/Solvers/Iterative/BiCgStabTest.cs
@@ -91,7 +91,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Solvers.Iterat
             var matrix = SparseMatrix.CreateIdentity(100);
 
             // Create the y vector
-            var y = DenseVector.Create(matrix.RowCount, i => Complex32.One);
+            var y = DenseVector.Create(matrix.RowCount, Complex32.One);
 
             // Create an iteration monitor which will keep track of iterative convergence
             var monitor = new Iterator<Complex32>(
@@ -135,7 +135,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Solvers.Iterat
             matrix.Multiply(new Complex32((float)Math.PI, (float)Math.PI), matrix);
 
             // Create the y vector
-            var y = DenseVector.Create(matrix.RowCount, i => Complex32.One);
+            var y = DenseVector.Create(matrix.RowCount, Complex32.One);
 
             // Create an iteration monitor which will keep track of iterative convergence
             var monitor = new Iterator<Complex32>(
@@ -212,7 +212,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Solvers.Iterat
             }
 
             // Create the y vector
-            var y = DenseVector.Create(matrix.RowCount, i => Complex32.One);
+            var y = DenseVector.Create(matrix.RowCount, Complex32.One);
 
             // Create an iteration monitor which will keep track of iterative convergence
             var monitor = new Iterator<Complex32>(
diff --git a/src/UnitTests/LinearAlgebraTests/Complex32/Solvers/StopCriterion/ResidualStopCriteriumTest.cs b/src/UnitTests/LinearAlgebraTests/Complex32/Solvers/StopCriterion/ResidualStopCriteriumTest.cs
index 1719e02c..1a89290f 100644
--- a/src/UnitTests/LinearAlgebraTests/Complex32/Solvers/StopCriterion/ResidualStopCriteriumTest.cs
+++ b/src/UnitTests/LinearAlgebraTests/Complex32/Solvers/StopCriterion/ResidualStopCriteriumTest.cs
@@ -131,7 +131,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Solvers.StopCr
                 1,
                 Vector<Complex32>.Build.Dense(3, 4),
                 Vector<Complex32>.Build.Dense(3, 4),
-                DenseVector.Create(4, i => 4)), Throws.ArgumentException);
+                DenseVector.Create(4, 4)), Throws.ArgumentException);
         }
 
         /// <summary>
diff --git a/src/UnitTests/LinearAlgebraTests/Double/DenseMatrixTests.cs b/src/UnitTests/LinearAlgebraTests/Double/DenseMatrixTests.cs
index 9193e926..830b7861 100644
--- a/src/UnitTests/LinearAlgebraTests/Double/DenseMatrixTests.cs
+++ b/src/UnitTests/LinearAlgebraTests/Double/DenseMatrixTests.cs
@@ -136,7 +136,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
         [Test]
         public void CanCreateMatrixWithUniformValues()
         {
-            var matrix = Matrix<double>.Build.Dense(10, 10, (i, j) => 10.0);
+            var matrix = Matrix<double>.Build.Dense(10, 10, 10.0);
             Assert.That(matrix, Is.TypeOf<DenseMatrix>());
             for (var i = 0; i < matrix.RowCount; i++)
             {
diff --git a/src/UnitTests/LinearAlgebraTests/Single/DenseVectorTests.cs b/src/UnitTests/LinearAlgebraTests/Single/DenseVectorTests.cs
index 9fd4d679..8adc053b 100644
--- a/src/UnitTests/LinearAlgebraTests/Single/DenseVectorTests.cs
+++ b/src/UnitTests/LinearAlgebraTests/Single/DenseVectorTests.cs
@@ -139,7 +139,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single
         [Test]
         public void CanCreateDenseVectorWithConstantValues()
         {
-            var vector = DenseVector.Create(5, i => 5);
+            var vector = DenseVector.Create(5, 5);
             foreach (var t in vector)
             {
                 Assert.AreEqual(t, 5);