Численно решить уравнение с 4-мя неизвестными

[Kevsh]

Мне нужно найти минимум данной функции:
$S(a,b,c,d)=\sum\limits_n(y_n-a\sin{(2\pi bx_n+c)-d)^2}$
Набор точек $y_n$ и $x_n$ известен.
Я уже делал тему с похожей проблемой, мне посоветовали использовать метод Марквардта. С помощью него можно минимизировать каждый член суммы, но ведь тогда мы найдём $n$ наборов искомых коэффициентов, равзе не так? Можно, конечно, потом найти сумму для каждого из этих наборов и найти такой, для которого она будет минимальной, но мне кажется, что это не гарантирует правильный ответ. Какой метод тут можно использовать?

[Amw]

Kevsh в сообщении #1592364 писал(а):

Мне нужно найти минимум данной функции... ...Какой метод тут можно использовать?

Любой численный метод нахождения минимума функции

[Евгений Машеров]

Kevsh в сообщении #1592364 писал(а):

С помощью него можно минимизировать каждый член суммы, но ведь тогда мы найдём $n$ наборов искомых коэффициентов, равзе не так?

Это неверно. Минимизируется сумма квадратов, получаем единственный вектор коэффициентов модели.

[Amw]

Kevsh в сообщении #1592364 писал(а):

Какой метод тут можно использовать?

1. Покоординатный спуск
2. Метод наискорейшего спуска
3. Метод деформированного многогранника
4. Монте-Карло
5. Генетический метод
Кто больше?
Забыл - ещё прямое сканирование...

[Kevsh]

Amw
Евгений Машеров
Я нашёл в интернете алгоритм для решения такой задачи методом Марквардта.

У меня что-то получилось, но явно не то, что я хотел.
Вот так я нахожу результат выполнения функции (тут $x$ и $y$ – известные векторы одинакового размера и $vars$ – вектор с искомыми параметрами):

function sum = SquareSum(x, y, vars)
    A = vars(1); B = vars(2); C = vars(3); D = vars(4);
    sum = 0;
    for n = 1 : length(x)
        sum = sum + (y(n) - A * sin(2 * pi * B * x(n) + C) - D)^2;
    end
end

Вот так я нахожу матрицу Гессе, на функции $f(x)=x^2+y^2+z^2+d^2$ она выдаёт правильный ответ, думаю проблема кроется не тут:

function HMatrix = HessianMatrix(x, y, vars, delta)
    A = vars(1); B = vars(2); C = vars(3); D = vars(4);
    n = size(vars);
    n = n(2);
    temp = zeros(1, n);
    temp2 = temp;
    temp3 = temp;
    temp4 = temp;
    HMatrix = zeros(n, n);

    for i = 1 : n
        for j = i : n
            if i == j
                for k = 1 : n
                    if k == j
                        temp(k) = vars(k) + delta;
                        temp2(k) = vars(k) - delta;
                    else
                        temp(k) = vars(k);
                        temp2(k) = vars(k);
                    end
                end
                HMatrix(i, j) = (SquareSum(x, y, temp) - 2 * SquareSum(x, y, vars) ...
                    + SquareSum(x, y, temp2)) / (delta^2);
            else
                for k = 1 : n
                    if k == i
                        temp(k) = vars(k) + delta;
                        temp2(k) = vars(k) + delta;
                        temp3(k) = vars(k) - delta;
                        temp4(k) = vars(k) - delta;
                    elseif k == j
                        temp(k) = vars(k) + delta;
                        temp2(k) = vars(k) - delta;
                        temp3(k) = vars(k) + delta;
                        temp4(k) = vars(k) - delta;
                    else
                        temp(k) = vars(k);
                        temp2(k) = vars(k);
                        temp3(k) = vars(k);
                        temp4(k) = vars(k);
                    end
                end
                HMatrix(i, j) = (SquareSum(x, y, temp) - SquareSum(x, y, temp2) ...
                    - SquareSum(x, y, temp3) + SquareSum(x, y, temp4)) / (4 * delta^2);
                HMatrix(j, i) = HMatrix(i, j);
            end
        end
    end    
end

Вот так я нахожу градиент (на той же $f(x)=x^2+y^2+z^2+d^2$ выдаёт всё верно, так что с этой функцией скорее всего тоже всё хорошо):

function grad = Grad(x, y, vars, delta)
    A = vars(1); B = vars(2); C = vars(3); D = vars(4);
    gradA = (SquareSum(x, y, [A + delta B C D]) - ...
    SquareSum(x, y, [A - delta B C D])) / (2 * delta);

    gradB = (SquareSum(x, y, [A B + delta C D]) - ...
    SquareSum(x, y, [A B - delta C D])) / (2 * delta);

    gradC = (SquareSum(x, y, [A B C + delta D]) - ...
    SquareSum(x, y, [A B C - delta D])) / (2 * delta);

    gradD = (SquareSum(x, y, [A B C D + delta]) - ...
    SquareSum(x, y, [A B C D - delta])) / (2 * delta);

    grad = [gradA; gradB; gradC; gradD];
end

И сам алгоритм:

%МЕТОД МАРКВАРДТА
A = (max(f_swam3) - min(f_swam3)) / 2;
B = max(f_swam3) - A;
C = 0;
D = B;
VARS = [A B C D];
M = 50; %Максимальное допустимое количество итераций
eps = 0.0001; %Параметр сходимости
k = 0; %Счётчик итераций
lambda = 10^4;
x = 1 : length(f_swam3);
delta = 0.0001; %Шаг вычисления частных производных
flag = true;

