Matlab R2015b, алгоритм sqp.

[vladimir-2013]

Есть функция, зависящая от нескольких переменных. Если бы была одна переменная, то с точностью до множителей она бы выглядела так: $ \left(1 - \frac{x^2}{1600}\right) \cdot \left(0.6\cos(\pi x^2) + 0.1\cos(2\pi x^2) + 0.3\cos(3\pi x^2)\right) \cdot \cos(2x)$

График выглядит, как упавшая ёлочка, изображение лучшего качества вставить не получилось, но оно здесь особой роли не играет.


Минимум функции на сегменте [-40, 40] достигается в нуле. Точку старта поиска я специально выбрал далеко, около $30$, чтобы между ней и глобальным минимумом было много глубоких минимумов. В случае одной переменной минимумов несколько сотен, а в многомерном пространстве их число возрастает по степенному закону. В частности, в случае 6 измерений их должно быть свыше триллиона, $10^{12}$. Почему алгоритм за такое небольшое ($10^4$) число итераций так хорошо перепрыгивает локальные минимумы и находит глобальный? Вопрос возник в связи с тем, что на некоторых других функциях алгоритм работает гораздо хуже, их код из-за громоздкости приводить не буду, производные проверял с помощью сверки с разностными, они считаются правильно. Связано ли это с тем, что у первые 2 множителя у этой функции обладают центральной симметрией, а также с тем, что третья часть (произведение косинусов) также с точностью до множителя обладает симметрией относительно перестановки переменных?

Используется синтаксис Matlab M

% command to run min search
% tic; [fmin, xmin] = DummyClass.findOptimumGrad(2000, 1e4, 6); toc;

classdef DummyClass
    properties
    end;

    properties ( Constant )
        M2 = 1600;
        USE_PROD_COS = 1;
    end;

    methods ( Static )
        function [f, g] = optfunc(x)
            x2 = x*x';
            par_part = 1 - (x2/(numel(x) * DummyClass.M2));
            cos_part = 0.6*cos(1*pi*x2) + 0.1*cos(2*pi*x2) + 0.3*cos(3*pi*x2);

            pro_part = 1;
            m = (1 : numel(x)) / (0.5*numel(x));
            if (DummyClass.USE_PROD_COS > 0)
                x_cos = cos(m.*x);
                pro_part = prod(x_cos);
            end;

            f = -(par_part*cos_part*pro_part);
            if (nargout == 2)
                d_par_part = -2*x/(numel(x) * DummyClass.M2);
                d_cos_part = -2*pi*x*(0.6*sin(1*pi*x2) + 0.2*sin(2*pi*x2) + 0.9*sin(3*pi*x2));

                d_pro_part = 0;
                if (DummyClass.USE_PROD_COS > 0)
                    d_pro_part = m .* sin(m .* x) .* pro_part;
                    for i = 1 : numel(x)
                        t_cos = x_cos(i);
                        if (abs(t_cos) > 0.05)
                            d_pro_part(i) = -d_pro_part(i)/t_cos;
                        else
                            d_pro_part(i) = -prod(x_cos(setdiff(1 : numel(x), i)));
                        end;
                    end;
                end;

                g = -(d_par_part*cos_part*pro_part + par_part*d_cos_part*pro_part + par_part*cos_part*d_pro_part);
            end;
        end;

        function [fmin, xmin] = findOptimumGrad(MaxIter, N, len_x)
            objFun = @(x) DummyClass.optfunc(x);
            nonlcon = [];

            lb = -2 * ones(1, len_x);
            ub = +sqrt(DummyClass.M2) * ones(1, len_x);

            options = optimset('Algorithm', 'sqp', 'GradObj', 'on', 'Hessian', 'off', 'MaxIter', MaxIter, 'MaxFunEvals', 2000, 'TolFun', 1e-12, 'TolCon', 1e-12, 'FinDiffRelStep', 1e-9, 'Display', 'off');

            xmin = lb;
            fmin = 1e10;
            for iter = 1 : N
                % try to choose bad start point
                xstart = 30 + rand(size(lb));

                [x0,fval,exitflag,output] = fmincon(objFun, xstart, [], [], [], [], lb, ub, nonlcon, options);
                if (fval < fmin)
                    fmin = fval;
                    xmin = x0;
                end;
            end;
        end;
    end;
