Problem with using spmd in the objective function of fmincon

I am experiencing trouble with using spmd in the objective function of fmincon. My code is as follows:

load('outervariables.mat');
theta1 = zeros(size(outstruct.x1, 2),1);
%% load data to workers here. 
spmd
    data = load(['dataworker' num2str(labindex) '.mat']);
end
theta20 = sig.x;

%objective function:
fun = @(theta2)gmmobj(theta2, outstruct, data);
options = optimoptions('fmincon', ...
    'SpecifyObjectiveGradient',true, 'Display', 'off', ...
    'MaxFunctionEvaluations', 2, ...
     'OptimalityTolerance', 1e-14, 'FunctionTolerance', 1e-14);
[x, fval] = fmincon(fun, theta20, [], [], [], [], [], [], [], options);

The objective function gmmobj() is coded as

function [f, grad] = gmmobj(theta2, outstruct, data)
spmd
    [mktmval, jac] = meanutil_jacEQN(theta2, data.struct);
end

delta = cell2mat(mktmval(:));
jacobian = cell2mat(jac(:));
clear mktmval jac;

% the following deals with cases were the min algorithm drifts into region where the objective is not defined
if max(isnan(delta)) == 1
	f = 1e+10	   
else
	temp1 = outstruct.x1'*outstruct.IVs;
	temp2 = delta'*outstruct.IVs;
   theta1 = (temp1/outstruct.W*temp1')\(temp1/outstruct.W*temp2');
   %theta1 = (temp1*temp1')\(temp1*temp2');
   clear temp1 temp2 
	gmmresid = delta - outstruct.x1*theta1;
	temp1 = gmmresid'*outstruct.IVs;
	f = temp1/outstruct.W*temp1';
    %f = temp1*temp1';
   clear temp1
   %save gmmresid gmmresid
end

if nargout > 1
  grad = 2*jacobian'*outstruct.IVs/outstruct.W*outstruct.IVs'*gmmresid;
end

disp(['GMM objective:  ' num2str(f)])
end

When I call the gmmobj() function on its own, it runs fine. But once I start running the optimization with "fmincon", I get the following output and error:

GMM objective:  36532.679
Error using gmmobj (line 13)
Error detected on workers 1 2 3 4 5 6.

Error in @(theta2)gmmobj(theta2,outstruct,data)

Error in barrier

