如何加速基于Numpy数组的大规模敏感性分析？

Hey there! Your current loop-based approach is going to crawl when dealing with 1e9 data points—let's fix that with vectorized NumPy operations and some mathematical shortcuts. The key insight here is that we can derive an analytical formula for each sensitivity value instead of recalculating the function f for every single element tweak.

Step 1: Derive Sensitivity Formulas #

First, let's formalize your sensitivity definition: for each element $a_{mn}$, the sensitivity $e_{mn}$ is the absolute value of $f(\text{tweaked } a_{mn}) - b_m$, where the tweak is $a_{mn} \rightarrow a_{mn}*(1-i)$.

Given your function:
$$f(a_{m1},a_{m2},a_{m3},a_{m4}) = a_{m1} - a_{m2} - (a_{m3} * a_{m4})$$

We can compute each $e_{mn}$ directly:

• For $n=1$ (first attribute): Tweaking $a_{m1}$ reduces the function value by $i*a_{m1}$, so $e_{m1} = i * |a_{m1}|$
• For $n=2$ (second attribute): Tweaking $a_{m2}$ increases the function value by $i*a_{m2}$, so $e_{m2} = i * |a_{m2}|$
• For $n=3$ (third attribute): Tweaking $a_{m3}$ increases the function value by $i*a_{m3}*a_{m4}$, so $e_{m3} = i * |a_{m3} * a_{m4}|$
• For $n=4$ (fourth attribute): Tweaking $a_{m4}$ has the same effect as tweaking $a_{m3}$, so $e_{m4} = i * |a_{m3} * a_{m4}|$

This matches exactly with your sample output—check it:

• For the first row ($a=[2,2,100,0.02]$): $0.052=0.1$, $0.052=0.1$, $0.05*(1000.02)=0.1$, $0.05(100*0.02)=0.1$
• For the second row ($a=[4,2,100,0.02]$): $0.054=0.2$, $0.052=0.1$, $0.05*(1000.02)=0.1$, $0.05(100*0.02)=0.1$

Step 2: Vectorized Implementation #

Now we can compute E in a single pass with no loops. Here's the code:

import numpy as np

# Sample input
A = np.array([ [2., 2., 100., 0.02], [4., 2., 100., 0.02] ])
i = 0.05

# Compute base function values B (vectorized, same result as your original)
B = (A[:, 0] - A[:, 1] - A[:, 2] * A[:, 3]).reshape(-1, 1)

# Compute sensitivity matrix E
E = np.zeros_like(A)
# Column 1 sensitivity
E[:, 0] = i * np.abs(A[:, 0])
# Column 2 sensitivity
E[:, 1] = i * np.abs(A[:, 1])
# Columns 3 and 4 sensitivity
prod_34 = A[:, 2] * A[:, 3]
E[:, 2] = i * np.abs(prod_34)
E[:, 3] = i * np.abs(prod_34)

# Verify output matches your sample
print(E)
# Output:
# [[0.1 0.1 0.1 0.1]
#  [0.2 0.1 0.1 0.1]]

Step 3: Handling 1e9-Scale Data #

For M=1e9, you might run into memory issues if you try to store the entire array in RAM at once. Here's how to handle it:

• Chunked Processing: Split A into smaller chunks (e.g., 1e6 rows per chunk), compute E for each chunk, and write the results to disk (using np.savez or a library like zarr) instead of keeping everything in memory.
• Out-of-Core Libraries: Tools like Dask can handle larger-than-memory arrays by automatically chunking and parallelizing operations, while maintaining a NumPy-like API.

Why This Is Way Faster #

• Your original loop-based approach is $O(M*N)$ with high constant overhead (array copies, function calls per element).
• The vectorized approach is $O(M)$ with minimal overhead, leveraging NumPy's optimized C backend. For large M, this will be orders of magnitude faster—we're talking seconds instead of hours/days.

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