end
 

[vladimir-2013]

Ни у кого нет никаких идей, почему при выборе далёкой от глобального минимума точки старта алгоритм легко перепрыгивает через минимумы, близкие к глобальному по значению?

[worm2]

Только по случайности.
Возможно, точка 0 в первую очередь проверяется, и запоминается, что ни в одной другой точке меньше значения не встретилось.
Насколько я знаю, SQP сам по себе — градиентный метод и не умеет искать глобальный минимум.

[vladimir-2013]

worm2 в сообщении #1661420 писал(а):

Возможно, точка 0 в первую очередь проверяется, и запоминается, что ни в одной другой точке меньше значения не встретилось.

Интересное предположение, что-то в этом есть, но оно объясняет не до конца. Я 5 раз сделал старый запуск, минимум в точке 0 нашёлся 5 раз из 5. Первый столбец это значение минимума, остальные - это координаты, при которых достигается минимум.

Используется синтаксис Matlab M

          fmin            x1            x2            x3            x4            x5            x6
   -1.00000000    0.00000013   -0.00000012   -0.00000001   -0.00000005    0.00000009    0.00000000
   -1.00000000    0.00000001   -0.00000038    0.00000014    0.00000007    0.00000020    0.00000002
   -1.00000000    0.00000089    0.00000091   -0.00000176   -0.00000007   -0.00000209    0.00000020
   -1.00000000    0.00000075   -0.00000053    0.00000052    0.00000007   -0.00000024   -0.00000032
   -1.00000000    0.00000022    0.00000015    0.00000002    0.00000015   -0.00000006   -0.00000019
 

А потом переместил минимум в точку [1, 2, 3, 4, 5, 6] и запустил ещё 5 раз, на этот раз минимум нашёлся 3 раза из 5

Используется синтаксис Matlab M

          fmin            x1            x2            x3            x4            x5            x6
   -0.99937500    2.38943238    0.53690395    3.58468955    4.97886423    4.54213486    6.64741579
   -0.99979167    0.71618837    1.72213225    2.46175129    4.65797176    4.60918504    6.98329181
   -1.00000000    1.00008110    2.00007829    3.00007779    3.99997902    4.99987654    6.00004164
   -1.00000000    1.00002309    1.99996780    3.00000154    4.00001327    5.00005072    5.99998861
   -1.00000000    1.00007335    1.99993667    2.99999294    4.00013785    5.00013431    6.00010122

worm2 в сообщении #1661420 писал(а):

Насколько я знаю, SQP сам по себе — градиентный метод и не умеет искать глобальный минимум.

А какие алгоритмы Matlab лучше подходят для поиска глобального минимума?

[worm2]

Я не знаю Матлаба. Но здравый смысл подсказывает, что такие методы в любом случае должны основываться на каком-то переборе, чтобы найти наиболее глубокую яму, а потом в ней градиентными методами резвиться. И чем больше размерность, тем больше (экспоненциально) должен быть перебор.

[vladimir-2013]

worm2 в сообщении #1661448 писал(а):

Но здравый смысл подсказывает, что такие методы в любом случае должны основываться на каком-то переборе, чтобы найти наиболее глубокую яму, а потом в ней градиентными методами резвиться.

Я делал примерно так же: много раз случайно выбирал начальную точку и перезапускал алгоритм.

worm2 в сообщении #1661448 писал(а):

И чем больше размерность, тем больше (экспоненциально) должен быть перебор.

Согласен, звучит для таких методов логично, но в описанном случае поиск почему-то происходит значительно быстрее, для 6 измерений достаточно $10^4$, а не $10^{12}$ итераций, стало интересно, почему.