Finding the input variable to which output is most sensitive I have a physical system which takes a number of inputs $x_i$ and produces an output $error$.
$$
 Y = f(x_1, x_2, x_3, .. x_{1000})
$$
The function $f()$ can be evaluated by running a compute-intensive simulation of a model.
I want to find the $x_i$ to which $Y$ is most sensitive. In practice, I want to optimize the values for the few input variable which would give maximum return (in terms improving the system performance).
I can think of randomly changing each of the $x_i$ by a small amount around the existing value and record the output. Repeat the experiment by few hundred times and compute the correlation between $x_i$ to $Y$ and pick the inputs with high correlation.
I am wondering if there is a more formal method to achieve this.
One important constraint in my particular problem is that each model evaluation requires a computationally intensive simulation of about 10 minutes and $x_i$ is of size $1000$ to $2000$.
 A: The function is basically a surface in your 100 dimensional space or so. And at any given point in the space, it's more sensitive to the inputs for which it has the highest value on the partial derivatives for that input.
That means, your function may be sensitive to one input at one point in your hight dimensional space, and sensitive to another input in another place in your space.
Mathematically speaking, you should calculate the partial derivatives of your function for each input, and then analyze those derivatives to see which one is usually higher than the others, if that's what you're interested in.
In other words, it may be the case that your function is very sensitive to $x_1$ if $x_2 = 5$, and extremely sensitive to $x_{10}$ if $x_2 = 100$.
You can calculate those derivatives either analytically if you have the closed form of the function, or estimate those derivatives at randomly given input points using numerical algorithms.
In practice, you usually don't have the analytical form of the function. What you can do is to randomly choose from your data set (an iid sample), and numerically estimate the derivatives for each input feature. With a large enough sample, you get a nice representation of the distribution of those partial derivatives, and then you can decide what to do with them depending on your calculated value sets (one set per input feature).
It may very well be the case that for some you get a nice gaussian distribution, for some you may see white noise, and for some you see a rather uniform distribution, etc. It really depends on your data AND your function.
A: If your output is linearly related to the inputs then, you can try the Regression Analysis in context of sensitivity. We can calculate the standardized beta coefficient as the direct measure of sensitivity. because in the calculation of standardized beta coefficients, every input is firstly standardized w.r.t mean and standard deviation. and it is very simple to interpret.
for reference, 
https://www.listendata.com/2015/04/standardized-vs-unstandardized.html
https://en.wikipedia.org/wiki/Standardized_coefficient
A: I think that you could look into the field of Sensitivity Analysis : https://en.wikipedia.org/wiki/Sensitivity_analysis.
In your case, I would advise to compute the Sobol' indices (https://en.wikipedia.org/wiki/Variance-based_sensitivity_analysis).
These indices represent the fraction of variance carried by a variable and/or a set of variables. Several R packages exist in order to compute first and second order indices quite efficiently, by using specific designs. 
In your case, as the number of model evaluation is pretty small and the number of inputs is large, you could try to look into surrogate based sensitivty analysis (see for instance https://doi.org/10.1016/j.apm.2013.01.019):
Take a well behaved initial design (Latin Hypersquare or other space filling designs), and based on these evaluations, build a surrogate model (using Kriging). This surrogate will then be used for intensive computations, and can give some insightful results.
Be aware however that due to the high number of inputs, an accurate surrogate will probably need a lot of initial runs to be generated. A usual rule of thumb is to take $10d$ initial design points, where $d$ is the input dimension.
