6

I am wondering if there is a C/C++ library or Matlab code technique to determine real and complex numbers using a minimization solver. Here is a code snippet showing what I would like to do. For example, suppose that I know Utilde, but not x and U variables. I want to use optimization (fminsearch) to determine x and U, given Utilde. Note that Utilde is a complex number. 

x = 1.5;
U = 50 + 1i*25;
x0 = [1 20];  % starting values
Utilde = U * (1 / exp(2 * x)) * exp( 1i * 2 * x);
xout = fminsearch(@(v)optim(v, Utilde), x0);

function diff = optim(v, Utilde)
x = v(1);
U = v(2);
diff =  abs( -(Utilde/U) + (1 / exp(2 * x)) * exp( 1i * 2 * x  ) );

The code above does not converge to the proper values, and xout = 1.7318 88.8760. However, if U = 50, which is not a complex number, then xout = 1.5000 50.0000, which are the proper values.

Is there a way in Matlab or C/C++ to ensure proper convergence, given Utilde as a complex number? Maybe I have to change the code above?

  • If there isn't a way to do this natively in Matlab, then perhaps one gist of the question is this: Is there a multivariate (i.e. Nelder-Mead or similar algorithm) optimization library that is able to work with real and complex inputs and outputs?

  • Yet another question is whether the function is convergent or not. I don't know if it is the algorithm or the function. Might I need to change something in the Utilde = U * (1 / exp(2 * x)) * exp( 1i * 2 * x) expression to make it convergent?

Nicholas Kinar
  • 1,440
  • 5
  • 24
  • 36
  • From my experience: using these built-in minimization procedures often gives you more headaches than it helps. If you definitely need to do it this way, I would stick to Python - with MATLAB it will probably not be better. – Michael Schlottke-Lakemper Jul 13 '12 at 17:16
  • Sure - what is the best way in Python to set up this optimization problem? Is there a tool for Python that can optimize using complex numbers? – Nicholas Kinar Jul 13 '12 at 23:06
  • @NicholasKinar I am a little bit uncertain about the arithmetic rules for complex numbers right now, but if you only want to retrieve the topmost `x` and `U` values in the optimization, would it not be more appropriate to specify `diff` as `diff = abs( Utilde - U * (1 / exp(2 * x)) * exp( 1i * 2 * x ) )`? Or better still from a differentiation point of view the difference squared instead of the absolute difference? – Anders Gustafsson Jul 24 '12 at 21:07
  • @AndersGustafsson: Thanks for your comment. Hmm...I tried this, and I still can't reach convergence for all `x` and `U`. For example, `x = 7` and `U = 10`. Maybe I am doing something wrong. – Nicholas Kinar Jul 25 '12 at 00:32
  • @NicholasKinar I have looked a little further at this problem and provided an answer. Please have a look to see if you agree. – Anders Gustafsson Jul 31 '12 at 07:33

3 Answers3

2

The documentation for fminsearch says how to deal with complex numbers in the limitations section:

fminsearch only minimizes over the real numbers, that is, x must only consist of real numbers and f(x) must only return real numbers. When x has complex variables, they must be split into real and imaginary parts.

You can use the functions real and imag to extract the real and imaginary parts, respectively.

zroth
  • 438
  • 4
  • 9
  • Actually, this probably doesn't help here since your function returns a complex output for a real input. If your function were real-valued, though, then I think that this method would work. – zroth Jul 13 '12 at 17:40
  • Isn't the `diff` variable always real since this is the absolute value of a complex number? If this is the case, then how do I split into real and imaginary parts? – Nicholas Kinar Jul 13 '12 at 23:04
  • @NicholasKinar You're right. Your function seems to take real inputs and be real-valued. Why do you say that `fminsearch` isn't converging properly with complex `Utilde`? Did you find the minimum analytically? – zroth Jul 20 '12 at 18:32
  • Hmm, I don't know if this is a problem with the function, or if there is another issue here. I will investigate. – Nicholas Kinar Jul 21 '12 at 18:54
2

The main problem here is that there is no unique solution to this optimization or parameter fitting problem. For example, looking at the expected and actual results above, Utilde is equivalent (ignoring round-off differences) for the two (x, U) pairs, i.e.

Utilde(x = 1.5, U = 50 + 25i) = Utilde(x = 1.7318, U = 88.8760)

Although I have not examined it in depth, I even suspect that for any value of x, you can find an U that computes to Utilde(x, U) = Utilde(x = 1.5, U = 50 + 25i).

The solution here would thus be to further constrain the parameter fitting problem so that the solver yields any solution that can be considered acceptable. Alternatively, reformulate Utilde to have a unique value for any (x, U) pair.

UPDATE, AUG 1

Given reasonable starting values, it actually seems like it is sufficient to restrict x to be real-valued. Performing unconstrained non-linear optimization using the diff function formulated above, I get the following result:

x = 1.50462926953244
U = 50.6977768845879 + 24.7676554234729i
diff = 3.18731710515855E-06

However, changing the starting guess to values more distant from the desired values does yield different solutions, so restricting x to be real-values does not alone provide a unique solution to the problem.

I have implemented this in C#, using the BOBYQA optimizer, but the numerics should be the same as above. If you want to try outside of Matlab, it should also be relatively simple to turn the C# code below into C++ code using the std::complex class and an (unconstrained) nonlinear C++ optimizer of your own choice. You could find some C++ compatible codes that do not require gradient computation here, and there is also various implementations available in Numerical Recipes. For example, you could access the C version of NR online here.

