基于Numpy/Scipy的任意维度输入输出向量值多元函数近似问询

Great question! When dealing with expensive vector-valued functions (n-dimensional input, m-dimensional output) in NumPy/SciPy, you have a few solid options for approximation or interpolation that work for arbitrary n and m. Let's walk through them:

1. Taylor Series Approximation (Local) #

SciPy doesn't have a built-in universal Taylor expansion function, but you can easily construct one using numerical derivative tools like scipy.optimize.approx_fprime to compute gradients/Jacobians. This is ideal for local approximation near a specific point (e.g., x0).

For a first-order Taylor approximation (linearization), the formula is:

( \hat{f}(x) = f(x_0) + J(x_0) \cdot (x - x_0) )
where ( J(x_0) ) is the Jacobian matrix of ( f ) at ( x_0 )

Here's a concrete implementation:

import numpy as np
from scipy.optimize import approx_fprime

def your_vector_func(x):
    # Example: 2D input → 3D output; replace with your actual function
    return np.array([x[0]**2 + x[1], np.sin(x[0]+x[1]), np.exp(x[0])])

def taylor_1st_order(x, x0, func):
    f0 = func(x0)
    # Compute Jacobian: each row is the gradient of one output component
    jacobian = np.array([
        approx_fprime(x0, lambda x: func(x)[i], epsilon=1e-6) 
        for i in range(len(f0))
    ])
    return f0 + jacobian @ (x - x0)

# Usage example
x0 = np.array([1.0, 0.5])  # Reference point
x_query = np.array([1.1, 0.6])  # Point to approximate
approx_result = taylor_1st_order(x_query, x0, your_vector_func)
print("First-order Taylor approximation:", approx_result)

• Adjust epsilon to balance numerical precision and computation cost.
• Extend to higher-order approximations (e.g., quadratic) by computing Hessian tensors, though this gets more complex for high n/m.

2. Radial Basis Function (RBF) Interpolation (Global/Local) #

If you have a set of sampled points ((x_i, f(x_i))), scipy.interpolate.RBFInterpolator is your best bet—it's a direct generalization of interpolation for arbitrary input/output dimensions, and behaves like a more flexible version of scipy.interpolate.approximate.

It works by fitting a radial basis function to your sampled data, allowing you to approximate the function at new points across your domain.

from scipy.interpolate import RBFInterpolator
import numpy as np

# Generate sample data (replace with your actual sampled points)
np.random.seed(42)
n_input = 2
m_output = 3
sample_x = np.random.uniform(-2, 2, (50, n_input))
sample_f = np.array([
    sample_x[:,0]**2 + sample_x[:,1], 
    np.sin(sample_x[:,0]+sample_x[:,1]), 
    np.exp(sample_x[:,0])
]).T

# Build the interpolator
rbf_interp = RBFInterpolator(
    sample_x, 
    sample_f, 
    kernel='gaussian',  # Try 'linear', 'cubic', or 'multiquadric' too
    epsilon=1.0  # Controls smoothness; adjust based on your data
)

# Approximate at new points
query_x = np.array([[1.1, 0.6], [0.2, -0.3]])
approx_results = rbf_interp(query_x)
print("RBF interpolation results:\n", approx_results)

• Works for any n (input dim) and m (output dim) out of the box.
• Choose the kernel based on your function's smoothness: linear for piecewise linear fits, gaussian for smooth approximations.

3. Gaussian Process Regression (GPR) #

For scenarios where you need uncertainty estimates or have limited sampling data, Gaussian Process Regression (via scipy.stats.gaussian_process) is a powerful option. While SciPy's GPR implementation is single-output by default, you can easily extend it to m-dimensional outputs by training a separate model for each output component.

from scipy.stats import gaussian_process
import numpy as np

def your_vector_func(x):
    return np.array([x[0]**2 + x[1], np.sin(x[0]+x[1]), np.exp(x[0])])

# Generate sample data
np.random.seed(42)
sample_x = np.random.uniform(-2, 2, (30, 2))
sample_f = your_vector_func(sample_x.T).T  # Shape: (30, 3)

# Train a GPR model for each output component
gpr_models = []
for output_idx in range(sample_f.shape[1]):
    gp = gaussian_process.GaussianProcess(
        theta0=1e-2, 
        thetaL=1e-4, 
        thetaU=1e-1
    )
    gp.fit(sample_x, sample_f[:, output_idx])
    gpr_models.append(gp)

# Predict at a query point
query_x = np.array([1.1, 0.6])
approx_result = np.array([
    model.predict([query_x])[0] 
    for model in gpr_models
])
print("GPR approximation result:", approx_result)

• Provides uncertainty bounds alongside predictions (use gp.predict(..., eval_MSE=True) to get variance).
• Computation cost grows with the number of samples, so it's best for smaller datasets.

• Local approximation near a point: Use Taylor series with numerical derivatives.
• General interpolation across a domain: Use RBFInterpolator (most similar to the generalized approximate you're looking for).
• Uncertainty-aware approximation: Use Gaussian Process Regression.

内容的提问来源于stack exchange，提问作者Niklas M.