diff --git a/src/Numerics/Fit.cs b/src/Numerics/Fit.cs
index f91695d2..2a829c5e 100644
--- a/src/Numerics/Fit.cs
+++ b/src/Numerics/Fit.cs
@@ -66,9 +66,20 @@ namespace MathNet.Numerics
         /// returning its best fitting parameters as [p0, p1, p2, ..., pk] array.
         /// If an intercept is added, its coefficient will be prepended to the resulting parameters.
         /// </summary>
-        public static double[] MultiDim(double[][] x, double[] y, bool intercept = false)
+        public static double[] MultiDim(double[][] x, double[] y, bool intercept = false, DirectRegressionMethod method = DirectRegressionMethod.NormalEquations)
         {
-            return MultipleRegression.NormalEquations(x, y, intercept);
+            return MultipleRegression.DirectMethod(x, y, intercept, method);
+        }
+
+        /// <summary>
+        /// Least-Squares fitting the points (X,y) = ((x0,x1,..,xk),y) to a linear surface y : X -> p0*x0 + p1*x1 + ... + pk*xk,
+        /// returning a function y' for the best fitting combination.
+        /// If an intercept is added, its coefficient will be prepended to the resulting parameters.
+        /// </summary>
+        public static Func<double[], double> MultiDimFunc(double[][] x, double[] y, bool intercept = false, DirectRegressionMethod method = DirectRegressionMethod.NormalEquations)
+        {
+            var parameters = MultipleRegression.DirectMethod(x, y, intercept, method);
+            return z => Control.LinearAlgebraProvider.DotProduct(parameters, z);
         }
 
         /// <summary>
@@ -81,23 +92,23 @@ namespace MathNet.Numerics
         }
 
         /// <summary>
-        /// Least-Squares fitting the points (X,y) = ((x0,x1,..,xk),y) to a linear surface y : X -> p0*x0 + p1*x1 + ... + pk*xk,
-        /// returning a function y' for the best fitting combination.
+        /// Least-Squares fitting the points (x,y) to a k-order polynomial y : x -> p0 + p1*x + p2*x^2 + ... + pk*x^k,
+        /// returning its best fitting parameters as [p0, p1, p2, ..., pk] array, compatible with Evaluate.Polynomial.
         /// </summary>
-        public static Func<double[], double> MultiDimFunc(double[][] x, double[] y)
+        public static double[] Polynomial(double[] x, double[] y, int order, DirectRegressionMethod method = DirectRegressionMethod.QR)
         {
-            var parameters = MultipleRegression.NormalEquations(x, y);
-            return z => Control.LinearAlgebraProvider.DotProduct(parameters, z);
+            var design = Matrix<double>.Build.Dense(x.Length, order + 1, (i, j) => Math.Pow(x[i], j));
+            return MultipleRegression.DirectMethod(design, Vector<double>.Build.Dense(y), method).ToArray();
         }
 
         /// <summary>
         /// Least-Squares fitting the points (x,y) to a k-order polynomial y : x -> p0 + p1*x + p2*x^2 + ... + pk*x^k,
-        /// returning its best fitting parameters as [p0, p1, p2, ..., pk] array, compatible with Evaluate.Polynomial.
+        /// returning a function y' for the best fitting polynomial.
         /// </summary>
-        public static double[] Polynomial(double[] x, double[] y, int order)
+        public static Func<double, double> PolynomialFunc(double[] x, double[] y, int order, DirectRegressionMethod method = DirectRegressionMethod.QR)
         {
-            var design = Matrix<double>.Build.Dense(x.Length, order + 1, (i, j) => Math.Pow(x[i], j));
-            return MultipleRegression.QR(design, Vector<double>.Build.Dense(y)).ToArray();
+            var parameters = Polynomial(x, y, order, method);
+            return z => Evaluate.Polynomial(z, parameters);
         }
 
         /// <summary>
@@ -111,32 +122,42 @@ namespace MathNet.Numerics
         }
 
         /// <summary>
-        /// Least-Squares fitting the points (x,y) to a k-order polynomial y : x -> p0 + p1*x + p2*x^2 + ... + pk*x^k,
-        /// returning a function y' for the best fitting polynomial.
+        /// Least-Squares fitting the points (x,y) to an arbitrary linear combination y : x -> p0*f0(x) + p1*f1(x) + ... + pk*fk(x),
+        /// returning its best fitting parameters as [p0, p1, p2, ..., pk] array.
         /// </summary>
-        public static Func<double, double> PolynomialFunc(double[] x, double[] y, int order)
+        public static double[] LinearCombination(double[] x, double[] y, params Func<double,double>[] functions)
         {
-            var parameters = Polynomial(x, y, order);
-            return z => Evaluate.Polynomial(z, parameters);
+            var design = Matrix<double>.Build.Dense(x.Length, functions.Length, (i, j) => functions[j](x[i]));
+            return MultipleRegression.QR(design, Vector<double>.Build.Dense(y)).ToArray();
+        }
+
+        /// <summary>
+        /// Least-Squares fitting the points (x,y) to an arbitrary linear combination y : x -> p0*f0(x) + p1*f1(x) + ... + pk*fk(x),
+        /// returning a function y' for the best fitting combination.
+        /// </summary>
+        public static Func<double, double> LinearCombinationFunc(double[] x, double[] y, params Func<double, double>[] functions)
+        {
+            var parameters = LinearCombination(x, y, functions);
+            return z => functions.Zip(parameters, (f, p) => p*f(z)).Sum();
         }
 
         /// <summary>
         /// Least-Squares fitting the points (x,y) to an arbitrary linear combination y : x -> p0*f0(x) + p1*f1(x) + ... + pk*fk(x),
         /// returning its best fitting parameters as [p0, p1, p2, ..., pk] array.
         /// </summary>
