diff --git a/src/FSharp/FSharp.fsproj b/src/FSharp/FSharp.fsproj
index 8c68a2e8..859175de 100644
--- a/src/FSharp/FSharp.fsproj
+++ b/src/FSharp/FSharp.fsproj
@@ -62,6 +62,8 @@
     <Compile Include="AssemblyInfo.fs" />
     <Compile Include="Random.fs" />
     <Compile Include="Distributions.fs" />
+    <Compile Include="LinearAlgebra.Vector.fs" />
+    <Compile Include="LinearAlgebra.Matrix.fs" />
     <Compile Include="LinearAlgebra.Double.Vector.fs" />
     <Compile Include="LinearAlgebra.Double.Matrix.fs" />
     <Compile Include="LinearAlgebra.Double.fs" />
@@ -73,8 +75,6 @@
     <Compile Include="Fit.fs" />
     <Compile Include="FindRoots.fs" />
     <Compile Include="RandomVariable.fs" />
-    <Compile Include="LinearAlgebra.Vector.fs" />
-    <Compile Include="LinearAlgebra.Matrix.fs" />
   </ItemGroup>
   <ItemGroup>
     <Reference Include="FSharp.Core" />
diff --git a/src/FSharp/LinearAlgebra.Double.Matrix.fs b/src/FSharp/LinearAlgebra.Double.Matrix.fs
index 966e0bf2..5508935c 100644
--- a/src/FSharp/LinearAlgebra.Double.Matrix.fs
+++ b/src/FSharp/LinearAlgebra.Double.Matrix.fs
@@ -37,15 +37,10 @@ open MathNet.Numerics.LinearAlgebra
 module Matrix =
 
     /// Returns the sum of all elements of a matrix.
