Matlab's GlobalSearch function : error message at execution
Show older comments
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.
Accepted Answer
Walter Roberson
on 12 Dec 2020
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]})
6 Comments
petit
on 12 Dec 2020
Walter Roberson
on 12 Dec 2020
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
on 12 Dec 2020
petit
on 12 Dec 2020
Walter Roberson
on 12 Dec 2020
Edited: Walter Roberson
on 13 Dec 2020
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(:)].'});
More Answers (1)
petit
on 13 Dec 2020
32 Comments
Walter Roberson
on 13 Dec 2020
Edited: Walter Roberson
on 13 Dec 2020
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
on 13 Dec 2020
Ah, sorry, I did not see you had posted the array sizes.
Walter Roberson
on 13 Dec 2020
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
on 13 Dec 2020
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
on 13 Dec 2020
Walter Roberson
on 13 Dec 2020
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
on 13 Dec 2020
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
on 13 Dec 2020
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
on 13 Dec 2020
Walter Roberson
on 13 Dec 2020
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
on 13 Dec 2020
Walter Roberson
on 13 Dec 2020
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);
Walter Roberson
on 13 Dec 2020
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
on 13 Dec 2020
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
on 13 Dec 2020
petit
on 13 Dec 2020
Walter Roberson
on 13 Dec 2020
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
on 21 Dec 2020
Walter Roberson
on 21 Dec 2020
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
on 21 Dec 2020
petit
on 31 Dec 2020
Walter Roberson
on 31 Dec 2020
I received my new drive yesterday; I will be trying to install it today. (It might take a while.)
petit
on 31 Dec 2020
petit
on 6 Jan 2021
Walter Roberson
on 6 Jan 2021
I have been testing my new hard drive system... which has some quirks.
petit
on 16 Jan 2021
Walter Roberson
on 16 Jan 2021
Edited: Walter Roberson
on 16 Jan 2021
You would be amazed how complicated setting up hard drives can be. Seriously.
petit
on 21 Jan 2021
Categories
Find more on Surrogate Optimization in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
