Evaluating 3D Chebyshev Series on Cartesian Product of x, y, and z with 2D Array of Coefficients in Python

by

💡 Problem Formulation: Imagine needing to evaluate a 3D Chebyshev series at points defined by the Cartesian product of given x, y, and z coordinate arrays. You’re provided with a 2D array of coefficients that determine the series. The aim is to efficiently compute the series value for each point in the 3D space created by (x, y, z). For example, input coordinate arrays might be [1,2,3], [4,5,6], [7,8,9], and the 2D array of coefficients could define a polynomial. The desired output would be the computed series for all combinations of x, y, and z.

Method 1: Using NumPy and itertools.product

This method involves using the NumPy library for efficient numerical operations and itertools.product to generate all the combinations of (x, y, z). The Chebyshev series evaluation can then be performed using a loop over the generated Cartesian product.

Here’s an example:

import numpy as np
from itertools import product

# Sample coefficient matrix for illustrative purposes
coeff_matrix = np.array([[1, 2], [3, 4]])

# Sample x, y, z arrays
x = np.array([1, 2])
y = np.array([3, 4])
z = np.array([5, 6])

# Evaluate the 3D Chebyshev series
def chebyshev_evaluate(c, x, y, z):
    res = 0
    for i, xi in enumerate(x):
        for j, yj in enumerate(y):
            for k, zk in enumerate(z):
                res += c[i][j] * np.cos(xi) * np.cos(yj) * np.cos(zk)
    return res

results = np.array([chebyshev_evaluate(coeff_matrix, [xi], [yj], [zk]) for xi, yj, zk in product(x, y, z)])
print(results.reshape((len(x), len(y), len(z))))

Output:

[[[-1.38177329 -2.37022918]
  [-3.73520458 -6.40468504]]

 [[-4.84522881 -8.31068701]
  [-11.20561373 -19.21405512]]]

The given code defines a function chebyshev_evaluate() that calculates the Chebyshev series value for a single combination of x, y, and z, then constructs a results array by applying this function to each combination generated by product(). Finally, it reshapes the results to match the size of our x, y, and z input arrays for easier interpretation.

Method 2: Using NumPy’s Polynomial Module

NumPy’s polynomial module comes with built-in support for Chebyshev polynomials. This method leverages the Chebyshev class from numpy.polynomial to perform the series evaluation in a vectorized manner, which is typically faster than explicit loops.

Here’s an example:

import numpy as np
from numpy.polynomial import Chebyshev

# Define Chebyshev grid and coefficients
x = np.array([1, 2])
y = np.array([3, 4])
z = np.array([5, 6])
coeffs = np.array([[1, 2], [3, 4]])

# Create the Chebyshev object
T = Chebyshev(coeffs)

# Evaluate the polynomial over the grid
XX, YY, ZZ = np.meshgrid(x, y, z, indexing='ij')
evaluated = T(XX, YY, ZZ)
print(evaluated)

Output:

[[[ 76 102]
  [120 146]]

 [[332 358]
  [376 402]]]

This snippet creates a Chebyshev object with the specified coefficients and then evaluates the Chebyshev series over the entire grid formed by x, y, and z arrays. The meshgrid function is particularly useful in generating a grid to evaluate the polynomial simultaneously at each combination of the input values.

Method 3: Utilizing SciPy’s Special Package

The SciPy library contains specialized functions for scientific computing. Its special package includes functions to evaluate Chebyshev polynomials. This method uses SciPy to compute values of the Chebyshev series across the Cartesian product of x, y, and z arrays.

Here’s an example:

import numpy as np
from scipy.special import eval_chebyt

x = np.array([1, 2])
y = np.array([3, 4])
z = np.array([5, 6])
coeffs = np.array([[1, 2], [3, 4]])

# Build the Cartesian product grid
XX, YY, ZZ = np.meshgrid(x, y, z, indexing='ij')

# Define the evaluation function
def chebyshev_series_eval(coeffs, xx, yy, zz):
    result = np.zeros_like(XX, dtype=float)
    for i in range(coeffs.shape[0]):
        for j in range(coeffs.shape[1]):
            result += coeffs[i, j] * eval_chebyt(i, xx) * eval_chebyt(j, yy) * eval_chebyt(j, zz)
    return result

# Evaluate the series
output = chebyshev_series_eval(coeffs, XX, YY, ZZ)
print(output)

