added function svd_rank() to calculate rank via svd. Added tests for the rank in math_spec.rb

Author: pgtgrly

ext/nmatrix/nmatrix.cpp

@@ -32,6 +32,7 @@
  */
 
 #include <ruby.h>
+#include <cfloat>
 #include <algorithm> // std::min
 #include <fstream>
 

ext/nmatrix/ruby_nmatrix.c

@@ -364,6 +364,12 @@ void Init_nmatrix() {
   rb_define_alias(cNMatrix, "effective_dim", "effective_dimensions");
   rb_define_alias(cNMatrix, "equal?", "eql?");
 
+  ////////////
+  //Epsilons//
+  ////////////
+  rb_define_const(cNMatrix, "FLOAT64_EPSILON", rb_const_get(rb_cFloat, rb_intern("EPSILON")));
+  rb_define_const(cNMatrix, "FLOAT32_EPSILON", DBL2NUM(FLT_EPSILON));
+
   ///////////////////////
   // Symbol Generation //
   ///////////////////////

lib/nmatrix/math.rb

@@ -698,6 +698,38 @@ def positive_definite?
     true
   end
 
+  #
+  # call-seq:
+  #   svd_rank() -> int
+  #   svd_rank(tolerence) ->int
+  # Gives rank of the matrix based on the singular value decomposition.
+  # The rank of a matrix  is computed as the number of diagonal elements in Sigma that are larger than a tolerance
+  #
+  #* *Returns* :
+  # - An integer equal to the rank of the matrix
+  #* *Raises* :
+  #  - +ShapeError+ -> Is only computable on 2-D matrices
+  #
+  def svd_rank(tolerence="default")
+    raise(ShapeError, "rank calculated only for 2-D matrices") unless
+      self.dim == 2 
+
+    sigmas = self.gesvd[1].to_a.flatten
+    eps = NMatrix::FLOAT64_EPSILON
+
+    # epsilon depends on the width of the number
+    if (self.dtype == :float32 || self.dtype == :complex64) 
+      eps = NMatrix::FLOAT32_EPSILON
+    end
+    case tolerence
+      when "default"
+        tolerence = self.shape.max * sigmas.max * eps # tolerence of a Matrix A is max(size(A))*eps(norm(A)). norm(A) is nearly equal to max(sigma of A)
+    end
+    return sigmas.map { |x| x > tolerence ? 1 : 0 }.reduce(:+)
+  end
+
+
+
 protected
   # Define the element-wise operations for lists. Note that the __list_map_merged_stored__ iterator returns a Ruby Object
   # matrix, which we then cast back to the appropriate type. If you don't want that, you can redefine these functions in

spec/math_spec.rb

@@ -1276,4 +1276,88 @@
         expect(n.positive_definite?).to be_truthy
       end
   end
+  
+  context "#svd_rank" do 
+    FLOAT_DTYPES.each do |dtype|
+      context dtype do
+        #examples from https://www.cliffsnotes.com/study-guides/algebra/linear-algebra/real-euclidean-vector-spaces/the-rank-of-a-matrix
+        it "calculates the rank of matrix using singular value decomposition with NMatrix on rectangular matrix without tolerence" do
+          pending("not yet implemented for NMatrix-JRuby") if jruby?
+          a = NMatrix.new([4,3],[2,-1,3, 1,0,1, 0,2,-1, 1,1,4], dtype: dtype)
+
+          begin
+            rank = a.svd_rank()
+          
+            rank_true = 3
+            expect(rank).to eq (rank_true)
+
+          rescue NotImplementedError
+            pending "Suppressing a NotImplementedError when the lapacke plugin is not available" 
+          end         
+        end
+
+        it "calculates the rank of matrix using singular value decomposition with NMatrix on rectangular matrix with tolerence" do
+        
+          a = NMatrix.new([4,3],[2,-1,3, 1,0,1, 0,2,-1, 1,1,4], dtype: dtype)
+          pending("not yet implemented for NMatrix-JRuby") if jruby?
+          begin
+            rank = a.svd_rank(4)
+  
+            rank_true = 1
+            expect(rank).to eq (rank_true)
+
+          rescue NotImplementedError
+             pending "Suppressing a NotImplementedError when the lapacke plugin is not available" 
+          end
+        end
+
+        it "calculates the rank of matrix using singular value decomposition with NMatrix on square matrix without tolerence" do
+        
+          a = NMatrix.new([4,4],[1,-1,1,-1, -1,1,-1,1, 1,-1,1,-1, -1,1,-1,1], dtype: dtype)
+          pending("not yet implemented for NMatrix-JRuby") if jruby?
+          begin
+            rank = a.svd_rank()
+          
+            rank_true = 1
+            expect(rank).to eq (rank_true)
+
+          rescue NotImplementedError
+             pending "Suppressing a NotImplementedError when the lapacke plugin is not available" 
+          end
+        end
+
+        it "calculates the rank of matrix using singular value decomposition with NMatrix on square matrix with very small tolerence(for float32)" do
+          pending("not yet implemented for NMatrix-JRuby") if jruby?
+          a = NMatrix.new([4,4],[1,-1,1,-1, -1,1,-1,1, 1,-1,1,-1, -1,1,-1,1], dtype: :float32)
+
+          begin
+            rank = a.svd_rank(1.7881389169360773e-08)
+          
+            rank_true = 2
+            expect(rank).to eq (rank_true)
+
+          rescue NotImplementedError
+             pending "Suppressing a NotImplementedError when the lapacke plugin is not available" 
+          end
+        end
+
+        it "calculates the rank of matrix using singular value decomposition with NMatrix on square matrix with very small tolerence(for float64)" do
+          pending("not yet implemented for NMatrix-JRuby") if jruby?
+          a = NMatrix.new([4,4],[1,-1,1,-1, -1,1,-1,1, 1,-1,1,-1, -1,1,-1,1], dtype: :float64)
+
+          begin
+            rank = a.svd_rank(1.7881389169360773e-08)
+          
+            rank_true = 1
+            expect(rank).to eq (rank_true)
+
+          rescue NotImplementedError
+             pending "Suppressing a NotImplementedError when the lapacke plugin is not available" 
+          end
+        end
+
+      end
+    end 
+  end 
+
 end