Resolver problemas no lineales con muchas variables
Este ejemplo muestra cómo manejar un gran número de variables en un problema no lineal. En general, para este tipo de problema:
Utilice la aproximación de matriz hessiana BFGS de memoria baja (LBFGS). Esta opción está disponible en los algoritmos predeterminados
fminunc
yfmincon
.Si tiene un gradiente explícito, también puede utilizar una matriz hessiana por diferencias finitas y el algoritmo del subproblema
'cg'
.Si tiene una matriz hessiana explícita, formule la matriz hessiana como dispersa.
Aunque no forme parte de este ejemplo, si tiene un problema estructurado y puede evaluar el producto de la matriz hessiana con un vector arbitrario, pruebe a utilizar una función de multiplicación de matriz hessiana. Consulte Minimization with Dense Structured Hessian, Linear Equalities.
El ejemplo utiliza la función auxiliar hfminunc0obj
que aparece al final de este ejemplo para los solvers no lineales generales fminunc
y fmincon
. Esta función es una generalización N-dimensional de la función de Rosenbrock, una función difícil de minimizar numéricamente. El valor mínimo de 0 se alcanza en el punto único x = ones(N,1)
.
La función es una suma de cuadrados explícita. Por lo tanto, el ejemplo también muestra la eficiencia de utilizar un solver de mínimos cuadrados. Para el solver de mínimos cuadrados lsqnonlin
, el ejemplo utiliza la función auxiliar hlsqnonlin0obj
que aparece al final de este ejemplo como una función objetivo de vector que es equivalente a la función hfminunc0obj
.
Configuración de problemas
Establezca el problema de modo que utilice la función objetivo hfminunc0obj
para 1000 variables. Establezca el punto inicial x0
en -2 para cada variable.
fun = @hfminunc0obj; N = 1e3; x0 = -2*ones(N,1);
Para la opción inicial, especifique que no se muestre ni se limite el número de evaluaciones o iteraciones de funciones.
options = optimoptions("fminunc",Display="none",... MaxFunctionEvaluations=Inf,MaxIterations=Inf);
Configure una tabla para almacenar los datos. Especifique etiquetas para las ejecuciones de ocho solvers y defina lugares para recopilar el tiempo de ejecución, el valor de función devuelto, el indicador de salida, el número de iteraciones y el tiempo por iteración.
ExperimentLabels = ["BFGS_NoGrad", "LBFGS_NoGrad",... "BFGS_Grad", "LBFGS_Grad", "Analytic", "fin-diff-grads",... "LSQ_NoJacob", "LSQ_Jacob"]; timetable = table('Size',[8 5],'VariableTypes',["double" "double" "double" "double" "double"],... 'VariableNames',["Time" "Fval" "Eflag" "Iters" "TimePerIter"],... 'RowNames',ExperimentLabels);
Las secciones de código siguientes crean las opciones apropiadas para cada ejecución de solver y recopilan la salida directamente en la tabla cuando sea posible.
Aproximación de matriz hessiana BFGS, sin gradiente
opts = options; opts.HessianApproximation = 'bfgs'; opts.SpecifyObjectiveGradient = false; overallTime = tic; [~,timetable.Fval("BFGS_NoGrad"),timetable.Eflag("BFGS_NoGrad"),output] =... fminunc(fun, x0, opts); timetable.Time("BFGS_NoGrad") = toc(overallTime); timetable.Iters("BFGS_NoGrad") = output.iterations; timetable.TimePerIter("BFGS_NoGrad") =... timetable.Time("BFGS_NoGrad")/timetable.Iters("BFGS_NoGrad");
Aproximación de matriz hessiana LBFGS, sin gradiente
opts = options; opts.HessianApproximation = 'lbfgs'; opts.SpecifyObjectiveGradient = false; overallTime = tic; [~,timetable.Fval("LBFGS_NoGrad"),timetable.Eflag("LBFGS_NoGrad"),output] =... fminunc(fun, x0, opts); timetable.Time("LBFGS_NoGrad") = toc(overallTime); timetable.Iters("LBFGS_NoGrad") = output.iterations; timetable.TimePerIter("LBFGS_NoGrad") =... timetable.Time("LBFGS_NoGrad")/timetable.Iters("LBFGS_NoGrad");
BFGS con gradiente
opts = options; opts.HessianApproximation = 'bfgs'; opts.SpecifyObjectiveGradient = true; overallTime = tic; [~,timetable.Fval("BFGS_Grad"),timetable.Eflag("BFGS_Grad"),output] =... fminunc(fun, x0, opts); timetable.Time("BFGS_Grad") = toc(overallTime); timetable.Iters("BFGS_Grad") = output.iterations; timetable.TimePerIter("BFGS_Grad") =... timetable.Time("BFGS_Grad")/timetable.Iters("BFGS_Grad");
LBFGS con gradiente
opts = options; opts.HessianApproximation = 'lbfgs'; opts.SpecifyObjectiveGradient = true; overallTime = tic; [~,timetable.Fval("LBFGS_Grad"),timetable.Eflag("LBFGS_Grad"),output] =... fminunc(fun, x0, opts); timetable.Time("LBFGS_Grad") = toc(overallTime); timetable.Iters("LBFGS_Grad") = output.iterations; timetable.TimePerIter("LBFGS_Grad") =... timetable.Time("LBFGS_Grad")/timetable.Iters("LBFGS_Grad");
Matriz hessiana analítica, algoritmo 'trust-region'
opts = options; opts.Algorithm = 'trust-region'; opts.SpecifyObjectiveGradient = true; opts.HessianFcn = "objective"; overallTime = tic; [~,timetable.Fval("Analytic"),timetable.Eflag("Analytic"),output] =... fminunc(fun, x0, opts); timetable.Time("Analytic") = toc(overallTime); timetable.Iters("Analytic") = output.iterations; timetable.TimePerIter("Analytic") =... timetable.Time("Analytic")/timetable.Iters("Analytic");
Matriz hessiana por diferencias finitas con gradiente, solver fmincon
opts = optimoptions("fmincon","SpecifyObjectiveGradient",true,... "Display","none","HessianApproximation","finite-difference",... SubproblemAlgorithm="cg",MaxFunctionEvaluations=Inf,MaxIterations=Inf); overallTime = tic; [~,timetable.Fval("fin-diff-grads"),timetable.Eflag("fin-diff-grads"),output] =... fmincon(fun, x0, [],[],[],[],[],[],[],opts); timetable.Time("fin-diff-grads") = toc(overallTime); timetable.Iters("fin-diff-grads") = output.iterations; timetable.TimePerIter("fin-diff-grads") =... timetable.Time("fin-diff-grads")/timetable.Iters("fin-diff-grads");
Mínimos cuadrados, sin matriz jacobiana
lsqopts = optimoptions("lsqnonlin","Display","none",... "MaxFunctionEvaluations",Inf,"MaxIterations",Inf); fun = @hlsqnonlin0obj; overallTime = tic; [~,timetable.Fval("LSQ_NoJacob"),~,timetable.Eflag("LSQ_NoJacob"),output] =... lsqnonlin(fun, x0, [],[],lsqopts); timetable.Time("LSQ_NoJacob") = toc(overallTime); timetable.Iters("LSQ_NoJacob") = output.iterations; timetable.TimePerIter("LSQ_NoJacob") =... timetable.Time("LSQ_NoJacob")/timetable.Iters("LSQ_NoJacob");
Mínimos cuadrados con matriz jacobiana
lsqopts.SpecifyObjectiveGradient = true; overallTime = tic; [~,timetable.Fval("LSQ_Jacob"),~,timetable.Eflag("LSQ_Jacob"),output] =... lsqnonlin(fun, x0, [],[],lsqopts); timetable.Time("LSQ_Jacob") = toc(overallTime); timetable.Iters("LSQ_Jacob") = output.iterations; timetable.TimePerIter("LSQ_Jacob") =... timetable.Time("LSQ_Jacob")/timetable.Iters("LSQ_Jacob");
Estudiar los resultados
disp(timetable)
Time Fval Eflag Iters TimePerIter ______ __________ _____ _____ ___________ BFGS_NoGrad 110.44 5.0083e-08 1 7137 0.015475 LBFGS_NoGrad 53.143 2.476e-07 1 4902 0.010841 BFGS_Grad 35.491 2.9865e-08 1 7105 0.0049952 LBFGS_Grad 1.2056 9.7505e-08 1 4907 0.0002457 Analytic 7.0991 1.671e-10 3 2301 0.0030852 fin-diff-grads 5.217 1.1422e-15 1 1382 0.003775 LSQ_NoJacob 94.708 3.7969e-25 1 1664 0.056916 LSQ_Jacob 6.5225 3.0056e-25 1 1664 0.0039197
Los resultados de tiempo muestran lo siguiente:
Para este problema, la aproximación de matriz hessiana LBFGS con gradiente es la más rápida hasta ahora.
Las siguientes ejecuciones de solver más rápidas son
fmincon
con una diferencia finita de gradientes de matriz hessiana, región de confianzafminunc
con matriz hessiana y gradiente analítico ylsqnonlin
con jacobiana analítica.El algoritmo BFGS
lsqnonlin
sin gradiente tiene una velocidad similar al solverfminunc
sin jacobiana. Tenga en cuenta quelsqnonlin
requiere muchas menos iteraciones quefminunc
para este problema, pero cada iteración dura mucho más.Las derivadas suponen una gran diferencia de velocidad para todos los solvers.
Funciones auxiliares
El siguiente código crea la función auxiliar hfminunc0obj
.
function [f,G,H] = hfminunc0obj(x) % Rosenbrock function in N dimensions N = numel(x); xx = x(1:N-1); xx_plus = x(2:N); f_vec = 100*(xx.^2 - xx_plus).^2 + (xx - 1).^2; f = sum(f_vec); if (nargout >= 2) % Gradient G = zeros(N,1); for k = 1:N if (k == 1) G(k) = 2*(x(k)-1) + 400*x(k)*(x(k)^2-x(k+1)); elseif (k == N) G(k) = 200*x(k) - 200*x(k-1)^2; else G(k) = 202*x(k) - 200*x(k-1)^2 - 400*x(k)*(x(k+1) - x(k)^2) - 2; end end if nargout >= 3 % Hessian H = spalloc(N,N,3*N); for i = 2:(N-1) H(i,i) = 202 + 1200*x(i)^2 - 400*x(i+1); H(i,i-1) = -400*x(i-1); H(i,i+1) = -400*x(i); end H(1,1) = 2 + 1200*x(1)^2 - 400*x(2); H(1,2) = -400*x(1); H(N,N) = 200; H(N,N-1) = -400*x(N-1); end end end
El siguiente código crea la función auxiliar hlsqnonlin0obj
.
function [f,G] = hlsqnonlin0obj(x) % Rosenbrock function in N dimensions N = numel(x); xx = x(1:N-1); xx_plus = x(2:N); f_vec = [10*(xx.^2 - xx_plus), (xx - 1)]; f = reshape(f_vec',[],1); % Vector of length 2*(N-1) % Jacobian if (nargout >= 2) G = spalloc(2*(N-1),N,3*N); % Jacobian size 2*(N-1)-by-N with 3*N nonzeros for k = 1:(N-1) G(2*k-1,k) = 10*2*x(k); G(2*k-1,k+1) = -10; G(2*k,k) = 1; end end end