Merging improvements from rank-constrained branch

.gitignore

@@ -1,3 +1,4 @@
+*.DS_Store
 *.swp
 *.pyc
 __pycache__

polyrat/aaa.py

@@ -4,8 +4,7 @@
 import numpy as np
 import scipy.linalg
 from iterprinter import IterationPrinter
-from .rational import RationalBarycentric
-
+from .rational import RationalFunction
 #def eval_barcentric(xeval, x, y, I, a, b):
 #	"""
 #	Parameters
@@ -58,7 +57,11 @@ def eval_aaa(xeval, x, y, I, b):
 
 
 def _build_cauchy(x,y):
-	return 1./(np.tile(x.reshape(-1,1), (1,len(y))) - np.tile(y.reshape(1,-1), (len(x),1)))
+	with np.errstate(divide = 'ignore'):
+		# FIXME: Properly handle divide by zeros in the Cauchy matrix
+		C = 1./(np.tile(x.reshape(-1,1), (1,len(y))) - np.tile(y.reshape(1,-1), (len(x),1)))
+		C[~np.isfinite(C)] = 1
+	return C
 
 def aaa(x, y, degree = None, tol = None, verbose = True):
 	r""" A vector-valued Adaptive Anderson-Antoulas implementation
@@ -123,6 +126,27 @@ def aaa(x, y, degree = None, tol = None, verbose = True):
 
 	return I, b 	
 
+class RationalBarycentric(RationalFunction):
+	r"""
+	"""
+	def __init__(self, degree):
+		self.degree = int(degree)
+		assert self.degree >= 0, "Degree must be non-negative"
+
+	@property
+	def num_degree(self):
+		return self.degree
+
+	@property
+	def denom_degree(self):
+		return self.degree
+
+	def poles(self):
+		raise NotImplementedError
+
+	def residues(self):
+		raise NotImplementedError
+
 class AAARationalApproximation(RationalBarycentric):
 
 	def __init__(self, degree = None, tol = None, verbose = True):

polyrat/arnoldi.py

@@ -246,10 +246,10 @@ def vandermonde_derivative(self, X, weight = None):
 	def roots(self, coef, *args, **kwargs):
 		from .basis import LegendrePolynomialBasis
 		from .polynomial import PolynomialApproximation
-		y = self.vandermonde_X @ coef
+		y = (self.vandermonde_X @ coef).flatten()
 		poly = PolynomialApproximation(self.degree, Basis = LegendrePolynomialBasis)
 		poly.fit(self.X, y)
-		roots = poly.roots()
+		roots = poly.roots().flatten()
 
 		return roots
 

polyrat/lagrange.py

@@ -173,6 +173,9 @@ def __init__(self, nodes):
 			np.prod(self.nodes[k] - self.nodes[0:k]) * np.prod(self.nodes[k] - self.nodes[k+1:]) for k in range(len(self.nodes))
 			])
 
+	@property
+	def dim(self):
+		return 1 
 
 	@cached_property
 	def vandermonde_X(self):
@@ -188,6 +191,7 @@ def roots(self, coef, deflation = True):
 		return lagrange_roots(self.nodes, self.weights, coef, deflation = deflation)
 
 
+
 class LagrangePolynomialInterpolant(Polynomial):
 	def __init__(self, X, y):
 		self.basis = LagrangePolynomialBasis(X)

polyrat/paaa.py

@@ -8,7 +8,7 @@
 import scipy.linalg
 from iterprinter import IterationPrinter
 
-from .rational import RationalBarycentric
+from .aaa import RationalBarycentric
 
 def _build_cauchy(x,y):
 	with np.errstate(divide = 'ignore', invalid = 'ignore'):

polyrat/polynomial.py

@@ -4,7 +4,7 @@
 from copy import deepcopy
 import scipy.linalg
 import cvxpy as cp
-
+from .util import _zeros
 
 
 class Polynomial:
@@ -32,13 +32,32 @@ def eval(self, X):
 		#return self.basis.vandermonde(X) @ self.coef
 		return np.einsum('ij,j...->i...', self.basis.vandermonde(X), self.coef)
 
-	def roots(self, *args, **kwargs): 
-		return self.basis.roots(self.coef, *args, **kwargs)	
+	def derivative(self, X):
+		r""" Compute the derivative 
+		"""
+		return np.einsum('ijk,j...->i...k', self.basis.vandermonde_derivative(X), self.coef)
+
+	def roots(self, *args, **kwargs):
+		assert self.basis.dim == 1, "Must have a single variable as input"
+		assert np.prod(self.coef[0].shape) == 1, "Must have one dimensional output"
+		return self.basis.roots(self.coef.reshape(-1), *args, **kwargs)	
 
 
-def _polynomial_fit_least_squares(P, y):
-	coef, _, _, _ = scipy.linalg.lstsq(P, y)
-	return coef.flatten()
+def _polynomial_fit_least_squares(P, Y, P_orth = False):
+	M, m = P.shape
+
+	coef = _zeros((P.shape[1], *Y.shape[1:]), P, Y)
+
+	if P_orth:
+		for j, idx in enumerate(np.ndindex(Y.shape[1:])):
+			coef[(slice(m), *idx)] = P.T.conj() @ Y[(slice(M), *idx)]
+	else:
+		Q, R = scipy.linalg.qr(P, mode = 'economic')
+		for j, idx in enumerate(np.ndindex(Y.shape[1:])):
+			a = scipy.linalg.solve_triangular(R, Q.T.conj() @ Y[(slice(M), *idx)])
+			coef[(slice(m), *idx)] = a
+
+	return coef
 
 def _polynomial_fit_pnorm(P, y, norm, **kwargs):
 	if np.iscomplexobj(P) or np.iscomplexobj(y):
@@ -65,10 +84,11 @@ def degree(self):
 		return self._degree
 
 	def fit(self, X, y, **kwargs):
+		from .arnoldi import ArnoldiPolynomialBasis
 		self.basis = self.Basis(X, self.degree)
 		P = self.basis.vandermonde_X
 		if self.norm == 2 or self.norm == 2.:
-			self.coef = _polynomial_fit_least_squares(P, y)
+			self.coef = _polynomial_fit_least_squares(P, y, isinstance(self.Basis, ArnoldiPolynomialBasis))
 		else:
 			self.coef = _polynomial_fit_pnorm(P, y, self.norm, **kwargs)
 

polyrat/rational.py

@@ -24,6 +24,12 @@ def eval(self, X):
 		"""
 		return self.__call__(X)
 
