Adding more documentation

doc/basis.rst

@@ -128,9 +128,18 @@ Arnoldi Polynomial Basis
 .. autoclass:: polyrat.ArnoldiPolynomialBasis
 
 
+For low-level access, the following functions are available. 
+
+.. autofunction:: polyrat.vandermonde_arnoldi_CGS 
+
+.. autofunction:: polyrat.vandermonde_arnoldi_eval
+
+.. autofunction:: polyrat.vandermonde_arnoldi_eval_der
 
 
 Lagrange Polynomial Basis
 =========================
 
 .. autoclass:: polyrat.LagrangePolynomialBasis
+
+.. autofunction:: polyrat.lagrange_roots

doc/references.rst

@@ -1,6 +1,21 @@
 References
 ==========
 
-.. [AKL+19x] Anthony P. Austin, Mohan Krishnamoorthy, Sven Leyffer, Stephen Mrenna, Juliane Müller, Holger Schulz.
-	Multivariate Rational Approximation.
+.. [AKL+20] Anthony P. Austin, Mohan Krishnamoorthy, Sven Leyffer, Stephen Mrenna, Juliane Müller, Holger Schulz.
+	Practical algorithms for multivariate rational approximation.
+	https://doi.org/10.1016/j.cpc.2020.107663
 	https://arxiv.org/abs/1912.02272
+
+.. [BNT19x] Pablo D. Brubeck, Yuji Nakatsukasa, and Lloyd N. Trefethen.
+	Vandermonde with Arnoldi.
+	https://arxiv.org/abs/1911.09988
+
+.. [BT04] Jean-Paul Berrut and Llyod N. Trefethen.
+	Barycentric Lagrange Interpolation.
+	SIAM Review, Vol 46, No 3, pp. 501--517, 2004.
+	https://doi.org/10.1137/S0036144502417715
+
+.. [LC14] Peirs W. Lawrence and Robert M. Corless.
+	Stability of rootfinding for barycentric Lagrange interpolants.
+	Numer. Algor. Vol. 65, pp. 447--464 (2014).
+	https://doi.org/10.1007/s11075-013-9770-3

polyrat/aaa.py

@@ -33,7 +33,7 @@
 
 
 def eval_aaa(xeval, x, y, I, b):
-	"""
+	""" Evaluate an AAA rational approximation
 	Parameters
 	----------
 	xeval: np.array

polyrat/arnoldi.py

@@ -24,6 +24,15 @@ def _update_rule(indices, idx):
 def vandermonde_arnoldi_CGS(X, degree, weight = None, mode = None):
 	r""" Multivariate Vandermode with Arnoldi using classical Gram-Schmidt with reorthogonalization
 
+
+	This function uses the Arnoldi proceedure to construct an orthogonal
+	basis for polynomials.  In the univariate case, this corresponds to
+	constructing a Krylov subspace for the starting vector :math:`\mathbf{b}`
+	and matrix :math:`\text{diag}(\mathbf{x})`; see [BNT19x]_. The multivariate
+	case uses a similar approach, but care is taken to use the right
+	combination of input coordinates.; see [AKL+20]_.
+
+
 	Notes
 	-----
 	* The use of Classical Gram-Schmidt with reorthogonalization
@@ -35,8 +44,8 @@ def vandermonde_arnoldi_CGS(X, degree, weight = None, mode = None):
 
 	Parameters
 	----------
-	X: np.array (M, dim)
-		Input coordinates
+	X: np.array
+		Input coordinates of size (M, dim)
 	degree: int or list of ints
 		Polynomial degree.  If an int, a total degree polynomial is constructed;
 		if a list, the list must be length m and a maximum degree polynomial is
@@ -45,7 +54,17 @@ def vandermonde_arnoldi_CGS(X, degree, weight = None, mode = None):
 		Initial vector in the Arnoldi iteration
 	mode: None or ['total', 'max']
 		What type of polynomial basis to construct; only matters if dim>1.
-		If None, the type of basis will be automatically detected. 
+		If None, the type of basis will be automatically detected.
+
+	
+	Returns
+	-------
+	Q: :class:`~numpy:numpy.ndarray`
+		Orthonormal basis for the desired polynomials on the input coordinates.
+	R: :class:`~numpy:numpy.ndarray`
+		Matrix containing orthogonalization information; needed to evaluate polynomial at new coordinates.
+	indices: :class:`~numpy:numpy.ndarray`
+		List of indices showing the order in which the basis was constructed.	 
 	"""
 	M, dim = X.shape
 
@@ -98,8 +117,27 @@ def vandermonde_arnoldi_CGS(X, degree, weight = None, mode = None):
 
 
 
