@@ -112,6 +112,71 @@ def invert
112112 end
113113 alias :inverse :invert
114114
115+
116+ #
117+ # call-seq:
118+ # pinv -> NMatrix
119+ #
120+ # Compute the Moore-Penrose pseudo-inverse of a matrix using its
121+ # singular value decomposition (SVD).
122+ #
123+ # This function requires the nmatrix-atlas gem installed.
124+ #
125+ # * *Arguments* :
126+ # - +tolerance(optional)+ -> Cutoff for small singular values.
127+ #
128+ # * *Returns* :
129+ # - Pseudo-inverse matrix.
130+ #
131+ # * *Raises* :
132+ # - +NotImplementedError+ -> If called without nmatrix-atlas or nmatrix-lapacke gem.
133+ # - +TypeError+ -> If called without float or complex data type.
134+ #
135+ # * *Examples* :
136+ #
137+ # a = NMatrix.new([2,2],[1,2,
138+ # 3,4], dtype: :float64)
139+ # a.pinv # => [ [-2.0000000000000018, 1.0000000000000007]
140+ # [1.5000000000000016, -0.5000000000000008] ]
141+ #
142+ # b = NMatrix.new([4,1],[1,2,3,4], dtype: :float64)
143+ # b.pinv # => [ [ 0.03333333, 0.06666667, 0.99999999, 0.13333333] ]
144+ #
145+ # == References
146+ #
147+ # * https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_pseudoinverse
148+ # * G. Strang, Linear Algebra and Its Applications, 2nd Ed., Orlando, FL, Academic Press
149+ #
150+ def pinv ( tolerance = 1e-15 )
151+ raise DataTypeError , "pinv works only with matrices of float or complex data type" unless
152+ [ :float32 , :float64 , :complex64 , :complex128 ] . include? ( dtype )
153+ if self . complex_dtype?
154+ u , s , vt = self . complex_conjugate . gesvd # singular value decomposition
155+ else
156+ u , s , vt = self . gesvd
157+ end
158+ rows = self . shape [ 0 ]
159+ cols = self . shape [ 1 ]
160+ if rows < cols
161+ u_reduced = u
162+ vt_reduced = vt [ 0 ..rows - 1 , 0 ..cols - 1 ] . transpose
163+ else
164+ u_reduced = u [ 0 ..rows - 1 , 0 ..cols - 1 ]
165+ vt_reduced = vt . transpose
166+ end
167+ largest_singular_value = s . max . to_f
168+ cutoff = tolerance * largest_singular_value
169+ ( 0 ...[ rows , cols ] . min ) . each do |i |
170+ s [ i ] = 1 / s [ i ] if s [ i ] > cutoff
171+ s [ i ] = 0 if s [ i ] <= cutoff
172+ end
173+ multiplier = u_reduced . dot ( NMatrix . diagonal ( s . to_a ) ) . transpose
174+ vt_reduced . dot ( multiplier )
175+ end
176+ alias :pseudo_inverse :pinv
177+ alias :pseudoinverse :pinv
178+
179+
115180 #
116181 # call-seq:
117182 # adjugate! -> NMatrix
0 commit comments