Norms (cond.)

lib/nmatrix/math.rb

@@ -803,6 +803,48 @@ def laswp(ary, opts={})
     self.clone.laswp!(ary, opts)
   end
 
+  #
+  # call-seq:
+  #   matrix_norm -> Numeric
+  #
+  #  Calculates the selected norm (defaults to 2-norm) of a 2D matrix.
+  #
+  #  This should be used for small or medium sized matrices.
+  #  For greater matrices, there should be a separate implementation where
+  #  the norm is estimated rather than computed, for the sake of computation speed.
+  #
+  #  Currently implemented norms are 1-norm, 2-norm, Frobenius, Infinity.
+  #  A minus on the 1, 2 and inf norms returns the minimum instead of the maximum value.
+  #
+  #  Tested mainly with dense matrices. Further checks and modifications might
+  #  be necessary for sparse matrices.
+  #
+  # * *Returns* :
+  # - The selected norm of the matrix.
+  # * *Raises* :
+  # - +NotImplementedError+ -> norm can be calculated only for 2D matrices
+  # - +ArgumentError+ -> unrecognized norm
+  #
+  def matrix_norm type = 2
+    raise(NotImplementedError, "norm can be calculated only for 2D matrices") unless self.dim == 2
+    raise(NotImplementedError, "norm only implemented for dense storage") unless self.stype == :dense
+    raise(ArgumentError, "norm not defined for byte dtype")if self.dtype == :byte
+    case type
+    when nil, 2, -2
+      return self.two_matrix_norm (type == -2)
+    when 1, -1
+      return self.one_matrix_norm (type == -1)
+    when :frobenius, :fro
+      return self.fro_matrix_norm
+    when :infinity, :inf, :'-inf', :'-infinity'
+      return self.inf_matrix_norm  (type == :'-inf' || type == :'-infinity')
+    else
+      raise ArgumentError.new("argument must be a valid integer or symbol")
+    end
+  end
+
+
+
   #
   # call-seq:
   #     det -> determinant
@@ -1205,6 +1247,55 @@ def nrm2 incx=1, n=nil
   end
   alias :norm2 :nrm2
 
+  # Norm calculation methods
+  # Frobenius norm: the Euclidean norm of the matrix, treated as if it were a vector
+  def fro_matrix_norm
+    #float64 has to be used in any case, since nrm2 will not yield correct result for float32
+    self_cast = self.cast(:dtype => :float64)
+
+    column_vector = self_cast.reshape([self.size, 1])
+
+    return column_vector.nrm2
+  end
+
+  # 2-norm: the largest/smallest singular value of the matrix
+  def two_matrix_norm minus = false
+
+    self_cast = self.cast(:dtype => :float64)
+
+    #TODO: confirm if this is the desired svd calculation
+    svd = self_cast.gesvd
+    return svd[1][0, 0] unless minus
+    return svd[1][svd[1].rows-1, svd[1].cols-1]
+  end
+
+  # 1-norm: the maximum/minimum absolute column sum of the matrix
+  def one_matrix_norm minus = false
+    #TODO: change traversing method for sparse matrices
+    number_of_columns = self.cols
+    col_sums = []
+
+    number_of_columns.times do |i|
+      col_sums << self.col(i).inject(0) { |sum, number| sum += number.abs}
+    end
+
+    return col_sums.max unless minus
+    return col_sums.min
+  end
+
+  # Infinity norm: the maximum/minimum absolute row sum of the matrix
+  def inf_matrix_norm minus = false
+    number_of_rows = self.rows
+    row_sums = []
+
+    number_of_rows.times do |i|
+      row_sums << self.row(i).inject(0) { |sum, number| sum += number.abs}
+    end
+
+    return row_sums.max unless minus
+    return row_sums.min
+  end
+
   #
   # call-seq:
   #     scale! -> NMatrix
@@ -1458,4 +1549,4 @@ def dtype_for_floor_or_ceil
       self.__dense_map__ { |l| l.send(op,rhs) }
     end
   end
-end
+end

spec/math_spec.rb

@@ -462,9 +462,9 @@
 
       it "should correctly invert a matrix in place (bang)" do
         pending("not yet implemented for :object dtype") if dtype == :object
