diff --git a/lib/nmatrix/math.rb b/lib/nmatrix/math.rb index caab679d..9fa8d97d 100644 --- a/lib/nmatrix/math.rb +++ b/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 diff --git a/spec/math_spec.rb b/spec/math_spec.rb index dd74ce44..b344332b 100644 --- a/spec/math_spec.rb +++ b/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|