Matlab's GlobalSearch function : error message at execution

I have a system of matricial equations where I want to find 2 matrices of 7x7 (so I am working with (1x98) vectors).

I have an issue when I want to use GlobalSearch Matlab function. Here my code :

    Eq1 = (P1.')*(a.')*a*P1 + (P1.')*(a.')*b*P2 + (P2.')*(b.')*a*P1 + (P2.')*(b.')*b*P2 - eye(7);
    Eq2 = F1*a*P1 + F1*b*P2 + F2*a*P1 + F2*b*P2 - (a*P1 + b*P2)*D;
    Eqs = reshape([Eq1,Eq2],2*7*7,[]);
    Fun = matlabFunction(Eqs);
    
    F = @(x) Fun(...
        x(  1   ),  x(  2   ),  x(  3   ),  x(  4   ),  x(  5   ),  x(  6   ),  x(  7   ),...
        x(  8   ),  x(  9   ),  x(  10  ),  x(  11  ),  x(  12  ),  x(  13  ),  x(  14  ),...
        x(  15  ),  x(  16  ),  x(  17  ),  x(  18  ),  x(  19  ),  x(  20  ),  x(  21  ),...
        x(  22  ),  x(  23  ),  x(  24  ),  x(  25  ),  x(  26  ),  x(  27  ),  x(  28  ),...
        x(  29  ),  x(  30  ),  x(  31  ),  x(  32  ),  x(  33  ),  x(  34  ),  x(  35  ),...
        x(  36  ),  x(  37  ),  x(  38  ),  x(  39  ),  x(  40  ),  x(  41  ),  x(  42  ),...
        x(  43  ),  x(  44  ),  x(  45  ),  x(  46  ),  x(  47  ),  x(  48  ),  x(  49  ),...
        x(  50  ),  x(  51  ),  x(  52  ),  x(  53  ),  x(  54  ),  x(  55  ),  x(  56  ),...
        x(  57  ),  x(  58  ),  x(  59  ),  x(  60  ),  x(  61  ),  x(  62  ),  x(  63  ),...
        x(  64  ),  x(  65  ),  x(  66  ),  x(  67  ),  x(  68  ),  x(  69  ),  x(  70  ),...
        x(  71  ),  x(  72  ),  x(  73  ),  x(  74  ),  x(  75  ),  x(  76  ),  x(  77  ),...
        x(  78  ),  x(  79  ),  x(  80  ),  x(  81  ),  x(  82  ),  x(  83  ),  x(  84  ),...
        x(  85  ),  x(  86  ),  x(  87  ),  x(  88  ),  x(  89  ),  x(  90  ),  x(  91  ),...
        x(  92  ),  x(  93  ),  x(  94  ),  x(  95  ),  x(  96  ),  x(  97  ),  x(  98  ));
    
    x0 =  ones(2*7*7,1);
    gs = GlobalSearch;
    options = optimoptions('fmincon','Algorithm', 'interior-point',...
        'UseParallel',true, 'Display','off','MaxFunctionEvaluations',6000, 'TolCon', 1e-4, 'TolFun', 1e-4);
    problem = createOptimProblem('fmincon', ...
    'objective',F, ...
    'x0',x0, ...
    'lb',[-1*ones(98)], ...
    'ub',[1*ones(98)], ...
    'options',options)
    [x,fval] = run(gs,problem)

But I get unfortunately the following error at execution :

    problem =
    
      struct with fields:
    
        objective: [function_handle]
               x0: [98x1 double]
            Aineq: []
            bineq: []
              Aeq: []
              beq: []
               lb: [98x98 double]
               ub: [98x98 double]
          nonlcon: []
           solver: 'fmincon'
          options: [1x1 optim.options.Fmincon]
    
    Warning: Length of lower bounds is > length(x); ignoring extra bounds.
    > In checkbounds (line 27)
      In checkglobalsearchnlpinputs (line 36)
      In globalsearchnlp
      In GlobalSearch/run (line 340)
      In compute_solving_Matricial_Global (line 96)
    Warning: Length of upper bounds is > length(x); ignoring extra bounds.
    > In checkbounds (line 41)
      In checkglobalsearchnlpinputs (line 36)
      In globalsearchnlp
      In GlobalSearch/run (line 340)
      In compute_solving_Matricial_Global (line 96)
    Error using fmincon (line 635)
    Supplied objective function must return a scalar value.
    
    Error in globalsearchnlp
    
    Error in GlobalSearch/run (line 340)
                    globalsearchnlp(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON,options,localOptions);
    
    Error in compute_solving_Matricial_Global (line 96)
    [x,fval] = run(gs,problem)
    
    Caused by:
        Failure in initial call to fmincon with user-supplied problem structure.