-        a = NMatrix.new(:dense, 5, [1, 8,-9, 7, 5, 
-                                    0, 1, 0, 4, 4, 
-                                    0, 0, 1, 2, 5, 
+        a = NMatrix.new(:dense, 5, [1, 8,-9, 7, 5,
+                                    0, 1, 0, 4, 4,
+                                    0, 0, 1, 2, 5,
                                     0, 0, 0, 1,-5,
                                     0, 0, 0, 0, 1 ], dtype)
         b = NMatrix.new(:dense, 5, [1,-8, 9, 7, 17,
@@ -1148,7 +1148,7 @@
                                       3, 4, 5,
                                       6, 7, 8], stype: stype, dtype: dtype)
           end
-          
+
           it "scales the matrix by a given factor and return the result" do
             pending("not yet implemented for :object dtype") if dtype == :object
             if integer_dtype? dtype
@@ -1160,7 +1160,7 @@
                                                              12, 14, 16], stype: stype, dtype: dtype))
             end
           end
-          
+
           it "scales the matrix in place by a given factor" do
             pending("not yet implemented for :object dtype") if dtype == :object
             if dtype == :int8
@@ -1177,4 +1177,78 @@
       end
     end
   end
+  context "matrix_norm" do
+    ALL_DTYPES.each do |dtype|
+      context dtype do
+        pending("not yet implemented for :object dtype") if dtype == :object
+        before do
+          @n = NMatrix.new([3,3], [-4,-3,-2,
+                                   -1, 0, 1,
+                                    2, 3, 4], dtype: dtype)
+
+          @matrix_norm_TOLERANCE = 1.0e-10
+        end
+
+        it "should default to 2-matrix_norm" do
+          pending("not yet implemented for NMatrix-JRuby") if jruby?
+          if(dtype == :byte)
+            expect{@n.matrix_norm}.to raise_error(ArgumentError)
+          else
+            begin
+              expect(@n.matrix_norm).to be_within(@matrix_norm_TOLERANCE).of(7.348469228349535)
+
+              rescue NotImplementedError
+                pending "Suppressing a NotImplementedError when the lapacke plugin is not available"
+            end
+          end
+        end
+
+        it "should reject invalid arguments" do
+          pending("not yet implemented for NMatrix-JRuby") if jruby?
+
+          expect{@n.matrix_norm(0.5)}.to raise_error(ArgumentError)
+        end
+
+        it "should calculate 1 and 2(minus) matrix_norms correctly" do
+          pending("not yet implemented for NMatrix-JRuby") if jruby?
+          if(dtype == :byte)
+              expect{@n.matrix_norm(1)}.to raise_error(ArgumentError)
+              expect{@n.matrix_norm(-2)}.to raise_error(ArgumentError)
+              expect{@n.matrix_norm(-1)}.to raise_error(ArgumentError)
+          else
+            expect(@n.matrix_norm(1)).to eq(7)
+            begin
+
+              #FIXME: change to the correct value when overflow issue is resolved
+              #expect(@n.matrix_norm(-2)).to eq(1.8628605857884395e-07)
+              expect(@n.matrix_norm(-2)).to be_within(@matrix_norm_TOLERANCE).of(0.0)
+              rescue NotImplementedError
+                pending "Suppressing a NotImplementedError when the lapacke plugin is not available"
+            end
+            expect(@n.matrix_norm(-1)).to eq(6)
+          end
+        end
+
+        it "should calculate infinity matrix_norms correctly" do
+          pending("not yet implemented for NMatrix-JRuby") if jruby?
+          if(dtype == :byte)
+            expect{@n.matrix_norm(:inf)}.to raise_error(ArgumentError)
+            expect{@n.matrix_norm(:'-inf')}.to raise_error(ArgumentError)
+          else
+            expect(@n.matrix_norm(:inf)).to eq(9)
+            expect(@n.matrix_norm(:'-inf')).to eq(2)
+          end
+        end
+
+        it "should calculate frobenius matrix_norms correctly" do
+          pending("not yet implemented for NMatrix-JRuby") if jruby?
+          if(dtype == :byte)
+            expect{@n.matrix_norm(:fro)}.to raise_error(ArgumentError)
+          else
+            expect(@n.matrix_norm(:fro)).to be_within(@matrix_norm_TOLERANCE).of(7.745966692414834)
+          end
+        end
+      end
+    end
+  end
 end