Implemented a provisional AAA

Reproducibility/ISK/fig_refine.py

@@ -32,10 +32,10 @@
 		sk = SKRationalApproximation(m,n, verbose = True, rebase = True, refine = False, maxiter = 10)
 		sk.fit(X,y)
 
-		err_isk[k] = np.linalg.norm(sk(X) - y, 2)
+		err_isk[k] = np.linalg.norm(sk(X) - y, 2)/np.linalg.norm(y, 2)
 
 		sk.refine(X, y)
-		err_iskr[k] = np.linalg.norm(sk(X) - y, 2)
+		err_iskr[k] = np.linalg.norm(sk(X) - y, 2)/np.linalg.norm(y, 2)
 
 		print(f"original {err_isk[k]:10.3e}; improved {err_iskr[k]:10.3e}")
 

Reproducibility/ISK/fig_scalar.py

@@ -0,0 +1,86 @@
+import numpy as np
+from polyrat import *
+from pgf import PGF
+from scipy.special import jv
+
+data = []
+
+# Two real-valued input problems
+
+# Absolute value test problem
+M = int(2e3)
+X = np.linspace(-1,1,M).reshape(-1,1)
+y = np.abs(X).flatten()
+data += [['abs', X, y]]
+
+# Exponential (NST18: Sec. 6.8)
+M = 4000
+X = -np.logspace(-3, 4, M).reshape(-1,1)
+y = np.exp(X.flatten())
+data += [['exp', X, y]]
+
+# Two complex-valued input problems
+
+# Tangent
+M = 1000
+X = np.exp(1j*np.linspace(0, 2*np.pi, M, endpoint = False)).reshape(-1,1)
+y = np.tan(256*X.flatten())
+data += [['tan256', X, y]]
+
+# Bessel function (NST18 Fig 6.5)
+# This converges too fast
+M = 2000
+np.random.seed(0)
+X = 10*np.random.rand(M) + 2j*(np.random.rand(M) - 0.5)
+y = 1./jv(0, X)
+data += [['bessel', X.reshape(-1,1), y]]
+
+# Log example: NST18 Fig 6.2;
+# note rate is independent of offset, so shinking does not harm convergence
+M = int(2e3)
+X = np.exp(1j*np.linspace(0, 2*np.pi, M, endpoint = False)).reshape(-1,1)
+y = np.log(1.1 - X.flatten())
+data += [['log11', X, y]]
+
+
+data = data[2:3]
+
+# Range of parameters
+mns = [(k,k) for k in range(2, 51, 1)]
+
+
+
+for name, X, y in data:
+
+	err_isk = np.zeros(len(mns))
+	err_iskr = np.zeros(len(mns))
+	err_aaa = np.zeros(len(mns))
+
+	for k, (m, n) in enumerate(mns):
+
+		print(f'\n======== AAA ({m},{m}) | {name} =======\n')
+		aaa = AAARationalApproximation(degree = m)
+		aaa.fit(X, y)
+		err_aaa[k] = np.linalg.norm(aaa(X) - y, 2)/np.linalg.norm(y,2)
+
+
+		print(f'\n======== SK Rebase ({m},{n}) | {name} =======\n')
+		sk = SKRationalApproximation(m,n, verbose = True, rebase = True, refine = False, maxiter = 10)
+		sk.fit(X,y)
+
+		err_isk[k] = np.linalg.norm(sk(X) - y, 2)/np.linalg.norm(y, 2)
+
+		sk.refine(X, y)
+		err_iskr[k] = np.linalg.norm(sk(X) - y, 2)/np.linalg.norm(y, 2)
+
+		print(f"original {err_isk[k]:10.3e}; improved {err_iskr[k]:10.3e}")
+
+
+	pgf = PGF()
+	pgf.add('m', [mn[0] for mn in mns])
+	pgf.add('n', [mn[1] for mn in mns])
+	pgf.add('err_isk', err_isk)
+	pgf.add('err_iskr', err_iskr)
+	pgf.add('err_aaa', err_aaa)
+
+	pgf.write(f'data/fig_scalar_{name}.dat')

