Hopefully this question is on topic for the statistics community (rather than the pure math) since the optimization method I am using relies on a statistical model. Well I am building an algorithm for constrained optimization, so a problem like, $$\min_x f(x)$$ such that $$c(x)\leq0$$ where $f(x)$ is scalar valued. Why this is a statistics problem, is that I don't know the analytical forms of $f(\cdot)$ or $c(\cdot)$ but I do have some inputs $x$ and their corresponding outputs $f(x)$ and $c(x)$. So you could kind of think of this as a regression problem where I know some inputs and outputs and I want to fit a response surface for $f(\cdot)$ and $c(\cdot)$ and try to minimize that true function based on the predictive response surface I have built (hopefully I am not omitting too many details). So its an iterative procedure where I gather my data (the inputs and outputs) in time and so far empirically the procedure works. However, now I want to prove theoretically that my algorithm will converge to the (or possibly a local) solution. And there is where I am stuck. I don't really have too much intuition about how to prove this so any suggestions or citations would be very much appreciated.
In one of your comments you state that you not only have values of $f(x)$ and $c(x)$ at various $x$, but can also evaluate $f(x)$ and $c(x)$ at new points. This means you are essentially doing ordinary numerical optimization. The lack of an analytical form for $f(x)$ and $c(x)$ is not an issue for most optimization algorithms, because they wouldn't use it even if they had it. They generally construct an easily optimized surrogate based on values and derivatives from the last few points, optimize the surrogate to get a new candidate optimum, and iterate. The only requirement for such a procedure is the ability to evaluate $f$ and $c$ and perhaps some derivatives.
All that is needed for convergence to a local optimum is some mild smoothness/regularity conditions on $f$ and $c$. Global convergence is hard to prove unless you're doing (quasi-)convex optimization.***
If you are able to compute gradients or Hessians of $f$ and $c$, there is an extensive classical literature on optimization algorithms that provably converge to local optima. Numerical Optimization by Nocedal and Wright is a good introduction. This literature is the foundation of all mathematical optimization theory, even if gradients and Hessians are not available in your problem.
If you cannot compute gradients or Hessians, you are doing derivative-free optimization (DFO). There are convergence results for unconstrained DFO, and for DFO with linear and bound constraints. For DFO with general nonlinear constraints I do not know of any proof results. There is the COBYLA algorithm by Michael Powell, one of the most productive researchers in DFO, but no proofs to my knowledge. This paper and this book are written by top DFO researchers and should help you get acquainted with the literature on what has been proven and how.
***It's always worth thinking for a moment about whether your problem has a quasiconvex flavor, because the benefits are potentially enormous. If you have two points $x_1,x_2$ and their costs $f(x_1), f(x_2)$, and you make a weighted average of the points $x = \alpha x_1 + (1-\alpha) x_2$ for $\alpha \in [0,1]$. Is the cost $f(x) \leq \max(f(x_1),f(x_2))$? Put another way, is a weighted average at least as good as the worst component that went into it?
You are not going to get very far without some assumptions about $f$ and $c$. The standard assumption is Lipschitz continuity, which guarantees that the function cannot be too badly behaved between evaluation points. Without Lipschitz continuity you cannot say much, because even in the one-dimensional case over an interval you cannot rule out the possibility of an arbitrarily deep hole between any two points where the function has been evaluated.
Can you sample the function values iteratively? You do not say so explicitly, and it does make a difference. If you can sample, there are a number of derivative-free optimization algorithms, such as COBYLA.
If you cannot sample, the best you can hope for, even with Lipschitz continuity, is a confidence interval for the optimal objective value, with very limited information about the location of the optimum.
Possibly useful things to think about: 1: Can you prove your estimate of the parameters in f and c is consistent? Do they always converge to the true values (assuming the form is correct) as the sample size gets bigger? For example, maximum likelihood estimators are consistent. 2: Can you prove the control variable x converges to the the optimal value? Something like the Karush-Kuhn-Tucker conditions might be useful.
There is a paper reviewing your kind of problem: A Taxonomy of Global optimization Methods Based on Response Surfaces. I skimmed it once but don't remember exactly what they have and whether or not they consider constrained optimization, but it's an OK place to start. Basically, the method is to model the function using Gaussian Process Regression and then take an optimization step based on the interpolant. This may or may not fit your line of thought. Let me know if you don't have access to the paper and I'll send it to you.
More broadly, what you might be interested in is derivative-free optimization (or so it seems). For this there are many methods that I am not familiar with and even software packages (I'm pretty sure pyton's scipy has some flags in its "minimize" funciton that allow you to do derivative free optimization). You can google "derivative-free optimization" and see what you find. Also, one specific method that I like (though never used) is adaptive coordinate descent.
If you have a method that you like and you want to use it for constrained optimization, one of the first things I would try is to take optimization steps according to the method but after that, project to the subspace / manifold that satisfies the constraint.