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

[Kevsh]

Amw
Про локальные минимумы очень понятно рассказали, спасибо. К сожалению, перебрать столько вариантов не получится – решать это всё планируется на микроконтроллере, он от такого заплачет и расплавится, наверное . Стоило изначально сказать, что ограничения по скорости работы алгоритма есть (просто до этого я с таким не сталкивался – все методы, о которых я знал раньше, сходились в худшем случае за несколько сотен итераций).

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

Ну, если микроконтроллер - есть резон повыжимать скорость.

[Kevsh]

svv
Предложенный вами код находит $d$ с очень высокой точностью, спасибо. Сейчас постараюсь разобраться.

[Geen]

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

К сожалению, перебрать столько вариантов не получится

Что известно о данных? Есть ли гарантия, что представлено, хотя бы, несколько периодов? Сколько точек на период ожидается?

[Kevsh]

Geen
Гарантия, что представлено несколько периодов есть. Сколько точек на период – сложно сказать. По сути, сколько надо – столько и будет (в пределах разумного). Сейчас я отталкиваюсь от того, что Теорема Котельникова выполняется с запасом в несколько раз, то есть точек 10 будет минимум.

[Kevsh]

svv
У меня возникли вопросы по предложенному вами коду:
1)Строка 26:

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

S=zeros(N);
 

Функция zeros от константы выдаёт просто ноль. Вы точно это хотели сделать? Или имелось в виду так:

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

S=zeros(1, N);
 

2)Рассмотрим цикл:

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

for j=1:N
   A=[ones(n,1), cos(2*pi*b(j)*x), sin(2*pi*b(j)*x)];
   c=pinv(A)*y;
   S(j)=norm(y-A*c);
   cc(:,j)=c;
end
 

-В первой строке цикла вы создаёте матрицу: в первом столбце единицы, во втором столбце значения косинуса во всех временных точках для $j$-го коэффициента $b$, в третьем – то же самое, но для синуса. Для чего первый столбец – не понимаю. Остальные два, видимо, выглядят так, потому что вы использовали описанное выше преобразование, чтобы сделать задачу линейной. Но куда делись $p$ и $q$?
-Что происходит дальше мне вообще не ясно. Что это за метод и откуда эти выражения появляются? То есть я понимаю, что вы ищете какую-то псевдообратную матрицу, ищете норму, но не понимаю, зачем. Какой смысл у вектора $c$ (вернее, почему в нём в итоге оказываются те самые $p$, $q$ и даже $d$? Как я понял, вектор $S$ определяет отклонение функции от заданного значения при $j$-том $b$.
3)Тут вы находите минимальное отклонение и его номер в массиве, откуда можно вытащить $b$, которое соответствует этому отклонению:

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

[Smin,Ind]=min(S);
 

Тут вы вытаскиваете искомые коэффициенты и соответствующий им коэффициент $b$:

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

b_=b(Ind(1))
c_=cc(:,Ind(1))
 

Концовку, как мне кажется, я правильно понял. Но вот какой принцип был в цикле – не ясно.
Большое вам спасибо, ваш код и короче моего, и работает, причём поиск выполняется не за такое уж большое количество операций.

[Geen]

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

Функция zeros от константы выдаёт просто ноль.

Нет. https://www.mathworks.com/help/matlab/r ... tid=doc_ta

[svv]

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

Вы точно это хотели сделать? Или имелось в виду так:

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

S=zeros(1, N);

Вы правы, имелся в виду Ваш вариант. Плохо ещё знаю язык. Интересно, что скрипт работает и без этих двух операторов с zeros. Но я решил, что так яснее: сразу видна размерность матриц S и cc. Обе эти матрицы содержат $N$ столбцов, по числу значений $b$.

Матрица S хранит невязки, а матрица cc — коэффициенты, найденные методом наименьших квадратов. При данном $b$ невязка — одно число, а коэффициентов три. Поэтому S — матрица $1\times N$, а cc — матрица $3\times N$.

Под "невязкой" (по-Вашему, "отклонение") понимается $\sum\limits_{i=1}^n (f(x_i)-y_i)^2$, где $f(x)$ — наилучшая аппроксимирующая функция для данного $b$. По техническим причинам (теория этого не требует) в данной программе из этого выражения ещё извлекается квадратный корень, подробности ниже.

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