polyrat/aaa.py

@@ -0,0 +1,124 @@
+import numpy as np
+import scipy.linalg
+from iterprinter import IterationPrinter
+
+
+def eval_barcentric(xeval, x, y, I, a, b):
+	"""
+	Parameters
+	----------
+	xeval: np.array
+		Locations to evaluate the function
+	x: np.array
+		input coordinates for fit
+	"""
+
+	with np.errstate(divide='ignore',invalid='ignore'):
+		A = np.hstack([ 1. /(xeval - xi) for xi in zip(x[I])])
+		# For the rows that have a divide by zero error, 
+		# we replace with value corresponding to the limit
+		for row in np.argwhere(~np.all(np.isfinite(A), axis = 1)):
+			A[row] = np.zeros(A.shape[1])
+			A[row, np.argmin(np.abs(xeval[row] - x[I])).flatten()] = 1 
+
+	num = np.einsum('ij,j...->i...', A, a)
+	denom = np.einsum('ij,j->i', A,b)
+	reval = np.einsum('i...,i->i...', num, 1./denom)
+
+	return reval	
+
+
+
+def eval_aaa(xeval, x, y, I, b):
+	"""
+	Parameters
+	----------
+	xeval: np.array
+		Locations to evaluate the function
+	x: np.array
+		input coordinates for fit
+	"""
+	xeval = xeval.reshape(-1,1)
+	with np.errstate(divide='ignore',invalid='ignore'):
+		A = np.hstack([ bi /(xeval - xi) for bi, xi in zip(b, x[I])])
+		# For the rows that have a divide by zero error, 
+		# we replace with value corresponding to the limit
+		for row in np.argwhere(~np.all(np.isfinite(A), axis = 1)):
+			A[row] = np.zeros(A.shape[1])
+			A[row, np.argmin(np.abs(xeval[row] - x[I])).flatten()] = 1 
+
+	num = np.einsum('ij,j...->i...', A, y[I])
+	denom = np.sum(A, axis = 1)
+	reval = np.einsum('i...,i->i...', num, 1./denom)
+
+	return reval	
+
+
+
+
+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)))
+
+
+def aaa(x, y, degree = None, tol = None, verbose = True):
+	r""" A vector-valued Adaptive Anderson-Antoulas implementation
+
+	Parameters
+	----------
+	x: np.array (M,)
+		input coordinates
+	y: np.array (M,...)
+		output coordinates
+	"""
+
+	assert not ((degree is None) and (tol is None)), "One or both of 'degree' and 'tol' must be specified" 
+
+	if degree is None:
+		degree = len(x)//2 - 1
+
+
+	mismatch = y
+	I = np.zeros(len(y), dtype = np.bool)	# Index of x values used as barycentric nodes 
+
+	if verbose:
+		printer = IterationPrinter(it = '4d', res = '20.16e', cond = '8.3e')
+		printer.print_header(it = 'iter', res = 'Frobenius norm residual', cond = 'condition number')
+
+	for it in range(degree+1):
+		# TODO: Is sum best? or how about maximum?
+		residual = np.sum(np.abs(mismatch), axis = tuple(range(1,len(y[0].shape)+1)))
+		residual[I] = 0 	# zero out residual at nodes that have already been selected
+
+		Inew = np.argmax(residual)
+		I[Inew] = True
+
+		# Construct the Cauchy matrix 
+		C = _build_cauchy(x[~I], x[I])
+
+		# Build the Loewner matrix associated with each input
+		Lten = []
+		for idx in np.ndindex(y[0].shape):
+			Lten.append( (C.T * y[(~I,*idx)]).T - C*y[(I,*idx)] )
+		L = np.vstack(Lten)
+
+		# Compute coefficients for denominator polynomial
+		U, s, VH = scipy.linalg.svd(L, full_matrices = False, compute_uv = True, overwrite_a = True)
+		b = VH.conj().T[:,-1] 
+		if len(s) >= 2:
+			with np.errstate(divide='ignore',invalid='ignore'):
+				cond = s[0]*np.sqrt(2)/(s[-2] - s[-1])
+		else:
+			cond = None
+
+		mismatch = y - eval_aaa(x, x, y, I, b) 	
+		res_norm = np.sqrt(np.sum(np.abs(mismatch)**2))
+
+		if verbose:
+			printer.print_iter(it=it, res = res_norm, cond = cond)	
+
+		if tol is not None:
+			if tol > res_norm:
+				if verbose: print("terminated due to small residual")
+				break
+
+	return I, b 	

