tessl/pypi-cupy-cuda11x
CuPy: NumPy & SciPy for GPU - CUDA 11.x optimized distribution providing GPU-accelerated computing with Python
—
Pending
polynomial-operations.mddocs/
Polynomial Operations
CuPy provides comprehensive polynomial operations through the cupy.polynomial module and legacy polynomial functions, enabling mathematical computation with polynomials including arithmetic, evaluation, fitting, root finding, and advanced polynomial manipulations optimized for GPU acceleration.
Capabilities
Legacy Polynomial Functions
Legacy polynomial functions providing basic polynomial operations and manipulations.
class poly1d:
"""
One-dimensional polynomial class.
A convenience class for dealing with polynomials, providing
methods for polynomial arithmetic, evaluation, and analysis.
"""
def __init__(self, c_or_r, r=False, variable=None):
"""
Parameters:
c_or_r: array_like - Polynomial coefficients or roots
r: bool, optional - If True, c_or_r contains roots instead of coefficients
variable: str, optional - Variable name for string representation
"""
def __call__(self, val):
"""Evaluate polynomial at given values."""
def __add__(self, other):
"""Add polynomials."""
def __sub__(self, other):
"""Subtract polynomials."""
def __mul__(self, other):
"""Multiply polynomials."""
def __div__(self, other):
"""Divide polynomials."""
def __pow__(self, val):
"""Raise polynomial to a power."""
def deriv(self, m=1):
"""Return derivative of polynomial."""
def integ(self, m=1, k=None):
"""Return integral of polynomial."""
@property
def coeffs(self):
"""Polynomial coefficients."""
@property
def order(self):
"""Order (degree) of polynomial."""
@property
def roots(self):
"""Roots of polynomial."""
def poly(seq_of_zeros):
"""
Find the coefficients of a polynomial with the given sequence of roots.
Parameters:
seq_of_zeros: array_like - Sequence of polynomial roots
Returns:
ndarray: Polynomial coefficients
"""
def polyfit(x, y, deg, rcond=None, full=False, w=None, cov=False):
"""
Least squares polynomial fit.
Parameters:
x: array_like - x-coordinates of data points
y: array_like - y-coordinates of data points
deg: int - Degree of fitting polynomial
rcond: float, optional - Relative condition number for fit
full: bool, optional - Switch determining nature of return value
w: array_like, optional - Weights to apply to y-coordinates
cov: bool or str, optional - Return covariance matrix
Returns:
ndarray or tuple: Polynomial coefficients (and additional info if requested)
"""
def polyval(p, x):
"""
Evaluate a polynomial at specific values.
Parameters:
p: array_like - Polynomial coefficients in decreasing powers
x: array_like - Values at which to evaluate polynomial
Returns:
ndarray: Polynomial evaluated at x
"""
def polyadd(a1, a2):
"""
Add two polynomials.
Parameters:
a1, a2: array_like - Polynomial coefficients
Returns:
ndarray: Coefficients of sum polynomial
"""
def polysub(a1, a2):
"""
Subtract second polynomial from first.
Parameters:
a1, a2: array_like - Polynomial coefficients
Returns:
ndarray: Coefficients of difference polynomial
"""
def polymul(a1, a2):
"""
Multiply two polynomials.
Parameters:
a1, a2: array_like - Polynomial coefficients
Returns:
ndarray: Coefficients of product polynomial
"""
def polydiv(u, v):
"""
Divide first polynomial by second.
Parameters:
u, v: array_like - Polynomial coefficients (dividend and divisor)
Returns:
tuple: (quotient, remainder) polynomial coefficients
"""
def roots(p):
"""
Return the roots of a polynomial with coefficients given in p.
Parameters:
p: array_like - Polynomial coefficients
Returns:
ndarray: Roots of the polynomial
"""Modern Polynomial Package
The cupy.polynomial package provides a comprehensive set of polynomial classes and functions.
# Polynomial base classes and utilities
class Polynomial:
"""
Base class for all polynomial classes.
Provides common interface and functionality for polynomial
operations across different polynomial bases.
"""
def __init__(self, coef, domain=None, window=None): ...
def __call__(self, x): ...
def __add__(self, other): ...
def __sub__(self, other): ...
def __mul__(self, other): ...
def __truediv__(self, other): ...
def __pow__(self, other): ...
def deriv(self, m=1): ...
def integ(self, m=1, k=[]): ...
def roots(self): ...
@classmethod
def fit(cls, x, y, deg, domain=None, rcond=None, full=False, w=None, window=None): ...
@classmethod
def fromroots(cls, roots, domain=[], window=None): ...
# Power series polynomials (standard polynomials)
class Polynomial:
"""Power series polynomial in standard form."""
def polyval(x, c):
"""Evaluate polynomial at x."""
def polyval2d(x, y, c):
"""Evaluate 2-D polynomial at x, y."""
def polyval3d(x, y, z, c):
"""Evaluate 3-D polynomial at x, y, z."""
def polyvander(x, deg):
"""Generate Vandermonde matrix."""
def polyvander2d(x, y, deg):
"""Generate 2-D Vandermonde matrix."""
def polyvander3d(x, y, z, deg):
"""Generate 3-D Vandermonde matrix."""
def polyfit(x, y, deg, rcond=None, full=False, w=None):
"""Least squares fit of polynomial to data."""
def polycompanion(c):
"""Return companion matrix."""
def polyroots(c):
"""Return roots of polynomial."""
def polyfromroots(roots):
"""Generate polynomial from roots."""
# Chebyshev polynomials
class Chebyshev:
"""Chebyshev polynomial class."""
def chebval(x, c):
"""Evaluate Chebyshev series at x."""
def chebval2d(x, y, c):
"""Evaluate 2-D Chebyshev series."""
def chebval3d(x, y, z, c):
"""Evaluate 3-D Chebyshev series."""
def chebvander(x, deg):
"""Generate Chebyshev Vandermonde matrix."""
def chebfit(x, y, deg, rcond=None, full=False, w=None):
"""Fit Chebyshev series to data."""
def chebroots(c):
"""Return roots of Chebyshev polynomial."""
def chebfromroots(roots):
"""Generate Chebyshev polynomial from roots."""
# Legendre polynomials
class Legendre:
"""Legendre polynomial class."""
def legval(x, c):
"""Evaluate Legendre series at x."""
def legval2d(x, y, c):
"""Evaluate 2-D Legendre series."""
def legval3d(x, y, z, c):
"""Evaluate 3-D Legendre series."""
def legvander(x, deg):
"""Generate Legendre Vandermonde matrix."""
def legfit(x, y, deg, rcond=None, full=False, w=None):
"""Fit Legendre series to data."""
def legroots(c):
"""Return roots of Legendre polynomial."""
def legfromroots(roots):
"""Generate Legendre polynomial from roots."""
# Laguerre polynomials
class Laguerre:
"""Laguerre polynomial class."""
def lagval(x, c):
"""Evaluate Laguerre series at x."""
def lagval2d(x, y, c):
"""Evaluate 2-D Laguerre series."""
def lagval3d(x, y, z, c):
"""Evaluate 3-D Laguerre series."""
def lagvander(x, deg):
"""Generate Laguerre Vandermonde matrix."""
def lagfit(x, y, deg, rcond=None, full=False, w=None):
"""Fit Laguerre series to data."""
def lagroots(c):
"""Return roots of Laguerre polynomial."""
def lagfromroots(roots):
"""Generate Laguerre polynomial from roots."""
# Hermite polynomials (physicist's)
class HermiteE:
"""Hermite polynomial class (probabilist's)."""
def hermeval(x, c):
"""Evaluate Hermite series at x."""
def hermeval2d(x, y, c):
"""Evaluate 2-D Hermite series."""
def hermeval3d(x, y, z, c):
"""Evaluate 3-D Hermite series."""
def hermevander(x, deg):
"""Generate Hermite Vandermonde matrix."""
def hermefit(x, y, deg, rcond=None, full=False, w=None):
"""Fit Hermite series to data."""
def hermeroots(c):
"""Return roots of Hermite polynomial."""
def hermefromroots(roots):
"""Generate Hermite polynomial from roots."""
# Hermite polynomials (probabilist's)
class Hermite:
"""Hermite polynomial class (physicist's)."""
def hermval(x, c):
"""Evaluate Hermite series at x."""
def hermval2d(x, y, c):
"""Evaluate 2-D Hermite series."""
def hermval3d(x, y, z, c):
"""Evaluate 3-D Hermite series."""
def hermvander(x, deg):
"""Generate Hermite Vandermonde matrix."""
def hermfit(x, y, deg, rcond=None, full=False, w=None):
"""Fit Hermite series to data."""
def hermroots(c):
"""Return roots of Hermite polynomial."""
def hermfromroots(roots):
"""Generate Hermite polynomial from roots."""Usage Examples
Basic Polynomial Operations
import cupy as cp
import cupy.polynomial as poly
# Create and manipulate polynomials using poly1d
p1 = cp.poly1d([1, 2, 3]) # x^2 + 2x + 3
p2 = cp.poly1d([2, 1]) # 2x + 1
print("p1:", p1)
print("p2:", p2)
# Polynomial arithmetic
sum_poly = p1 + p2
diff_poly = p1 - p2
prod_poly = p1 * p2
print("Sum:", sum_poly)
print("Product:", prod_poly)
# Evaluate polynomials
x_values = cp.linspace(-5, 5, 100)
y1 = p1(x_values)
y2 = p2(x_values)
print("p1 evaluated at x=2:", p1(2))
print("p2 evaluated at x=2:", p2(2))
# Find roots
roots_p1 = p1.roots
print("Roots of p1:", roots_p1)
# Derivatives and integrals
p1_deriv = p1.deriv()
p1_integ = p1.integ()
print("Derivative of p1:", p1_deriv)
print("Integral of p1:", p1_integ)Polynomial Fitting
# Generate noisy data
x_data = cp.linspace(0, 10, 100)
true_coeffs = [0.5, -2, 1, 3] # 0.5x^3 - 2x^2 + x + 3
y_true = cp.polyval(true_coeffs, x_data)
noise = cp.random.normal(0, 0.5, x_data.shape)
y_noisy = y_true + noise
# Fit polynomial of different degrees
degrees = [2, 3, 4, 5]
fitted_polys = {}
for degree in degrees:
coeffs = cp.polyfit(x_data, y_noisy, degree)
fitted_polys[degree] = coeffs
# Evaluate fitted polynomial
y_fitted = cp.polyval(coeffs, x_data)
mse = cp.mean((y_noisy - y_fitted) ** 2)
print(f"Degree {degree}: MSE = {mse:.4f}")
print(f" Coefficients: {coeffs}")
# Best fit analysis
best_degree = min(fitted_polys.keys(),
key=lambda d: cp.mean((y_noisy - cp.polyval(fitted_polys[d], x_data)) ** 2))
print(f"Best degree: {best_degree}")
# Weighted fitting for heteroscedastic noise
weights = 1.0 / (1.0 + cp.abs(x_data)) # Less weight for larger x
weighted_coeffs = cp.polyfit(x_data, y_noisy, 3, w=weights)
print("Weighted fit coefficients:", weighted_coeffs)
# Robust fitting with full output
coeffs, residuals, rank, singular_vals, rcond = cp.polyfit(
x_data, y_noisy, 3, full=True
)
print("Fit diagnostics:")
print(f" Residuals: {residuals}")
print(f" Rank: {rank}")
print(f" Condition number: {1/rcond:.2e}")Advanced Polynomial Types
# Chebyshev polynomials for better numerical properties
x_cheb = cp.linspace(-1, 1, 100)
y_data = cp.exp(x_cheb) + 0.1 * cp.random.normal(size=x_cheb.shape)
# Fit using Chebyshev polynomials
cheb_coeffs = poly.chebyshev.chebfit(x_cheb, y_data, 5)
y_cheb_fitted = poly.chebyshev.chebval(x_cheb, cheb_coeffs)
print("Chebyshev coefficients:", cheb_coeffs)
print("Chebyshev fit MSE:", cp.mean((y_data - y_cheb_fitted) ** 2))
# Compare with standard polynomial fit
std_coeffs = cp.polyfit(x_cheb, y_data, 5)
y_std_fitted = cp.polyval(std_coeffs, x_cheb)
print("Standard polynomial MSE:", cp.mean((y_data - y_std_fitted) ** 2))
# Legendre polynomial approximation
leg_coeffs = poly.legendre.legfit(x_cheb, y_data, 5)
y_leg_fitted = poly.legendre.legval(x_cheb, leg_coeffs)
print("Legendre fit MSE:", cp.mean((y_data - y_leg_fitted) ** 2))
# Hermite polynomials for Gaussian-weighted data
x_herm = cp.linspace(-3, 3, 100)
gaussian_weight = cp.exp(-x_herm**2 / 2)
y_weighted = cp.sin(x_herm) * gaussian_weight + 0.05 * cp.random.normal(size=x_herm.shape)
herm_coeffs = poly.hermite_e.hermefit(x_herm, y_weighted, 6)
y_herm_fitted = poly.hermite_e.hermeval(x_herm, herm_coeffs)
print("Hermite fit MSE:", cp.mean((y_weighted - y_herm_fitted) ** 2))Multi-dimensional Polynomial Operations
# 2D polynomial fitting
x2d = cp.random.rand(200) * 4 - 2
y2d = cp.random.rand(200) * 4 - 2
z_true = 2*x2d**2 + 3*y2d**2 - x2d*y2d + x2d + 2*y2d + 1
z_noisy = z_true + 0.1 * cp.random.normal(size=z_true.shape)
# Fit 2D polynomial surface
degree_2d = [2, 2] # Degree in x and y directions
coeffs_2d = poly.polynomial.polyfit2d(x2d, y2d, z_noisy, degree_2d)
# Evaluate on regular grid
x_grid = cp.linspace(-2, 2, 50)
y_grid = cp.linspace(-2, 2, 50)
X_grid, Y_grid = cp.meshgrid(x_grid, y_grid)
Z_fitted = poly.polynomial.polyval2d(X_grid, Y_grid, coeffs_2d)
print("2D polynomial coefficients shape:", coeffs_2d.shape)
print("Fitted surface shape:", Z_fitted.shape)
# 3D polynomial evaluation
x3d = cp.linspace(-1, 1, 20)
y3d = cp.linspace(-1, 1, 20)
z3d = cp.linspace(-1, 1, 20)
X3d, Y3d, Z3d = cp.meshgrid(x3d, y3d, z3d, indexing='ij')
# Simple 3D polynomial coefficients
coeffs_3d = cp.zeros((3, 3, 3))
coeffs_3d[1, 0, 0] = 1 # x term
coeffs_3d[0, 1, 0] = 1 # y term
coeffs_3d[0, 0, 1] = 1 # z term
coeffs_3d[2, 0, 0] = 0.5 # x^2 term
# Evaluate 3D polynomial
result_3d = poly.polynomial.polyval3d(X3d, Y3d, Z3d, coeffs_3d)
print("3D polynomial result shape:", result_3d.shape)Root Finding and Analysis
# Complex polynomial root finding
# Polynomial: x^4 - 2x^3 + 3x^2 - 2x + 1
coeffs = [1, -2, 3, -2, 1]
roots = cp.roots(coeffs)
print("Polynomial coefficients:", coeffs)
print("Roots:", roots)
print("Root types:", [type(r).__name__ for r in roots])
# Verify roots by evaluation
for i, root in enumerate(roots):
value = cp.polyval(coeffs, root)
print(f"Root {i}: {root}, P(root) = {value}")
# Create polynomial from roots
reconstructed_coeffs = cp.poly(roots)
print("Reconstructed coefficients:", reconstructed_coeffs)
print("Original coefficients (normalized):", cp.array(coeffs) / coeffs[0])
# Companion matrix approach for root finding
companion = poly.polynomial.polycompanion([1, -2, 3, -2, 1])
eigenvals = cp.linalg.eigvals(companion)
print("Roots from companion matrix:", eigenvals)
# Analysis of polynomial behavior
def analyze_polynomial(coeffs, x_range=(-5, 5), n_points=1000):
"""Analyze polynomial behavior over a range."""
x = cp.linspace(x_range[0], x_range[1], n_points)
y = cp.polyval(coeffs, x)
# Find extrema (approximate)
dy = cp.gradient(y, x[1] - x[0])
extrema_indices = cp.where(cp.abs(dy) < cp.std(dy) * 0.1)[0]
analysis = {
'min_value': cp.min(y),
'max_value': cp.max(y),
'range': cp.max(y) - cp.min(y),
'approximate_extrema': len(extrema_indices),
'zeros_in_range': cp.sum(cp.abs(y) < cp.std(y) * 0.01)
}
return analysis
# Analyze different polynomials
test_polynomials = [
[1, 0, -1], # x^2 - 1
[1, -3, 3, -1], # (x-1)^3
[1, 0, 0, 0, -1], # x^4 - 1
[1, -2, 1] # (x-1)^2
]
for i, coeffs in enumerate(test_polynomials):
analysis = analyze_polynomial(coeffs)
print(f"\nPolynomial {i+1}: {coeffs}")
for key, value in analysis.items():
print(f" {key}: {value}")Specialized Polynomial Applications
# Interpolation using different polynomial bases
def compare_interpolation_methods(x_data, y_data, eval_points):
"""Compare different polynomial interpolation methods."""
n_points = len(x_data)
degree = n_points - 1
methods = {}
# Standard polynomial interpolation
std_coeffs = cp.polyfit(x_data, y_data, degree)
methods['Standard'] = cp.polyval(std_coeffs, eval_points)
# Chebyshev interpolation
cheb_coeffs = poly.chebyshev.chebfit(x_data, y_data, degree)
methods['Chebyshev'] = poly.chebyshev.chebval(eval_points, cheb_coeffs)
# Legendre interpolation
leg_coeffs = poly.legendre.legfit(x_data, y_data, degree)
methods['Legendre'] = poly.legendre.legval(eval_points, leg_coeffs)
return methods
# Test interpolation with Runge function (challenging for high-degree polynomials)
def runge_function(x):
return 1 / (1 + 25 * x**2)
# Interpolation points
n_interp = 10
x_interp = cp.linspace(-1, 1, n_interp)
y_interp = runge_function(x_interp)
# Evaluation points
x_eval = cp.linspace(-1, 1, 200)
y_true = runge_function(x_eval)
# Compare methods
methods = compare_interpolation_methods(x_interp, y_interp, x_eval)
print("Interpolation comparison (MSE vs true function):")
for method_name, y_approx in methods.items():
mse = cp.mean((y_true - y_approx) ** 2)
print(f" {method_name}: {mse:.6f}")
# Orthogonal polynomial series expansion
def orthogonal_series_approximation(func, x_range, n_terms, poly_type='chebyshev'):
"""Approximate a function using orthogonal polynomial series."""
x = cp.linspace(x_range[0], x_range[1], 1000)
y = func(x)
if poly_type == 'chebyshev':
coeffs = poly.chebyshev.chebfit(x, y, n_terms-1)
y_approx = poly.chebyshev.chebval(x, coeffs)
elif poly_type == 'legendre':
coeffs = poly.legendre.legfit(x, y, n_terms-1)
y_approx = poly.legendre.legval(x, coeffs)
elif poly_type == 'hermite':
coeffs = poly.hermite.hermfit(x, y, n_terms-1)
y_approx = poly.hermite.hermval(x, coeffs)
mse = cp.mean((y - y_approx) ** 2)
return coeffs, y_approx, mse
# Test with exponential function
exponential = lambda x: cp.exp(x)
x_range = (-2, 2)
for n_terms in [5, 10, 15, 20]:
for poly_type in ['chebyshev', 'legendre']:
coeffs, y_approx, mse = orthogonal_series_approximation(
exponential, x_range, n_terms, poly_type
)
print(f"{poly_type.capitalize()} ({n_terms} terms): MSE = {mse:.8f}")
# Polynomial differentiation and integration
p = cp.poly1d([1, -2, 1, 3]) # x^3 - 2x^2 + x + 3
print(f"Original polynomial: {p}")
# Multiple derivatives
for k in range(1, 4):
p_deriv = p.deriv(k)
print(f"Derivative {k}: {p_deriv}")
# Multiple integrals
for k in range(1, 3):
p_integ = p.integ(k)
print(f"Integral {k}: {p_integ}")
# Definite integration
x_vals = cp.linspace(0, 2, 100)
y_vals = p(x_vals)
numerical_integral = cp.trapz(y_vals, x_vals)
analytical_integral = p.integ()(2) - p.integ()(0)
print(f"Numerical integral (0 to 2): {numerical_integral}")
print(f"Analytical integral (0 to 2): {analytical_integral}")Performance Optimization
# Efficient polynomial evaluation using Horner's method
def horner_evaluation(coeffs, x):
"""Evaluate polynomial using Horner's method for numerical stability."""
result = cp.zeros_like(x, dtype=cp.result_type(coeffs.dtype, x.dtype))
for coeff in coeffs:
result = result * x + coeff
return result
# Compare with standard evaluation
coeffs_high_degree = cp.random.rand(50) # High-degree polynomial
x_test = cp.linspace(-10, 10, 10000)
# Time both methods
import time
start_time = time.time()
result_standard = cp.polyval(coeffs_high_degree, x_test)
standard_time = time.time() - start_time
start_time = time.time()
result_horner = horner_evaluation(coeffs_high_degree, x_test)
horner_time = time.time() - start_time
print(f"Standard evaluation time: {standard_time:.4f}s")
print(f"Horner evaluation time: {horner_time:.4f}s")
print(f"Results match: {cp.allclose(result_standard, result_horner)}")
# Vectorized polynomial operations for multiple polynomials
def batch_polynomial_evaluation(poly_coeffs_list, x_values):
"""Evaluate multiple polynomials efficiently."""
results = []
for coeffs in poly_coeffs_list:
results.append(cp.polyval(coeffs, x_values))
return cp.stack(results)
# Create batch of polynomials
n_polynomials = 100
poly_list = [cp.random.rand(cp.random.randint(3, 10)) for _ in range(n_polynomials)]
x_batch = cp.linspace(-5, 5, 1000)
# Batch evaluation
batch_results = batch_polynomial_evaluation(poly_list, x_batch)
print(f"Batch evaluation shape: {batch_results.shape}")
print(f"Memory usage: {batch_results.nbytes / 1024**2:.2f} MB")
# Memory-efficient polynomial fitting for large datasets
def streaming_polynomial_fit(data_generator, degree, max_batch_size=10000):
"""Fit polynomial to streaming data in batches."""
x_accumulator = []
y_accumulator = []
for x_batch, y_batch in data_generator:
x_accumulator.extend(x_batch)
y_accumulator.extend(y_batch)
# Process when batch is full
if len(x_accumulator) >= max_batch_size:
x_array = cp.array(x_accumulator)
y_array = cp.array(y_accumulator)
# Fit polynomial to current batch
coeffs = cp.polyfit(x_array, y_array, degree)
yield coeffs
# Clear accumulators
x_accumulator = []
y_accumulator = []
# Process remaining data
if x_accumulator:
x_array = cp.array(x_accumulator)
y_array = cp.array(y_accumulator)
coeffs = cp.polyfit(x_array, y_array, degree)
yield coeffs
# Example streaming data generator
def synthetic_data_generator(n_batches=10, batch_size=5000):
for i in range(n_batches):
x = cp.random.rand(batch_size) * 10 - 5
y = 2*x**3 - 3*x**2 + x + 1 + 0.1*cp.random.randn(batch_size)
yield cp.asnumpy(x), cp.asnumpy(y)
# Process streaming data
coeffs_list = list(streaming_polynomial_fit(synthetic_data_generator(), degree=3))
print(f"Processed {len(coeffs_list)} batches")
print("Average coefficients:", cp.mean(cp.stack(coeffs_list), axis=0))Polynomial operations in CuPy provide comprehensive mathematical tools for polynomial manipulation, fitting, evaluation, and analysis, supporting both legacy interfaces and modern orthogonal polynomial bases with GPU acceleration for high-performance computational mathematics and numerical analysis applications.
Install with Tessl CLI
npx tessl i tessl/pypi-cupy-cuda11x