В первой строке цикла вы создаёте матрицу: в первом столбце единицы, во втором столбце значения косинуса во всех временных точках для $j$-го коэффициента $b$, в третьем – то же самое, но для синуса. Для чего первый столбец – не понимаю. Остальные два, видимо, выглядят так, потому что вы использовали описанное выше преобразование, чтобы сделать задачу линейной. Но куда делись $p$ и $q$?

Рассмотрим выражение, которое надо минимизировать при заданном $b$:
$\sum\limits_{i=1}^n(d+p\cos 2\pi bx_i+q\sin 2\pi bx_i-y_i)^2$
Согласно МНК, аппроксимирующая функция $f(x)$ ищется в виде линейной комбинации $m$ заданных базисных функций $u_k(x)$ с неизвестными коэффициентами $c_k$:
$f(x)=\sum\limits_{k=1}^m c_k u_k(x)=c_1 u_1(x)+c_2 u_2(x)+...+c_m u_m(x)$

Возьмём $m=3$ известных базисных функции
$u_1(x)=1,\;u_2(x)=\cos 2\pi b x,\;u_3(x)=\sin 2\pi b x$
с неизвестными коэффициентами
$c_1=d,\;c_2=p,\;c_3=q$ (вот куда они делись!)
Тогда наша аппроксимирующая функция принимает нужную форму:
$d+p\cos 2\pi bx+q\sin 2\pi bx=c_1u_1(x)+c_2u_2(x)+c_3u_3(x)$
Функция $u_1(x)=1$ в действительности является константой. Но это особенно ничего не упрощает, и я её всё равно должен ввести, чтобы обеспечить форму $f(x)=\sum\limits_{k=1}^m c_k u_k(x)$.

svv в сообщении #1592824 писал(а):

Пусть $A$ — матрица $n\times m$ с элементами $a_{ik}=u_k(x_i)$, где $i=1...n, k=1...m$.

В нашем случае это матрица $74\times 3$. Каждая её строка соответствует одной из $n=74$ точек $x_i$. Каждый из её столбцов соответствует одной из $m=3$ базисных функций $u_k(x)$. А на их пересечении стоит значение этой базисной функции $u_k$ в этой точке $x_i$.

Поскольку $u_1(x_i)=1$ при любом $x_i$, в первом столбце матрицы $A$ стоят только единицы. Во втором и третьем столбцах стоят соответственно $u_2(x_i)=\cos 2\pi b x_i$ и $u_3(x_i)=\sin 2\pi b x_i$.

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

2)Рассмотрим цикл:

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

for j=1:N
   A=[ones(n,1), cos(2*pi*b(j)*x), sin(2*pi*b(j)*x)];
   c=pinv(A)*y;
   S(j)=norm(y-A*c);
   cc(:,j)=c;
end
 

Что происходит дальше мне вообще не ясно. Что это за метод и откуда эти выражения появляются? То есть я понимаю, что вы ищете какую-то псевдообратную матрицу, ищете норму, но не понимаю, зачем. Какой смысл у вектора $c$ (вернее, почему в нём в итоге оказываются те самые $p$, $q$ и даже $d$? Как я понял, вектор $S$ определяет отклонение функции от заданного значения при $j$-том $b$.

Про вектор $S$ Вы поняли правильно.
Про вектор $c$ я уже рассказал — это вектор $3\times 1$ из неизвестных коэффициентов:
$c=\begin{bmatrix}c_1\\c_2\\c_3\end{bmatrix}=\begin{bmatrix}d\\p\\q\end{bmatrix}$
Вычисленный вектор $c$ потом с помощью оператора cc(:,j)=c становится $j$-м столбцом матрицы cc, хранящей такие векторы для всех $b$.

Вообще, MATLAB разрабатывался так, чтобы в матричных вычислениях пользователю приходилось писать как можно меньше явных циклов. Составив из $d,p,q$ вектор $c$, я могу их находить или использовать "одним махом". В противном случае пришлось бы всюду, где идёт работа с этими коэффициентами, записать что-то с участием $d$, потом что-то аналогичное с $p$, а потом аналогичное с $q$.

Итак, матрица $A$ заполнена.

svv в сообщении #1592824 писал(а):

Согласно теории, минимальную невязку обеспечивает вектор $c=(A^\top A)^{-1}A^\top y$.

Здесь написано (в пересчёте на вещественный случай), что если столбцы $A$ линейно независимы (а у нас это почти всегда так), то псевдообратная матрица $A^+=(A^\top A)^{-1}A^\top$. Значит, $c=A^+ y$. Но для псевдообратной матрицы есть функция pinv. Получаем оператор
c=pinv(A)*y;

svv в сообщении #1592824 писал(а):

Сама невязка равна
$(Ac-y)^\top (Ac-y)=y^\top y-y^\top Ac$

Это можно закодировать так: (y-A*c)'*(y-A*c). Но мне хотелось записать это короче. Я искал в MATLAB функцию, которая бы для произвольного вектора $v$ возвращала $v^T v$. Не нашёл. Зато нашлась функция norm, которая выдаёт $\sqrt{v^T v}$. Но это тоже подходит: точки минимума $v^T v$ и $\sqrt{v^T v}$ совпадают. Итак,
S(j)=norm(y-A*c);

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

3)Тут вы находите минимальное отклонение и его номер в массиве, откуда можно вытащить $b$, которое соответствует этому отклонению:

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

