Resolver problemas no lineales con muchas variables

Abrir script en vivo

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 y fmincon.
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 confianza fminunc con matriz hessiana y gradiente analítico y lsqnonlin con jacobiana analítica.
El algoritmo BFGS lsqnonlin sin gradiente tiene una velocidad similar al solver fminunc sin jacobiana. Tenga en cuenta que lsqnonlin requiere muchas menos iteraciones que fminunc 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