Fixed a scaling issue when computing roots

polyrat/arnoldi.py

@@ -165,11 +165,6 @@ def roots(self, coef, *args, **kwargs):
 		poly = PolynomialApproximation(self.degree, Basis = LegendrePolynomialBasis)
 		poly.fit(self.X, y)
 		roots = poly.roots()
-#		def print_roots(r):
-#			fr = self.vandermonde(r.reshape(-1,1)) @ coef
-#			for ri, fri in zip(r, fr):
-#				print(f'root: {ri.real:+10.5e} {ri.imag:+10.5e} I \t abs fun value', np.abs(fri))
-#
-#		print_roots(roots)
+
 		return roots
 

polyrat/basis.py

@@ -49,6 +49,8 @@ def _scale(self, X):
 		"""
 		return 2*(X-self._lb[None,:])/(self._ub[None,:] - self._lb[None,:]) - 1
 
+	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
@@ -81,7 +83,7 @@ def _vander(self, *args, **kwargs):
 	def _der(self, *args, **kwargs):
 		return polyder(*args, **kwargs)
 	def roots(self, *args, **kwargs):
-		return polyroots(*args, **kwargs)
+		return self._inv_scale(polyroots(*args, **kwargs))
 
 
 class LegendrePolynomialBasis(TensorProductPolynomialBasis):
@@ -90,23 +92,23 @@ def _vander(self, *args, **kwargs):
 	def _der(self, *args, **kwargs):
 		return legder(*args, **kwargs)
 	def roots(self, *args, **kwargs):
-		return legroots(*args, **kwargs)
+		return self._inv_scale(legroots(*args, **kwargs))
 
 class ChebyshevPolynomialBasis(TensorProductPolynomialBasis):
 	def _vander(self, *args, **kwargs):
 		return chebvander(*args, **kwargs)
 	def _der(self, *args, **kwargs):
 		return chebder(*args, **kwargs)
 	def roots(self, *args, **kwargs):
-		return chebroots(*args, **kwargs)
+		return self._inv_scale(chebroots(*args, **kwargs))
 
 class HermitePolynomialBasis(TensorProductPolynomialBasis):
 	def _vander(self, *args, **kwargs):
 		return hermvander(*args, **kwargs)
 	def _der(self, *args, **kwargs):
 		return hermder(*args, **kwargs)
 	def roots(self, *args, **kwargs):
-		return hermroots(*args, **kwargs)
+		return self._inv_scale(hermroots(*args, **kwargs))
 
 	def _set_scale(self, X):
 		self._mean = np.mean(X, axis = 0)
@@ -115,14 +117,16 @@ def _set_scale(self, X):
 	def _scale(self, X):		
 		return (X - self._mean[None,:])/self._std[None,:]/np.sqrt(2)
 
+	def _inv_scale(self, X):
+		return np.sqrt(2)*self._std[None,:]*X + self._mean[None,:]
 
 class LaguerrePolynomialBasis(TensorProductPolynomialBasis):
 	def _vander(self, *args, **kwargs):
 		return lagvander(*args, **kwargs)
 	def _der(self, *args, **kwargs):
 		return lagder(*args, **kwargs)
 	def roots(self, *args, **kwargs):
-		return lagroots(*args, **kwargs)
+		return self._inv_scale(lagroots(*args, **kwargs))
 
 	def _set_scale(self, X):
 		r""" Laguerre polynomial expects x[i] to be distributed like exp[-lam*x] on [0,infty)
@@ -134,6 +138,8 @@ def _set_scale(self, X):
 	def _scale(self, X):
 		return self._a[None,:]*(X - self._lb[None,:])
 
+	def _inv_scale(self, X):
+		return X/self._a[None,:] + self._lb[None,:]
 
 
 

tests/test_arnoldi.py

@@ -21,12 +21,17 @@ def wilkinson(x):
 	y = wilkinson(X)
 	arn = PolynomialApproximation(n, Basis = ArnoldiPolynomialBasis)
 	arn.fit(X, y)