[Smin,Ind]=min(S);
 

Тут вы вытаскиваете искомые коэффициенты и соответствующий им коэффициент $b$:

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

b_=b(Ind(1))
c_=cc(:,Ind(1))
 

Концовку, как мне кажется, я правильно понял.

Да, тут Вы всё поняли правильно, отлично.

[Geen]

svv в сообщении #1593114 писал(а):

Но я решил, что так яснее: сразу видна размерность матриц

Так быстрее - не происходит перевыделения памяти и копирования старого в новое.

-- 09.05.2023, 07:38 --

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

Гарантия, что представлено несколько периодов есть. Сколько точек на период – сложно сказать. По сути, сколько надо – столько и будет (в пределах разумного). Сейчас я отталкиваюсь от того, что Теорема Котельникова выполняется с запасом в несколько раз, то есть точек 10 будет минимум.

Тогда, мне кажется, стоит найти грубое начальное приближение так: "ноль" найти как среднее, амплитуду - через разницу максимума и минимума, а частоту - через количество пересечений "нуля" (отбросив близкие к нему точки, скажем, процентов 10 от амплитуды; зависит от уровня шумов). А далее, по предложенному алгоритму (вместо сканирования минимизация по частоте).
Только придётся поискать алгоритм вычисления псевдообратной...

[Kevsh]

Geen
Действительно, я тестировал код на SciLab – там функция zeros работает чуть-чуть иначе. На Matlab будет квадратная матрица из нулей.

[svv]

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

Так быстрее - не происходит перевыделения памяти и копирования старого в новое.

Понятно. Я об этом тоже думал.

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

Только придётся поискать алгоритм вычисления псевдообратной...

Это не проблема. Чтобы избавиться от слишком высокоуровневых функций и приблизить код MATLAB к тому, что будет исполняться на микроконтроллере, перепишем цикл так:

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

for j=1:N
   A=[ones(n,1), cos(2*pi*b(j)*x), sin(2*pi*b(j)*x)];