-        public static double[] LinearCombination(double[] x, double[] y, params Func<double,double>[] functions)
+        public static double[] LinearCombination(double[] x, double[] y, DirectRegressionMethod method, params Func<double, double>[] functions)
         {
             var design = Matrix<double>.Build.Dense(x.Length, functions.Length, (i, j) => functions[j](x[i]));
-            return MultipleRegression.QR(design, Vector<double>.Build.Dense(y)).ToArray();
+            return MultipleRegression.DirectMethod(design, Vector<double>.Build.Dense(y), method).ToArray();
         }
 
         /// <summary>
         /// Least-Squares fitting the points (x,y) to an arbitrary linear combination y : x -> p0*f0(x) + p1*f1(x) + ... + pk*fk(x),
         /// returning a function y' for the best fitting combination.
         /// </summary>
-        public static Func<double, double> LinearCombinationFunc(double[] x, double[] y, params Func<double, double>[] functions)
+        public static Func<double, double> LinearCombinationFunc(double[] x, double[] y, DirectRegressionMethod method, params Func<double, double>[] functions)
         {
-            var parameters = LinearCombination(x, y, functions);
+            var parameters = LinearCombination(x, y, method, functions);
             return z => functions.Zip(parameters, (f, p) => p*f(z)).Sum();
         }
 
@@ -160,6 +181,26 @@ namespace MathNet.Numerics
             return z => functions.Zip(parameters, (f, p) => p * f(z)).Sum();
         }
 
+        /// <summary>
+        /// Least-Squares fitting the points (X,y) = ((x0,x1,..,xk),y) to an arbitrary linear combination y : X -> p0*f0(x) + p1*f1(x) + ... + pk*fk(x),
+        /// returning its best fitting parameters as [p0, p1, p2, ..., pk] array.
+        /// </summary>
+        public static double[] LinearMultiDim(double[][] x, double[] y, DirectRegressionMethod method, params Func<double[], double>[] functions)
+        {
+            var design = Matrix<double>.Build.Dense(x.Length, functions.Length, (i, j) => functions[j](x[i]));
+            return MultipleRegression.DirectMethod(design, Vector<double>.Build.Dense(y), method).ToArray();
+        }
+
+        /// <summary>
+        /// Least-Squares fitting the points (X,y) = ((x0,x1,..,xk),y) to an arbitrary linear combination y : X -> p0*f0(x) + p1*f1(x) + ... + pk*fk(x),
+        /// returning a function y' for the best fitting combination.
+        /// </summary>
+        public static Func<double[], double> LinearMultiDimFunc(double[][] x, double[] y, DirectRegressionMethod method, params Func<double[], double>[] functions)
+        {
+            var parameters = LinearMultiDim(x, y, method, functions);
+            return z => functions.Zip(parameters, (f, p) => p * f(z)).Sum();
+        }
+
         /// <summary>
         /// Least-Squares fitting the points (T,y) = (T,y) to an arbitrary linear combination y : X -> p0*f0(T) + p1*f1(T) + ... + pk*fk(T),
         /// returning its best fitting parameters as [p0, p1, p2, ..., pk] array.
@@ -179,5 +220,25 @@ namespace MathNet.Numerics
             var parameters = LinearGeneric(x, y, functions);
             return z => functions.Zip(parameters, (f, p) => p * f(z)).Sum();
         }
+
+        /// <summary>
+        /// Least-Squares fitting the points (T,y) = (T,y) to an arbitrary linear combination y : X -> p0*f0(T) + p1*f1(T) + ... + pk*fk(T),
+        /// returning its best fitting parameters as [p0, p1, p2, ..., pk] array.
+        /// </summary>
+        public static double[] LinearGeneric<T>(T[] x, double[] y, DirectRegressionMethod method, params Func<T, double>[] functions)
+        {
+            var design = Matrix<double>.Build.Dense(x.Length, functions.Length, (i, j) => functions[j](x[i]));
+            return MultipleRegression.DirectMethod(design, Vector<double>.Build.Dense(y), method).ToArray();
+        }
+
+        /// <summary>
+        /// Least-Squares fitting the points (T,y) = (T,y) to an arbitrary linear combination y : X -> p0*f0(T) + p1*f1(T) + ... + pk*fk(T),
+        /// returning a function y' for the best fitting combination.
+        /// </summary>
+        public static Func<T, double> LinearGenericFunc<T>(T[] x, double[] y, DirectRegressionMethod method, params Func<T, double>[] functions)
+        {
+            var parameters = LinearGeneric(x, y, method, functions);
+            return z => functions.Zip(parameters, (f, p) => p * f(z)).Sum();
+        }
     }
 }