-def vandermonde_arnoldi_eval(X, R, indices, mode, weight = None):
-	r"""
+def vandermonde_arnoldi_eval(X, R, indices, weight = None):
+	r""" Evaluate a Vandermonde with Arnoldi polynomial
+
+	Parameters
+	----------
+	X: array_like
+		Coordiantes at which to evaluate the polynomial 
+	R: :class:`~numpy:numpy.ndarray`
+		Upper triangular matrix containing orthogonalization information from 
+		:class:`~polyrat:polyrat.vandermonde_arnoldi_CGS`.
+	indices: :class:`~numpy:numpy.ndarray`
+		The ordering of polynomial powers returned from
+		:class:`~polyrat:polyrat.vandermonde_arnoldi_CGS`.
+	weight: None or array_like
+		Starting vector for Arnoldi; by default is the ones vector.
+		In general, this should not change without good reason.
+
+	Returns
+	-------
+	W: :class:`~numpy:numpy.ndarray`
+		Polynomial basis evaluated at the specified points.
 	"""
 
 	X = np.array(X)
@@ -129,9 +167,32 @@ 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):
+def vandermonde_arnoldi_eval_der(X, R, indices, weight = None, V = None):
+	r""" Evaluate the derivative of Vandermonde with Arnoldi polynomial basis
+
+	Parameters
+	----------
+	X: array_like
+		Coordiantes at which to evaluate the polynomial 
+	R: :class:`~numpy:numpy.ndarray`
+		Upper triangular matrix containing orthogonalization information from 
+		:class:`~polyrat:polyrat.vandermonde_arnoldi_CGS`.
+	indices: :class:`~numpy:numpy.ndarray`
+		The ordering of polynomial powers returned from
+		:class:`~polyrat:polyrat.vandermonde_arnoldi_CGS`.
+	weight: None or array_like
+		Starting vector for Arnoldi; by default is the ones vector.
+		In general, this should not change without good reason.
+	V: :class:`~numpy:numpy.ndarray`
+		Evaluation of the polynomial without the derivative at X.
+	
+	Returns
+	-------
+	W: :class:`~numpy:numpy.ndarray`
+		Polynomial basis deriative evaluated at the specified points.
+	"""
 	if V is None:
-		V = vandermonde_arnoldi_eval(X, R, indices, mode, weight = weight)
+		V = vandermonde_arnoldi_eval(X, R, indices, weight = weight)
 
 	M = X.shape[0]
 	N = R.shape[1]
@@ -154,7 +215,14 @@ def vandermonde_arnoldi_eval_der(X, R, indices, mode, weight = None, V = None):
 
 
 class ArnoldiPolynomialBasis(PolynomialBasis):
-	r""" A polynomial basis constructed using Vandermonde with Arnoldi
+	r""" A polynomial basis constructed using Vandermonde with Arnoldi.
+
+	A polynomial basis constructed with Vandermonde with Arnoldi;
+	see [BNT19x]_, [AKL+20]_ for discussion.
+
+	Notes
+	-----
+	* In order to compute the roots in the univariate case we convert to Legendre polynomial. 
 
 	"""
 	def __init__(self, X, degree, weight = None):
@@ -167,13 +235,13 @@ def vandermonde_X(self):
 		return self._Q   
 
 	def vandermonde(self, X, weight = None):
-		return vandermonde_arnoldi_eval(X, self._R, self._indices, self.mode, weight = weight)
+		return vandermonde_arnoldi_eval(X, self._R, self._indices, 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)
+			return vandermonde_arnoldi_eval_der(X, self._R, self._indices,  weight = weight, V = self._Q)
 		else:
-			return vandermonde_arnoldi_eval_der(X, self._R, self._indices, self.mode, weight = weight)
+			return vandermonde_arnoldi_eval_der(X, self._R, self._indices,  weight = weight)
 
 	def roots(self, coef, *args, **kwargs):
 		from .basis import LegendrePolynomialBasis

polyrat/lagrange.py

@@ -13,19 +13,32 @@
 def lagrange_roots(nodes, weights, coef, deflation = True):
 	r""" Compute the roots of a Lagrange polynomial
 
-	This implements the deflation algorithm from LC14
-
+	This implements the deflation algorithm from [LC14]_
+	to compute the roots of a Lagrange polynomial 
+	without changing basis and to high accuraccy.
 
 
 	Parameters
 	----------
 	nodes: numpy.array (n,)
 		Nodes :math:`x_j` where the polynomial value is defined.
 	weights: numpy.array (n,)
-		The barycentric weights :math:`prod_{k\ne j} (x_j - x_k)^{-1}`
+		The barycentric weights :math:`\prod_{k\ne j} (\xi_j - \xi_k)^{-1}`
 	coef: numpy.array (n,)
-		The coeffients :math:`c_j` defining the value of the polynomial at each point;
-		i.e., :math:`p(x_j) = c_j`.
+		The coeffients :math:`f_j` defining the value of the polynomial at each point;
+		i.e., :math:`p(\xi_j) = f_j`.
+	deflation: bool
+		In the standard formulation to compute these roots
+		we solve a generalized eigenvalue problem (GEP)
+		which has two infinite poles.
+		If True, we explicitly remove these infinite poles
+		by shrinking the matrices;
+		if False we do not.
+	
+	Returns
+	-------
+	roots: :class:`~numpy:numpy.ndarray`
+		Roots of the polynomial.
 	"""
 	n = len(nodes)
 	assert (n == len(weights)) and  (n == len(coef)), "Dimensions of nodes, weights, and coef should be the same"
@@ -131,9 +144,24 @@ def lagrange_vandermonde(nodes, weights, X):
 
 
 class LagrangePolynomialBasis(PolynomialBasis):
-	r""" Constructs a Lagrange polynomial basis in barycentric form
+	r""" Constructs a Lagrange polynomial basis in barycentric form.
+
+
+	Here we construct a univariate polynomial in barycentric form
+	as described in [BT04]_:
+	
+	.. math::
+
+		p(x) = \sum_{j=1}^d \frac{ f_j w_j (x - \xi_j)^{-1}}{w_j (x - \xi_j)^{-1}},
+		\quad w_j := \prod_{k\ne j} (\xi_j - \xi_k)^{-1}
+
+	such that :math:`p(\xi_j) = f_j`.
+
+	Parameters
+	----------
+	nodes: array_like
+		List of the nodes :math:`\xi_j` specifying the basis.
 
-	See: BT04
 	"""
 	def __init__(self, nodes):
 		self.nodes = np.array(nodes).flatten()
@@ -157,8 +185,6 @@ def vandermonde_derivative(self, X):
 		raise NotImplementedError
 
 	def roots(self, coef, deflation = True):
-		r"""
-		"""
 		return lagrange_roots(self.nodes, self.weights, coef, deflation = deflation)