If your optimization problem is black-box, then you can use optimization solvers such as COBYLA (https://handwiki.org/wiki/COBYLA), Nomad (https://www.gerad.ca/fr/software/nomad), EGO (https://smt.readthedocs.io/en/latest/_src_docs/applications/ego.html) or CMAES (https://cma-es.github.io). The regularization term can be then added to the objective function of your nonlinear optimization problem.

If your optimization problem is a least-squares one, then you can use the Levenberg-Marquardt (LM) algorithm (i.e., a regularized version of the Gauss-Newton method). The LM method is regularization optimization method typically used to solve inverse problems.

Youssef Diouane Thank you for your answer!

My optimization problem is not least-square. The objective function is defined, then we have not used the gradient method for solving it. we used Brent's Method. I mean we have such a problem

arg E(x)+Reg(x)

We solve arg E(x) by Brent's method. My question is now how I can solve it when I add regularization such as Tikhonov regularization. since the method that I used for solving it is not based on a gradient base. what approach should I follow?

I have attached a C code and also an Excel spreadsheet that both implement Brent's method. There is no derivative involved. What do you mean by "regularization"? Are you only considering discrete values of x? If so, that is easily done. For example if dx=0.1, x=0.0, 0.1, 0.2, 0.3,... then simply make x3=10*(int)(x/10) My book Numerical Methods will be free this Saturday (4/30). The software is always free. https://www.amazon.com/dp/B09HHFDRZB

Dudley J Benton Thank you so much for your answer!

I am working on image processing. Hence my data is discrete.

The very common method is to solve a convex optimization problem such as

minimizeX∈R^m×n ‖A(X)−B‖^2+γ‖Dx‖

where D is a discrete gradient operator (in which case, for example, we're using "total variation regularization")

For the problem that I am going to solve, I consider ssim as the norm instead of the l2 norm. When γ‖Dx‖ is not considered, it is solved by Brent's method.

I do not know how it is possible. I consider γ‖Dx‖ and solve the problem since we usually solve it by taking the derivative of the whole problem that we defined. Can you guide me regarding it?

Orthogonal functions are often used for image processing (e.g., FFT, Radon transform, sinogram, SVD used in transmitting images from deep space probes). There are lots of examples in my book, which will be free on 5/10. I have attached a zip file containing all of the examples. https://www.amazon.com/dp/B07GT8TLDV