-	roots = arn.roots()
+	roots = arn.roots().flatten()
+	I = hungarian_sort(true_roots, roots)
+	roots = roots[I]	
+	for tr, r, fr in zip(true_roots, roots, arn(roots.reshape(-1,1))):	
+		print(f'true root: {tr.real:+10.5e} {tr.imag:+10.5e}I \t root: {r.real:+10.5e} {r.imag:+10.5e} I \t abs fun value {np.abs(fr):10.5e}')
+
 	print("computed roots", roots)
 	print("true roots    ", true_roots)
 	err = sorted_norm(roots, true_roots, np.inf)
 	print("error", err)	
 	assert err < 1e-7, "Error too large"
 
 if __name__ == '__main__':
-	test_arnoldi_roots(15)
+	test_arnoldi_roots(20)

tests/test_basis.py

@@ -163,9 +163,29 @@ def test_subspace_angles(seed, degree, dim,  Basis):
 		assert np.max(ang) < tol
 
 
+@pytest.mark.parametrize("Basis", 
+	[MonomialPolynomialBasis,
+	 LegendrePolynomialBasis,
+	 ChebyshevPolynomialBasis,
+	 HermitePolynomialBasis, 
+	 LaguerrePolynomialBasis])
+@pytest.mark.parametrize("dim", [1,2,3])
+def test_scale(Basis, dim):
+	X = np.random.randn(1000, dim)
+	basis = Basis(X, 2)
+	Y = basis._scale(X)
+	Z = basis._inv_scale(Y)
+	err = np.linalg.norm(Z - X, np.inf)
+	print("error", err)
+	assert err < 1e-10, "Scaling/inverse scaling does not map to identity"
 
 
 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_subspace_angles(10, [1,4], 2, ArnoldiPolynomialBasis)
+	test_scale(MonomialPolynomialBasis, 1)
+
+
+
+

tests/test_polynomial.py

@@ -0,0 +1,46 @@
+import numpy as np
+import pytest
+from polyrat import *
+
+
+
+@pytest.mark.parametrize("Basis", 
+	[MonomialPolynomialBasis,
+	 LegendrePolynomialBasis,
+	 ChebyshevPolynomialBasis,
+	 HermitePolynomialBasis, 
+	 LaguerrePolynomialBasis,
+	 ArnoldiPolynomialBasis])
+@pytest.mark.parametrize("n", [2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 20])
+def test_wilkinson_roots(Basis, n):
+	r""" Check root computation in Arnoldi polynomials
+	"""
+
+	if Basis in [LaguerrePolynomialBasis, HermitePolynomialBasis] and n>= 8:
+		# These tests fail due to the ill-conditioning of this basis
+		return
+
+	true_roots = np.arange(1, n+1)
+
+	def wilkinson(x):
+		value = np.zeros(x.shape, dtype = np.complex)
+		for i, xi in enumerate(x):
+			value[i] = np.prod(xi - true_roots)
+		return value
+
+	# It is important that we sample at the roots to avoid the large
+	# values the Wilkinson polynomial takes away from these points.
+	X = np.arange(0, n+1, step = 0.1, dtype = np.float).reshape(-1,1)
+	y = wilkinson(X)
+	poly = PolynomialApproximation(n, Basis = Basis)
+	poly.fit(X, y)
+	roots = poly.roots().flatten()
+	I = hungarian_sort(true_roots, roots)
+	roots = roots[I]	
+	for tr, r, fr in zip(true_roots, roots, poly(roots.reshape(-1,1))):	
+		print(f'true root: {tr.real:+10.5e} {tr.imag:+10.5e}I \t root: {r.real:+10.5e} {r.imag:+10.5e} I \t abs fun value {np.abs(fr):10.5e}')
+
+	err = sorted_norm(roots, true_roots, np.inf)
+	print("error", err)	
+	assert err < 1e-7, "Error too large"
+