fminsearch error re/ matrix dimensions

Hi, I am new to Matlab, and have spent 2 days trying to understand why I get the below error message with fminsearch -- any help / hint / pointer would be GREATLY appreciated:

Error message:

Matrix dimensions must agree.

Error in mp_new23 (line 138)
            vhn_calc  = (- cs  + sum(gamma_bb  .* max(y2(:,5) - p  - vhn_m,0)) ...

Error in fminsearch (line 324)
    x(:) = xr; fxr = funfcn(x,varargin{:});

Error in mainsearch_2nd (line 429)
[solx1,fval1]  =
fminsearch('mp_new23',solx0,options,forsale,deal_inv,dtrans,rdtrans,avgage,active,A,md,...

mainsearch_2nd (line 429):

[solx1,fval1]  = fminsearch('mp_new23',solx0,options,forsale,deal_inv,dtrans,rdtrans,avgage,active,A,md,...
                            beta_vec,maxt,imax,k0,tf,rho,inv_vc,weights,opts);

mp_new23 (line 138):

            vhn_calc  = (- cs  + sum(gamma_bb  .* max(y2(:,5) - p  - vhn_m,0)) ...
                               + sum(gamma_bdd .* max(y2(:,5) - pa - vhn_m,0))) /(rho+lambda);

I also post the whole mp_new23 function:

            function [f,g] = mp_new23(b,forsale,deal_inv,dtrans,rdtrans,avgage,active,A,md,beta_vec,...
                 maxt,imax,k0,tf,rho,vc,weight,opts)

t0 = 0;

lambda      = exp(b(1));
gamma_s     = exp(b(2));
gamma_sd    = exp(b(3));
alpha_ds    = exp(b(4));
delta2      = exp(b(5)); 

zh          = exp(b(6));
zl          = exp(b(7));
cs          = exp(b(8));
thetas      = exp(b(9))/(1+exp(b(9)));
thetad      = exp(b(10))/(1+exp(b(10)));

tspan1 = t0:1:tf;
tspan1 = tspan1';

a1 = ones(length(tspan1),1);
a2 = a1;
a3 = a1;
a4 = a1;
a5 = a1;
a6 = a1;

flowx = A/tf;

mdo = zeros(tf+1,1); mlo = zeros(tf+1,1);

mho(2:tf+1,1) = flowx - mdo(2:tf+1) - mlo(2:tf+1);

mho(1) = flowx; mlo(1) = 0; mdo(1) = 0; 

mdn         = max(md - sum(mdo),0);

tspan2 = tf+t0-tspan1;
tspan2 = tspan2';

jdn_m    = 0;
vln_m    = 0;
pscrap   = 0;

delta1 = 1; delta0 = 0;

    d3 = 100; d4 = 100; d5 = 100;
    i  = 0; 

    options1 = [];
    %scraph2  = tf;

    while (d3>0.005 || d4>.01) && i<imax 

            aa1 = flipdim(a1,1); aa2 = flipdim(a2,1); aa3 = flipdim(a3,1); 
            aa4 = flipdim(a4,1); aa5 = flipdim(a5,1); aa6 = flipdim(a6,1); 

            %scraph = min(max(find(aa2==1)),tf);
            %scrapd = min(max(find(aa5==1)),tf);

            gammaprime  = gamma_sd / mdn;
            mhn         = alpha_ds / gammaprime;
            gamma       = gamma_s  / mhn;        

            mho2 = mho;
            mlo2 = mlo;

            X1 = [0 flowx]';
            
            [t1,y1]=ode45(@masses6,tspan1,X1,options1,mhn,gamma_s,gamma_sd,...
                          flowx,gamma,gammaprime,lambda,aa1,aa2,aa3,aa4,aa5,aa6,tf);

            mho = y1(:,2).*aa2; mlo = y1(:,1); mdo = max((flowx - mlo - mho).*aa5,0);

            scraph = min(max(find(mho>0)),tf);
            %scrapd = min(max(find(mdo>0)),tf);

            gamma_b  = gamma*mlo.*aa1;
            gamma_bd = gammaprime*mdo.*aa5;
            alpha_db = gammaprime*mlo.*aa1;

            mdn  = max(md - sum(mdo),0);

            gamma_bb  = flipdim(gamma_b,1);  gamma_bdd = flipdim(gamma_bd,1); 
            alpha_dbb = flipdim(alpha_db,1);

            d0 = 100; d4 = 100; 
            i2 = 0;

            vhn_calc = 0;

            while d0>0.005 && i2<imax  

            vhn_m = vhn_calc; k = k0;
            %%%%% terminal condition
            X0 = [0 0 0 0 vhn_m+pscrap pscrap jdn_m+pscrap]';

            [t2,y2]=ode45(@fun_my4,tspan2,X0,options1, lambda, gamma_s, gamma_sd, ...
                          alpha_ds, delta0, delta1, delta2, zh, zl, vhn_m, cs, pscrap, k0, rho, thetas, thetad);

            p  = (1-thetas)*(y2(:,6) - vln_m) + thetas*(y2(:,5) - vhn_m);
            pb = (1-thetad)*(y2(:,7) - jdn_m) + thetad*(y2(:,6) - vln_m);
            pa = (1-thetad)*(y2(:,7) - jdn_m) + thetad*(y2(:,5) - vhn_m);

            vhn_calc  = (- cs  + sum(gamma_bb  .* max(y2(:,5) - p  - vhn_m,0)) ...
                               + sum(gamma_bdd .* max(y2(:,5) - pa - vhn_m,0))) /(rho+lambda);

            k0        =   sum (alpha_dbb.* max(y2(:,7) - pb- jdn_m,0));

            d1 = [abs(vhn_calc./vhn_m-1),abs((k0+1)./(k+1)-1),any(isnan(p)),isnan(vhn_calc),isnan(k0)];
            
            d0 = max([d1,isnan(d1)]);
            
            d4 = min(d4,d0);
            
            if d4 == d0
                
            vhn_m1 = vhn_calc;
            
            a1 = (y2(:,2) >= y2(:,3));
            a2 = (y2(:,1) >= pscrap*ones(length(t2),1)+vhn_m1);
            a3 = (y2(:,5) - pa - vhn_m1>=0);
            a4 = (y2(:,7) - pb - jdn_m>=0);
            a5 = (y2(:,7) >= (pscrap+jdn_m)*ones(length(t2),1));
            a6 = (y2(:,5) - p - vhn_m1>=0);       
            
            p1  = p;
            pb1 = pb;
            k1  = k0;
            
            end
 
            i2=i2+1;

            end
            
            s1 = scraph;
            d2 = [mean(abs(mlo(2:s1)./mlo2(2:s1)-1)),mean(abs(mho(1:s1)./mho2(1:s1)-1)),...
                  sum(isnan(mlo)),sum(isnan(mho)),100*isempty(s1),100*isnan(s1)];
              
            d3 = min(100,max(d2));
            
            d5 = min(d5,d3);
            
            if d5 == d3 
                mho1 = mho; mhn1 = mhn; mlo1 = mlo; mdo1 = mdo; a11 = a1; k11 = k1; p11 = p1; pb11 = pb1;
            end

            i = i+1;
            
       end
        
mho = mho1; mhn = mhn1; mlo = mlo1; mdo = mdo1; p = p11; pb = pb11;  a1 = a11; k = k11;
                   
mu = lambda*(sum(mho) + mhn);

pp  = flipdim(p,1);
ppb = flipdim(pb,1);

aa1 = flipdim(a1,1);

mlo_fors = mlo.*aa1;

assets = sum(mlo+mdo+mho);     

%maxage_om = max(find(aa1>0));

%%%%%%%%%%%% simulate prices 

pr_sim = pp(1:maxt+1);
pw_sim = ppb(1:maxt+1);

p_sim = [pr_sim;pw_sim];

nl = length(p_sim);

dummy_sim = [zeros(nl/2,1);ones(nl/2,1)];

qagevec_sim = [(0:1:maxt)';(0:1:maxt)'];

x_sim = [dummy_sim,qagevec_sim];

beta0 = [beta_vec(1)-beta_vec(3); ...
         beta_vec(2) - beta_vec(1) - beta_vec(3) - beta_vec(4); ...
         beta_vec(3); ...
         - beta_vec(3) + beta_vec(4); ...
         beta_vec(5); ...
         -beta_vec(5)+beta_vec(6)];

options  = optimset('Display','off','MaxFunEvals',20000,'MaxIter',5000);
beta_sim = fminsearch('pnlfit',beta0,options,p_sim,x_sim);
beta_sim = fminsearch('pnlfit',beta_sim,options,p_sim,x_sim);
beta_sim = fminsearch('pnlfit',beta_sim,options,p_sim,x_sim);

modelfun = @(b,x) (b(1)+b(2)*x(:,1)) + (b(3)+b(4)*x(:,1)).*exp(-(b(5)+b(6)*x(:,1)).*x(:,2));

% opts                = statset('nlinfit');
% opts.Display        = 'off';
% opts.RobustWgtFun   = 'bisquare';
% opts.MaxIter        = 1000;

beta0 = beta_sim;
beta_sim = nlinfit(x_sim,p_sim,modelfun,beta0,opts);

beta_vec_sim = [beta_sim(1)+beta_sim(3);beta_sim(1)+beta_sim(2)+beta_sim(3)+beta_sim(4);beta_sim(3);beta_sim(3)+beta_sim(4);beta_sim(5);beta_sim(5)+beta_sim(6)];
 
g = zeros(1,11);
 
%%%%%%%
constr  = [(zh<zl);(vhn_calc<0);(cs>k);(mu/lambda<flowx*(tf+1)/1.1)];
constr2 = (cs>35000);

       if d5>0.025
       f = 10e30;
       else

    g(1)  = wmean((forsale./active) / (sum(mlo_fors+mdo)/assets)                            - 1,weight);
    g(2)  = wmean((deal_inv./active) / (sum(mdo)/assets)                                    - 1,weight); 
    g(3)  = wmean((rdtrans./active) / (gamma_s*sum(mlo_fors)/assets)                        - 1,weight); 
    g(4)  = wmean((dtrans./active) / (alpha_ds*sum(mdo)/assets)                             - 1,weight); 
    g(5)  = wmean(avgage./(sum(tspan1.*(mlo_fors+mdo))/sum(mlo_fors+mdo))                   - 1,weight);     
    g(6:11)  = beta_vec./beta_vec_sim                                                              - 1;

        f = g*vc*g' + 10e5*sum(constr) + 10e5*constr2;
       end

--- THANKS!!!!

[pepepe]

P.S.: Another reason why I'm puzzled is that I get the error message only after the 39th iteration. This means matrix dimensions are not an issue for the first 38 iterations - how can this be?

[dpb]

That's a whole lot of code to try to wade thru; but a quick perusal found

 scraph = min(max(find(mho>0)),tf);

I didn't see where that result is used, but the result of the find operation could lead to a variable number of results -- particularly an empty result. That or other similar logic is undoubtedly the culprit.

NB that intermediate function evaluations during search can be outside ranges of variables you normally expect.

Set a break on error in the debugger and poke around to see which variable(s) aren't the dimensions expected first. That can then help you figure out why/how they're not what are expected and to protect against the problem.

There's quite possibly/(probably?) other areas in the code a glance-over didn't find that have potential as well, but that's the kind of thing to be looking at -- a non-deterministic result.

[Matt J]

Just a note that 10 unknown variables to minimize over is a lot for fminsearch.

[pepepe]

@dpb and @Matt J: Thank you so much for your quick response.

@dpb: scraph ...

 scraph = min(max(find(mho>0)),tf);

...is actually needed to define s1. I didn't find any other redundant lines.

To use ...

dbstop

... to set a break on errors in the debugger was a great suggestion. The problem is with (and only with) the following equation in the mp_new23 function which fminsearch solves:

vhn_calc  = (- cs  + sum(gamma_bb  .* max(y2(:,5) - p  - vhn_m,0)) ...
              + sum(gamma_bdd .* max(y2(:,5) - pa - vhn_m,0))) /(rho+lambda);

I then checked each variable in this equation and traced it back to each input variable. The dimensions ARE correct. This leaves me puzzled.

I want to make to additional observations, which may help you @dpb and @Matt J to give additional suggestions.

1- When I evaluate the line "vhn_calc = ..." of the mp_new23 function in the command line ONLY, I get a result which seems correct.

2- When I run fminsearch for the function mp_new23, the first 38 iterations run without a problem, only the 39th iteration throws out the error.

I'M PUZZLED -- WHY CAN I EVALUATE THE LINE "vhn_calc = ..." IN THE COMMAND LINE BUT NOT WITHIN FMINSEARCH? WHY DOES FMINSEARCH WORK FOR THE FIRST 38 ITERATIONS? (There is nothing in the code that changes dimensions.)

Thanks again!!

[Walter Roberson]

At the time that you get the error generated, save the current workspace (the one with the assignment statement) to a .mat and attach the .mat for us to test.

[dpb]

"(There is nothing in the code that changes dimensions.)"

Well, nothing that is intended to do so, but the error pretty-much conclusively says there is -- now it may end up being something that happens inside some of the optimization code, but definitely something caused the result; that it apparently is reproducible indicates it must have something to do with the path the optimization is following.

Despite the frustration of inital debugging, that's good -- at least it is something that doesn't change behavior every time you run the code; perserverance will uncover the root cause.

"WHY CAN I EVALUATE THE LINE "vhn_calc = ..." IN THE COMMAND LINE...?"

Did you try that after the error occurred? One would suspect at that point the problem would still exist.

I don't have time to do anything extensive at the moment; one thought that might help in debugging would be to simplify the expressions to further isolate the area of the difficulty.

vhn_calc  = (-cs+sum(gamma_bb .*max(y2(:,5)-p -vhn_m,0)) + ...
                 sum(gamma_bdd.*max(y2(:,5)-pa-vhn_m,0))) /(rho+lambda);

could be written as

vhn_a=  sum(gamma_bb .*max(y2(:,5)-p -vhn_m,0));
vhn_b=  sum(gamma_bdd.*max(y2(:,5)-pa-vhn_m,0));
vhn_c=(rho+lambda);
vhn_calc=(-cs+vhn_a + vhn_b)/vhn_c;

and see which of the three pieces shows the mismatch error.

Also, is the "/" operator the correct one -- should it perhaps be "./" instead or are rho, lambda scalar?

That's another thing to inspect for; an inadvertent change there can silently change dimensions with a matrix operation where an element-by-element one was intended.

Which release are you using? There's also some enhanced automagic array expansion that occurs in some later releases (forget just when it was incorporated) that also can unexpectedly run but not be what was expected.

Still, the change with iterations makes me strongly suspect something like I suggested first -- there's a logic operation that is data dependent that is the culprit in not always returning the same number of elements...

I just ran into one in my own present problem of a missing separator so split() returned 4 cells of 2x1 and one of 1x1 and boom! when tried to put those together. Isn't anything in the code to change the dimensions, but an input data error caused the symptom. I'm betting it's something in your evaluation going on you're not expecting to be possible but happens because the optimization goes into places you didn't expect and code for (like I didn't test the input data on the presumption it would be correct).