NLopt SLSQP отказывается от хорошего решения в пользу более старого, худшего решения

Question

Я решаю стандартную задачу оптимизации из Финансов - оптимизация портфеля. Подавляющее большинство времени NLopt возвращает разумное решение. Однако в редких случаях алгоритм SLSQP, по-видимому, выполняет итерацию к правильному решению, а затем без всякой очевидной причины он решает вернуть решение примерно через одну треть итеративного процесса, который, очевидно, является неоптимальным. Интересно, что изменение вектора начальных параметров на очень небольшое количество может решить проблему.

Мне удалось выделить относительно простой рабочий пример поведения, о котором я говорю. Извиняюсь, что цифры немного грязные. Это было лучшее, что я мог сделать. Следующий код может быть вырезан и вставлен в REPL Julia и будет запускать и печатать значения целевой функции и параметров каждый раз NLopt вызывает целевую функцию. Я вызываю процедуру оптимизации дважды. Если вы прокрутите обратно выходные данные, напечатанные приведенным ниже кодом, вы заметите, что при первом вызове подпрограмма оптимизации превращается в хорошее решение с целевым значением функции 0.0022 но тогда без видимой причины возвращается к гораздо более раннему решению, где целевая функция 0.0007и возвращает его вместо. Во второй раз, когда я вызываю функцию оптимизации, я использую немного другой начальный вектор параметров. Опять же, подпрограмма оптимизации повторяет то же самое хорошее решение, но на этот раз она возвращает хорошее решение со значением целевой функции 0.0022,

Итак, вопрос: кто-нибудь знает, почему в первом случае SLSQP отказывается от хорошего решения в пользу гораздо более бедного - всего около трети пути итеративного процесса? Если так, есть ли способ, как я могу это исправить?

#-------------------------------------------
#Load NLopt package
using NLopt
#Define objective function for the portfolio optimisation problem (maximise expected return subject to variance constraint)
function obj_func!(param::Vector{Float64}, grad::Vector{Float64}, meanVec::Vector{Float64}, covMat::Matrix{Float64})
    if length(grad) > 0
        tempGrad = meanVec - covMat * param
        for j = 1:length(grad)
            grad[j] = tempGrad[j]
        end
        println("Gradient vector = " * string(grad))
    end
    println("Parameter vector = " * string(param))
    fOut = dot(param, meanVec) - (1/2)*dot(param, covMat*param)
    println("Objective function value = " * string(fOut))
    return(fOut)
end
#Define standard equality constraint for the portfolio optimisation problem
function eq_con!(param::Vector{Float64}, grad::Vector{Float64})
    if length(grad) > 0
        for j = 1:length(grad)
            grad[j] = 1.0
        end
    end
    return(sum(param) - 1.0)
end
#Function to call the optimisation process with appropriate input parameters
function do_opt(meanVec::Vector{Float64}, covMat::Matrix{Float64}, paramInit::Vector{Float64})
    opt1 = Opt(:LD_SLSQP, length(meanVec))
    lower_bounds!(opt1, [0.0, 0.0, 0.05, 0.0, 0.0, 0.0])
    upper_bounds!(opt1, [1.0, 1.0, 1.0, 1.0, 1.0, 1.0])
    equality_constraint!(opt1, eq_con!)
    ftol_rel!(opt1, 0.000001)
    fObj = ((param, grad) -> obj_func!(param, grad, meanVec, covMat))
    max_objective!(opt1, fObj)
    (fObjOpt, paramOpt, flag) = optimize(opt1, paramInit)
    println("Returned parameter vector = " * string(paramOpt))
    println("Return objective function = " * string(fObjOpt))
end
#-------------------------------------------
#Inputs to optimisation
meanVec = [0.00238374894628471,0.0006879970888824095,0.00015027322404371585,0.0008440624572209092,-0.004949409024535505,-0.0011493778903180567]
covMat = [8.448145928621056e-5 1.9555283947528615e-5 0.0 1.7716366331331983e-5 1.5054664977783003e-5 2.1496436765051825e-6;
          1.9555283947528615e-5 0.00017068536691928327 0.0 1.4272576023325365e-5 4.2993023110905543e-5 1.047156519965148e-5;
          0.0 0.0 0.0 0.0 0.0 0.0;
          1.7716366331331983e-5 1.4272576023325365e-5 0.0 6.577888700124854e-5 3.957059294420261e-6 7.365234067319808e-6
          1.5054664977783003e-5 4.2993023110905543e-5 0.0 3.957059294420261e-6 0.0001288060347757139 6.457128839875466e-6
          2.1496436765051825e-6 1.047156519965148e-5 0.0 7.365234067319808e-6 6.457128839875466e-6 0.00010385067478418426]
paramInit = [0.0,0.9496114216578236,0.050388578342176464,0.0,0.0,0.0]

#Call the optimisation function
do_opt(meanVec, covMat, paramInit)

#Re-define initial parameters to very similar numbers
paramInit = [0.0,0.95,0.05,0.0,0.0,0.0]

#Call the optimisation function again
do_opt(meanVec, covMat, paramInit)

Примечание: я знаю, что моя ковариационная матрица является положительно-полуопределенной, а не положительно определенной. Это не источник проблемы. Я подтвердил это, изменив диагональный элемент нулевой строки на небольшое, но существенно ненулевое значение. Проблема все еще присутствует в приведенном выше примере, а также в других, которые я могу генерировать случайным образом.

Answer 1

SLSQP - это ограниченный алгоритм оптимизации. Каждый раунд должен проверять наличие наилучшего объективного значения и удовлетворение ограничений. Окончательный результат является лучшим значением при выполнении ограничений.

Распечатка значения ограничения путем изменения eq_con! чтобы:

function eq_con!(param::Vector{Float64}, grad::Vector{Float64})
    if length(grad) > 0
        for j = 1:length(grad)
            grad[j] = 1.0
        end
    end
    @show sum(param)-1.0
    return(sum(param) - 1.0)
end

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

Objective function value = 0.0007628202546187453
sum(param) - 1.0 = 0.0

Во время второго запуска все точки оценки удовлетворяют ограничению. Это объясняет поведение и показывает, что это разумно.

ДОПОЛНЕНИЕ:

Существенной проблемой, приводящей к нестабильности параметров, является точная природа ограничения равенства. Цитирование из ссылки NLopt ( http://ab-initio.mit.edu/wiki/index.php/NLopt_Reference):

Для ограничений равенства настоятельно рекомендуется небольшая положительная допустимая погрешность, чтобы позволить NLopt сходиться, даже если ограничение равенства немного отлично от нуля.

Действительно, переключение equality_constraint! вызывать do_opt в

    equality_constraint!(opt1, eq_con!,0.00000001)

Дает 0,0022 решение для обоих начальных параметров.