Fix for #864 - Adding tolerance to pseudo inverse matrix (#865)

Author: raggaraluz

math/src/main/codegen/breeze/linalg/pinv.scala

@@ -14,18 +14,40 @@ import breeze.macros.expand
  * Solving A^T (AX-B) = 0 for X yields
  * A^T AX = A^T B
  * =>      X = (A^T A)^(-1) A^T B
+ *
+ * @param v: Matrix to be pseudo-inverted
+ * @param rcond: Cutoff for small singular values. Singular values less than or equal to rcond * largest_singular_value
+ *               are set to zero. Default: 1e-15. To deactivate this option, set rcond to zero.
  */
 object pinv extends UFunc with pinvLowPrio {
 
+  private val DEFAULT_RCOND = 1e-15f
+
   @expand
   @expand.valify
   implicit def pinvFromSVD[@expand.args(Float, Double) T]: Impl[DenseMatrix[T], DenseMatrix[T]] = {
     new Impl[DenseMatrix[T], DenseMatrix[T]] {
-      // http://en.wikipedia.org/wiki/Singular_value_decomposition#Applications_of_the_SVD
+      // No rcond passed as parameter, use default value
       override def apply(v: DenseMatrix[T]): DenseMatrix[T] = {
+          val rcond: T = DEFAULT_RCOND
+          pinv(v, rcond)
+      }
+    }
+  }
+
+
+  @expand
+  @expand.valify
+  implicit def pinvFromSVDRcond[@expand.args(Float, Double) T]: Impl2[DenseMatrix[T], T, DenseMatrix[T]] = {
+    new Impl2[DenseMatrix[T], T, DenseMatrix[T]] {
+      // http://en.wikipedia.org/wiki/Singular_value_decomposition#Applications_of_the_SVD
+      override def apply(v: DenseMatrix[T], rcond: T): DenseMatrix[T] = {
+        require(rcond >= 0, "rcond must be non-negative")
+
         val svd.SVD(s, svs, d) = svd(v)
+        val cutoff = max(svs) * rcond
         val vi = svs.map { v =>
-          if (v == 0.0) 0 else 1 / v
+          if (v <= cutoff) 0 else 1 / v
         }
 
         val svDiag = DenseMatrix.tabulate[T](s.cols, d.rows) { (i, j) =>
@@ -58,12 +80,12 @@ trait pinvLowPrio { self: pinv.type =>
    * @return
    */
   implicit def implFromTransposeAndSolve[T, TransT, MulRes, Result](
-                                                                     implicit numericT: T => NumericOps[T],
-                                                                     trans: CanTranspose[T, TransT],
-                                                                     numericTrans: TransT => NumericOps[TransT],
-                                                                     mul: OpMulMatrix.Impl2[TransT, T, MulRes],
-                                                                     numericMulRes: MulRes => NumericOps[MulRes],
-                                                                     solve: OpSolveMatrixBy.Impl2[MulRes, TransT, Result]): Impl[T, Result] = {
+      implicit numericT: T => NumericOps[T],
+      trans: CanTranspose[T, TransT],
+      numericTrans: TransT => NumericOps[TransT],
+      mul: OpMulMatrix.Impl2[TransT, T, MulRes],
+      numericMulRes: MulRes => NumericOps[MulRes],
+      solve: OpSolveMatrixBy.Impl2[MulRes, TransT, Result]): Impl[T, Result] = {
     new Impl[T, Result] {
       def apply(X: T): Result = {
         (X.t * X) \ X.t

math/src/test/scala/breeze/linalg/LinearAlgebraTest.scala

@@ -159,6 +159,27 @@ class LinearAlgebraTest extends AnyFunSuite with Checkers with DoubleImplicits {
     matricesNearlyEqual(mi * m, eye)
   }
 
+  test("pinv test singular values cutoff doubles: #864") {
+    val m = DenseMatrix((0.1d, 0.1d, 0d), (0.1d, 0.1d, 0d), (0d, 0d, 0d))
+    val mi = pinv(m)
+    val expected = DenseMatrix((2.5d, 2.5d, 0d), (2.5d, 2.5d, 0d), (0d, 0d, 0d))
+    matricesNearlyEqual(mi, expected)
+  }
+
+  test("pinv test singular values cutoff floats: #864") {
+    // Cutoff needs to be to a more strict value than for Doubles
+    val m = DenseMatrix((0.1f, 0.1f, 0f), (0.1f, 0.1f, 0f), (0f, 0f, 0f))
+    val mi = pinv(m, 1e-7f)
+    val expected = DenseMatrix((2.5f, 2.5f, 0f), (2.5f, 2.5f, 0f), (0f, 0f, 0f))
+    matricesNearlyEqual_Float(mi, expected)
+  }
+
+  test("pinv test bad conditioned matrix, no cutoff: #864") {
+    val m = DenseMatrix((0.1d, 0.1d, 0d), (0.1d, 0.1d, 0d), (0d, 0d, 0d))
+    val mi = pinv(m, 0.0)
+    mi(0,0) should be > 1e10
+  }
+
   test("cross") {
     // specific example; with prime elements
     val (v1, v2, r) = (DenseVector(13, 3, 7), DenseVector(5, 11, 17), DenseVector(-26, -186, 128))