From 93553febc0c922f44792fc6e146c645f08ad9ded Mon Sep 17 00:00:00 2001 From: Arafatk Date: Wed, 23 Mar 2016 13:37:19 +0530 Subject: [PATCH 1/2] Adding pseudo inverse Squashed commit of the following: commit a7730047ca8db153ccc8edb397dc05b6fa8a12f2 Author: Arafatk Date: Wed Mar 23 13:33:47 2016 +0530 removing loop commit 9c04728788b3a50a9ce082e7304d15f072c7fd4e Author: Arafatk Date: Fri Mar 11 10:32:24 2016 +0530 Adding Moore pseudo inverse Squashed commit of the following: commit 2d335f84de81032811b840cbabb1cb58adbbc294 Author: Arafatk Date: Fri Mar 11 10:31:07 2016 +0530 Improving tests commit f17061f1d1d6fb0aafe835d823e7c578d1a2424e Author: Arafatk Date: Sat Mar 5 12:25:59 2016 +0530 Removing un-necessary spaces commit fa908ac3347f399e2ebf1f4d0619458ff1333e23 Author: Arafatk Date: Fri Mar 4 10:33:23 2016 +0530 Corrections to spec tests commit 9d9ab352a192089c8880f8d11889f30f537a2d7a Author: Arafatk Date: Tue Mar 1 20:40:33 2016 +0530 Adding Moore-Penrose pseudo-inverse of a matrix Squashed commit of the following: commit 8c02601d5cd065071c4378c16a5a525477033b87 Author: Arafatk Date: Tue Mar 1 20:39:17 2016 +0530 correcting Docs commit b8b7f1f68962251853eb2fea1414212182e7fe77 Author: Arafatk Date: Mon Feb 29 23:55:14 2016 +0530 Squashed commit of the following: commit fd2cba085e44d5434ca324302d3e5ad4b7ce4426 Author: Arafatk Date: Mon Feb 29 23:52:33 2016 +0530 correcting docs commit 87eeaaed679e15ba0166b7b335bd874f31000f07 Author: Arafatk Date: Mon Feb 29 23:26:17 2016 +0530 Adding complex data type commit f3950b11b8189903f1d70c8012f1bd15e153cca4 Author: Arafatk Date: Sun Feb 28 09:39:05 2016 +0530 Adding Moore-Penrose pseudo-inverse of a matrix Squashed commit of the following: commit dbba0c17dc0b91c91984e19d2845819bdfeae0c9 Author: Arafatk Date: Sat Feb 27 23:20:31 2016 +0530 Adding Pseudo inverse commit 32a49e675c6db03787b7de4fc0ab321b08d5b4b5 Author: Arafatk Date: Sat Feb 27 22:24:19 2016 +0530 correcting error in docs commit 8031cc83764a310991ea55b12f2889b780bce529 Author: Arafatk Date: Sat Feb 27 21:06:55 2016 +0530 Correcting load error commit 4b0b37596552277cd6417fe9c824bc2be9b4f3ca Author: Arafatk Date: Sat Feb 27 18:59:25 2016 +0530 Correcting Load error in Spec Tests commit 2fffafca46fdf965f3000b299e88751e8c7b627a Author: Arafatk Date: Sat Feb 27 17:22:08 2016 +0530 Improving Docs Adding Tests commit 6e0ab741ad4f4f69eb02e07192a4ceb93fb09a73 Author: Arafatk Date: Wed Feb 24 22:28:21 2016 +0530 Adding Moore-Penrose pseudo-inverse of a matrix --- lib/nmatrix/math.rb | 65 +++++++++++++++++++++++++++++++++++++++++++++ spec/math_spec.rb | 57 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 122 insertions(+) diff --git a/lib/nmatrix/math.rb b/lib/nmatrix/math.rb index 05c9c981..5c667ca6 100644 --- a/lib/nmatrix/math.rb +++ b/lib/nmatrix/math.rb @@ -112,6 +112,71 @@ def invert end alias :inverse :invert + + # + # call-seq: + # pinv -> NMatrix + # + # Compute the Moore-Penrose pseudo-inverse of a matrix using its + # singular value decomposition (SVD). + # + # This function requires the nmatrix-atlas gem installed. + # + # * *Arguments* : + # - +tolerance(optional)+ -> Cutoff for small singular values. + # + # * *Returns* : + # - Pseudo-inverse matrix. + # + # * *Raises* : + # - +NotImplementedError+ -> If called without nmatrix-atlas or nmatrix-lapacke gem. + # - +TypeError+ -> If called without float or complex data type. + # + # * *Examples* : + # + # 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(tolerance = 1e-15) + raise DataTypeError, "pinv works only with matrices of float or complex data type" unless + [:float32, :float64, :complex64, :complex128].include?(dtype) + if self.complex_dtype? + u, s, vt = self.complex_conjugate.gesvd # singular value decomposition + else + u, s, vt = self.gesvd + end + 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.max.to_f + cutoff = tolerance * 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.dot(NMatrix.diagonal(s.to_a)).transpose + vt_reduced.dot(multiplier) + end + alias :pseudo_inverse :pinv + alias :pseudoinverse :pinv + + # # call-seq: # getrf! -> Array diff --git a/spec/math_spec.rb b/spec/math_spec.rb index 127c7b13..290f2459 100644 --- a/spec/math_spec.rb +++ b/spec/math_spec.rb @@ -332,6 +332,63 @@ end end + NON_INTEGER_DTYPES.each do |dtype| + next if dtype == :object + context dtype do + err = Complex(1e-3, 1e-3) + it "should correctly invert a 2x2 matrix" do + if dtype == :complex64 || dtype == :complex128 + a = NMatrix.new([2, 2], [Complex(16, 81), Complex(91, 51), \ + Complex(13, 54), Complex(71, 24)], dtype: dtype) + b = NMatrix.identity(2, dtype: dtype) + + begin + expect(a.dot(a.pinv)).to be_within(err).of(b) + rescue NotImplementedError + pending "Suppressing a NotImplementedError when the atlas plugin is not available" + end + + else + a = NMatrix.new([2, 2], [141, 612, 9123, 654], dtype: dtype) + b = NMatrix.identity(2, dtype: dtype) + + begin + expect(a.dot(a.pinv)).to be_within(err).of(b) + rescue NotImplementedError + pending "Suppressing a NotImplementedError when the atlas plugin is not available" + end + end + end + + it "should verify a.dot(b.dot(a)) == a and b.dot(a.dot(b)) == b" do + if dtype == :complex64 || dtype == :complex128 + a = NMatrix.new([3, 2], [Complex(94, 11), Complex(87, 51), Complex(82, 39), \ + Complex(45, 16), Complex(25, 32), Complex(91, 43) ], dtype: dtype) + + begin + b = a.pinv # pseudo inverse + expect(a.dot(b.dot(a))).to be_within(err).of(a) + expect(b.dot(a.dot(b))).to be_within(err).of(b) + rescue NotImplementedError + pending "Suppressing a NotImplementedError when the atlas plugin is not available" + end + + else + a = NMatrix.new([3, 3], [9, 4, 52, 12, 52, 1, 3, 55, 6], dtype: dtype) + + begin + b = a.pinv # pseudo inverse + expect(a.dot(b.dot(a))).to be_within(err).of(a) + expect(b.dot(a.dot(b))).to be_within(err).of(b) + rescue NotImplementedError + pending "Suppressing a NotImplementedError when the atlas plugin is not available" + end + 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| From ba0f2b44dd57f053249db84fd886b6f346dbb871 Mon Sep 17 00:00:00 2001 From: John Woods Date: Sat, 11 Mar 2017 13:04:56 -0600 Subject: [PATCH 2/2] Added pendings for jruby. Also made sure pendings were marked for object dtypes (instead of just using `next`) --- spec/math_spec.rb | 25 ++++++++++++++++++++----- 1 file changed, 20 insertions(+), 5 deletions(-) diff --git a/spec/math_spec.rb b/spec/math_spec.rb index 290f2459..6c229d3f 100644 --- a/spec/math_spec.rb +++ b/spec/math_spec.rb @@ -203,7 +203,6 @@ end NON_INTEGER_DTYPES.each do |dtype| - next if dtype == :object context dtype do before do @m = NMatrix.new([3,4], GETRF_EXAMPLE_ARRAY, dtype: dtype) @@ -217,11 +216,13 @@ #haven't check this spec yet. Also it doesn't check all the elements of the matrix. it "should correctly factorize a matrix" do + pending("not yet implemented for :object dtype") if dtype == :object a = @m.factorize_lu expect(a).to be_within(@err).of(NMatrix.new([3,4], GETRF_SOLUTION_ARRAY, dtype: dtype)) end it "also returns the permutation matrix" do + pending("not yet implemented for :object dtype") if dtype == :object a, p = @m.factorize_lu perm_matrix: true expect(a).to be_within(@err).of(NMatrix.new([3,4], GETRF_SOLUTION_ARRAY, dtype: dtype)) @@ -244,6 +245,7 @@ # 1 2 6 # which is symmetric and positive-definite as required, but # we need only store the lower-half of the matrix. + pending("not yet implemented for :object dtype") if dtype == :object a = NMatrix.new([3,3],[1,0,0, 1,2,0, 1,2,6], dtype: dtype) begin r = a.potrf!(:lower) @@ -257,6 +259,7 @@ end it "calculates cholesky decomposition using potrf (upper)" do + pending("not yet implemented for :object dtype") if dtype == :object a = NMatrix.new([3,3],[1,1,1, 0,2,2, 0,0,6], dtype: dtype) begin r = a.potrf!(:upper) @@ -270,6 +273,7 @@ end it "calculates cholesky decomposition using #factorize_cholesky" do + pending("not yet implemented for :object dtype") if dtype == :object a = NMatrix.new([3,3],[1,2,1, 2,13,5, 1,5,6], dtype: dtype) begin u,l = a.factorize_cholesky @@ -287,7 +291,6 @@ ALL_DTYPES.each do |dtype| next if dtype == :byte #doesn't work for unsigned types - next if dtype == :object context dtype do err = case dtype @@ -298,6 +301,7 @@ end 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, @@ -319,6 +323,7 @@ end it "should correctly invert a matrix out-of-place" do + pending("not yet implemented for :object dtype") if dtype == :object a = NMatrix.new(:dense, 3, [1,2,3,0,1,4,5,6,0], dtype) if a.integer_dtype? @@ -333,10 +338,11 @@ end NON_INTEGER_DTYPES.each do |dtype| - next if dtype == :object context dtype do err = Complex(1e-3, 1e-3) it "should correctly invert a 2x2 matrix" do + pending("not yet implemented for NMatrix-JRuby") if jruby? + pending("not yet implemented for :object dtype") if dtype == :object if dtype == :complex64 || dtype == :complex128 a = NMatrix.new([2, 2], [Complex(16, 81), Complex(91, 51), \ Complex(13, 54), Complex(71, 24)], dtype: dtype) @@ -361,6 +367,8 @@ end it "should verify a.dot(b.dot(a)) == a and b.dot(a.dot(b)) == b" do + pending("not yet implemented for NMatrix-JRuby") if jruby? + pending("not yet implemented for :object dtype") if dtype == :object if dtype == :complex64 || dtype == :complex128 a = NMatrix.new([3, 2], [Complex(94, 11), Complex(87, 51), Complex(82, 39), \ Complex(45, 16), Complex(25, 32), Complex(91, 43) ], dtype: dtype) @@ -609,9 +617,9 @@ context "#solve" do NON_INTEGER_DTYPES.each do |dtype| - next if dtype == :object # LU factorization doesnt work for :object yet it "solves linear equation for dtype #{dtype}" do + pending("not yet implemented for :object dtype") if dtype == :object a = NMatrix.new [2,2], [3,1,1,2], dtype: dtype b = NMatrix.new [2,1], [9,8], dtype: dtype @@ -619,6 +627,7 @@ end it "solves linear equation for #{dtype} (non-symmetric matrix)" do + pending("not yet implemented for :object dtype") if dtype == :object a = NMatrix.new [3,3], [1,1,1, -1,0,1, 3,4,6], dtype: dtype b = NMatrix.new [3,1], [6,2,29], dtype: dtype @@ -633,6 +642,7 @@ end it "solves linear equation for dtype #{dtype} (non-vector rhs)" do + pending("not yet implemented for :object dtype") if dtype == :object a = NMatrix.new [3,3], [1,0,0, -1,0,1, 2,1,1], dtype: dtype b = NMatrix.new [3,2], [1,0, 1,2, 4,2], dtype: dtype @@ -814,8 +824,8 @@ context "determinants" do ALL_DTYPES.each do |dtype| - next if dtype == :object context dtype do + pending("not yet implemented for :object dtype") if dtype == :object before do @a = NMatrix.new([2,2], [1,2, 3,4], dtype: dtype) @@ -836,15 +846,19 @@ end end it "computes the determinant of 2x2 matrix" do + pending("not yet implemented for :object dtype") if dtype == :object expect(@a.det).to be_within(@err).of(-2) end it "computes the determinant of 3x3 matrix" do + pending("not yet implemented for :object dtype") if dtype == :object expect(@b.det).to be_within(@err).of(-8) end it "computes the determinant of 4x4 matrix" do + pending("not yet implemented for :object dtype") if dtype == :object expect(@c.det).to be_within(@err).of(-18) end it "computes the exact determinant of 2x2 matrix" do + pending("not yet implemented for :object dtype") if dtype == :object if dtype == :byte expect{@a.det_exact}.to raise_error(DataTypeError) else @@ -852,6 +866,7 @@ end end it "computes the exact determinant of 3x3 matrix" do + pending("not yet implemented for :object dtype") if dtype == :objectx if dtype == :byte expect{@a.det_exact}.to raise_error(DataTypeError) else