diff --git a/clifford/_layout.py b/clifford/_layout.py index 4ccb8cd2..05e535c9 100644 --- a/clifford/_layout.py +++ b/clifford/_layout.py @@ -690,6 +690,29 @@ def leftLaInvJIT(value): return leftLaInvJIT + @_cached_property + def pinv_func(self): + """ + Get a function that returns left-pseudoinverse using a computational + linear algebra method proposed by Christian Perwass. Identical to + `inv_func` except the inverse has been replaced by a pseudo inverse. + """ + mult_table = self.gmt + k_list, l_list, m_list = mult_table.coords + mult_table_vals = mult_table.data + n_dims = mult_table.shape[1] + + identity = np.zeros((n_dims,)) + identity[self._basis_blade_order.bitmap_to_index[0]] = 1 + + @_numba_utils.njit + def leftLaPInvJIT(value): + intermed = _numba_val_get_left_gmt_matrix(value, k_list, l_list, m_list, mult_table_vals, n_dims) + sol = np.linalg.pinv(intermed) @ identity + return sol + + return leftLaPInvJIT + def get_left_gmt_matrix(self, x): """ This produces the matrix X that performs left multiplication with x diff --git a/clifford/_multivector.py b/clifford/_multivector.py index a394a9a4..4aaef1bf 100644 --- a/clifford/_multivector.py +++ b/clifford/_multivector.py @@ -773,6 +773,13 @@ def inv(self) -> 'MultiVector': leftInv = leftLaInv rightInv = leftLaInv + def pinv(self) -> 'MultiVector': + """Returns the pseudoinverse using a computational + linear algebra method proposed by Christian Perwass. Identical to + `inv` when an inverse exists. + """ + return self._newMV(self.layout.pinv_func(self.value)) + def dual(self, I=None) -> 'MultiVector': r"""The dual of the multivector against the given subspace I, :math:`\tilde M = MI^{-1}` diff --git a/clifford/test/test_clifford.py b/clifford/test/test_clifford.py index 8690f7a4..89c5ce4e 100644 --- a/clifford/test/test_clifford.py +++ b/clifford/test/test_clifford.py @@ -184,6 +184,17 @@ def test_indexing(self, g3): e12[1 + e12] assert e12[(2, 1)] == -1 + def test_pinv(self, g3): + layout, blades = g3, g3.blades + e2 = blades['e2'] + e3 = blades['e3'] + x = 0.5*(1-e3) # Singular + y = 2*e2 + 5*e3 # invertible + with pytest.raises(ValueError): + x.leftLaInv() + assert x == x.pinv() + assert y.pinv() == y.leftLaInv() + def test_add_float64(self, g3): ''' test array_wrap method to take control addition from numpy array