How to specify limits for lsqnonlin

I would like to specify that variables in my lsqnonlin fit are real. The other vectors in my problem are complex so I cannot use a fully real solver.

My function has the following form. r and theta are real-valued coordinate matrices. Intensity is a complex valued matrix. p is the real-valued and positive vector (at least that is what I want to specify). Hej is the function to be minimized.

Hej=  @(p)(p(1)*besselj(1,p(2)*r)+(p(4)*besselj(1,p(6)*r).*cos(theta+p(8))+p(9)*besselj(1,p(6)*r).*sin(theta+p(8)))*exp(1i*p(5))).*(r/a<=1) ...
                            + (p(1)*besselj(1, p(2)*a)/besselk(1, p(3)*a)*besselk(1,p(3)*r)+...
                               (p(4)*besselj(1, p(6)*a)/besselk(1, p(7)*a).*cos(theta+p(8))+p(9)*besselj(1, p(6)*a)/besselk(1, p(7)*a).*sin(theta+p(8))).*exp(1i*p(5))).*(r/a>1)...
                            -Intensity;

opts = optimoptions(@lsqnonlin,'DiffMaxChange', 0.1,'FinDiffType', 'central', 'Display','off','MaxFunEvals',2E7,'TolFun',1E-18,'TolX',1E-24,'MaxIter',4E3);
    x0 = st; % arbitrary initial guess
    lb = 0.0*ones(size(st));
    [p_estimated,resnorm,residuals,exitflag,output] = lsqnonlin(Hej,x0, lb,[], opts);

The initial guess is given as real-valued and positive vector.

So how do I specify that the only valid solution are positive and realvalued?

Answer 1

[Matt J]

Implement the objective function as follows, splitting Hej into real and imaginary parts,

 function out=objective(p,Intensity,r,theta)

    Hej = (p(1)*besselj(1,p(2)*r)+(p(4)*besselj(1,p(6)*r).*cos(theta+p(8))+p(9)*besselj(1,p(6)*r).*sin(theta+p(8)))*exp(1i*p(5))).*(r/a<=1)  + (p(1)*besselj(1, p(2)*a)/besselk(1, p(3)*a)*besselk(1,p(3)*r)+(p(4)*besselj(1, p(6)*a)/besselk(1, p(7)*a).*cos(theta+p(8))+p(9)*besselj(1, p(6)*a)/besselk(1, p(7)*a).*sin(theta+p(8))).*exp(1i*p(5))).*(r/a>1)-Intensity;

             out=[real(Hej); imag(Hej)]; %Split into real and imaginary parts

 end

and the minimization as

 lsqnonlin(@(p)objective(p,Intensity,r,theta) , x0, lb,[], opts);

[Roger Wohlwend]

I considered that solution, too, but it is basically the same as the one I proposed. In both cases you sum up the squared real and imaginary parts.

[Matt J]

Both implementations correspond to the same sum of squares. However, the underlying algorithms will try to compute the Jacobian of the residual vector that your objective function returns as output, and expects this vector to be differentiable. But note that abs(Hej) will be non-differentiable where Hej=0. Conversely, [real(Hej); imag(Hej)] shouldn't have that problem, assuming Hej(p) has differentiable real and imaginary parts.

[Stine Larsen]

Thank you for the suggestion. I just tried to implement the suggestion but it seems that Matlab does not accept using real to find the real part of a function handle (the same goes for imag). I get the error:

Undefined function 'real' for input arguments of type 'function_handle'.

[Matt J]

Note that Hej is no longer an anonymous function in my version of the code. The @(p) is gone.

[Stine Larsen]

I managed to find a way around the first problem reformulating the objective function to:

function out=objective(p, Intensity, theta, r)

       Hej=  @(p) (p(1)*besselj(0,p(2)*r)+besselj(1,p(6)*r).*(p(4)*cos(theta+p(8))+p(9)*sin(theta+p(8)))*exp(1i*p(5))).*(r/a<=1) ...
                              + (p(1)*besselj(0, p(2)*a)/besselk(0, p(3)*a)*besselk(0,p(3)*r)+...
                                 (p(4)*besselj(1, p(6)*a)/besselk(1, p(7)*a)*cos(theta+p(8))+p(9)*besselj(1, p(6)*a)/besselk(1, p(7)*a)*sin(theta+p(8))).*besselk(1, p(7)*r).*exp(1i*p(5))).*(r/a>1)...
                                 -Intensity;

                             Hej1= @(p) real((p(1)*besselj(0,p(2)*r)+besselj(1,p(6)*r).*(p(4)*cos(theta+p(8))+p(9)*sin(theta+p(8)))*exp(1i*p(5))).*(r/a<=1) ...
                              + (p(1)*besselj(0, p(2)*a)/besselk(0, p(3)*a)*besselk(0,p(3)*r)+...
                                 (p(4)*besselj(1, p(6)*a)/besselk(1, p(7)*a)*cos(theta+p(8))+p(9)*besselj(1, p(6)*a)/besselk(1, p(7)*a)*sin(theta+p(8))).*besselk(1, p(7)*r).*exp(1i*p(5))).*(r/a>1)...
                                 -Intensity);

                             Hej2= @(p) imag((p(1)*besselj(0,p(2)*r)+besselj(1,p(6)*r).*(p(4)*cos(theta+p(8))+p(9)*sin(theta+p(8)))*exp(1i*p(5))).*(r/a<=1) ...
                              + (p(1)*besselj(0, p(2)*a)/besselk(0, p(3)*a)*besselk(0,p(3)*r)+...
                                 (p(4)*besselj(1, p(6)*a)/besselk(1, p(7)*a)*cos(theta+p(8))+p(9)*besselj(1, p(6)*a)/besselk(1, p(7)*a)*sin(theta+p(8))).*besselk(1, p(7)*r).*exp(1i*p(5))).*(r/a>1)...
                                 -Intensity);

                             out= @(p) [Hej1; Hej2]; %Split into real and imaginary parts

But now I run into a second problem, lsqnonlin only accepts inputs in double format and I get the following error:

Error using lsqnonlin (line 233)
LSQNONLIN requires all values returned by user functions to be of data type double.

Error in Mode_decomposition_2modes_field_real_imag (line 191)
    [p_estimated,resnorm,residuals,exitflag,output] = lsqnonlin(@(p) objective(p,Intensity,r,theta), x0, lb,[], opts);

[Matt J]

See my last comment.

[Stine Larsen]

Using a cell array in 'out' this is solved.

[Matt J]

Way, way, way over-complicated. You have functions within functions within cells, whereas my original proposal was just two lines of numeric MATLAB operations...

Answer 2

[Roger Wohlwend]

I am afraid you cannot order the optimizing function to search for a real solution. If you want a real solution you have to make sure that your function Hej returns only real values. Since you use lsqnonlin as optimizing function your goal is to minimize the sum of the squared elements of the vector Hej. In my opinion, it does not make sense to minimize the sum of squared complex numbers. It would make more sense to minimize the sum of the absolute values of the complex numbers. So perhaps you should adjust Hej in a way that it returns the absolute values of the complex numbers instead of the complex numbers. If you do that, Hej returns only real values. As a consequence lsqnonlin should search for real solutions only.

[Stine Larsen]

Thank you for the suggestion. I will try it out and see if it gives a physical valid solution. I would still like to consider the full complex structure in case anyone has other suggestions.