For reference, here are the relevant parts of my C# code:

class Program
{
    private static readonly Complex Coeff = new Complex(-2.0, 2.0);
    private static readonly Complex UTilde0 = GetUTilde(1.5, new Complex(50.0, 25.0));

    static void Main(string[] args)
    {
        double[] vars = new[] {1.0, 25.0, 0.0}; // xstart = 1.0, Ustart = 25.0
        BobyqaExitStatus status = Bobyqa.FindMinimum(GetObjfnValue, vars.Length, vars);
    }

    public static Complex GetUTilde(double x, Complex U)
    {
        return U * Complex.Exp(Coeff * x);
    }

    public static double GetObjfnValue(int n, double[] vars)
    {
        double x = vars[0]; 
        Complex U = new Complex(vars[1], vars[2]);
        return Complex.Abs(-UTilde0 / U + Complex.Exp(Coeff * x));
    }
}
Anders Gustafsson
  • 15,837
  • 8
  • 56
  • 114
  • Thanks for your insightful answer. How might I reformulate `Utilde` to have a unique value for any `(x, U)` pair? – Nicholas Kinar Jul 31 '12 at 13:57
  • @NicholasKinar What is `Utilde` meant to represent? Can you link to some context where this is illustrated? – Anders Gustafsson Jul 31 '12 at 14:07
  • Utilde (as a complex number) is the output of a product of two numbers in the frequency domain. The `(1 / exp(2 * x)) * exp( 1i * 2 * x)` is the filter kernel, and the `U` is a Gabor transformed signal. This is a seismic Q filtering application, and there is a similar equation on pg. 128 of the book Seismic Inverse Q filtering: http://www.scribd.com/doc/45448335/SEISMIC-INVERSE-q-FILTERING. Could you suggest a good reference book on this type of optimization problem? Maybe if `U` is unknown, this is a blind deconvolution problem? – Nicholas Kinar Jul 31 '12 at 14:46
  • The real part of the expression can be found on pg. 67 of the book, as Equation (4.22) and Equation (4.20). Moreover, the `Utilde` expression is on pg. 128 as Equation (7.21). The `U` is a Gabor-transformed seismic input trace. – Nicholas Kinar Jul 31 '12 at 14:55
  • @NicholasKinar Oh boy! Those were quite _complex_ (pun intended :-) expressions. It is not immediately obvious to me how to deal with these expressions in a parameter fitting scenario (in fact I cannot even see straight away that your `Utilde` corresponds to those equations). In order to go forward I think you need to consider the following: in a real fitting scenario you would only have `Utilde`, right? Does it matter which pair of `x` and `U` that satisfies this `Utilde`? In that case, what are the restrictions on either of the `x` and `U` parameters? – Anders Gustafsson Jul 31 '12 at 19:25
  • Thanks for looking at this, Anders. I have a physical copy of that particular book, and some of the pages are getting very worn through with continued use. Yes, it does matter which pair of `x` and `U` satisfies the `Utilde`, and suppose that we place restrictions on the `x` and `U` parameters, such that the real and complex parts of `x` and `U` are within two values. How do I place these restrictions? Can you suggest a book detailing how to do this? – Nicholas Kinar Jul 31 '12 at 20:23
  • So perhaps placing limits on `x` and `U` would be the way to go? At least, `x` must be between a higher and lower limit, and `U` must be between a higher and lower limit. – Nicholas Kinar Jul 31 '12 at 20:41
  • Yes, I would start in this direction. Maybe there are even further restrictions, i.e. that `x` should be real-valued? – Anders Gustafsson Jul 31 '12 at 20:51
  • Yes, `x` should be real-valued as well. What numerical procedure or algorithm would be the best to do the optimization with these restrictions? – Nicholas Kinar Aug 01 '12 at 13:46
  • 1
    @NicholasKinar Please see my encouraging :-) update to this answer. – Anders Gustafsson Aug 01 '12 at 14:47
  • That's just lovely. Thank you so much for this great answer! Three questions: (1) I am wondering if the starting values for `x` and `U` have to be sufficiently "close" to the optimization values; (2) How does your code restrict `x` to be real-valued; and (3) Could you suggest any reasonably good "cross-platform" C++ optimization library that I can use for this problem as well? – Nicholas Kinar Aug 01 '12 at 15:36
  • Or perhaps I could use Mono with your BOBYQA optimizer. – Nicholas Kinar Aug 01 '12 at 15:51
  • See further update of the answer, in particular I have emphasized `x` to be real-valued (`double`). If you are running on Mac or Linux and it is otherwise no problem for you to develop in C#, then yes, go for Mono. – Anders Gustafsson Aug 01 '12 at 16:09
  • 1
    Thanks, Anders; this is very much appreciated. Thank you for giving me an excellent framework and code for this problem. I will experiment with the exponential function argument a little bit more and then get back to you. – Nicholas Kinar Aug 01 '12 at 16:19
  • After much experimentation, I've found that `x` and `U` need to be close to the actual values, and the framework that you've provided is the best answer to this question. Thank you. – Nicholas Kinar Aug 21 '12 at 14:19
0

It appears that there is no easy way to do this, even if both x and U are real numbers. The equation for Utilde is not well-posed for an optimization problem, and so it must be modified.

I've tried to code up my own version of the Nelder-Mead optimization algorithm, as well as tried Powell's method. Neither seem to work well for this problem, even when I attempted to modify these methods.

Nicholas Kinar
  • 1,440
  • 5
  • 24
  • 36