while(true)
    grad = Grad(x, f_swam3, VARS, delta);
    if (norm(grad) < eps && flag)
        disp(k)
        break
    end
    if (k >= M && flag)
        break
    end
    
    H = HessianMatrix(x, f_swam3, VARS, delta);
    s = -inv(H + lambda * eye(4)) * grad;
    temp = VARS;
    VARS = VARS + s';
    
    if SquareSum(x, f_swam3, VARS) < SquareSum(x, f_swam3, temp)
        flag = true;
        lambda = lambda / 2;
        k = k + 1;
    else
        flag = false;
        lambda = lambda * 2;
    end
end

disp(VARS)
figure
hold on
grid on
plot(x, VARS(1)*sin(2*pi*VARS(2)*x + VARS(3)) + VARS(4))
plot(x, f_swam3)
legend('МНК', 'То, что было изначально')

Полученные коэффициенты:

0.0000   70.5002    3.7404   70.1103

Вместо нормальной синусоиды я получил практически прямую линию из-за нулевой амплитуды. Хотя, кстати, постоянная составляющая исходной функции была 70, так что её получилось найти относительно точно. Но почему такая амплитуда – мне вообще не ясно. Если поставить вместо $VARS(1)$ изначальную $A$, то становится видно, что и частота подобрана неверно. Что могло пойти не так?

Вот мои векторы x и f_swam3:
x:

1   2   3   4   5   6   7   8   9   10   11   12   13   14   15   16   17   18   19   20   21   22   23   24   25   26   27   28   29   30   31   32   33   34   35   36   37   38   39   40   41   42   43   44   45   46   47   48   49   50   51   52   53   54   55   56   57   58   59   60   61   62   63   64   65   66   67   68   69   70   71   72   73   74

f_swam3:

61.8345252303583   64.2186685630058   68.5780907223303   73.8351943469310   78.1178746381488   79.4447599032848   77.1209584326325   72.3384843492439   67.1780874592201   63.3083342155163   61.5848686426315   62.2811015529467   65.2944676251564   70.0570903999628   75.2539157928855   78.8647097424086   79.1332309482167   75.9175990711235   70.8145100998371   65.8890541188997   62.5767534515232   61.5348068381206   62.9195242022718   66.5179434492526   71.5768723085928   76.5420755341518   79.3273502810222   78.5249447120497   74.5575769012328   69.3099594259272   64.7368823911916   62.0334551154671   61.6848623894562   63.7418912006921   67.8660696583989   73.0933449460809   77.6481344114599   79.4840538567364   77.6481344114599   73.0933449460811   67.8660696583989   63.7418912006923   61.6848623894562   62.0334551154669   64.7368823911916   69.3099594259272   74.5575769012327   78.5249447120497   79.3273502810222   76.5420755341517   71.5768723085929   66.5179434492527   62.9195242022719   61.5348068381207   62.5767534515231   65.8890541188995   70.8145100998372   75.9175990711235   79.1332309482166   78.8647097424087   75.2539157928858   70.0570903999627   65.2944676251561   62.2811015529469   61.5848686426318   63.3083342155164   67.1780874592199   72.3384843492438   77.1209584326324   79.4447599032847   78.1178746381488   73.8351943469311   68.5780907223305   64.2186685630060

[Geen]

Kevsh в сообщении #1592439 писал(а):

Что могло пойти не так?

Локальный минимум... очень большое число локальных минимумов...

[Евгений Машеров]

Какая-то у Вас лямбда... Грандиозная. Прибавление лямбды к диагонали имеет простой прикладной смысл. Страховка на случай, если матрица необратима (а если формально обратима, но близка к вырожденной, и результат оказывает заложником ошибок вычислений, то "тормозит" увеличение элементов обратной). Попробуйте не 10000, а 1/10000.

[Евгений Машеров]

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

[Kevsh]

Евгений Машеров
При лямбде $10^{-4}$ получилось так:

 0.0166   70.5135   -4.7439   70.1103

Лучше, конечно, но всё равно не то, мягко говоря

[Amw]

Geen в сообщении #1592442 писал(а):

Локальный минимум... очень большое число локальных минимумов...

Предварительно отсканировать $100^4$

[Евгений Машеров]

Что-то этот алгоритм не совсем походит на Марквардта.

[Amw]

Евгений Машеров в сообщении #1592521 писал(а):

Что-то этот алгоритм не совсем походит на Марквардта.

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

[Евгений Машеров]

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

[Kevsh]

Amw
Можете, пожалуйста, пояснить, что значит "отсканировать". У меня в университете был очень общий курс вычислительных методов, я даже примерно не понимаю, о чём идет речь. Кстати, что за локальные минимумы я тоже не очень понимаю (точнее, при чём они тут здесь).

[svv]

Kevsh в сообщении #1592364 писал(а):

Мне нужно найти минимум данной функции:
$S(a,b,c,d)=\sum\limits_n(y_n-a\sin{(2\pi bx_n+c)-d)^2}$

Пусть $p=a\sin c, q=a\cos c$, тогда функция равна
$\sum\limits_n(y_n-d-p\cos 2\pi bx_n-q\sin 2\pi bx_n)^2$
При фиксированном $b$ задача становится линейной, и значения $d,p,q$, обеспечивающие минимум, (ну, и сам минимум) находятся методом наименьших квадратов. Тут всего три базисных функций $1, \cos 2\pi bx, \sin 2\pi bx$, так что, наверное, и явные формулы можно выписать.

Так Ваша задача сводится к задаче минимизации функции одной переменной $b$. В каких-то случаях, вероятно, можно просто построить график и посмотреть, где минимум.

Kevsh, а Ваши $x_n$ равноотстоящие? И сколько их примерно?