If there are persons who were faced to the same problem, how could I fix this error message ?

I think that I did a correct using of GlobalSearch but maybe this issue comes from incorrect lower or upper bounds, or even with the size of vectors I use.

Any help would be welcome.

Answer 1

[Walter Roberson]

    Eqs = reshape([Eq1,Eq2],2*7*7,[]);

That tells us that Eqs is (2*7*7) by something

    Fun = matlabFunction(Eqs);

matlabFunction of an array is going to evaluate to an array the same size. Therefore Fun() with an appropriate argument is not going to return a scalar.

   problem = createOptimProblem('fmincon', ...
    'objective',F, ...
    'x0',x0, ...
    'lb',[-1*ones(98)], ...
    'ub',[1*ones(98)], ...
    'options',options)
    [x,fval] = run(gs,problem)

When you use fmincon, your function must return a scalar.

The scalar your Fun must return must be some measure of how "good" the input is, with smaller value (less positive, more negative if appropriate) being "better".

        F = @(x) Fun(...
        x(  1   ),  x(  2   ),  x(  3   ),  x(  4   ),  x(  5   ),  x(  6   ),  x(  7   ),...
        x(  8   ),  x(  9   ),  x(  10  ),  x(  11  ),  x(  12  ),  x(  13  ),  x(  14  ),...
        x(  15  ),  x(  16  ),  x(  17  ),  x(  18  ),  x(  19  ),  x(  20  ),  x(  21  ),...
        x(  22  ),  x(  23  ),  x(  24  ),  x(  25  ),  x(  26  ),  x(  27  ),  x(  28  ),...
        x(  29  ),  x(  30  ),  x(  31  ),  x(  32  ),  x(  33  ),  x(  34  ),  x(  35  ),...
        x(  36  ),  x(  37  ),  x(  38  ),  x(  39  ),  x(  40  ),  x(  41  ),  x(  42  ),...
        x(  43  ),  x(  44  ),  x(  45  ),  x(  46  ),  x(  47  ),  x(  48  ),  x(  49  ),...
        x(  50  ),  x(  51  ),  x(  52  ),  x(  53  ),  x(  54  ),  x(  55  ),  x(  56  ),...
        x(  57  ),  x(  58  ),  x(  59  ),  x(  60  ),  x(  61  ),  x(  62  ),  x(  63  ),...
        x(  64  ),  x(  65  ),  x(  66  ),  x(  67  ),  x(  68  ),  x(  69  ),  x(  70  ),...
        x(  71  ),  x(  72  ),  x(  73  ),  x(  74  ),  x(  75  ),  x(  76  ),  x(  77  ),...
        x(  78  ),  x(  79  ),  x(  80  ),  x(  81  ),  x(  82  ),  x(  83  ),  x(  84  ),...
        x(  85  ),  x(  86  ),  x(  87  ),  x(  88  ),  x(  89  ),  x(  90  ),  x(  91  ),...
        x(  92  ),  x(  93  ),  x(  94  ),  x(  95  ),  x(  96  ),  x(  97  ),  x(  98  ));
    

Oh, don't do that. Instead

       Fun = matlabFunction(EXPRESSION, 'vars', {[row vector of 98 variables]})

[petit]

@Walter Roberson

Thanks for your quick answer.

1. Unfortunately, I don't know how to deal with your solution :

Fun = matlabFunction(EXPRESSION, 'vars', {[row vector of 98 variables]})

the EXPRESSION corresponds to what ?

Eqs ?

when I do :

Fun = matlabFunction(Eqs);

?

2. And what about the "row vector of 98 variables" ? How could I define it in practice ? : is it an array of symbolic variables "x" of size (1x98) ? How to declare it please ?

3. Once I have defined "Fun" , how can I made the link with "F" function :

F = @(x) Fun

is good ?

Best Regards

[Walter Roberson]