%    c=pinv(A)*y;      
%    S(j)=norm(y-A*c);
   v=A'*y;
   c=inv(A'*A)*v;
   S(j)=y'*y-v'*c;
   cc(:,j)=c;
end

Пусть $v=A^\top y$. Новый код соответствует формулам
$\begin{array}{ll}c=(A^\top A)^{-1}A^\top y&=(A^\top A)^{-1}v\\S=y^\top (y-Ac)&=y^\top y-v^\top c\end{array}$
($y^\top y$ не зависит от $b$ и его можно вычислить заранее)

Матрица $A^\top A$ имеет детский размер $3\times 3$, поэтому обратную к ней матрицу можно записать в явном виде (для обращения потребуется около 20 умножений и одно деление). Нужно только следить, чтобы $A^\top A$ не была вырожденной. Это случится, если $b$ будет полуцелым или целым.

[Kevsh]

svv
Получается, что в предложенном вами алгоритме начальное приближение вообще ни на что не влияет? Как я понял, тут даже каких-то последовательных приближений нет. Тут, видимо, стоит выбирать не правильное начальное приближение, а шаг для $b$.

Также всё-таки не очень понятно, откуда такая формула для невязки. Это стандартный МНК? Я просто первый раз в жизни увидел псевдообратную матрицу, в МНК (по крайней мере у нас в университете) такого точно не было. Где можно почитать про это? Допустим, вы в одном сообщении писали, что $(Ac-y)^\top (Ac-y)=y^\top y-y^\top Ac$, а в коде $Ac$ и $y$ уже поменялись местами (в скобках). Хотелось бы понять, почему так можно делать.

А так – всё понятно, ещё раз большое спасибо. В скором времени попытаюсь оценить, сколько в худшем случае операций нужно будет сделать, чтобы применить этот алгоритм.

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

Начальное приближение для перебора по сетке это задание максимальной и минимальной точки поиска и шага сетки.

[svv]

Тут MATLAB подсказывает мне, что вместо обращения матрицы проще и эффективнее находить $c$ как решение системы уравнений
$(A^\top A)c=v$
Кодируется это так:
c=(A'*A)\v;
Ну, не знаю... матрица $A^\top A$ симметричная, обратная к ней тоже, и при обращении вручную это можно использовать. Оба варианта очень просты для реализации на микроконтроллере из-за крохотного размера матрицы.

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

Получается, что в предложенном вами алгоритме начальное приближение вообще ни на что не влияет? Как я понял, тут даже каких-то последовательных приближений нет. Тут, видимо, стоит выбирать не правильное начальное приближение, а шаг для $b$.

Да, в МНК нет последовательных приближений. Действительно, надо правильно выбрать шаг для $b$, а также интервал изменения $b$, на котором ищется минимум.

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

Также всё-таки не очень понятно, откуда такая формула для невязки. Это стандартный МНК?
...
Допустим, вы в одном сообщении писали, что $(Ac-y)^\top (Ac-y)=y^\top y-y^\top Ac$, а в коде $Ac$ и $y$ уже поменялись местами (в скобках). Хотелось бы понять, почему так можно делать.

Да, это стандартный МНК. Из системы
$A^\top A c=A^\top y$
получается формула для наилучших коэффициентов (см. Вики)
$c=(A^\top A)^{-1}A^\top y=(A^\top A)^{-1}v$ (где $v=A^\top y$)
Систему также можно записать в виде (это пригодится ниже)
$A^\top (y-Ac)=0\qquad(*)$

Теперь найдём невязку. Писать ли $S=(y-Ac)^\top(y-Ac)$, или $S=(Ac-y)^\top (Ac-y)$ — всё равно: $Ac-y=-(y-Ac)$, поэтому во второй формуле два минуса компенсируются и результат остаётся тем же. Но Вы правы: я довольно беспорядочно чередую первый и второй вариант.
Остановимся на первом варианте и раскроем скобки:
$S=y^\top (y-Ac) - c^\top A^\top(y-Ac)$
В силу $(*)$ второе слагаемое равно нулю, и остаётся
$S=y^\top (y-Ac)=y^\top y-y^\top Ac = y^\top y-v^\top c$
Здесь удобно то, что $y^\top y$ достаточно вычислить один раз, а $v$ и $c$ уже найдены.

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

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

Это понятие полезное, но пока читать про это не надо, это в вычислении никак не поможет, только отвлечёт. Сейчас достаточно знать, что
1) выражение $(A^\top A)^{-1}A^\top$ равно матрице, псевдообратной для $A$, если только $A^\top A$ обратима;
2) для вычисления псевдообратной матрицы в MATLAB есть готовая функция.

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

В скором времени попытаюсь оценить, сколько в худшем случае операций нужно будет сделать, чтобы применить этот алгоритм.

Можете опираться на мои оценки, сделанные здесь. У Вас может получиться и меньше операций, но не должно быть больше.

[Amw]

svv в сообщении #1593196 писал(а):

для вычисления псевдообратной матрицы в MATLAB есть готовая функция.

Не знаю как в MATLAB, не пользовался, а в MATHCAD есть готовая функция "нахождение минимума функции" - Minimize называется.
(... у нас управдом - друг человека.)