-    let inline sum (A: #Matrix<_>) =
-        let mutable f = 0.0
-        for i=0 to A.RowCount-1 do
-            for j=0 to A.ColumnCount-1 do
-                f <- f + A.At(i,j)
-        f
+    let inline sum (A: #Matrix<float>) = A |> Matrix.foldnz (+) 0.0
 
     /// Fold all columns into one row vector.
-    let inline foldByCol (f: float -> float -> float) acc (A: #Matrix<float>) =
+    let inline foldByCol f acc (A: #Matrix<float>) =
         let v = new DenseVector(A.ColumnCount)
         for k=0 to A.ColumnCount-1 do
             let mutable macc = acc
@@ -55,7 +50,7 @@ module Matrix =
         v :> _ Vector
 
     /// Fold all rows into one column vector.
-    let inline foldByRow (f: float -> float -> float) acc (A: #Matrix<float>) =
+    let inline foldByRow f acc (A: #Matrix<float>) =
         let v = new DenseVector(A.RowCount)
         for k=0 to A.RowCount-1 do
             let mutable macc = acc
diff --git a/src/FSharp/LinearAlgebra.Matrix.fs b/src/FSharp/LinearAlgebra.Matrix.fs
index 98f49acf..f81408b2 100644
--- a/src/FSharp/LinearAlgebra.Matrix.fs
+++ b/src/FSharp/LinearAlgebra.Matrix.fs
@@ -99,29 +99,167 @@ module Matrix =
 
 
     /// Transform a matrix into a sequence.
-    let inline toSeq (v: #Matrix<_>) = v.Enumerate()
+    let inline toSeq (m: #Matrix<_>) = m.Enumerate()
 
     /// Transform a matrix into an indexed sequence.
-    let inline toSeqi (v: #Matrix<_>) = v.EnumerateIndexed()
+    let inline toSeqi (m: #Matrix<_>) = m.EnumerateIndexed()
 
-    /// Transform a matrix into a sequence where zero-values are skipped.
-    let inline toSeqnz (v: #Matrix<_>) = v.EnumerateNonZero()
-
-    /// Transform a matrix into an indexed sequence where zero-values are skipped.
-    let inline toSeqinz (v: #Matrix<_>) = v.EnumerateNonZeroIndexed()
+    /// Transform a matrix into a sequence where zero-values are skipped. Skipping zeros is efficient on sparse data.
+    let inline toSeqnz (m: #Matrix<_>) = m.EnumerateNonZero()
 
+    /// Transform a matrix into an indexed sequence where zero-values are skipped. Skipping zeros is efficient on sparse data.
+    let inline toSeqinz (m: #Matrix<_>) = m.EnumerateNonZeroIndexed()
 
     /// Transform a matrix into a column sequence.
-    let inline toColSeq (v: #Matrix<_>) = v.EnumerateColumns()
+    let inline toColSeq (m: #Matrix<_>) = m.EnumerateColumns()
 
     /// Transform a matrix into an indexed column sequence.
-    let inline toColSeqi (v: #Matrix<_>) = v.EnumerateColumnsIndexed()
+    let inline toColSeqi (m: #Matrix<_>) = m.EnumerateColumnsIndexed()
 
     /// Transform a matrix into a row sequence.
-    let inline toRowSeq (v: #Matrix<_>) = v.EnumerateRows()
+    let inline toRowSeq (m: #Matrix<_>) = m.EnumerateRows()
 
     /// Transform a matrix into an indexed row sequence.
-    let inline toRowSeqi (v: #Matrix<_>) = v.EnumerateRowsIndexed()
+    let inline toRowSeqi (m: #Matrix<_>) = m.EnumerateRowsIndexed()
+
+
+    /// Applies a function to all elements of the matrix.
+    let inline iter f (m: #Matrix<_>) = m.Enumerate() |> Seq.iter f
+
+    /// Applies a function to all indexed elements of the matrix.
+    let inline iteri f (m: #Matrix<_>) = m.EnumerateIndexed() |> Seq.iter (fun (i, j, x) -> f i j x)
+
+    /// Applies a function to all non-zero elements of the matrix. Skipping zeros is efficient on sparse data.
+    let inline iternz f (m: #Matrix<_>) = m.EnumerateNonZero() |> Seq.iter f
+
+    /// Applies a function to all non-zero indexed elements of the matrix. Skipping zeros is efficient on sparse data.
+    let inline iterinz f (m: #Matrix<_>) = m.EnumerateNonZeroIndexed() |> Seq.iter (fun (i, j, x) -> f i j x)
+
+    /// Applies a function to all columns of the matrix.
+    let inline iterCols f (m: #Matrix<_>) = m.EnumerateColumns() |> Seq.iter f
+
+    /// Applies a function to all indexed columns of the matrix.
+    let inline iteriCols f (m: #Matrix<_>) = m.EnumerateColumns() |> Seq.iteri f
+
+    /// Applies a function to all rows of the matrix.
+    let inline iterRows f (m: #Matrix<_>) = m.EnumerateRows() |> Seq.iter f
+
+    /// Applies a function to all indexed rows of the matrix.
+    let inline iteriRows f (m: #Matrix<_>) = m.EnumerateRows() |> Seq.iteri f
+
+
+    /// Fold all entries of a matrix.
+    let inline fold f state (m: #Matrix<_>) = m.Enumerate() |> Seq.fold f state
+
+    /// Fold all entries of a matrix with an indexed folding function.
+    let inline foldi f state (m: #Matrix<_>) = m.EnumerateIndexed() |> Seq.fold (fun s (i,j,x) -> f i j s x) state
+
+    /// Fold all non-zero entries of a matrix. Skipping zeros is efficient on sparse data.
+    let inline foldnz f state (m: #Matrix<_>) = m.EnumerateNonZero() |> Seq.fold f state
+
+    /// Fold all non-zero entries of a matrix with an indexed folding function. Skipping zeros is efficient on sparse data.
+    let inline foldinz f state (m: #Matrix<_>) = m.EnumerateNonZeroIndexed() |> Seq.fold (fun s (i,j,x) -> f i j s x) state
+
+    /// Fold all columns of a matrix.
+    let inline foldCols f state (m: #Matrix<_>) = m.EnumerateColumns() |> Seq.fold f state
+
+    /// Fold all columns of a matrix with an indexed folding function.
+    let inline foldiCols f state (m: #Matrix<_>) = m.EnumerateColumnsIndexed() |> Seq.fold (fun s (j,x) -> f j s x) state
+
+    /// Fold all rows of a matrix.
+    let inline foldRows f state (m: #Matrix<_>) = m.EnumerateRows() |> Seq.fold f state
+
+    /// Fold all rows of a matrix with an indexed folding function.
+    let inline foldiRows f state (m: #Matrix<_>) = m.EnumerateRowsIndexed() |> Seq.fold (fun s (i,x) -> f i s x) state
+
+
+    /// Scan all entries of a matrix.
+    let inline scan f state (m: #Matrix<_>) = m.Enumerate() |> Seq.scan f state
+
+    /// Scan all entries of a matrix with an indexed folding function.
+    let inline scani f state (m: #Matrix<_>) = m.EnumerateIndexed() |> Seq.scan (fun s (i,j,x) -> f i j s x) state
+
+    /// Scan all non-zero entries of a matrix. Skipping zeros is efficient on sparse data.
+    let inline scannz f state (m: #Matrix<_>) = m.EnumerateNonZero() |> Seq.scan f state
+
+    /// Scan all non-zero entries of a matrix with an indexed folding function. Skipping zeros is efficient on sparse data.
+    let inline scaninz f state (m: #Matrix<_>) = m.EnumerateNonZeroIndexed() |> Seq.scan (fun s (i,j,x) -> f i j s x) state
+
+    /// Scan all columns of a matrix.
+    let inline scanCols f state (m: #Matrix<_>) = m.EnumerateColumns() |> Seq.scan f state
+
+    /// Scan all columns of a matrix with an indexed folding function.
+    let inline scaniCols f state (m: #Matrix<_>) = m.EnumerateColumnsIndexed() |> Seq.scan (fun s (j,x) -> f j s x) state
+
+    /// Scan all rows of a matrix.
+    let inline scanRows f state (m: #Matrix<_>) = m.EnumerateRows() |> Seq.scan f state
+
+    /// Scan all rows of a matrix with an indexed folding function.
+    let inline scaniRows f state (m: #Matrix<_>) = m.EnumerateRowsIndexed() |> Seq.scan (fun s (i,x) -> f i s x) state
+
+
+    /// Reduce all entries of a matrix.
+    let inline reduce f (m: #Matrix<_>) = m.Enumerate() |> Seq.reduce f
+
+    /// Reduce all non-zero entries of a matrix. Skipping zeros is efficient on sparse data.
+    let inline reducenz f (m: #Matrix<_>) = m.EnumerateNonZero() |> Seq.reduce f
+
+    /// Reduce all columns of a matrix.
+    let inline reduceCols f (m: #Matrix<_>) = m.EnumerateColumns() |> Seq.reduce f
+
+    /// Reduce all rows of a matrix.
+    let inline reduceRows f (m: #Matrix<_>) = m.EnumerateColumns() |> Seq.reduce f
+
+
+    /// Checks whether there is an entry in the matrix that satisfies a predicate.
+    let inline exists p (m: #Matrix<_>) = m.Enumerate() |> Seq.exists p
+
+    /// Checks whether there is an entry in the matrix that satisfies a position dependent predicate.
+    let inline existsi p (m: #Matrix<_>) = m.EnumerateIndexed() |> Seq.exists (fun (i,j,x) -> p i j x)
+
+    /// Checks whether there is a non-zero entry in the matrix that satisfies a predicate. Skipping zeros is efficient on sparse data.
+    let inline existsnz p (m: #Matrix<_>) = m.EnumerateNonZero() |> Seq.exists p
+
+    /// Checks whether there is a non-zero entry in the matrix that satisfies a position dependent predicate. Skipping zeros is efficient on sparse data.
+    let inline existsinz p (m: #Matrix<_>) = m.EnumerateNonZeroIndexed() |> Seq.exists (fun (i,j,x) -> p i j x)
+
+    /// Checks whether there is a column in the matrix that satisfies a predicate.
+    let inline existsCol p (m: #Matrix<_>) = m.EnumerateColumns() |> Seq.exists p
+    
+    /// Checks whether there is a column in the matrix that satisfies a position dependent predicate.
+    let inline existsiCol p (m: #Matrix<_>) = m.EnumerateColumnsIndexed() |> Seq.exists (fun (j,x) -> p j x)
+    
+    /// Checks whether there is a row in the matrix that satisfies a predicate.
+    let inline existsRow p (m: #Matrix<_>) = m.EnumerateRows() |> Seq.exists p
+    
+    /// Checks whether there is a row in the matrix that satisfies a position dependent predicate.
+    let inline existsiRow p (m: #Matrix<_>) = m.EnumerateRowsIndexed() |> Seq.exists (fun (i,x) -> p i x)
+
+
+    /// Checks whether all entries in the matrix that satisfies a given predicate.
+    let inline forall p (m: #Matrix<_>) = m.Enumerate() |> Seq.forall p
+
+    /// Checks whether all entries in the matrix that satisfies a given position dependent predicate.
+    let inline foralli p (m: #Matrix<_>) = m.EnumerateIndexed() |> Seq.forall (fun (i,j,x) -> p i j x)
+
+    /// Checks whether all non-zero entries in the matrix that satisfies a given predicate. Skipping zeros is efficient on sparse data.
+    let inline forallnz p (m: #Matrix<_>) = m.EnumerateNonZero() |> Seq.forall p
+
+    /// Checks whether all non-zero entries in the matrix that satisfies a given position dependent predicate. Skipping zeros is efficient on sparse data.
+    let inline forallinz p (m: #Matrix<_>) = m.EnumerateNonZeroIndexed() |> Seq.forall (fun (i,j,x) -> p i j x)
+
+    /// Checks whether all columns in the matrix that satisfy a predicate.
+    let inline forallCols p (m: #Matrix<_>) = m.EnumerateColumns() |> Seq.forall p
+    
+    /// Checks whether all columns in the matrix that satisfy a position dependent predicate.
+    let inline foralliCols p (m: #Matrix<_>) = m.EnumerateColumnsIndexed() |> Seq.forall (fun (j,x) -> p j x)
+    
+    /// Checks whether all rows in the matrix that satisfy a predicate.
+    let inline forallRows p (m: #Matrix<_>) = m.EnumerateRows() |> Seq.forall p
+    
+    /// Checks whether all rows in the matrix that satisfy a position dependent predicate.
+    let inline foralliRows p (m: #Matrix<_>) = m.EnumerateRowsIndexed() |> Seq.forall (fun (i,x) -> p i x)
+
 
 
     /// In-place map of every matrix element using a function.
@@ -142,7 +280,6 @@ module Matrix =
     let inline mapinzInPlace f (A: #Matrix<_>) =
         A.MapIndexedInplace((fun i j x -> f i j x), false)
 
-
     /// In-place map every matrix column using the given position dependent function.
     let inline mapColsInPlace (f: int -> Vector<'a> -> Vector<'a>) (A: #Matrix<_>) =
         for j = 0 to A.ColumnCount-1 do
@@ -153,6 +290,7 @@ module Matrix =
         for i = 0 to A.RowCount-1 do
             A.SetRow(i, f i (A.Row(i)))
 
+
     /// Map every matrix element using the given function.
     let inline map f (A: #Matrix<_>) =
         let A = A.Clone()
@@ -191,110 +329,36 @@ module Matrix =
         mapRowsInPlace f A
         A
 
-    /// Fold a function over all matrix elements.
-    let inline fold f acc0 (A: #Matrix<_>) =
-        let n = A.RowCount
-        let m = A.ColumnCount
-        let mutable acc = acc0
-        for i=0 to n-1 do
-            for j=0 to m-1 do
-                acc <- f acc (A.At(i,j))
+
+
+    /// Fold one column.
+    let inline foldCol f state (A: #Matrix<_>) k =
+        let mutable acc = state
+        for i=0 to A.RowCount-1 do
+            acc <- f acc (A.Item(i,k))
+        acc
+
+    /// Fold one row.
+    let inline foldRow f state (A: #Matrix<_>) k =
+        let mutable acc = state
+        for i=0 to A.ColumnCount-1 do
+            acc <- f acc (A.Item(k,i))
         acc
 
     /// Fold a function over all matrix elements in reverse order.
-    let inline foldBack f acc0 (A: #Matrix<_>) =
+    let inline foldBack f state (A: #Matrix<_>) =
         let n = A.RowCount
         let m = A.ColumnCount
-        let mutable acc = acc0
+        let mutable acc = state
         for i in n-1 .. -1 .. 0 do
             for j in m-1 .. -1 .. 0 do
                 acc <- f (A.At(i,j)) acc
         acc
 
-    /// Fold a matrix by applying a given function to all matrix elements.
-    let inline foldi f acc0 (A: #Matrix<_>) =
-        let n = A.RowCount
-        let m = A.ColumnCount
-        let mutable acc = acc0
-        for i=0 to n-1 do
-            for j=0 to m-1 do
-                acc <- f i j acc (A.At(i,j))
-        acc
-
-    /// Checks whether a predicate holds for all elements of a matrix.
-    let inline forall p (A: #Matrix<_>) =
-        let mutable b = true
-        let mutable i = 0
-        let mutable j = 0
-        while b && i < A.RowCount do
-            b <- b && (p (A.At(i,j)))
-            j <- j+1
-            if j = A.ColumnCount then i <- i+1; j <- 0
-        b
-
-    /// Chechks whether a predicate holds for at least one element of a matrix.
-    let inline exists p (A: #Matrix<_>) =
-        let mutable b = false
-        let mutable i = 0
-        let mutable j = 0
-        while not(b) && i < A.RowCount do
-            b <- b || (p (A.At(i,j)))
-            j <- j+1
-            if j = A.ColumnCount then i <- i+1; j <- 0
-        b
-
-    /// Checks whether a position dependent predicate holds for all elements of a matrix.
-    let inline foralli p (A: #Matrix<_>) =
-        let mutable b = true
-        let mutable i = 0
-        let mutable j = 0
-        while b && i < A.RowCount do
-            b <- b && (p i j (A.At(i,j)))
-            j <- j+1
-            if j = A.ColumnCount then i <- i+1; j <- 0
-        b
-
-    /// Checks whether a position dependent predicate holds for at least one element of a matrix.
-    let inline existsi p (A: #Matrix<_>) =
-        let mutable b = false
-        let mutable i = 0
-        let mutable j = 0
-        while not(b) && i < A.RowCount do
-            b <- b || (p i j (A.At(i,j)))
-            j <- j+1
-            if j = A.ColumnCount then i <- i+1; j <- 0
-        b
-
     /// In-place assignment.
     let inline inplaceAssign f (A: #Matrix<_>) =
         A.MapIndexedInplace((fun i j x -> f i j), true)
 
-    /// Iterates over all elements of a matrix.
-    let inline iter f (A: #Matrix<_>) =
-        for i=0 to A.RowCount-1 do
-            for j=0 to A.ColumnCount-1 do
-                f (A.At(i,j))
-
-    /// Iterates over all elements of a matrix using the element indices.
-    let inline iteri f (A: #Matrix<_>) =
-        for i=0 to A.RowCount-1 do
-            for j=0 to A.ColumnCount-1 do
-                f i j (A.At(i,j))
-
-    /// Fold one column.
-    let inline foldCol f acc (A: #Matrix<_>) k =
-        let mutable macc = acc
-        for i=0 to A.RowCount-1 do
-            macc <- f macc (A.Item(i,k))
-        macc
-
-    /// Fold one row.
-    let inline foldRow f acc (A: #Matrix<_>) k =
-        let mutable macc = acc
-        for i=0 to A.ColumnCount-1 do
-            macc <- f macc (A.Item(k,i))
-        macc
-
     /// Returns the sum of the results generated by applying a position dependent function to each column of the matrix.
     let inline sumColsBy f (A: #Matrix<_>) =
         A.EnumerateColumnsIndexed() |> Seq.map (fun (j,col) -> f j col) |> Seq.reduce (+)
@@ -302,6 +366,3 @@ module Matrix =
     /// Returns the sum of the results generated by applying a position dependent function to each row of the matrix.
     let inline sumRowsBy f (A: #Matrix<_>) =
         A.EnumerateRowsIndexed() |> Seq.map (fun (i,row) -> f i row)  |> Seq.reduce (+)
-
-    /// Creates a sequence that iterates the non-zero entries in the matrix.
-    let nonZeroEntries (A: #Matrix<_>) = A.EnumerateNonZero()
diff --git a/src/FSharp/LinearAlgebra.Vector.fs b/src/FSharp/LinearAlgebra.Vector.fs
index 8ad605e1..10b77f7f 100644
--- a/src/FSharp/LinearAlgebra.Vector.fs
+++ b/src/FSharp/LinearAlgebra.Vector.fs
@@ -75,13 +75,86 @@ module Vector =
     /// Transform a vector into an indexed sequence.
     let inline toSeqi (v: #Vector<_>) = v.EnumerateIndexed()
 
-    /// Transform a vector into a sequence where zero-values are skipped.
+    /// Transform a vector into a sequence where zero-values are skipped. Skipping zeros is efficient on sparse data.
     let inline toSeqnz (v: #Vector<_>) = v.EnumerateNonZero()
 
-    /// Transform a vector into an indexed sequence where zero-values are skipped.
+    /// Transform a vector into an indexed sequence where zero-values are skipped. Skipping zeros is efficient on sparse data.
     let inline toSeqinz (v: #Vector<_>) = v.EnumerateNonZeroIndexed()
 
 
+    /// Applies a function to all elements of the vector.
+    let inline iter f (v: #Vector<_>) = v.Enumerate() |> Seq.iter f
+
+    /// Applies a function to all indexed elements of the vector.
+    let inline iteri f (v: #Vector<_>) = v.Enumerate() |> Seq.iteri f
+
+    /// Applies a function to all non-zero elements of the vector. Skipping zeros is efficient on sparse data.
+    let inline iternz f (v: #Vector<_>) = v.EnumerateNonZero() |> Seq.iter f
+
+    /// Applies a function to all non-zero indexed elements of the vector. Skipping zeros is efficient on sparse data.
+    let inline iterinz f (v: #Vector<_>) = v.EnumerateNonZeroIndexed() |> Seq.iter (fun (i,x) -> f i x)
+
+
+    /// Fold all entries of a vector.
+    let inline fold f state (v: #Vector<_>) = v.Enumerate() |> Seq.fold f state
+
+    /// Fold all entries of a vector using a position dependent folding function.
+    let inline foldi f state (v: #Vector<_>) = v.EnumerateIndexed() |> Seq.fold (fun s (i,x) -> f i s x) state
+
+    /// Fold all non-zero entries of a vector. Skipping zeros is efficient on sparse data.
+    let inline foldnz f state (v: #Vector<_>) = v.EnumerateNonZero() |> Seq.fold f state
+
+    /// Fold all non-zero entries of a vector using a position dependent folding function. Skipping zeros is efficient on sparse data.
+    let inline foldinz f state (v: #Vector<_>) = v.EnumerateNonZeroIndexed() |> Seq.fold (fun s (i,x) -> f i s x) state
+
+
+    /// Scan all entries of a vector.
+    let inline scan f state (v: #Vector<_>) = v.Enumerate() |> Seq.scan f state
+
+    /// Scan all entries of a vector using a position dependent folding function.
+    let inline scani f state (v: #Vector<_>) = v.EnumerateIndexed() |> Seq.scan (fun s (i,x) -> f i s x) state
+
+    /// Scan all non-zero entries of a vector. Skipping zeros is efficient on sparse data.
+    let inline scannz f state (v: #Vector<_>) = v.EnumerateNonZero() |> Seq.scan f state
+
+    /// Scan all non-zero entries of a vector using a position dependent folding function. Skipping zeros is efficient on sparse data.
+    let inline scaninz f state (v: #Vector<_>) = v.EnumerateNonZeroIndexed() |> Seq.scan (fun s (i,x) -> f i s x) state
+
+
+    /// Reduce all entries of a vector.
+    let inline reduce f (v: #Vector<_>) = v.Enumerate() |> Seq.reduce f
+
+    /// Reduce all non-zero entries of a vector. Skipping zeros is efficient on sparse data.
+    let inline reducenz f (v: #Vector<_>) = v.EnumerateNonZero() |> Seq.reduce f
+
+
+    /// Checks whether there is an entry in the vector that satisfies a predicate.
+    let inline exists p (v: #Vector<_>) = v.Enumerate() |> Seq.exists p
+
+    /// Checks whether there is an entry in the vector that satisfies a position dependent predicate.
+    let inline existsi p (v: #Vector<_>) = v.EnumerateIndexed() |> Seq.exists (fun (i,x) -> p i x)
+
+    /// Checks whether there is a non-zero entry in the vector that satisfies a predicate. Skipping zeros is efficient on sparse data.
+    let inline existsnz p (v: #Vector<_>) = v.EnumerateNonZero() |> Seq.exists p
+
+    /// Checks whether there is a non-zero entry in the vector that satisfies a position dependent predicate. Skipping zeros is efficient on sparse data.
+    let inline existsinz p (v: #Vector<_>) = v.EnumerateNonZeroIndexed() |> Seq.exists (fun (i,x) -> p i x)
+
+
+    /// Checks whether all entries in the vector that satisfies a given predicate.
+    let inline forall p (v: #Vector<_>) = v.Enumerate() |> Seq.forall p
+
+    /// Checks whether all entries in the vector that satisfies a given position dependent predicate.
+    let inline foralli p (v: #Vector<_>) = v.EnumerateIndexed() |> Seq.forall (fun (i,x) -> p i x)
+
+    /// Checks whether all non-zero entries in the vector that satisfies a given predicate. Skipping zeros is efficient on sparse data.
+    let inline forallnz p (v: #Vector<_>) = v.EnumerateNonZero() |> Seq.forall p
+
+    /// Checks whether all non-zero entries in the vector that satisfies a given position dependent predicate. Skipping zeros is efficient on sparse data.
+    let inline forallinz p (v: #Vector<_>) = v.EnumerateNonZeroIndexed() |> Seq.forall (fun (i,x) -> p i x)
+
+
+
     /// In-place mutation by applying a function to every element of the vector.
     let inline mapInPlace f (v: #Vector<_>) =
         v.MapInplace((fun x -> f x), true)
@@ -128,110 +201,33 @@ module Vector =
         w
 
 
+
     /// In-place vector addition.
     let inline addInPlace (v: #Vector<_>) (w: #Vector<_>) = v.Add(w, v)
 
     /// In place vector subtraction.
     let inline subInPlace (v: #Vector<_>) (w: #Vector<_>) = v.Subtract(w, v)
 
-
-    /// Applies a function to all elements of the vector.
-    let inline iter f (v: #Vector<_>) = v.Enumerate() |> Seq.iter f
-
-    /// Applies a function to all indexed elements of the vector.
-    let inline iteri f (v: #Vector<_>) = v.Enumerate() |> Seq.iteri f
-
-    /// Applies a function to all non-zero elements of the vector.
-    let inline iternz f (v: #Vector<_>) = v.EnumerateNonZero() |> Seq.iter f
-
-    /// Applies a function to all non-zero indexed elements of the vector.
-    let inline iterinz f (v: #Vector<_>) = v.EnumerateNonZeroIndexed() |> Seq.iter (fun (i,v) -> f i v)
-
-
-    /// Fold all entries of a vector.
-    let inline fold f state (v: #Vector<_>) = v.Enumerate() |> Seq.fold f state
-
-    /// Fold all entries of a vector using a position dependent folding function.
-    let inline foldi f state (v: #Vector<_>) = v.EnumerateIndexed() |> Seq.fold (fun s (i,x) -> f i s x) state
-
-    /// Fold all non-zero entries of a vector.
-    let inline foldnz f state (v: #Vector<_>) = v.EnumerateNonZero() |> Seq.fold f state
-
-    /// Fold all non-zero entries of a vector using a position dependent folding function.
-    let inline foldinz f state (v: #Vector<_>) = v.EnumerateNonZeroIndexed() |> Seq.fold (fun s (i,x) -> f i s x) state
-
-
     /// Fold all entries of a vector in reverse order.
-    let inline foldBack f acc0 (v: #Vector<_>) =
-        let mutable acc = acc0
+    let inline foldBack f state (v: #Vector<_>) =
+        let mutable acc = state
         for i=2 to v.Count do
             acc <- f (v.At (v.Count - i)) acc
         acc
 
-    /// Checks whether a predicate is satisfied for every element in the vector.
-    let inline forall p (v: #Vector<_>) =
-        let mutable b = true
-        let mutable i = 0
-        while b && i < v.Count do
-            b <- b && (p (v.At i))
-            i <- i+1
-        b
-
-    /// Checks whether a predicate is true for all entries in a vector.
-    let inline foralli p (v: #Vector<_>) =
-        let mutable b = true
-        let mutable i = 0
-        while b && i < v.Count do
-            b <- b && (p i (v.At i))
-            i <- i+1
-        b
-
-    /// Checks whether there is an entry in the vector that satisfies a given predicate.
-    let inline exists p (v: #Vector<_>) =
-        let mutable b = false
-        let mutable i = 0
-        while not(b) && i < v.Count do
-            b <- b || (p (v.At i))
-            i <- i+1
-        b
-
-    /// Checks whether there is an entry in the vector that satisfies a given position dependent predicate.
-    let inline existsi p (v: #Vector<_>) =
-        let mutable b = false
-        let mutable i = 0
-        while not(b) && i < v.Count do
-            b <- b || (p i (v.At i))
-            i <- i+1
-        b
-
-    /// Scans a vector; like fold but returns the intermediate result.
-    let inline scan f (v: #Vector<_>) =
-        let w = v.Clone()
-        let mutable p = v.Item(0)
-        for i=1 to v.Count-1 do
-            p <- f p (v.At i)
-            w.At(i, p)
-        w
-
-    /// Scans a vector in reverse order; like foldBack but returns the intermediate result.
-    let inline scanBack f (v: #Vector<_>) =
-        let w = v.Clone()
-        let mutable p = v.At (v.Count-1)
-        for i=2 to v.Count do
-            p <- f (v.At (v.Count - i)) p
-            w.At(v.Count - i, p)
-        w
-
-    /// Reduces a vector: the result of this function will be f(...f(f(v[0],v[1]), v[2]),..., v[n]).
-    let inline reduce f (v: #Vector<_>) =
-        let mutable p = v.Item(0)
-        for i=1 to v.Count-1 do
-            p <- f p (v.At i)
-        p
-
     /// Reduces a vector in reverse order: the result of this function will be f(v[1], ..., f(v[n-2], f(v[n-1],v[n]))...).
     let inline reduceBack f (v: #Vector<_>) =
         let mutable p = v.Item(v.Count-1)
         for i=2 to v.Count do
             p <- f (v.At (v.Count - i)) p
         p
+
+    /// Scans a vector in reverse order; like foldBack but returns the intermediate result.
+    let inline scanBack f state (v: #Vector<_>) =
+        seq {
+            let rstate = ref state
+            yield !rstate
+            for i in v.Count-1..-1..0 do
+                rstate := f (v.At(i)) !rstate
+                yield !rstate
+        }
diff --git a/src/FSharpUnitTests/MatrixTests.fs b/src/FSharpUnitTests/MatrixTests.fs
index 2132d3bd..4e24259c 100644
--- a/src/FSharpUnitTests/MatrixTests.fs
+++ b/src/FSharpUnitTests/MatrixTests.fs
@@ -11,14 +11,14 @@ module MatrixTests =
     let approximately_equal tolerance = equalWithin (10.0 ** (float -tolerance))
 
     /// A small uniform matrix.
-    let smallM = DenseMatrix.OfArray( Array2D.create 2 2 0.3 )
+    let smallM = DenseMatrix.ofArray2 (Array2D.create 2 2 0.3)
     let failingFoldBackM = DenseMatrix.init 2 3 (fun i j -> 1.0)
 
     /// A small sparse matrix.
     let sparseM = SparseMatrix.ofListi 2 3 [(1,0,0.3)]
 
     /// A large matrix with increasingly large entries
-    let largeM = DenseMatrix.OfArray( Array2D.init 100 100 (fun i j -> float i * 100.0 + float j) )
+    let largeM = DenseMatrix.init 100 100 (fun i j -> float i * 100.0 + float j)
 
     [<Test>]
     let ``Matrix.GetSlice`` () =
@@ -149,8 +149,8 @@ module MatrixTests =
         N |> should equal (0.0 * smallM)
 
     [<Test>]
-    let ``Matrix.nonZeroEntries`` () =
-        Seq.length (Matrix.nonZeroEntries smallM) |> should equal 4
+    let ``Matrix.toSeqnz`` () =
+        Seq.length (Matrix.toSeqnz smallM) |> should equal 4
 
     [<Test>]
     let ``Matrix.sum`` () =
diff --git a/src/FSharpUnitTests/VectorTests.fs b/src/FSharpUnitTests/VectorTests.fs
index baed54b0..258eec1d 100644
--- a/src/FSharpUnitTests/VectorTests.fs
+++ b/src/FSharpUnitTests/VectorTests.fs
@@ -11,13 +11,13 @@ module VectorTests =
     let approximately_equal tolerance = equalWithin (10.0 ** (float -tolerance))
 
     /// A small uniform vector.
-    let smallv = DenseVector([|0.3;0.3;0.3;0.3;0.3|])
+    let smallv = DenseVector.raw [|0.3;0.3;0.3;0.3;0.3|]
 
     /// A small sparse vector.
     let sparsev = SparseVector.ofListi 5 [(1,0.3)]
 
     /// A large vector with increasingly large entries
-    let largev = DenseVector(Array.init 100 (fun i -> float i / 100.0))
+    let largev = DenseVector.init 100 (fun i -> float i / 100.0)
 
     [<Test>]
     let ``Vector.GetSlice`` () =
@@ -147,11 +147,11 @@ module VectorTests =
 
     [<Test>]
     let ``Vector.scan`` () =
-        Vector.scan (fun acc x -> acc + x) smallv |> should (approximately_equal 14) (new DenseVector( [|0.3;0.6;0.9;1.2;1.5|] ) :> Vector<float>)
+        Vector.scan (fun acc x -> acc + x) 0.0 smallv |> should (approximately_equal 14) (DenseVector.raw [|0.0;0.3;0.6;0.9;1.2;1.5|])
 
     [<Test>]
     let ``Vector.scanBack`` () =
-        Vector.scanBack (fun x acc -> acc + x) smallv |> should (approximately_equal 14) (new DenseVector( [|1.5;1.2;0.9;0.6;0.3|] ) :> Vector<float>)
+        Vector.scanBack (fun x acc -> acc + x) 0.0 smallv |> should (approximately_equal 14) (DenseVector.raw [|0.0;0.3;0.6;0.9;1.2;1.5|])
 
     [<Test>]
     let ``Vector.reduce`` () =