Minimización con hessian estructurado denso, Equalidades lineales

Función de multiplicar de hessian para memoria inferior

Los algoritmos y el algoritmo pueden resolver problemas donde el hessian es denso pero estructurado.fminconinterior-pointtrust-region-reflectivefminunctrust-region Para estos problemas, y no computar con el hessian directamente, porque la formación sería la memoria intensiva.fminconfminuncH*YHH En su lugar, debe proporcionar o con una función que, dada una matriz e información sobre, calcula =.fminconfminuncYHWH*Y

En este ejemplo, la función objetiva es no lineal y existen ecualidades lineales por lo que se utiliza.fmincon La descripción se aplica al algoritmo de reflexión de la región de confianza; el algoritmo es similar.fminunctrust-region Para el algoritmo de punto interior, vea la opción en.HessianMultiplyFcnFunción de multiplicar de hessian La función objetiva tiene la estructura

$f (x) = \hat{f} (x) - \frac{1}{2} x^{T} V V^{T} x,$

donde hay una matriz 1000-by-2.V El Hessiano de es denso, pero el Hessiano def $\hat{f}$ es escaso. Si el Hessiano de Es $\hat{H}$ , entonces, el Hessiano de, esHf

$H = \hat{H} - V V^{T} .$

Para evitar el uso excesivo de memoria que podría suceder trabajando directamente, el ejemplo proporciona una función de multiplicar de hessian,.Hhmfleq1 Esta función, cuando se pasa una matriz, utiliza matrices dispersas, que corresponde aYHinfo , y para computar el producto matriz hessianV