Eqs_vars = num2cell(symvar(Eqs));
F = matlabFunction(Eqs(:).', 'var', Eqs_vars);

No Fun; this F can be used directly in

    problem = createOptimProblem('fmincon', ...
    'objective',F, ...
    'x0',x0, ...
    'lb',[-1*ones(1,98)], ...
    'ub',[1*ones(1,98)], ...
    'options',options)

Note: I corrected your lb and ub here.

[petit]

@Walter Roberson

Thanks for your patience. Here the part of code concerned :

%% unknown variables
a = sym('a_',[7,7]);
b = sym('b_',[7,7]);
%% generate function from equations
D = (eye(7).*diag(a).^2)*D1 + (eye(7).*diag(b).^2)*D2;
Eq1 = (P1.')*(a.')*a*P1 + (P1.')*(a.')*b*P2 + (P2.')*(b.')*a*P1 + (P2.')*(b.')*b*P2 - eye(7);
Eq2 = F1*a*P1 + F1*b*P2 + F2*a*P1 + F2*b*P2 - (a*P1 + b*P2)*D;
Eqs = reshape([Eq1,Eq2],2*7*7,[]);
%% Your solution
Eqs_vars = num2cell(symvar(Eqs));
F = matlabFunction(Eqs(:).', 'var', Eqs_vars);
%% Global search 
x0 =  ones(1,98);
gs = GlobalSearch;
options = optimoptions('fmincon','Algorithm', 'interior-point',...
    'UseParallel',true, 'Display','off','MaxFunctionEvaluations',6000, 'TolCon', 1e-10, 'TolFun', 1e-10);
problem = createOptimProblem('fmincon', ...
'objective',F, ...
'x0',x0, ...
'lb',[-10*ones(1,98)], ...
'ub',[10*ones(1,98)], ...
'options',options);
run(gs,problem);

Unfortunately, there is a first error at execution :

Not enough input arguments.

Error in
symengine>@(a_1_1,a_1_2,a_1_3,a_1_4,a_1_5,a_1_6,a_1_7,a_2_1,a_2_2,a_2_3,a_2_4,a_2_5,a_2_6,a_2_7,a_3_1,a_3_2,a_3_3,a_3_4,a_3_5,a_3_6,a_3_7,a_4_1,a_4_2,a_4_3,a_4_4,a_4_5,a_4_6,a_4_7,a_5_1,a_5_2,a_5_3,a_5_4,a_5_5,a_5_6,a_5_7,a_6_1,a_6_2,a_6_3,a_6_4,a_6_5,a_6_6,a_6_7,a_7_1,a_7_2,a_7_3,a_7_4,a_7_5,a_7_6,a_7_7,b_1_1,b_1_2,b_1_3,b_1_4,b_1_5,b_1_6,b_1_7,b_2_1,b_2_2,b_2_3,b_2_4,b_2_5,b_2_6,b_2_7,b_3_1,b_3_2,b_3_3,b_3_4,b_3_5,b_3_6,b_3_7,b_4_1,b_4_2,b_4_3,b_4_4,b_4_5,b_4_6,b_4_7,b_5_1,b_5_2,b_5_3,b_5_4,b_5_5,b_5_6,b_5_7,b_6_1,b_6_2,b_6_3,b_6_4,b_6_5,b_6_6,b_6_7,b_7_1,b_7_2,b_7_3,b_7_4,b_7_5,b_7_6,b_7_7)[a_1_1.*(b_1_1.*1.325090276395953e-2+b_1_2.*2.088183766129291e-3+b_1_3.*2.470923064794312e-1-b_1_4.*9.687105021306056e-1-b_1_5.*1.326763761399711e-2+b_1_6.*2.014497696984806e-3-b_1_7.*1.361322137542385e-2).*(-3.203162853943744e-2)-a_1_2.*(b_1_1.*1.325090276395953e-2+b_1_2.*2.088183766129291e-3+b_1_3.*2.470923064794312e-1-b_1_4.*9.687105021306056e-1-b_1_5.*1.326763761399711e-2+b_1_6.*2.014497696984806e-3-b_1_7.*1.361322137542385e-2).*6.462760798949253e-3-a_1_3.*(b_1_1.*1.325090276395953e-2+b_1_2.*2.088183766129291e-3+b_1_3.*2.470923064794312e-1-b_1_4.*9.687105021306056e-
...
...

First, it shows the 98 unknown corresponding to flattened a and b wanted matrices and after a mix between symbolic system and numerical values.

Which is the reason of this symengine error ?

Best regards

[petit]

the issue seems to be pretty due to not enough input arguments

[petit]

Is really no one who could give me a little help about this last error ? Regards

[Walter Roberson]

I was asleep, and I still need some more sleep.

You have not posted enough for me to be able to test your code. I need to know D1, D2, the F variables, the P variables, or at least all their sizes.

In the meantime,

Eqs_vars = {[a(:); b(:)].'});

Answer 2

[petit]

@Walter Roberson

thanks, I provide all the code below. D1 and D2 are 7x7 diagonal matrices, F1, F1 and P1, P2 are 7x7 matrices. The system of equations is formulated with Eq1 and Eq2 ( 2 matricial equations actually : solution is a vector 1x98 elements).

I don't understand the suggestion in your last post :

Eqs_vars = {[a(:), b(:)].'};

Indeed, previously, you suggested me to do :

Eqs_vars = num2cell(symvar(Eqs));
F = matlabFunction(Eqs(:).', 'var', Eqs_vars);

Here the full code :

clear
clc

% Load
FISH_GCsp = load('Fisher_GCsp_flat');
FISH_XC = load('Fisher_XC_GCph_WL_flat');

% Marginalizing over uncommon parameters between the two matrices
COV_GCsp_first = inv(FISH_GCsp);
COV_XC_first = inv(FISH_XC);
COV_GCsp = COV_GCsp_first(1:7,1:7);
COV_XC = COV_XC_first(1:7,1:7);

% Invert to get Fisher matrix
FISH_sp = inv(COV_GCsp);
FISH_xc = inv(COV_XC);

% Get eigen values (marginalized error since we handle Covariance error)
% and eigen vectors giving matrix "P" with convention : F = P D P^-1
[eigenv_sp, eigen_sp] = eig(COV_GCsp);
[eigenv_xc, eigen_xc] = eig(COV_XC);

% Fisher eigen scalar vectors
FISH_eigen_sp = inv(eigen_sp);
FISH_eigen_xc = inv(eigen_xc);

% Eigen sum for  Fisher matrices
FISH_eigen_sum = FISH_eigen_sp + FISH_eigen_xc;
% Eigen sum matrix for  Fisher matrices
%FISH_eigenv_sum = FISH_eigenv_sp + FISH_eigenv_xc

% Fisher matrices
F1 = FISH_sp;
F2 = FISH_xc;
P1 = eigenv_sp;
P2 = eigenv_xc;
D1 = FISH_eigen_sp;
D2 = FISH_eigen_xc;

%% unknown variables
a = sym('a_',[7,7]);
b = sym('b_',[7,7]);

%% generate function from equations
D = (eye(7).*diag(a).^2)*D1 + (eye(7).*diag(b).^2)*D2;
Eq1 = (P1.')*(a.')*a*P1 + (P1.')*(a.')*b*P2 + (P2.')*(b.')*a*P1 + (P2.')*(b.')*b*P2 - eye(7);
Eq2 = F1*a*P1 + F1*b*P2 + F2*a*P1 + F2*b*P2 - (a*P1 + b*P2)*D;
Eqs = reshape([Eq1,Eq2],2*7*7,[]);
%Fun = matlabFunction(Eqs);

% Solution
Eqs_vars = num2cell(symvar(Eqs));
F = matlabFunction(Eqs(:), 'var', Eqs_vars);

options = optimoptions('fmincon','Algorithm', 'interior-point',...
    'UseParallel',true, 'Display','off','MaxFunctionEvaluations',6000, 'TolCon', 1e-4, 'TolFun', 1e-4);
problem = createOptimProblem('fmincon', ...
'objective',F, ...
'x0',x0, ...
'lb',[-1e3*ones(1,98)], ...
'ub',[1e3*ones(1,98)], ...
'options',options);
gs = GlobalSearch;
[x, fval] = run(gs,problem);

%% output
a = zeros(7,7);
b = zeros(7,7);
for ii=1:7
    a(ii,:) = x(1+7*(ii-1):7*ii);
    b(ii,:) = x(50+7*(ii-1):49+7*ii);
end

disp('a*transpose(a) = ')
a*a'
disp('b*transpose(b) = ')
b*b'

dlmwrite('a_normal.txt',a,'delimiter',' ')
dlmwrite('b_normal.txt',b,'delimiter',' ')

% Sum of passing matrix P = aX + bY with X = eigenv_sp, Y = eigenv_xc 
% a and b to infer
eigenv_final = a*eigenv_sp + b*eigenv_xc

% Fisher final
FISH_final = eigenv_final*FISH_eigen_sum*(eigenv_final.');

% Save Fisher_final
dlmwrite('Fisher_final.txt',FISH_final,'delimiter',' ')

Hoping you will understand, Regards

[Walter Roberson]

Replace

Eqs_vars = num2cell(symvar(Eqs));
F = matlabFunction(Eqs(:), 'var', Eqs_vars);

with

Eqs_vars = {[a(:); b(:)].'};
F = matlabFunction(Eqs(:), 'var', Eqs_vars);

Note:

FISH_GCsp = load('Fisher_GCsp_flat');
FISH_XC = load('Fisher_XC_GCph_WL_flat');

You did not attach those files, and you did not tell me how big the resulting arrays are, so I cannot substitute random data.

[Walter Roberson]

Ah, sorry, I did not see you had posted the array sizes.

[petit]

@Walter Roberson

Excuse me, the code and the files are available on :

code

Fisher_GCsp_flat

Fisher_XC_GCph_WL_flat

Normally, you should be able to run the code.

[petit]

@Walter Roberson

with your correction, I get the new error messages

Error using fmincon (line 635)
Supplied objective function must return a scalar value.

Error in globalsearchnlp

Error in GlobalSearch/run (line 340)
                globalsearchnlp(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON,options,localOptions);

Error in compute_solving_Matricial_Global (line 103)
[x, fval] = run(gs,problem);

Caused by:
    Failure in initial call to fmincon with user-supplied problem structure.

[Walter Roberson]

Here is repaired and generalized code, with all warnings cleaned up.

The function is generated from the equations, and the function is properly configured to expect a vector of 2*N^2 values.

Now what you have to do is concentrate on the fact that your function is returning 2*N^2 values, which is totally unsuited for any of the optimizers other than gamultiobj() (though you could also use it with paretosearch() )

You need to figure out what you want to optimize. Currently you are asking to optimize 98 different things simultaneously.

N = 7;
Nsq2 = 2*N*N;

if ispc()
    FISH_GCsp = load('Fisher_GCsp_flat');
    FISH_XC = load('Fisher_XC_GCph_WL_flat');
else
    
    FISH_GCsp = randn(N,N);
    FISH_XC = randn(N,N);
end


% Marginalizing over uncommon parameters between the two matrices
COV_GCsp_first = inv(FISH_GCsp);
COV_XC_first = inv(FISH_XC);
COV_GCsp = COV_GCsp_first(1:N,1:N);
COV_XC = COV_XC_first(1:N,1:N);

% Invert to get Fisher matrix
FISH_sp = inv(COV_GCsp);
FISH_xc = inv(COV_XC);

% Get eigen values (marginalized error since we handle Covariance error)
% and eigen vectors giving matrix "P" with convention : F = P D P^-1
[eigenv_sp, eigen_sp] = eig(COV_GCsp);
[eigenv_xc, eigen_xc] = eig(COV_XC);

% Fisher eigen scalar vectors
FISH_eigen_sp = inv(eigen_sp);
FISH_eigen_xc = inv(eigen_xc);

% Eigen sum for  Fisher matrices
FISH_eigen_sum = FISH_eigen_sp + FISH_eigen_xc;
% Eigen sum matrix for  Fisher matrices
%FISH_eigenv_sum = FISH_eigenv_sp + FISH_eigenv_xc

% Fisher matrices
F1 = FISH_sp;
F2 = FISH_xc;
P1 = eigenv_sp;
P2 = eigenv_xc;
D1 = FISH_eigen_sp;
D2 = FISH_eigen_xc;

%% unknown variables
a = sym('a_',[N,N]);
b = sym('b_',[N,N]);

%% generate function from equations
D = (eye(N).*diag(a).^2)*D1 + (eye(N).*diag(b).^2)*D2;
Eq1 = (P1.')*(a.')*a*P1 + (P1.')*(a.')*b*P2 + (P2.')*(b.')*a*P1 + (P2.')*(b.')*b*P2 - eye(N);
Eq2 = F1*a*P1 + F1*b*P2 + F2*a*P1 + F2*b*P2 - (a*P1 + b*P2)*D;
Eqs = reshape([Eq1,Eq2], 1, []);

% Solution
Eqs_vars = {[a(:); b(:)].'};
F = matlabFunction(Eqs(:), 'var', Eqs_vars);

x0 =  ones(1, Nsq2);

options = optimoptions('fmincon','Algorithm', 'interior-point',...
    'UseParallel',true, 'Display','off','MaxFunctionEvaluations',6000, 'TolCon', 1e-4, 'TolFun', 1e-4);
problem = createOptimProblem('fmincon', ...
    'objective', F, ...
    'x0', x0, ...
    'lb', -1e3*ones(1,Nsq2), ...
    'ub', 1e3*ones(1,Nsq2), ...
    'options', options);
gs = GlobalSearch;
[x, fval] = run(gs,problem);

%% output
a = zeros(N,N);
b = zeros(N,N);
for ii=1:N
    a(ii,:) = x(1+N*(ii-1):N*ii);
    b(ii,:) = x(N*N+1+N*(ii-1):N*N+N*ii);
end

disp('a*transpose(a) = ')
disp(a*a');
disp('b*transpose(b) = ')
disp(b*b');

dlmwrite('a_normal.txt',a,'delimiter',' ')
dlmwrite('b_normal.txt',b,'delimiter',' ')

% Sum of passing matrix P = aX + bY with X = eigenv_sp, Y = eigenv_xc
% a and b to infer
eigenv_final = a*eigenv_sp + b*eigenv_xc;
disp(eigenv_final)

% Fisher final
FISH_final = eigenv_final*FISH_eigen_sum*(eigenv_final.');

% Save Fisher_final
dlmwrite('Fisher_final.txt',FISH_final,'delimiter',' ')

[Walter Roberson]

My tests suggest that your equations are likely to be complex-valued -- the eig() will be complex valued unless the matrices have some narrow properties. Not impossible, but at the very least the system is fragile.

Remember that using inv() is numerically fragile.

[petit]

thanks for your quick answer.

Are you sure that the code of your last answer is working well ? If I execute it with Matlab 2020a macOS, I get the following errors :

Error using fmincon (line 635)
Supplied objective function must return a scalar value.

Error in globalsearchnlp

Error in GlobalSearch/run (line 340)
                globalsearchnlp(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON,options,localOptions);

Error in compute_solving_Matricial_last_dev_bis (line 58)
[x, fval] = run(gs,problem);

Caused by:
    Failure in initial call to fmincon with user-supplied problem structure.

However, we had already solved this problem with your first suggestion :

Eqs_vars = num2cell(symvar(Eqs));
F = matlabFunction(Eqs(:).', 'var', Eqs_vars);

which has been replaced after by :

Eqs_vars = {[a(:); b(:)].'};
F = matlabFunction(Eqs(:), 'var', Eqs_vars);

I am surprised that you have no errors with your code.

I have also noticed that you wrote :

Eqs = reshape([Eq1,Eq2], 1, []);

instead of before :

Eqs = reshape([Eq1,Eq2],2*7*7,[]);

If you could check again if there is no bug in the code you provided in your last answer, this would be fine from your part, I don't want to be boring.

Regards

[Walter Roberson]

I am sure there is no bug in my code.

Your Eq1 and Eq2 are both 7 x 7 systems, so [Eq1,Eq2] is a 7 x 14 system, and reshaping that to 2*7*7, [] is reshaping to 98 x 1 . reshape([Eq1, Eq2], 1, []) reshapes to 1 x whatever needed, which is 1 x 98 . So the difference between the old code and the new code is whether you are creating Eqs as a 98 x 1 vector or as a 1 x 98 vector.

The vector of either 98 x 1 (old) or 1 x 98 (new) equations is passed to matlabFunction() . In the old case, that would cause matlabFunction() to generate an anonymous function that returns an column vector of values; in the new code, that would cause matlabFunction() to generate an anonymous function that returns a row vector of the same values as in the other case. So the two cases are just transposes of the other.

fmincon() does not care whether the function returns a row vector or a column vector.

With fmincon() not caring whether the function returns a row vector or a column vector, the "right" code is the simpler code -- and reshape(expression, 1, []) is simpler than reshape(expression, 2*7*7, 1) because the new version does not have to build-in the sizes.

[Walter Roberson]

If I execute it with Matlab 2020a macOS

Heh. I needed random data to test with, and I figured you were probably on pc, so instead of adding some code to test the existance of the files and use random data if not found, I used the shortcut of testing ispc()... since I am on Mac myself and did not realize you are too. You will find that with that code, random data is being created instead of your files being loaded.

[Walter Roberson]

Supplied objective function must return a scalar value.

However, we had already solved this problem with your first suggestion :

No, you misunderstood the purpose of that code with creating vectors of symbolic variables. The purpose of that code with num2cell() (before) and with a() and b() now, is to generate a vector of symbolic variable names inside a cell array, to pass to the 'vars' option of matlabFunction . Passing a cell array containing a vector of variable names, causes matlabFunction to extract all of the variables from a single parameter instead of needing multiple parameters.

I was not dealing at all with the objective function needing to return a scalar value: I was dealing with cleaning up the

        F = @(x) Fun(...
        x(  1   ),  x(  2   ),  x(  3   ),  x(  4   ),  x(  5   ),  x(  6   ),  x(  7   ),...
        x(  8   ),  x(  9   ),  x(  10  ),  x(  11  ),  x(  12  ),  x(  13  ),  x(  14  ),...

code, having it generate code that could use a vector of values directly instead of having to have that layer of splitting the vector into individual values.

I did that to get you past the error of about not enough input arguments, because until that problem was solved, you were refusing to pay attention to what I was telling you about your function having to return a scalar. As I posted early on:

When you use fmincon, your function must return a scalar.

The scalar your Fun must return must be some measure of how "good" the input is, with smaller value (less positive, more negative if appropriate) being "better".

and more recently I posted

Now what you have to do is concentrate on the fact that your function is returning 2*N^2 values, which is totally unsuited for any of the optimizers other than gamultiobj() (though you could also use it with paretosearch() )

You need to figure out what you want to optimize. Currently you are asking to optimize 98 different things simultaneously.

You needed the rest of the code cleaned up before you would pay attention to the fact you are asking to optimized the wrong thing .

A moment ago, I posted that fmincon does not care whether the function returns a row vector or a column vector. The reason that fmincon does not care about that, is that fmincon will generate an error either way . fmincon will generate an error unless the function returns a scalar.

You have 98 equations. What are you searching for? How would you know if you found it?

Are you possibly trying to find the zeros of the equations? If so then use fsolve() instead of fmincon(), if you have reason to believe that there is a vector of 98 input values that will return simultaneous zeros for the 98 different equations. If you are not certain that there are exact solutions, and you are just hoping to find as close as you can get to solving all 98 simultaneously, then minimize the L2-norm of the set of equations.

I said L2-norm instead of sum of squares of the values because, as I posted earlier, it looks to me as if your equations typically generate complex values. There might be some input matrices that lead to real-valued equations being generated... but as I pointed out, inv() is numerically fragile, so you can accidentally introduce matrices with complex eigenvalues because of numeric round-off even if mathematically you would not expect complex.

[petit]

Thanks for all these detailled explanations.

At the beginning, I have the 2 matricial equations where matrices "a" and "b" are unknown (so I have to find 2*49 values) :

what I am looking for since the beginning of this personal work is :

1) diagnalize Fisher F1 and F2 matrices

2) search an endomorphism which has D=D1+D2 as eigen values

3) search this endomorphism with passing matrix under the form : P = a P1 + b 2 with a and b unknown matrices (so 98 unknown)

4) once a and b known, coming back to start space, I can build the endomorphism by : F = P D P^-1

Initially, I used fsolve but is very sensitive to the guess values, so I thought that I had to looking for a global minimum and not just different local minimum.

That's why I am redirected myself to GlobalSearch which I think looks for global minimum.

Do you understand my approach ? I have to explore a parameters space of 98 unknown, but maybe with a lower/upper bound, I can find this global minimum.

Are there any solver for this kind of problem ?

Regards

[Walter Roberson]

When you were trying fsolve(), were you trying to fsolve() a single equation, or where you trying to solve 98 different equations simulteneously?

so I thought that I had to looking for a global minimum

Global minimum of what though?

[petit]

When I am trying fsolve(), I am looking for fill "a" and "b" matrices, so I have 98 unknown. But as I said, this is very guess dependent.

So I conclude that instead of searching different local minimum (which are very guess dependent I recall) , it would deserve to use a function that explore a space of parameters more extended and find a global minimum.

Does all this approach make sense for you ? or Isn't just wrong to have this reasoning ?

Regards

[Walter Roberson]

Are you looking for the a and b such that the equations work out to zero ?

If so then minimizing the same equations is the wrong operation. If you want to find out how close to zero you can get, then you need to minimize the squares of the functions. But then you are stuck because you are only asking how close to zero each individual equation can get, not how close combinations can get.

If you want to minimize the squares of the equations simultaneously then you need a relative weighting system for how important it is to satisfy each equation relative to each other. For example, support one equation only has a range of +/- 0.1, and a change of 0.03 for it might be huge, but for another equation that has a range of 10^7, a change of 0.03 for it might be round-off error.

I suggest you try with:

objective = sum(Eqs .* conj(Eqs));
F = matlabFunction(objective, 'var', Eqs_vars);

[petit]

@Walter Roberson

Yes, I am looking for the 2 matricial equations (with 98 unknown for "a" and "b" matrices) which must be equal to null matrix, i.e each equations must be equal to zero.

I tried with your suggestion, below the snippet code that I run :

D = D1 + D2
Eq1 = (P1.')*(a.')*a*P1 + (P1.')*(a.')*b*P2 + (P2.')*(b.')*a*P1 + (P2.')*(b.')*b*P2 - eye(N);
Eq2 = F1*a*P1 + F1*b*P2 + F2*a*P1 + F2*b*P2 - (a*P1 + b*P2)*D;
Eqs = reshape([Eq1,Eq2], 98, []);
% Solution
Eqs_vars = {[a(:); b(:)].'};
objective = sum(Eqs .* conj(Eqs));
F = matlabFunction(objective, 'var', Eqs_vars);

x0 =  ones(1, Nsq2);
options = optimoptions('fmincon','Algorithm', 'interior-point',...
    'UseParallel',true, 'Display','off','MaxFunctionEvaluations',100000, 'TolCon', 1e-6, 'TolFun', 1e-6);
problem = createOptimProblem('fmincon', ...
    'objective', F, ...
    'x0', x0, ...
    'lb', -1e10*ones(1,Nsq2), ...
    'ub', 1e10*ones(1,Nsq2), ...
    'options', options);
gs = GlobalSearch;
[x, fval] = run(gs,problem);
% output
a = zeros(N,N);
b = zeros(N,N);
for ii=1:N
    a(ii,:) = x(1+N*(ii-1):N*ii);
    b(ii,:) = x(N*N+1+N*(ii-1):N*N+N*ii);

1) Do you think it is correct ?

2) For a reason I ignore, after few minutes, the following message appears :

GlobalSearch stopped because it analyzed all the trial points.

1 out of 2 local solver runs converged with a positive local solver exit flag.

and the elements of matrices "a" and "b" are all equal to 0.

3) From this message, it seems that 1 solver converged but how to explain the zero values for matrix ?

But how to deal with this exit flag ?

Thanks again

[Walter Roberson]

GlobalSearch uses multiple starting positions. In your case it tried two configurations. One of the ones it tried did not converge; the other one did. This is not generally a problem.

However, my testing finds that even though mathematically the x.*conj(x) process should only result in non-negative real-only values, that in practice it is resulting in values with imaginary components as well, probably due to round-off error. I am testing now with a correction for that.

[Walter Roberson]

If you change

objective = sum(Eqs .* conj(Eqs));

to

objective = real(sum(Eqs .* conj(Eqs)));

then it should run notably longer.

Mathematically it should be the same...

[petit]

I cite you :

"GlobalSearch uses multiple starting positions. In your case it tried two configurations. One of the ones it tried did not converge; the other one did. This is not generally a problem."

I cite myself :

"3) From this message, it seems that 1 solver converged but how to explain the zero values for matrix ?

But how to deal with this exit flag ?"

I would like to get access to the solution that has converged and not the solution that gives all zeros values for "a" and "b" matrices, i.e null matrices. How to get this converged solution ?

[petit]

and why do you say that "GlobalSearch" tries 2 configurations, you talk about the lower/upper bounds ?

[Walter Roberson]

The all-zero solution was the converged version. Evaluating the objective on the all-zero vector gave a real-valued result, and evaluating everywhere else gave a complex-value result, so it decided the minimum was at 0.

Or at least that is my guess, based upon the rand() that I put in the code to generate test data as I still do not have actual data to test with. :(

GlobalSearch tries multiple start points; for each start point, it runs fmincon, which take however long it takes. WIth the options you used, if you have the Parallel Computing Toolbox, GlobalSearch runs several different fmincon in parallel, each from a different start point.

Some of the start points might not produce a useful result; some of them might give complex answers for example. Even if only due to round-off error in the way the calculations were turned into code.

In your case, the all-one start point produced all complex values, and the all-zero start point produced a non-complex value for that one location only.

The suggestion I made to use real() gets rid of the tiny imaginary component that mathematically does not exist, allowing the real computations to proceed.

[petit]

@Walter Roberson

I have provided to you the 2 Matrices FISH_GCsp and FISH_XC.

Caution : remove the .txt extension for Fisher_GCsp_flat.txt and Fisher_XC_GCph_WL_flat.txt

It would be fine to compare our results on the solution found.

Best regards

[petit]

@Walter Roberson

Hi Walter, did you try my code above with real data that I provided in attachment in my previous answer ( Fisher_GCsp_flat.txt and Fisher_XC_GCph_WL_flat.txt ?

I would be grateful if you could, that would allow me to compare with my results.

Best regards

[Walter Roberson]

I did try that. It filled up my hard drive (I have a suspicion as to why).

Unfortunately it turns out that with this new operating system that I am using, when you fill up the hard drive, you literally cannot delete any files. Unless you get lucky when you reboot. Unfortunately when I rebooted, my computer lost track of my drive, and it took me a day to recover it...

I won't be trying that again until I install more hard disk space (I'm really short on drive space at the moment as I had a major disk failure recently.)

[petit]

@Walter Roberson

Thanks to keep me informed, on my computer, it takes a long time before producing results (unfortunately, the previous results that returns the 2 arrays with full zeros).

I hope this will get better for your hardware. Don't hesitate to leave a comment or suggestion if you think where is my error to not get 2 null arrays as final solution.

Best regards

[petit]

Isn't there anyone whc could to run the code in Matlab ? this would be fine.

@Walter Robsertson can't run it since a possible issue of disk space or leak memory.

I am using Matlab 2018a on Linux Debian 10.

Thanks in advance.

[Walter Roberson]

I received my new drive yesterday; I will be trying to install it today. (It might take a while.)

[petit]

Cool ! you will tell me ...

Best wishes

[petit]

Helo Walter !

how are you ? Did you manage to make run the code above with the input Fisher matrices in attachment ?

Regards

[Walter Roberson]

I have been testing my new hard drive system... which has some quirks.

[petit]

Hi Walter ! always in trouble with your hard drive ? If you need further informations about my issue, don't hesitate to ask ...

Best regards

[Walter Roberson]

You would be amazed how complicated setting up hard drives can be. Seriously.

[petit]

Hi Walter, you can tell me directly that you haven't enough time to testing by execute my code, I will understand , I don't want to be boring and I won't be annoyed. Best regards.