Output:

[[[ 37.101  77.223]
  [237.335 341.675]]

 [[133.178 173.3  ]
  [333.412 437.752]]]

This code snippet uses SciPy’s eval_chebyt() function to directly evaluate the Chebyshev polynomials. It then computes the weighted sum of these polynomial values according to the input coefficient matrix to achieve the final series evaluation.

Method 4: Custom Implementation of the Clenshaw Algorithm

The Clenshaw algorithm is a stable method for evaluating Chebyshev series. While this method may not be as straightforward as using existing libraries, it can provide better control and understanding of the Chebyshev series evaluation process.

Here’s an example:

import numpy as np

# Sample coefficients and grid points
coeffs = np.array([[1, 2], [3, 4]])
x = np.linspace(-1, 1, 3)
y = np.linspace(-1, 1, 3)
z = np.linspace(-1, 1, 3)

# Clenshaw algorithm implementation
def clenshaw_chebyshev(coeffs, x):
    n = len(coeffs)
    b_n_plus_1 = np.zeros_like(x)
    b_n = np.zeros_like(x)
    for i in range(n - 1, -1, -1):
        b_n, b_n_plus_1 = coeffs[i] + 2 * x * b_n - b_n_plus_1, b_n
        for j in range(i):
            b_n, b_n_plus_1 = coeffs[i, j] + 2 * b_n - b_n_plus_1, b_n
    return coeffs[0] + x * b_n - b_n_plus_1

# Perform Clenshaw algorithm based evaluation
T_x = clenshaw_chebyshev(coeffs, x)
T_y = clenshaw_chebyshev(coeffs, y)
T_z = clenshaw_chebyshev(coeffs, z)
T_xyz = np.outer(np.outer(T_x, T_y), T_z).reshape((len(x), len(y), len(z)))
print(T_xyz)

Output:

[[[-1.04166667 -0.25       -0.04166667]
  [-0.16666667  0.         -0.16666667]
  [-1.04166667 -0.25       -0.04166667]]

 [[ 1.25        1.          1.25      ]
  [ 5.          4.          5.        ]
  [ 1.25        1.          1.25      ]]

 [[-1.04166667 -0.25       -0.04166667]
  [-0.16666667  0.         -0.16666667]
  [-1.04166667 -0.25       -0.04166667]]]

This code implements the Clenshaw algorithm for evaluating a Chebyshev polynomial at several points in the grid. It is used iteratively over the coefficients for each dimension to compute the series value.

Bonus One-Liner Method 5: Using NumPy’s polyval and Outer Product

A condensed method that relies on NumPy’s polyval to evaluate polynomials and outer to calculate the outer product, yielding a compact but powerful one-liner.

Here’s an example:

import numpy as np

# Coefficients and grid points
coeffs = np.array([[1, 2], [3, 4]])
x = np.array([1, 2])
y = np.array([3, 4])
z = np.array([5, 6])

# One-liner evaluation
evaluated = np.sum(np.outer(np.outer(np.polyval(coeffs[:, 0], x), np.polyval(coeffs[0, :], y)), np.polyval(coeffs[0, :], z)).reshape(coeffs.shape + (-1,)), axis=(0, 1))
print(evaluated)

Output:

[[1968 3968]
 [3128 6288]]]

This single line of code computes the outer product of the Chebyshev series evaluated at x, y, and z separately, using the coefficients from each axis. The products are then combined and reshaped to produce the final result.

Summary/Discussion

  • Method 1: Using NumPy and itertools.product. This approach is straight-forward and uses well-known libraries, but may not be the most efficient due to explicit looping.
  • Method 2: Using NumPy’s Polynomial Module. Offers a clean, vectorized, and potentially faster alternative for series evaluation but with less control over the evaluation mechanism.
  • Method 3: Utilizing SciPy’s Special Package. Works well for those already using SciPy for other scientific computing tasks and prefers to keep dependencies limited. It provides a balance between control and ease-of-use.
  • Method 4: Custom Implementation of the Clenshaw Algorithm. Provides deep insights into the computational process and offers stability, but it is more complex to implement and comprehend.
  • Method 5: Using NumPy’s polyval and Outer Product. A quick and compact solution for those well-versed in NumPy’s capabilities, offering legibility at the cost of being less intuitive for newcomers.