diff --git a/math/src/main/scala/breeze/stats/regression/Lasso.scala b/math/src/main/scala/breeze/stats/regression/Lasso.scala
index 66f40b521..29bd11666 100644
--- a/math/src/main/scala/breeze/stats/regression/Lasso.scala
+++ b/math/src/main/scala/breeze/stats/regression/Lasso.scala
@@ -91,7 +91,7 @@ private case class LassoCalculator(
     r2
   }
 
-  private def estimateOneColumn(column: Int): LeastSquaresRegressionResult = {
+  private def estimateOneColumn(column: Int): LeastSquaresRegressionResult[Double] = {
     /*
      * Goal of this routine is to use the specified column to explain as much of the residual
      * as possible, after using the already specified values in other columns.
diff --git a/math/src/main/scala/breeze/stats/regression/LeastSquares.scala b/math/src/main/scala/breeze/stats/regression/LeastSquares.scala
index 6f1817bd0..cb6c45989 100644
--- a/math/src/main/scala/breeze/stats/regression/LeastSquares.scala
+++ b/math/src/main/scala/breeze/stats/regression/LeastSquares.scala
@@ -2,15 +2,19 @@ package breeze.stats.regression
 
 import breeze.generic.UFunc
 import breeze.linalg._
+import breeze.linalg.operators.{OpMulInner, OpMulMatrix}
 import org.netlib.util.intW
 import dev.ludovic.netlib.lapack.LAPACK.{getInstance => lapack}
+
 import java.util.Arrays
 
 private object leastSquaresImplementation {
-  def doLeastSquares(
+
+
+  def doLeastSquaresDouble(
       data: DenseMatrix[Double],
       outputs: DenseVector[Double],
-      workArray: Array[Double]): LeastSquaresRegressionResult = {
+      workArray: Array[Double]): LeastSquaresRegressionResult[Double] = {
     require(data.rows == outputs.size)
     require(data.rows > data.cols + 1)
     require(workArray.length >= 2 * data.rows * data.cols)
@@ -37,69 +41,183 @@ private object leastSquaresImplementation {
     for (i <- 0 until (data.rows - data.cols)) {
       r2 = r2 + math.pow(outputs.data(data.cols + i), 2)
     }
-    LeastSquaresRegressionResult(coefficients, r2)
+    LeastSquaresRegressionResult[Double](coefficients, r2)
+  }
+
+  def doLeastSquaresFloat(
+                            data: DenseMatrix[Float],
+                            outputs: DenseVector[Float],
+                            workArray: Array[Float]): LeastSquaresRegressionResult[Float] = {
+    require(data.rows == outputs.size)
+    require(data.rows > data.cols + 1)
+    require(workArray.length >= 2 * data.rows * data.cols)
+
+    val info = new intW(0)
+    lapack.sgels(
+      "N",
+      data.rows,
+      data.cols,
+      1,
+      data.data,
+      data.rows,
+      outputs.data,
+      data.rows,
+      workArray,
+      workArray.length,
+      info)
+    if (info.`val` < 0) {
+      throw new ArithmeticException("Least squares did not converge.")
+    }
+
+    val coefficients = new DenseVector[Float](Arrays.copyOf(outputs.data, data.cols))
+    var r2 = 0.0
+    for (i <- 0 until (data.rows - data.cols)) {
+      r2 = r2 + math.pow(outputs.data(data.cols + i), 2)
+    }
+    LeastSquaresRegressionResult[Float](coefficients, r2.toFloat)
   }
 }
 
-case class LeastSquaresRegressionResult(coefficients: DenseVector[Double], rSquared: Double)
-    extends RegressionResult[DenseVector[Double], Double] {
-  def apply(x: DenseVector[Double]): Double = coefficients.dot(x)
 
-  def apply(X: DenseMatrix[Double]): DenseVector[Double] = X * coefficients
+
+
+
+case class LeastSquaresRegressionResult[T](coefficients: DenseVector[T], rSquared: T)(implicit can_dot: OpMulInner.Impl2[DenseVector[T], DenseVector[T], T])
+  extends RegressionResult[DenseVector[T], T] {
+
+  def apply(x: DenseVector[T]): T = coefficients.dot(x)
+
+
+  def apply(X: DenseMatrix[T])(implicit can_dot: OpMulMatrix.Impl2[DenseMatrix[T], DenseVector[T], DenseVector[T]]): DenseVector[T]  = {
+       X * coefficients
+    }
 }
 
+
 object leastSquares extends UFunc {
-  implicit val matrixVectorWithWorkArray
-    : Impl3[DenseMatrix[Double], DenseVector[Double], Array[Double], LeastSquaresRegressionResult] =
-    new Impl3[DenseMatrix[Double], DenseVector[Double], Array[Double], LeastSquaresRegressionResult] {
+  implicit val matrixVectorWithWorkArrayDouble
+    : Impl3[DenseMatrix[Double], DenseVector[Double], Array[Double], LeastSquaresRegressionResult[Double]] =
+    new Impl3[DenseMatrix[Double], DenseVector[Double], Array[Double], LeastSquaresRegressionResult[Double]] {
       def apply(
           data: DenseMatrix[Double],
           outputs: DenseVector[Double],
-          workArray: Array[Double]): LeastSquaresRegressionResult =
-        leastSquaresImplementation.doLeastSquares(data.copy, outputs.copy, workArray)
+          workArray: Array[Double]): LeastSquaresRegressionResult[Double] =
+        leastSquaresImplementation.doLeastSquaresDouble(data.copy, outputs.copy, workArray)
+    }
+
+  implicit val matrixVectorWithWorkArrayFloat
+  : Impl3[DenseMatrix[Float], DenseVector[Float], Array[Float], LeastSquaresRegressionResult[Float]] =
+    new Impl3[DenseMatrix[Float], DenseVector[Float], Array[Float], LeastSquaresRegressionResult[Float]] {
+      def apply(
+                 data: DenseMatrix[Float],
+                 outputs: DenseVector[Float],
+                 workArray: Array[Float]): LeastSquaresRegressionResult[Float] =
+        leastSquaresImplementation.doLeastSquaresFloat(data.copy, outputs.copy, workArray)
     }
 
-  implicit val matrixVectorSpecifiedWork
-    : Impl3[DenseMatrix[Double], DenseVector[Double], Int, LeastSquaresRegressionResult] =
-    new Impl3[DenseMatrix[Double], DenseVector[Double], Int, LeastSquaresRegressionResult] {
-      def apply(data: DenseMatrix[Double], outputs: DenseVector[Double], workSize: Int): LeastSquaresRegressionResult =
-        leastSquaresImplementation.doLeastSquares(data.copy, outputs.copy, new Array[Double](workSize))
+  implicit val matrixVectorSpecifiedWorkDouble
+    : Impl3[DenseMatrix[Double], DenseVector[Double], Int, LeastSquaresRegressionResult[Double]] =
+    new Impl3[DenseMatrix[Double], DenseVector[Double], Int, LeastSquaresRegressionResult[Double]] {
+      def apply(
+                 data: DenseMatrix[Double],
+                 outputs: DenseVector[Double],
+                 workSize: Int): LeastSquaresRegressionResult[Double] =
+        leastSquaresImplementation.doLeastSquaresDouble(data.copy, outputs.copy, new Array[Double](workSize))
+    }
+
+  implicit val matrixVectorSpecifiedWorkFloat
+  : Impl3[DenseMatrix[Float], DenseVector[Float], Int, LeastSquaresRegressionResult[Float]] =
+    new Impl3[DenseMatrix[Float], DenseVector[Float], Int, LeastSquaresRegressionResult[Float]] {
+      def apply(
+                 data: DenseMatrix[Float],
+                 outputs: DenseVector[Float],
+                 workSize: Int): LeastSquaresRegressionResult[Float] =
+        leastSquaresImplementation.doLeastSquaresFloat(data.copy, outputs.copy, new Array[Float](workSize))
     }
 
-  implicit val matrixVector: Impl2[DenseMatrix[Double], DenseVector[Double], LeastSquaresRegressionResult] =
-    new Impl2[DenseMatrix[Double], DenseVector[Double], LeastSquaresRegressionResult] {
-      def apply(data: DenseMatrix[Double], outputs: DenseVector[Double]): LeastSquaresRegressionResult =
-        leastSquaresImplementation.doLeastSquares(
+  implicit val matrixVectorDouble: Impl2[DenseMatrix[Double], DenseVector[Double], LeastSquaresRegressionResult[Double]] =
+    new Impl2[DenseMatrix[Double], DenseVector[Double], LeastSquaresRegressionResult[Double]] {
+      def apply(
+                 data: DenseMatrix[Double],
+                 outputs: DenseVector[Double]): LeastSquaresRegressionResult[Double] =
+        leastSquaresImplementation.doLeastSquaresDouble(
           data.copy,
           outputs.copy,
           new Array[Double](math.max(1, data.rows * data.cols * 2)))
     }
+
+  implicit val matrixVectorFloat: Impl2[DenseMatrix[Float], DenseVector[Float], LeastSquaresRegressionResult[Float]] =
+    new Impl2[DenseMatrix[Float], DenseVector[Float], LeastSquaresRegressionResult[Float]] {
+      def apply(
+                 data: DenseMatrix[Float],
+                 outputs: DenseVector[Float]): LeastSquaresRegressionResult[Float] =
+        leastSquaresImplementation.doLeastSquaresFloat(
+          data,
+          outputs,
+          new Array[Float](math.max(1, data.rows * data.cols * 2)))
+    }
 }
 
 object leastSquaresDestructive extends UFunc {
-  implicit val matrixVectorWithWorkArray
-    : Impl3[DenseMatrix[Double], DenseVector[Double], Array[Double], LeastSquaresRegressionResult] =
-    new Impl3[DenseMatrix[Double], DenseVector[Double], Array[Double], LeastSquaresRegressionResult] {
+  implicit val matrixVectorWithWorkArrayDouble
+    : Impl3[DenseMatrix[Double], DenseVector[Double], Array[Double], LeastSquaresRegressionResult[Double]] =
+    new Impl3[DenseMatrix[Double], DenseVector[Double], Array[Double], LeastSquaresRegressionResult[Double]] {
       def apply(
           data: DenseMatrix[Double],
           outputs: DenseVector[Double],
-          workArray: Array[Double]): LeastSquaresRegressionResult =
-        leastSquaresImplementation.doLeastSquares(data, outputs, workArray)
+          workArray: Array[Double]): LeastSquaresRegressionResult[Double] =
+        leastSquaresImplementation.doLeastSquaresDouble(data, outputs, workArray)
+    }
+
+  implicit val matrixVectorWithWorkArrayFloat
+  : Impl3[DenseMatrix[Float], DenseVector[Float], Array[Float], LeastSquaresRegressionResult[Float]] =
+    new Impl3[DenseMatrix[Float], DenseVector[Float], Array[Float], LeastSquaresRegressionResult[Float]] {
+      def apply(
+                 data: DenseMatrix[Float],
+                 outputs: DenseVector[Float],
+                 workArray: Array[Float]): LeastSquaresRegressionResult[Float] =
+        leastSquaresImplementation.doLeastSquaresFloat(data, outputs, workArray)
     }
 
-  implicit val matrixVectorSpecifiedWork
-    : Impl3[DenseMatrix[Double], DenseVector[Double], Int, LeastSquaresRegressionResult] =
-    new Impl3[DenseMatrix[Double], DenseVector[Double], Int, LeastSquaresRegressionResult] {
-      def apply(data: DenseMatrix[Double], outputs: DenseVector[Double], workSize: Int): LeastSquaresRegressionResult =
-        leastSquaresImplementation.doLeastSquares(data, outputs, new Array[Double](workSize))
+  implicit val matrixVectorSpecifiedWorkDouble
+    : Impl3[DenseMatrix[Double], DenseVector[Double], Int, LeastSquaresRegressionResult[Double]] =
+    new Impl3[DenseMatrix[Double], DenseVector[Double], Int, LeastSquaresRegressionResult[Double]] {
+      def apply(
+                 data: DenseMatrix[Double],
+                 outputs: DenseVector[Double],
+                 workSize: Int): LeastSquaresRegressionResult[Double] =
+        leastSquaresImplementation.doLeastSquaresDouble(data, outputs, new Array[Double](workSize))
     }
 
-  implicit val matrixVector: Impl2[DenseMatrix[Double], DenseVector[Double], LeastSquaresRegressionResult] =
-    new Impl2[DenseMatrix[Double], DenseVector[Double], LeastSquaresRegressionResult] {
-      def apply(data: DenseMatrix[Double], outputs: DenseVector[Double]): LeastSquaresRegressionResult =
-        leastSquaresImplementation.doLeastSquares(
+  implicit val matrixVectorSpecifiedWorkFloat
+  : Impl3[DenseMatrix[Float], DenseVector[Float], Int, LeastSquaresRegressionResult[Float]] =
+    new Impl3[DenseMatrix[Float], DenseVector[Float], Int, LeastSquaresRegressionResult[Float]] {
+      def apply(
+                 data: DenseMatrix[Float],
+                outputs: DenseVector[Float],
+                workSize: Int): LeastSquaresRegressionResult[Float] =
+        leastSquaresImplementation.doLeastSquaresFloat(data, outputs, new Array[Float](workSize))
+    }
+
+  implicit val matrixVectorDouble: Impl2[DenseMatrix[Double], DenseVector[Double], LeastSquaresRegressionResult[Double]] =
+    new Impl2[DenseMatrix[Double], DenseVector[Double], LeastSquaresRegressionResult[Double]] {
+      def apply(
+                 data: DenseMatrix[Double],
+                outputs: DenseVector[Double]): LeastSquaresRegressionResult[Double] =
+        leastSquaresImplementation.doLeastSquaresDouble(
           data,
           outputs,
           new Array[Double](math.max(1, data.rows * data.cols * 2)))
     }
+
+  implicit val matrixVectorFloat: Impl2[DenseMatrix[Float], DenseVector[Float], LeastSquaresRegressionResult[Float]] =
+    new Impl2[DenseMatrix[Float], DenseVector[Float], LeastSquaresRegressionResult[Float]] {
+      def apply(
+                 data: DenseMatrix[Float],
+                 outputs: DenseVector[Float]): LeastSquaresRegressionResult[Float] =
+        leastSquaresImplementation.doLeastSquaresFloat(
+          data,
+          outputs,
+          new Array[Float](math.max(1, data.rows * data.cols * 2)))
+    }
 }