Error in fmincon (line 834)
    [X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] = barrier(funfcn,X,A,B,Aeq,Beq,l,u,confcn,options.HessFcn, ...

Caused by:
    Error using levenbergMarquardt (line 16)
    Objective function is returning undefined values at initial point. fsolve cannot continue.

Note that line 13 of the "gmmobj" function is where the "spmd" block starts. The gmmobj() function was obviously evaluated once, but on the second iteration, fmincon doesn't seem to recognize the spmd block. Can someone help me with making this work?

[Walter Roberson]

By the way,

if max(isnan(delta)) == 1

can be written as

if any(isnan(delta))

[Walter Roberson]

Please post meanutil_jacEQN and your .mat for testing.

[Claire Fang]

Thank you! any chance you know the source of the errors that I posted above? Walter, I have read many of your answers and comments, you seem really knowlegeable about parallel computing and optimizaiton.

[Walter Roberson]

Could you explain more why you are doing

spmd
    [mktmval, jac] = meanutil_jacEQN(theta2, data.struct);
end

You are not (in any obvious way) using different data depending upon iterations or spmd lab number, so the same calculation would seem to be carried out on all of the workers.

[Claire Fang]

Reason why I am doing this is: I need to back out a set of intermediate values based on my data and the parameter values of theta2. This set of intermediate values has over 77,000X17 values. Doing this sequentially takes a lot of time. So I broke down my data into 6 chunks, and distributed one chunk to each worker (I have 6 cores in total). That was outside of the gmmobj() function in the initial statement:

%% load data to workers here. 
spmd
    data = load(['dataworker' num2str(labindex) '.mat']);
end

These data on each work do not change once they have been distributed to the workers. However, the intermediate values at each iteration will change because of theta2, the parameters that I am searching for to minimize the objective function gmmobj(). Is that clear?

Answer 1

[Walter Roberson]

spmd
    data = load(['dataworker' num2str(labindex) '.mat']);
end

You are loading data into each worker, not into the client.

fun = @(theta2)gmmobj(theta2, outstruct, data);

There you make a reference to the composite on the client.

function [f, grad] = gmmobj(theta2, outstruct, data)
spmd
    [mktmval, jac] = meanutil_jacEQN(theta2, data.struct);
end

There you expect data to resolve to the per-worker value. data is present as a parameter because it matches the positional parameter data wired into the fun function handle.

However, what got wired in to the function handle is the composite as a whole, not a magic reference to the appropriate member of the composite for the lab. And when you have a reference to the composite as a whole, what gets received on the client is an empty composite.

Solution: you can evalin('base') to get at the local instance of the composite.

I also notice that outstruct is being imported into each worker every time you call. I suggest instead adding it using parallel.pool.Constant

load('outervariables.mat');
theta1 = zeros(size(outstruct.x1, 2),1);
%% load data to workers here. 
spmd
    data = load(['dataworker' num2str(labindex) '.mat']);
end
OS = parallel.pool.Constant(outstruct);

theta20 = sig.x;

%objective function:
fun = @gmmobj;
options = optimoptions('fmincon', ...
    'SpecifyObjectiveGradient',true, 'Display', 'off', ...
    'MaxFunctionEvaluations', 2, ...
     'OptimalityTolerance', 1e-14, 'FunctionTolerance', 1e-14);
[x, fval] = fmincon(fun, theta20, [], [], [], [], [], [], [], options);

function [f, grad] = gmmobj(theta2)
outstruct = evalin('base', 'OS.value');  %grab the parallel.pool.Constant
data = evalin('base', 'data');           %grab the local instance of the composite object

spmd
    [mktmval, jac] = meanutil_jacEQN(theta2, data.struct);
end
%etc

I tested the overall idea successfully.

[Claire Fang]

Walter, it still doesn't work on mine. I still have the following errors:

Error using gmmobj (line 14)
Error detected on workers 1 2 3 5 6.

Error in @(theta2)gmmobj(theta2,outstruct)

Error in barrier

Error in fmincon (line 834)
    [X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] = barrier(funfcn,X,A,B,Aeq,Beq,l,u,confcn,options.HessFcn, ...

Caused by:
    Error using levenbergMarquardt (line 16)
    Objective function is returning undefined values at initial point. fsolve cannot continue.

Per your request, I can post the data and the code for other functions here tomorrow. I can't do it today as the web says that I exceeded the limit of 10 daily uploads.

[Claire Fang]

Walter, here is a link to the data and complete set of code for the above program to run: https://www.dropbox.com/sh/eutyjvthycbpdyi/AAAFQwnwRuDbmVmC_CDTvxL7a?dl=0

Note once again that the gmmobj() function runs on its own. I did the 'CheckGradients' in the fmincon option, and the gradient checks ran just fine. But once fmincon starts the optimization process, the errors occur on the second function evaluation of gmmobj after having successfully evaluated it the first time.

[Claire Fang]

It turns out that there is nothing wrong with the spmd block. The fmincon() travelled to parameter values that render the nested objective function for fsolve invalid. Once I added lower and upper bounds to fmincon(), the program ran very smoothly and successfully.

Answer 2

[Matt J]

Clearly, the evaluation of this line was unsuccessful due to an error within meanutil_jacEQN

[mktmval, jac] = meanutil_jacEQN(theta2, data.struct);

My recommendation would be to replace, for debugging purposes, the spmd line with an equivalent for-loop implementation

for i=1:numel(data.struct)
    [mktmval{i}, jac{i}] = meanutil_jacEQN(theta2, data{i}.struct);
end

That way you can use the debugger to observe the cause of the error.

[Claire Fang]

Matt, the solution you suggested does not work. I get the following errors:

Error using Composite/subsref (line 17)
Composite objects only support simple subscripting.

Error in gmmobj (line 20)
    [mktmval{i}, jac{i}] = meanutil_jacEQN(theta2, data{i}.struct);

Thanks for the suggestion nonetheless.

[Matt J]

The idea is to replace the parallel spmd with an equivalent calculation that can be run serially. Please make the fixes necesary to get this working.

data=cell2mat(data(:));
for i=1:numel(data)
    [mktmval{i}, jac{i}] = meanutil_jacEQN(theta2, data(i).struct);
end

[Claire Fang]

Matt, the code runs fine serially, but I want the calculation to be parallelized because meanutil_jacEQN() takes time. Distributing the tasks to different workers saves a lot of time.

[Matt J]

If it runs fine in a serial for- loop, I would try replacing the ordinary for-loop with a parfor loop. This should avoid whatever issues that your meanutil_jacEQN function is having with spmd composite data types.

[Claire Fang]

I believe parfor sends data to each worker at every iteration. I want to preload the data onto the workers, so that when fmincon does optimization iterations, I don't need to send data over everytime. How do I do that?

[Matt J]

Using parallel.pool.Constant

https://www.mathworks.com/help/parallel-computing/parallel.pool.constant.html

[Claire Fang]

Matt, see my answer below. It turns out that there is nothing wrong with the spmd block. I just needed to set the upper and lower bounds for parameters in fmincon(). After that, the function ran very smoothly and successfully.

[Matt J]

Very strange then, that you didn'st see that in the serial version...

[Claire Fang]

yeah, you are right. I tested out the serial version by using a different sample, so I didn't see that...

Answer 3

[Claire Fang]

It turns out that there is nothing wrong with the spmd block. The fmincon() travelled to parameter values that render the nested objective function for fsolve invalid. Once I added lower and upper bounds to fmincon(), the program ran very smoothly and successfully.