From 93553febc0c922f44792fc6e146c645f08ad9ded Mon Sep 17 00:00:00 2001 From: Arafatk Date: Wed, 23 Mar 2016 13:37:19 +0530 Subject: [PATCH] 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|