+	@abc.abstractmethod
+	def poles(self, *args, **kwargs):
+		r""" Poles of this rational function
+		"""
+		raise NotImplementedError
+
 class RationalApproximation(RationalFunction):
 	def __init__(self, num_degree, denom_degree):
 		self.num_degree = num_degree
@@ -74,8 +80,6 @@ def eval(self, X):
 		q = self.denominator(X)
 		return np.multiply(1./q.reshape(-1, *[1 for i in p.shape[1:]]), p)
 
-
-
 	def refine(self, X, y, norm = 2, verbose = False, **kwargs):
 		r""" Refine the rational approximation using optimization
 
@@ -92,25 +96,23 @@ def refine(self, X, y, norm = 2, verbose = False, **kwargs):
 		#	res_norm = np.linalg.norm( (self.P @ a)/(self.Q @ b) - y, norm)
 		#	print(f"final residual norm {res_norm:21.15e}")
 
+	def pole_residue(self, *args, **kwargs):
+		poles = self.denominator.roots(*args, **kwargs).flatten()
+
+		residues = []
+		for k, x in enumerate(poles):
+			R = self.numerator(x.reshape(-1,1))/self.denominator.derivative(x.reshape(-1,1))
+			residues.append(R.reshape(self.a.shape[1:]))
 
+		residues = np.array(residues)
+
+		return poles, residues	
+
 
-class RationalBarycentric(RationalFunction):
-	r"""
-	"""
-	def __init__(self, degree):
-		self.degree = int(degree)
-		assert self.degree >= 0, "Degree must be non-negative"
-
-	@property
-	def num_degree(self):
-		return self.degree
-
-	@property
-	def denom_degree(self):
-		return self.degree
-
-
+	def poles(self, *args, **kwargs):
+		return self.denominator.roots(*args, **kwargs)
 
+
 
 
 

polyrat/util.py

@@ -6,6 +6,14 @@
 import scipy.linalg
 
 
+def _zeros(size, *args):
+	r""" allocate a zeros matrix of the given size matching the type of the arguments
+	"""
+	if all([np.isrealobj(a) for a in args]):
+		return np.zeros(size, dtype = np.float)
+	else:
+		return np.zeros(size, dtype = np.complex)
+
 def linearized_ratfit_operator_dense(P, Q, Y):
 	r""" Dense analog of LinearizedRatfitOperator
 	"""

polyrat/vecfit.py

@@ -162,6 +162,8 @@ def eval(self, X):
 		nout_dim = len(self._a.shape[1:])
 		return np.multiply(1./denom.reshape(-1, *([1,]*nout_dim)), num)
 
+	def pole_residue(self):
+		return self.poles, self._a
 
 
 class VectorFittingRationalApproximation(VectorFittingRationalFunction):

tests/test_arnoldi.py

@@ -18,12 +18,16 @@ def wilkinson(x):
 	# 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)
+	y = wilkinson(X).flatten()
+
 	arn = PolynomialApproximation(n, Basis = ArnoldiPolynomialBasis)
 	arn.fit(X, y)
 	roots = arn.roots().flatten()
 	I = hungarian_sort(true_roots, roots)
 	roots = roots[I]	
+	print("true_roots", true_roots)
+	print("roots", roots)
+	print("value", arn(roots.reshape(-1,1)))
 	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}')
 

tests/test_data.py

@@ -33,13 +33,13 @@ def array_absolute_value(M, output_dim = ()):
 	return X.reshape(-1,1), Y
 
 
-def random_data(M, dim, complex_, seed ):
+def random_data(M, dim, complex_, seed, output_dim = ()):
 	np.random.seed(seed)
 	X = np.random.randn(M, dim)
-	y = np.random.randn(M)
+	y = np.random.randn(M, *output_dim)
 	if complex_:
 		X = X + 1j*np.random.randn(M, dim)
-		y = y + 1j*np.random.randn(M)
+		y = y + 1j*np.random.randn(M, *output_dim)
 
 	return X, y
 

tests/test_polynomial.py

@@ -31,7 +31,7 @@ def wilkinson(x):
 	# 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)
+	y = wilkinson(X).flatten()
 	poly = PolynomialApproximation(n, Basis = Basis)
 	poly.fit(X, y)
 	roots = poly.roots().flatten()
@@ -66,6 +66,55 @@ def test_approx(complex_, Basis, norm):
 	assert np.all(np.isclose(p(X), y2))	
 
 
+@pytest.mark.parametrize("Basis", 
+	[MonomialPolynomialBasis,
+	 LegendrePolynomialBasis,
+	 ChebyshevPolynomialBasis,
+	 HermitePolynomialBasis, 
+	 LaguerrePolynomialBasis,
+	 ArnoldiPolynomialBasis])
+@pytest.mark.parametrize("dim", [1,2])
+@pytest.mark.parametrize("output_dim", [None, 1, 3, (3,2)]) 
+def test_derivative(Basis, dim, output_dim):
+	np.random.seed(0)
+
+	N = 1000
+	degree = 5
+	X = np.random.randn(N,dim)
+	if output_dim is None:
+		y = np.random.randn(N)
+	else:
+		try:
+			y = np.random.randn(N, output_dim)
+		except TypeError:
+			y = np.random.randn(N, *output_dim)
+
+	poly = PolynomialApproximation(degree, Basis = Basis)
+	poly.fit(X, y)
+
+	Xhat = np.random.randn(5,dim)
+
+	D = poly.derivative(Xhat)
+
+	h = 1e-6
+	for i in range(len(Xhat)):
+		for k in range(dim):
+			ek = np.eye(dim)[k]
+			x1 = (Xhat[i] + h*ek).reshape(1,-1)
+			x2 = (Xhat[i] - h*ek).reshape(1,-1)
+			dest = (poly(x1) - poly(x2))/(2*h)
+			print(f"point {i}, direction {k}")
+			print("finite difference")
+			print(dest)
+			print("nomial value")
+			print(D[i,...,k])
+			print('difference')
+			print(D[i,...,k] - dest)
+			assert np.all(np.isclose(dest, D[i,...,k], atol = 1e-4, rtol = 1e-4))
+
+
+
 if __name__ == '__main__':
-	test_approx(True, None, 1)
+	#test_approx(True, None, 1)
+	test_derivative(MonomialPolynomialBasis,2, (2,1))