Квадратичная минимизация с плотным структурированным Гессианом

Воспользуйтесь структурированным Hessian

quadprog метод, отражающий доверительную область, может решить большие проблемы, где Гессиан плотен, но структурирован. Для этих задач quadprog не вычисляет H*Y с Гессианской H непосредственно, как это делает для проблем, отражающих доверительные области с разреженными H, потому что формирование H было бы интенсивным для памяти. Вместо этого вы должны предоставить quadprog с функцией, которая, задавая матричное Y и информацию о H, вычисляет W = H*Y.

В этом примере матрица Гессия H имеет структуру H = B + A*A' где B является разреженной симметричной матрицей 512 на 512 и A является разреженной матрицей 512 на 10, составленной из ряда плотных столбцов. Чтобы избежать чрезмерного использования памяти, которое может произойти при работе с H непосредственно потому, что H плотный, пример обеспечивает функцию умножения Гессиана, qpbox4mult. Эта функция, когда передана матрица Y, использует разреженные матрицы A и B для вычисления матричного продукта Гессия       W = H*Y = (B + A*A')*Y.

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

Вторая часть примера показывает, как затянуть TolPCG допуск для компенсации приблизительного предварительного условия вместо точного H матрица.

Шаг 1: Решите, какую часть H передать квадрограмме в качестве первого аргумента.

Либо A или B может быть передан как первый аргумент в quadprog. Пример выбирает пройти B как первый аргумент, поскольку это приводит к лучшему предварительному условию (см. Предварительное кондиционирование).

quadprog(B,f,[],[],[],[],l,u,xstart,options)

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

Теперь задайте функцию runqpbox4 это

  • Содержит вложенную функцию qpbox4mult который использует A и B для вычисления матричного продукта Гессия W, где       W = H*Y = (B + A*A')*Y. Вложенная функция должна иметь вид

    W = qpbox4mult(Hinfo,Y,...)

    Первые два аргумента Hinfo и Y требуются.

  • Загружает параметры задачи из qpbox4.mat.

  • Использование optimoptions для установки HessianMultiplyFcn опция указателя на функцию, которая указывает на qpbox4mult.

  • Вызовы quadprog с B как первый аргумент.

Первый аргумент вложенной функции qpbox4mult должен быть таким же, как и первый переданный аргумент quadprog, которая в данном случае является матрицей B.

Второй аргумент в qpbox4mult является матрицей Y (из   W = H*Y). Потому что quadprog ожидает Y для формирования матричного продукта Гессия, Y всегда является матрицей с n строки, где n - количество размерностей в задаче. Количество столбцов в Y может варьироваться. Функция qpbox4mult вложено так, что значение матрицы A происходит от внешней функции. Программное обеспечение Optimization Toolbox™ включает в себя runqpbox4.m файл.

function [fval, exitflag, output, x] = runqpbox4
%RUNQPBOX4 demonstrates 'HessianMultiplyFcn' option for QUADPROG with bounds.

problem = load('qpbox4'); % Get xstart, u, l, B, A, f
xstart = problem.xstart; u = problem.u; l = problem.l;
B = problem.B; A = problem.A; f = problem.f;
mtxmpy = @qpbox4mult; % function handle to qpbox4mult nested function

% Choose algorithm and the HessianMultiplyFcn option
options = optimoptions(@quadprog,'Algorithm','trust-region-reflective','HessianMultiplyFcn',mtxmpy);

% Pass B to qpbox4mult via the H argument. Also, B will be used in
% computing a preconditioner for PCG.
[x, fval, exitflag, output] = quadprog(B,f,[],[],[],[],l,u,xstart,options);

    function W = qpbox4mult(B,Y)
        %QPBOX4MULT Hessian matrix product with dense structured Hessian.
        %   W = qpbox4mult(B,Y) computes W = (B + A*A')*Y where
        %   INPUT:
        %       B - sparse square matrix (512 by 512)
        %       Y - vector (or matrix) to be multiplied by B + A'*A.
        %   VARIABLES from outer function runqpbox4:
        %       A - sparse matrix with 512 rows and 10 columns.
        %
        %   OUTPUT:
        %       W - The product (B + A*A')*Y.
        %

        % Order multiplies to avoid forming A*A',
        %   which is large and dense
        W = B*Y + A*(A'*Y);
    end

end

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

Чтобы вызвать квадратичную стандартную программу минимизации, содержащуюся в runqpbox4, введите

[fval,exitflag,output] = runqpbox4;

чтобы запустить предыдущий код. Затем отобразите значения для fval, exitflag, output.iterations, и output.cgiterations.

fval,exitflag,output.iterations, output.cgiterations

fval =

  -1.0538e+03


exitflag =

     3


ans =

    18


ans =

    30

После 18 итераций с общей суммой 30 итераций PCG, значение функции уменьшается до

fval
fval =
 -1.0538e+003

и оптимальность первого порядка

output.firstorderopt
ans =
    0.0043

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

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

Потому что предварительный кондиционер более аппроксимирован, чем когда H доступен явно, настраивая TolPCG может потребоваться параметр в несколько меньшем значении. Этот пример аналогичен предыдущему, но уменьшает TolPCG от 0,1 до 0,01 по умолчанию.

function [fval, exitflag, output, x] = runqpbox4prec
%RUNQPBOX4PREC demonstrates 'HessianMultiplyFcn' option for QUADPROG with bounds.

problem = load('qpbox4'); % Get xstart, u, l, B, A, f
xstart = problem.xstart; u = problem.u; l = problem.l;
B = problem.B; A = problem.A; f = problem.f;
mtxmpy = @qpbox4mult; % function handle to qpbox4mult nested function

% Choose algorithm, the HessianMultiplyFcn option, and override the TolPCG option
options = optimoptions(@quadprog,'Algorithm','trust-region-reflective',...
    'HessianMultiplyFcn',mtxmpy,'TolPCG',0.01);

% Pass B to qpbox4mult via the H argument. Also, B will be used in
% computing a preconditioner for PCG.
% A is passed as an additional argument after 'options'
[x, fval, exitflag, output] = quadprog(B,f,[],[],[],[],l,u,xstart,options);

    function W = qpbox4mult(B,Y)
        %QPBOX4MULT Hessian matrix product with dense structured Hessian.
        %   W = qpbox4mult(B,Y) computes W = (B + A*A')*Y where
        %   INPUT:
        %       B - sparse square matrix (512 by 512)
        %       Y - vector (or matrix) to be multiplied by B + A'*A.
        %   VARIABLES from outer function runqpbox4prec:
        %       A - sparse matrix with 512 rows and 10 columns.
        %
        %   OUTPUT:
        %       W - The product (B + A*A')*Y.
        %

        % Order multiplies to avoid forming A*A',
        %   which is large and dense
        W = B*Y + A*(A'*Y);
    end

end

Теперь введите

[fval,exitflag,output] = runqpbox4prec; 

чтобы запустить предыдущий код. После 18 итераций и 50 итераций PCG значение функции равняется пяти значащим цифрам

fval
fval = 
-1.0538e+003

и оптимальность первого порядка по существу одинаковая.

output.firstorderopt
ans =
    0.0043

Примечание

Уменьшение TolPCG слишком много может существенно увеличить количество итераций PCG.

Похожие темы