W = H*Y = (Hinfo - V*V')*Y

En este ejemplo, la función de multiplicar de hessian necesita y calcular el producto de matriz hessian. es una constante, por lo que puede capturar en un identificador de función para una función anónima.VVV

Sin embargo no es una constante y debe calcularse en la corriente.x Puede hacerlo calculando en la función objetiva y devolver como en el tercer argumento de salida.Hinfo Mediante el uso de establecer las opciones a, sabe para obtener el valor de la función objetiva y pasarlo a la función de multiplicar hessian.optimoptions'Hessian''on'fminconHinfohmfleq1

Paso 1: escriba un archivo brownvv. m que calcule la función objetiva, el degradado y la parte dispersa del hessian.

El ejemplo se pasa como la función objetiva.brownvvfmincon el brownvv.m archivo es largo y no se incluye aquí. Puede ver el código con el comando

type brownvv

Dado que calcula el degradado, así como la función objetiva, el ejemplo () utiliza para establecer la opción.brownvvPaso 3 optimoptionsSpecifyObjectiveGradienttrue

Paso 2: escribir una función para calcular los productos de matriz de hessian para H dada una matriz Y.

Ahora, defina una función que utiliza, que se calcula en y, que se puede capturar en un identificador de función para una función anónima, para calcular el producto de matriz de hessian donde.hmfleq1HinfobrownvvVWW = H*Y = (Hinfo - V*V')*Y Esta función debe tener el formato

W = hmfleq1(Hinfo,Y)

El primer argumento debe ser el mismo que el tercer argumento devuelto por la función objetiva.brownvv El segundo argumento de la función de multiplicar de Hessian es la matriz (de).YW = H*Y

Porque espera que el segundo argumento que se utiliza para formar el producto de matriz de hessian, es siempre una matriz con filas donde es el número de dimensiones en el problema.fminconYYnn El número de columnas en puede variar.Y Por último, puede usar un identificador de función para una función anónima para capturar V, por lo que V puede ser el tercer argumento.'hmfleqq'

function W = hmfleq1(Hinfo,Y,V); %HMFLEQ1 Hessian-matrix product function for BROWNVV objective. %   W = hmfleq1(Hinfo,Y,V) computes W = (Hinfo-V*V')*Y %   where Hinfo is a sparse matrix computed by BROWNVV  %   and V is a 2 column matrix. W = Hinfo*Y - V*(V'*Y);

Nota

La función está disponible en la carpeta como.hmfleq1optimdemoshmfleq1.m

Paso 3: llame a una rutina de minimización no lineal con un punto inicial y restricciones de igualdad lineales.

Cargue el parámetro del problema y las matrices de restricción de igualdad dispersa y, desde, que está disponible en la carpeta.VAeqbeqfleq1.matoptimdemos Se usa para establecer la opción y establecer la opción en un identificador de función al que apunta.optimoptionsSpecifyObjectiveGradienttrueHessianMultiplyFcnhmfleq1 Llamada con función objetiva y con como parámetro adicional:fminconbrownvvV

function [fval,exitflag,output,x] = runfleq1 % RUNFLEQ1 demonstrates 'HessMult' option for FMINCON with linear % equalities.  problem = load('fleq1');   % Get V, Aeq, beq V = problem.V; Aeq = problem.Aeq; beq = problem.beq; n = 1000;             % problem dimension xstart = -ones(n,1); xstart(2:2:n,1) = ones(length(2:2:n),1); % starting point options = optimoptions(@fmincon,...     'Algorithm','trust-region-reflective',...     'SpecifyObjectiveGradient',true, ...     'HessianMultiplyFcn',@(Hinfo,Y)hmfleq1(Hinfo,Y,V),...     'Display','iter',...     'OptimalityTolerance',1e-9,...     'FunctionTolerance',1e-9); [x,fval,exitflag,output] = fmincon(@(x)brownvv(x,V),xstart,[],[],Aeq,beq,[],[], ...                                     [],options);

Para ejecutar el código anterior, escriba

[fval,exitflag,output,x] = runfleq1;

Dado que la visualización iterativa se estableció mediante, este comando genera la siguiente pantalla iterativa:optimoptions

                                Norm of      First-order   Iteration        f(x)          step          optimality   CG-iterations      0            2297.63                      1.41e+03                      1            1084.59         6.3903            578           1      2            1084.59            100            578           3      3            1084.59             25            578           0      4            1084.59           6.25            578           0      5            1047.61         1.5625            240           0      6            761.592          3.125           62.4           2      7            761.592           6.25           62.4           4      8            746.478         1.5625            163           0      9            546.578          3.125           84.1           2     10            274.311           6.25           26.9           2     11            55.6193        11.6597             40           2     12            55.6193             25             40           3     13            22.2964           6.25           26.3           0     14            -49.516           6.25             78           1     15           -93.2772         1.5625             68           1     16           -207.204          3.125           86.5           1     17           -434.162           6.25           70.7           1     18           -681.359           6.25           43.7           2     19           -681.359           6.25           43.7           4     20           -698.041         1.5625            191           0     21           -723.959          3.125            256           7     22            -751.33        0.78125            154           3     23           -793.974         1.5625           24.4           3     24           -820.831        2.51937           6.11           3     25           -823.069       0.562132           2.87           3     26           -823.237       0.196753          0.486           3     27           -823.245      0.0621202          0.386           3     28           -823.246      0.0199951           0.11           6     29           -823.246     0.00731333         0.0404           7     30           -823.246     0.00505883         0.0185           8     31           -823.246     0.00126471        0.00268           9     32           -823.246     0.00149326        0.00521           9     33           -823.246    0.000373314        0.00091           9  Local minimum possible.  fmincon stopped because the final change in function value relative to  its initial value is less than the value of the function tolerance.

La convergencia es rápida para un problema de este tamaño, ya que el costo de la iteración PCG aumenta modestamente a medida que avanza la optimización. La viabilidad de las restricciones de igualdad se mantiene en la solución.

problem = load('fleq1');   % Get V, Aeq, beq V = problem.V; Aeq = problem.Aeq; beq = problem.beq; norm(Aeq*x-beq,inf)
ans =     1.8874e-14

Preacondicionamiento

En este ejemplo, no se puede utilizar para calcular un preacondicionador porque sólo existe implícitamente.fminconHH En lugar de, utiliza, el tercer argumento devuelto por, para calcular un preacondicionador. es una buena opción porque es del mismo tamaño que y se aproxima a algún grado.HfminconHinfobrownvvHinfoHH Si no fueran del mismo tamaño que, calcularían un preacondicionador basado en algunas matrices de escalado diagonal determinadas a partir del algoritmo.HinfoHfmincon Típicamente, esto no llevaría a cabo también.