Solving for different types of variables on large dataset with lsqnonlin

Hello!

I am currently working on some code (see below) in order to solve for two unknown variables: T and a. T and a are derived from three n x 1 datasets: Data_1, Data_2 and Data_3, where T is n x 1 and a is a scalar. Basically, I want to find the T for each row across the three datasets so that f(Data_1(i),Data_2(i))=Data_3(i) based either on measured Data_1 and Data_2, or calculated values for them (how I intend to find T and a). a is a scaling factor that is unknown, but is the same across all points of the measured values of Data_1 and Data_2. The code below works in theory, but in practise, the data vectors are approx 3M x 1 large, and so lsqnonlin wants to make a transformation matrix approx 100TB in size... In an older version of the script where I wasn´t also trying to solve for a, I got around this by feeding smaller sections of the data into the solver, but since a needs to be a single optimal solution that holds for all datapoints, this is no longer an option. I need to use a solver (or options for lsqnonlin) that are more memory friendly.

I know solvers like fmincon are less memory intensive, but after trying for a bit, I don't think it actually is appropriate for this problem, even if rewritten as an optimisation problem rather than an equations problem. At the very least, it does not seem to handle multiple unknown variables of different sizes well. Meanwhile, lsqnonlin should be able to do what I need, but I am unsure of it is technically feasible.

Is what I am trying to do feasible, and if so, any suggestions on how to do it?

Let me know if you need any more info. Thanks!

function [T_values] = calc_T(Data_1,Data_2,Data_3,PARAMS,M_0,B_1)

T = optimvar('T',length(Data_1),LowerBound=0,UpperBound=4);

x0.T=ones(length(Data_1),1).*2;

a = optimvar('a',LowerBound=0,UpperBound=1);

x0.a=sind(max(PARAMS.alpha));

% Define mathematical functions

N_slices_pf =192*[PARAMS.PF-0.5 0.5];

N_echoes_bef_T_E = N_slices_pf(1);

N_echoes_aft_T_E = N_slices_pf(2);

N_echoes_tot = sum(N_slices_pf);

T_acq = N_echoes_tot * PARAMS.TE;

T_acq_bef_T_E = N_echoes_bef_T_E * PARAMS.TE;

T_acq_aft_T_E = N_echoes_aft_T_E * PARAMS.TE;

E_R = exp(-PARAMS.TE./T);

ca1 = cos(PARAMS.alpha(1)*pi/180);

ca2 = cos(PARAMS.alpha(2)*pi/180);

ca1E_R = ca1.*E_R;

ca2E_R = ca2.*E_R;

T_A = PARAMS.TI(1) - T_acq_bef_T_E;

T_C = PARAMS.TR - PARAMS.TI(2) - T_acq_aft_T_E;

T_B = PARAMS.TR - 2*T_acq - T_A - T_C;

E_A = exp(-T_A./T);

E_B = exp(-T_B./T);

E_C = exp(-T_C./T);

X = M_0.*(1-E_R).*(1-ca2E_R.^N_echoes_tot)./(1-ca2E_R).*E_C + M_0.*(1-E_C);

Y = M_0.*(1-E_R).*(1-ca1E_R.^N_echoes_tot)./(1-ca1E_R).*E_B.*ca2E_R.^N_echoes_tot.*E_C + M_0.*(1-E_B).*ca2E_R.^N_echoes_tot.*E_C;

Z = M_0.*(1-E_A).*ca1E_R.^N_echoes_tot.*E_B.*ca2E_R.^N_echoes_tot.*E_C;

M_z_ss_denom = 1+PARAMS.Eff_Inv*ca1E_R.^N_echoes_tot.*ca2E_R.^N_echoes_tot.*E_A.*E_B.*E_C;

M_z_ss_numer = X+Y+Z;

M_z_ss = M_z_ss_numer./M_z_ss_denom;

Mz1 = -PARAMS.Eff_Inv.*M_z_ss;

Mz2 = Mz1.*E_A + M_0.*(1-E_A);

Mz3_1 = Mz2.*ca1E_R.^N_echoes_bef_T_E + M_0.*(1-E_R).*(1-ca1E_R.^N_echoes_bef_T_E)./(1-ca1E_R);

Mz4 = ((Data_1.*a)./sin(PARAMS.alpha(1)*pi/180)).*ca1E_R.^N_echoes_aft_T_E + M_0.*(1-E_R).*(1-ca1E_R.^N_echoes_aft_T_E)./(1-ca1E_R);

Mz5 = Mz4.*E_B+M_0.*(1-E_B);

Mz6_2 = Mz5.*ca2E_R.^N_echoes_bef_T_E + M_0.*(1-E_R).*(1-ca2E_R.^N_echoes_bef_T_E)./(1-ca2E_R);

% Define equasions

eq1 = Data_3 == (Mz3_1.*sin(PARAMS.alpha(1).*pi./180).*Mz6_2.*sin(PARAMS.alpha(2).*pi./180))./((Mz3_1.*sin(PARAMS.alpha(1).*pi./180)).^2+(Mz6_2.*sin(PARAMS.alpha(2).*pi./180)).^2);

eq2 = Data_3 == (Data_1.*(((Mz6_2.*sin(PARAMS.alpha(2).*pi./180)))./a)) ./ ((Data_1).^2+((Mz6_2.*sin(PARAMS.alpha(2).*pi./180))./a).^2);

eq3 = Data_3 == (((Mz3_1.*sin(PARAMS.alpha(1).*pi./180))./a).*Data_2)./((((Mz3_1.*sin(PARAMS.alpha(1).*pi./180)))./a).^2+(Data_2).^2);

prob = eqnproblem;

prob.Equations.eq1 = eq1;

prob.Equations.eq2 = eq2;

prob.Equations.eq3 = eq3;

options=optimoptions('lsqnonlin','Display','none');

sol=solve(prob,x0,'Options',options);

Tvalues=sol.T;

end

Answer 1

[Matt J]

You could imoplement a nested optimization, with the outer optimization over a done with fminsearch while that over T can be done as you did in the old days before the unknown a was introduced,

a_optimal = fminsearch(@(a) optimize_over_T(Data,a)  ,1)

[Lenno]

Thanks, I will give it a go! While it is churning away though, how will this perform in terms of time? I imagine like fmincon it runs a certain amount of iterations (3000 for fmincon, iirc). The old script would run for about 1.5 hours at a time, so 3000 iterations of that would take over half a year ^^'.

[Matt J]

Come to think of it, fminbnd would probably be better, since a is a scalar parameter on a bounded interval

a_optimal = fminbnd(@(a) optimize_over_T(Data,a)  ,0,1,options)

As for the time, that would depenend on what you'd set for the TolX option parameter. I imagine fminbnd would be similar to bisection in the number of iterations required, maybe 10 for an uncertainy of 2^-10 on a.

You might also consider running the inner optimzations with some subset of the n Data points.

a_optimal = fminbnd(@(a) optimize_over_T(Data(subset,:),a)  ,0,1)

I doubt you need so much data just to accurately determine a. It's just a single scalar parameter. After a_optimal is determined, you can go back and recover all the T parameters for the full Data set with a fixed at a_optimal.

[Lenno]

Alright! It'll run over the weekend, if it isn't done by the time I get back in, I will see if using a subset of the data gives good results.

[Lenno]

It seems to work perfectly, thanks! And by changing the original problem-based equations problem into a solver-based optimisation problem, I also reduced the runtime of that sevenfold ^^