polyrat/rational.py

@@ -5,14 +5,15 @@
 from .polynomial import *
 from .skiter import *
 from .rational_ratio import *
+from .aaa import *
 from copy import deepcopy
 
 
 class RationalFunction:
 	pass
 
 
-class RationalApproximation:
+class RationalApproximation(RationalFunction):
 	def __init__(self, num_degree, denom_degree):
 		self.num_degree = num_degree
 		self.denom_degree = denom_degree
@@ -38,6 +39,11 @@ def a(self):
 	@property
 	def b(self):
 		return self.denominator.coef
+
+	def __call__(self, X):
+		p = self.numerator(X)
+		q = self.denominator(X)
+		return p/q	
 
 	def refine(self, X, y, **kwargs):
 		a, b = rational_ratio_optimize(y, self.P, self.Q, self.a, self.b, norm = self.norm, **kwargs)
@@ -50,6 +56,12 @@ def refine(self, X, y, **kwargs):
 			print(f"final residual norm {res_norm:21.15e}")
 
 
+class LinearizedRationalApproximation(RationalApproximation, RationalRatio):
+	def __init__(self, num_degree, denom__degree):
+		RationalApproximation.__init__(self, num_degree, denom_degree)
+
+	def fit(self, X, y):
+		self.numerator, self.denominator = linearized_ratfit(X, y, self.num_degree, self.denom_degree)
 
 class SKRationalApproximation(RationalApproximation, RationalRatio):
 	r"""
@@ -106,11 +118,50 @@ def fit(self, X, y, denom0 = None):
 			self.refine(X, y)
 
 
-
+
+
+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
+
+
+class AAARationalApproximation(RationalBarycentric):
+
+	def __init__(self, degree = None, tol = None, verbose = True):
+		self.degree = degree
+		self.tol = tol
+		self.verbose = verbose
+
+	def fit(self, X, y):
+		X = np.array(X)
+		self.y = np.array(y)
+		assert len(X) == len(y), "Length of X and y do not match"
+		self.x = X.flatten()
+		assert len(self.x) == len(y), "AAA only supports scalar-valued inputs"
+
+		self.I, self.b = aaa(self.x, self.y, degree = self.degree, tol = self.tol, verbose = self.verbose)
 
 	def __call__(self, X):
-		p = self.numerator(X)
-		q = self.denominator(X)
-		return p/q	
-
+		x = np.array(X).flatten().reshape(-1,1)
+		assert len(x) == len(X), "X must be a scalar-valued input"
+		return eval_aaa(x, self.x, self.y, self.I, self.b)
+
+
+class AAALawsonRationalApproximation(RationalBarycentric):
+	pass
+
+
+
 
+

polyrat/skiter.py

@@ -78,6 +78,24 @@ def _minimize_1_norm(A):
 #	return x, s
 
 
+def linearized_ratfit(X, y, num_degree, denom_degree):
+	r""" This solves the linearized rational approximation problem
+
+	See: AKL+19x
+	"""
+	num_basis = ArnoldiPolynomialBasis(X, num_degree)
+	denom_basis = ArnoldiPolynomialBasis(X, denom_degree)
+
+
+	A = np.hstack([P, np.multiply(-y[:,None], Q) ])
+	x, cond = linearized_solution(A)
+	a = x[:P.shape[1]]
+	b = x[-Q.shape[1]:]
+
+	numerator = Polynomial(num_basis, a)
+	denominator = Polynomial(denom_basis, b)
+
+	return numerator, denominator
 
 def skfit(y, P, Q, maxiter = 20, verbose = True, history = False, denom0 = None, norm = 2, xtol = 1e-7):
 	r"""