exponenta event banner

Минимизация с помощью плотных структурированных гессенских линейных уравнений

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

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 ^, то Н, гессен 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 в качестве целевой функции. 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 вычисляет устройство предварительного кондиционирования на основе некоторых матриц диагонального масштабирования, определенных из алгоритма. Обычно это не так хорошо.