迭代式支配曲线求解性能劣化原因及优化方案问询

Alright, let's break down why your evaluation time is spiking with each iteration and fix that recursion depth error at the same time.

Root Causes of the Performance Degradation #

1. Nested Piecewise Function Bloat
   Your early recursion depth error is a dead giveaway: each iteration's dominating curve is being built as a nested piecewise function that wraps the previous iteration's curve. When you evaluate a point on this nested structure, it has to traverse all previous layers of piecewise checks to get to the actual function calculation. By iteration 4, that's 4 layers of nested conditionals; by iteration N, it's N layers—leading to exponential growth in evaluation time.

2. Redundant Integral Calculations in evalv
   If your evalv method recalculates Gauss-Hermite nodes/weights every time it's called, or doesn't cache results for repeated m/s parameter pairs, you're wasting cycles on identical computations. This redundancy gets amplified as the nested dominating curve forces more calls to evalv.

3. Unoptimized Dominating Curve Evaluation
   When starting from the rightmost x-axis point each iteration, if you're re-evaluating all three curves (red/blue/green) from scratch every time without caching, you're doing far more work than necessary—especially as those curves become more complex with each iteration.

Optimized Implementation Plan #

Let's fix these issues step by step with concrete code changes:

1. Flatten the Dominating Curve Representation #

Ditch nested piecewise functions entirely. Instead, represent the dominating curve as a flat list of intervals and their corresponding functions. This lets you evaluate points with a simple binary search (no recursion, no nested checks).

Implement a custom class to handle this cleanly:

import bisect

class DominatingCurve:
    def __init__(self, breakpoints, functions):
        # breakpoints: sorted list of x-values (from right to left, e.g., [10, 5, 0])
        # functions: list of functions corresponding to intervals [breakpoints[i], breakpoints[i+1]]
        self.breakpoints = breakpoints
        self.functions = functions
    
    def __call__(self, x):
        # Find the index of the first breakpoint <= x (using bisect for speed)
        idx = bisect.bisect_right(self.breakpoints, x) - 1
        # Handle edge cases (x beyond the rightmost/leftmost breakpoint)
        if idx < 0:
            return self.functions[0](x)
        elif idx >= len(self.functions):
            return self.functions[-1](x)
        return self.functions[idx](x)

When building the curve in find_all_intersections, instead of nesting piecewise functions, collect the breakpoints and their associated functions directly into this class.

2. Cache & Precompute Integral Components in evalv #

• Precompute Gauss-Hermite nodes/weights once: Don't regenerate these every time evalv runs—they're fixed for a given number of nodes.
• Cache evalv results: Use functools.lru_cache to store results for repeated m/s pairs (ensure parameters are hashable, e.g., use floats instead of numpy scalars if needed).

Example:

from functools import lru_cache
import numpy as np

# Precompute Gauss-Hermite nodes/weights once (adjust n for your precision needs)
GAUSS_HERMITE_NODES, GAUSS_HERMITE_WEIGHTS = np.polynomial.hermite.hermgauss(20)

@lru_cache(maxsize=None)
def evalv(m: float, s: float) -> float:
    # Scale nodes to match your normal distribution parameters
    scaled_nodes = m + s * GAUSS_HERMITE_NODES
    # Compute your integrand for all nodes at once (vectorized for speed)
    integrand = ...  # Replace with your actual integrand logic, using scaled_nodes
    # Calculate integral using precomputed weights
    integral = np.sum(GAUSS_HERMITE_WEIGHTS * integrand)
    # Apply Gauss-Hermite scaling factor
    return integral * np.sqrt(np.pi)

3. Optimize Intersection Finding & Dominance Checks #

• Track current dominant curve from right to left: Instead of re-evaluating all three curves at every point, start at the rightmost x, note which curve is dominant, then move left only until you hit an intersection with another curve. Update the dominant curve at that breakpoint and repeat.
• Avoid redundant curve evaluations: Cache the values of each curve at breakpoints so you don't recalculate them when checking dominance.

4. Eliminate Recursion in Curve Generation #

Instead of passing the previous dominating curve as a nested function into the next iteration's _function_template, pass the DominatingCurve instance directly. Your template function can call the instance like any other function (thanks to the __call__ method), without nesting.

Why This Works #

• No more recursion depth errors: The flat DominatingCurve class uses no recursion—evaluation is a simple index lookup.
• Linear evaluation time: Each point evaluation now takes O(log N) time (for binary search over breakpoints), instead of O(N) time for nested conditionals.
• Reduced redundant computation: Precomputing nodes/weights and caching evalv results cuts down on repeated work that was eating up your runtime.

内容的提问来源于stack exchange，提问作者math