Минимизация с плотным структурированным Гессианом, линейными равенствами

Функция умножения Гессиана для меньшей памяти

fmincon interior-point и trust-region-reflective алгоритмы и fminunc trust-region алгоритм может решить задачи, где Гессиан плотен, но структурирован. Для этих задач fmincon и fminunc не вычисляйте H*Y с Hessian H непосредственно, потому что формирование H было бы интенсивным. Вместо этого вы должны предоставить fmincon или fminunc с функцией, которая, задавая матричное Y и информацию о H, вычисляет W = H*Y.

В этом примере целевая функция является нелинейной, и существуют линейные равенства, поэтому fmincon используется. Описание относится к отражающему алгоритму доверительной области; fminunc trust-region алгоритм аналогичен. Для алгоритма внутренней точки смотрите HessianMultiplyFcn опция в Функция умножения Гессиана. Целевая функция имеет структуру

f(x)=f^(x)12xTVVTx,

где V - матрица 1000 на 2. Гессиан f плотен, но Гессиан из f^ является разреженным. Если Гессиан из f^ является H^, тогда H, Гессиан f,

H=H^VVT.

Чтобы избежать чрезмерного использования памяти, которое может произойти при непосредственной работе с H, в примере предусмотрена функция умножения Гессиана, hmfleq1. Эта функция, когда передана матрица Y, использует разреженные матрицы Hinfo, что соответствует H^, и V для вычисления матричного продукта Гессия

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

В этом примере функции умножения Гессиана требуется H^ и V для вычисления матричного продукта Гессия. V является константой, поэтому можно захватывать V в указателе на функцию в анонимную функцию.

Однако, H^ не является константой и должна вычисляться в текущей x. Вы можете сделать это, вычисляя H^ в целевой функции и возврат H^ как Hinfo в третьем выходном аргументе. При помощи optimoptions для установки 'Hessian' опции для 'on', fmincon знает, чтобы получить Hinfo значение от целевой функции и передать его в функцию умножения Гессиана hmfleq1.

Шаг 1: Напишите файл brownvv.m, который вычисляет целевую функцию, градиент и разреженную часть Гессиана.

Пример проходит brownvv на fmincon как целевая функция. The brownvv.m файл длинный и не включен сюда. Просмотреть код можно с помощью команды

type brownvv

Потому что brownvv вычисляет градиент, а также целевую функцию, пример (Шаг 3) использует optimoptions для установки SpecifyObjectiveGradient опция для true.

Шаг 2: Напишите функцию, чтобы вычислить Гессиано-матричные продуктов для H, заданной матрицы Y.

Теперь задайте функцию hmfleq1 который использует Hinfo, который вычисляется в brownvv, и V, который можно захватить в указатель на функцию в анонимную функцию, чтобы вычислить матрицу Гессиана продукта W где       W = H*Y = (Hinfo - V*V')*Y. Эта функция должна иметь вид

W = hmfleq1(Hinfo,Y)

Первый аргумент должен быть таким же, как и третий аргумент, возвращенный целевой функцией brownvv. Вторым аргументом функции умножения Гессиана является матрица Y (из W = H*Y).

Потому что fmincon ожидает второго аргумента Y для формирования матричного продукта Гессия, Y всегда является матрицей с n строки, где n - количество размерностей в задаче. Количество столбцов в Y может варьироваться. Наконец, можно использовать указатель на функцию для анонимной функции для захвата V, поэтому V может быть третьим аргументом для '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);

Примечание

Функция hmfleq1 доступно в optimdemos папка как hmfleq1.m.

Шаг 3: Вызовите нелинейную стандартную программу минимизации с начальной точкой и линейными ограничениями равенства.

Загрузите параметр задачи, Vи разреженные матрицы ограничений равенства, Aeq и beq, от fleq1.mat, который доступен в optimdemos папка. Использование optimoptions для установки SpecifyObjectiveGradient опция для true и установить HessianMultiplyFcn опция указателя на функцию, которая указывает на hmfleq1. Функции fmincon с целевой функцией brownvv и с V как дополнительный параметр:

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);

Чтобы запустить предыдущий код, введите

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

Поскольку итеративное отображение было задано с помощью 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.

Сходимость является быстрой для задачи такого размера с незначительным увеличением стоимости итерации PCG по мере прогрессирования оптимизации. Выполнимость ограничений равенства сохраняется в решении.

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

Предварительное создание условий

В этом примере fmincon невозможно использовать H для вычисления средства предварительной подготовки из-за H существует только неявно. Вместо H, fmincon использует Hinfo, третий аргумент, возвращенный brownvv, для вычисления средства предварительной подготовки. Hinfo это хороший выбор, потому что он такого же размера, как и H и аппроксимирует H до некоторой степени. Если Hinfo не совпадали с размером H, fmincon вычисляет предварительный кондиционер на основе некоторых диагональных матриц масштабирования, определенных из алгоритма. Как правило, это не так хорошо.