Added support for vandermonde derivatives

polyrat/arnoldi.py

@@ -128,6 +128,31 @@ def vandermonde_arnoldi_eval(X, R, indices, mode, weight = None):
 
 	return W
 
+
+def vandermonde_arnoldi_eval_der(X, R, indices, mode, weight = None, V = None):
+	if V is None:
+		V = vandermonde_arnoldi_eval(X, R, indices, mode, weight = weight)
+
+	M = X.shape[0]
+	N = R.shape[1]
+	n = X.shape[1]
+	DV = np.zeros((M, N, n), dtype = (R[0,0] * X[0,0]).dtype)
+
+	for ell in range(n):
+		index_iterator = enumerate(indices)
+		next(index_iterator)
+		for k, ids in index_iterator:
+			i, j = _update_rule(indices[:k], ids)
+			# Q[:,k] = X[:,i] * Q[:,j] - sum_s Q[:,s] * R[s, k]
+			if i == ell:
+				DV[:,k,ell] = V[:,j] + X[:,i] * DV[:,j,ell] - DV[:,0:k,ell] @ R[0:k,k] 
+			else:
+				DV[:,k,ell] = X[:,i] * DV[:,j,ell] - DV[:,0:k,ell] @ R[0:k,k] 
+			DV[:,k,ell] /= R[k,k]	
+
+	return DV
+
+
 class ArnoldiPolynomialBasis(PolynomialBasis):
 	r""" A polynomial basis constructed using Vandermonde with Arnoldi
 
@@ -158,6 +183,12 @@ def basis(self):
 	def vandermonde(self, X, weight = None):
 		return vandermonde_arnoldi_eval(X, self._R, self._indices, self.mode, weight = weight)
 
+	def vandermonde_derivative(self, X, weight = None):
+		if np.array_equal(X, self.X):
+			return vandermonde_arnoldi_eval_der(X, self._R, self._indices, self.mode, weight = weight, V = self.Q)
+		else:
+			return vandermonde_arnoldi_eval_der(X, self._R, self._indices, self.mode, weight = weight)
+
 	def roots(self, coef, *args, **kwargs):
 		from .basis import LegendrePolynomialBasis
 		from .polynomial import PolynomialApproximation

polyrat/basis.py

@@ -14,7 +14,34 @@
 from .index import *
 
 class PolynomialBasis:
+
+	def basis(self):
+		r""" Alias for vandermonde(X) where X is the points where the basis was initialized
+
+		Returns
+		-------
+		V: np.array (M, N)
+			generalized Vandermonde matrix evaluated using the points the basis was initialized with
+		"""
+		raise NotImplementedError
+
+
 	def vandermonde(self, X):
+		r""" Construct the generalized Vandermonde matrix associated with the polynomial basis 
+
+		Parameters
+		----------
+		X: np.array (M, dim)
+			Coordinates at which to evaluate the basis at 
+
+		Returns
+		-------
+		V: np.array (M, N)
+			generalized Vandermonde matrix 
+		"""
+		raise NotImplementedError
+
+	def vandermonde_derivative(self, X):
 		raise NotImplementedError
 
 
@@ -53,8 +80,6 @@ def _inv_scale(self, X):
 		return X*(self._ub[None,:] - self._lb[None,:])/2.0 + (self._ub[None,:] + self._lb[None,:])/2.0
 
 	def vandermonde(self, X):
-		r""" Construct the Vandermonde matrix
-		"""
 		X = self._scale(X)
 
 		if self.mode == 'total':		 
@@ -68,15 +93,73 @@ def vandermonde(self, X):
 			for k in range(self.dim):
 				V[:,j] *= V_coordinate[k][:,alpha[k]]
 		return V
-
+
+	def vandermonde_derivative(self, X):
+		r""" Construct a column-wise derivative of the generalized Vandermonde matrix
+		"""
+		M, dim = X.shape
+		N = len(self._indices)
+		DV = np.ones((M, N, dim), dtype = X.dtype)
+
+		Y = self._scale(X) 
+
+		if self.mode == 'total':		 
+			V_coordinate = [self._vander(Y[:,k], self.degree) for k in range(self.dim)]
+		elif self.mode == 'max':
+			V_coordinate = [self._vander(Y[:,k], d) for k,d in enumerate(self.degree)]
+
+		for k in range(dim):
+			for j, alpha in enumerate(self._indices):
+				for q in range(self.dim):
+					if q == k:
+						if self.mode == 'total':
+							DV[:,j,k] *= V_coordinate[q][:,0:-1] @ self._Dmat[alpha[q],:] 
+						elif self.mode == 'max':
+							DV[:,j,k] *= V_coordinate[q][:,0:-1] @ self._Dmat[alpha[q],:self.degree[q]]
+					else:
+						DV[:,j,k] *= V_coordinate[q][:,alpha[q]]
+
+			DV[:,:,k] *= self._dscale[k] 			
+
+		return DV	
+
+
 
 	@lru_cache(maxsize = 1)
 	def basis(self):
 		r""" The basis for the input coordinates
 		""" 
 		return self.vandermonde(self.X)
+
+	@property
+	@lru_cache(maxsize = 1)
+	def _Dmat(self):
+		r""" The matrix specifying the action of the derivative operator in this basis
+		"""
+		if self.mode == 'total':
+			max_degree = self.degree
+		elif self.mode == 'max':
+			max_degree = max(self.degree)
+		Dmat = np.zeros( (max_degree+1, max_degree))
+		I = np.eye(max_degree + 1)
+		for j in range(max_degree + 1):
+			Dmat[j,:] = self._der(I[:,j])
+		return Dmat
 
 
+	@property
+	@lru_cache(maxsize = 1)
+	def _dscale(self):
+		r""" Derivative of the scaling applied to the coordinates
+		"""
+		# As we assume the transformation is linear, we simply compute the finite difference
+		# with a unit step size
+		XX = np.zeros((2, self.dim))
+		XX[1,:] = 1
+		sXX = self._scale(XX)
+		dscale = sXX[1] - sXX[0]
+		return dscale
+
 class MonomialPolynomialBasis(TensorProductPolynomialBasis):
 	def _vander(self, *args, **kwargs):
 		return polyvander(*args, **kwargs)

polyrat/version.py

@@ -1,2 +1,2 @@
-__version__ = '0.1.1'
+__version__ = '0.1.2'
 

tests/test_basis.py

@@ -180,12 +180,62 @@ def test_scale(Basis, dim):
 	assert err < 1e-10, "Scaling/inverse scaling does not map to identity"
 
 
+
+@pytest.mark.parametrize("Basis",
+	[MonomialPolynomialBasis,
+	 LegendrePolynomialBasis,
+	 ChebyshevPolynomialBasis,
+	 HermitePolynomialBasis, 
+	 LaguerrePolynomialBasis,
+	ArnoldiPolynomialBasis,
+	]
+)
+@pytest.mark.parametrize("dim", [1,2,3])
+@pytest.mark.parametrize("degree",
+	[3,
+	(3,2),
+	(2,3),
+	(3,0,2),
+	]
+)
+def test_vandermonde_derivative(Basis, dim, degree):
+
+	# Exit without error if testing a total degree problem
+	try:
+		degree = int(degree)
+	except (TypeError, ValueError):
+		if len(degree) != dim: return
+
+	np.random.seed(0)
+	X = np.random.randn(1000, dim) 
+	#X = np.linspace(-2,2, 100).reshape(-1,1)
+	basis = Basis(X, degree)
+
+	X = np.random.randn(100, dim)
+	#X = np.linspace(-1, 1, 11).reshape(-1,1)
+	V = basis.vandermonde(X)
+	DV = basis.vandermonde_derivative(X)
+	for k in range(dim):
+		tau = 1e-7
+		dX = np.zeros_like(X)
+		dX[:,k] = tau
+		dV = basis.vandermonde(X + dX)
+		DV_est = (dV - V)/tau
+		err = DV[:,:,k]  - DV_est
+		norm_err = np.linalg.norm(err)
+		print(f"err {norm_err:10.6e}")
+		if Basis == HermitePolynomialBasis:
+			assert norm_err < 5e-5
+		else:
+			assert norm_err < 1e-6
+
 if __name__ == '__main__':
 	#test_monomial_total(100, 10, 3, MonomialPolynomialBasis,0) 
 	#test_monomial_max(100, [4], MonomialPolynomialBasis,0) 
 	#test_subspace_angles(10, [1,4], 2, ArnoldiPolynomialBasis)
-	test_scale(MonomialPolynomialBasis, 1)
-
+	#test_scale(MonomialPolynomialBasis, 1)
+	test_vandermonde_derivative(ArnoldiPolynomialBasis, 2, (3,2))
+	pass