Matrix norms implementation

lib/nmatrix/math.rb

@@ -282,7 +282,7 @@ def laswp(ary)
   # * *Returns* :
   #   - The determinant of the matrix. It's the same type as the matrix's dtype.
   # * *Raises* :
-  #   - +NotImplementedError+ -> Must be used in 2D matrices.
+  #   - +NotImplementedError+ -> determinant can be calculated only for 2D matrices
   #
   def det
     raise(NotImplementedError, "determinant can be calculated only for 2D matrices") unless self.dim == 2
@@ -506,6 +506,47 @@ def abs
       self.__yale_map_stored__ { |v| v.abs }
     end.cast(self.stype, abs_dtype)
   end
+
+
+  #
+  # call-seq:
+  #   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 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 
+
+    case type
+    when nil, 2, -2 
+      return self.two_norm (type == -2)
+    when 1, -1 
+      return self.one_norm (type == -1)
+    when :frobenius, :fro 
+      return self.fro_norm
+    when :infinity, :inf, :'-inf', :'-infinity' 
+      return self.inf_norm  (type == :'-inf' || type == :'-infinity')
+    else
+      raise ArgumentError.new("argument must be a valid integer or symbol")
+    end
+  end
 
 
   #
@@ -669,4 +710,53 @@ def __dense_unary_log__(base)
       self.__dense_map__ { |l| l.send(op,rhs) }
     end
   end
+
+  # Norm calculation methods
+  # Frobenius norm: the Euclidean norm of the matrix, treated as if it were a vector
+  def fro_norm   
+    #float64 has to be used in any case, since nrm2 will not yield correct result for float32
+    self_cast = self.cast(:dtype => :float64) unless self.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_norm minus = false
+    #self.dtype == :int32 ? self_cast = self.cast(:dtype => :float32) : self_cast = self.cast(:dtype => :float64)
+    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_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.sort!.last unless minus
+    return col_sums.sort!.first
+  end
+
+  # Infinity norm: the maximum/minimum absolute row sum of the matrix  
+  def inf_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.sort!.last unless minus
+    return row_sums.sort!.first    
+  end
 end

spec/math_spec.rb

@@ -162,6 +162,40 @@
       end
     end
   end
+
+  #TODO: check validity of SVD values in 2-norm
+  context "#norm" do
+    it "should default to 2-norm" do
+      n = NMatrix.new([3, 3], [-4, -3, -2, -1, 0, 1, 2, 3, 4])
+      expect(n.norm).to eq(7.348469228349535)
+    end
+
+    it "should reject invalid arguments" do
+      n = NMatrix.new([3, 3], [-4, -3, -2, -1, 0, 1, 2, 3, 4])
+
+      expect{n.norm(0.5)}.to raise_error(ArgumentError)
+    end
+
+    it "should calculate 1 and 2 norms correctly" do
+      n = NMatrix.new([3, 3], [-4, -3, -2, -1, 0, 1, 2, 3, 4])
+      expect(n.norm(1)).to eq(7)
+      #FIXME: change to the correct value when overflow issue is resolved
+      #expect(n.norm(-2)).to eq(1.8628605857884395e-07)
+      expect(n.norm(-2)).to eq(0.0)
+      expect(n.norm(-1)).to eq(6)
+    end
+
+    it "should calculate infinity norms correctly" do
+      n = NMatrix.new([3, 3], [-4, -3, -2, -1, 0, 1, 2, 3, 4])
+      expect(n.norm(:inf)).to eq(9)
+      expect(n.norm(:'-inf')).to eq(2)
+    end
+
+    it "should calculate frobenius norms correctly" do
+      n = NMatrix.new([3, 3], [-4, -3, -2, -1, 0, 1, 2, 3, 4])
+      expect(n.norm(:fro)).to eq(7.745966692414834)
+    end
+  end
 
   # TODO: Get it working with ROBJ too
   [:byte,:int8,:int16,:int32,:int64,:float32,:float64,:rational64,:rational128].each do |left_dtype|