Adding Moore-Penrose pseudo-inverse of a matrix

lib/nmatrix/math.rb

@@ -112,6 +112,79 @@ def invert
   end
   alias :inverse :invert
 
+
+  #
+  # call-seq:
+  #     pinv -> NMatrix
+  #
+  #  Compute the (Moore-Penrose) pseudo-inverse of a matrix.
+  #  Works with ATLAS.
+  #
+  #  Calculate the generalized inverse of a matrix using its
+  #  singular-value decomposition (SVD).
+  #
+  # * *Arguments* :
+  #  - +a+ -> (M, N) array_like
+  #      Matrix to be pseudo-inverted.
+  #
+  #  - +rcond+ -> float(optional)
+  #      Cutoff for small singular values.Singular values smaller (in modulus) than
+  #      rcond * largest_singular_value (again, in modulus)
+  #      are set to zero.
+  #
+  # * *Returns* :
+  #   -  Pseudo-inverse matrix.
+  #
+  # * *Raises* :
+  #   - +NotImplementedError+ -> If called without nmatrix-atlas or nmatrix-lapacke gem
+  #   - +TypeError+ -> If called without float data type
+  # === Usage
+  #
+  #  a = NMatrix.new([2,2],[1,2,
+  #                         3,4], dtype: :float64)
+  #  a.pinv # => [ [-2.0000000000000018, 1.0000000000000007]
+  #                [1.5000000000000016, -0.5000000000000008] ]
+  #
+  #  b = NMatrix.new([4,1],[1,2,3,4], dtype: :float64)
+  #  b.pinv # => [ [ 0.03333333, 0.06666667, 0.99999999, 0.13333333] ]
+  #
+  # == References
+  #
+  # * https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_pseudoinverse
+  # * G. Strang, Linear Algebra and Its Applications, 2nd Ed., Orlando,
+  #   FL, Academic Press
+  def pinv(rcond = 1e-15)
+    raise TypeError, "Works with float matrices only" unless
+      [:float64,:float32].include?(dtype)
+    u, s, vt = self.gesvd # singular value decomposition
+    rows = self.shape[0]
+    cols = self.shape[1]
+    if rows < cols
+      u_reduced = u
+      vt_reduced = vt[0..rows - 1, 0..cols - 1].transpose
+    else
+      u_reduced = u[0..rows-1,0..cols-1]
+      vt_reduced = vt.transpose
+    end
+    largest_singular_value = s[0]
+    (0...s.shape[0]).each do |i|
+      largest_singular_value = s[i]   if largest_singular_value < s[i]
+    end
+    cutoff = rcond * largest_singular_value
+    (0...[rows, cols].min).each do |i|
+      s[i] = 1 / s[i] if s[i] > cutoff
+      s[i] = 0        if s[i] <= cutoff
+    end
+    multiplier = u_reduced.transpose
+    0.upto(multiplier.rows - 1) do |i|
+      0.upto(multiplier.cols - 1) do |j|
+        multiplier[i, j] = multiplier[i, j] * s[i]
+      end
+    end
+    vt_reduced.dot(multiplier)
+  end
+
+
   #
   # call-seq:
   #     getrf! -> Array

spec/math_spec.rb

@@ -332,6 +332,21 @@
     end
   end
 
+  [:float32,:float64].each do |dtype|
+    context dtype do
+      err = 1e-4
+      it "should correctly invert a" do
+        a = NMatrix.new([2,2],[1,2,3,4], dtype: dtype)
+        b = NMatrix.new([2,2],[-2,1,1.5,-0.5], dtype: dtype)
+        begin
+          expect(a.pinv).to be_within(err).of(b)
+        rescue NotImplementedError
+          pending "Suppressing a NotImplementedError when the lapacke or atlas plugin is not available"
+        end
+        end
+      end
+  end
+
   # TODO: Get it working with ROBJ too
   [:byte,:int8,:int16,:int32,:int64,:float32,:float64].each do |left_dtype|
     [:byte,:int8,:int16,:int32,:int64,:float32,:float64].each do |right_dtype|