Разложение Фурье, нахождение гармоник, прогноз

[svv]

Schrodinger's cat
Моя вина, простите!

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

$a_1e^{\lambda_1 n}+a_2e^{\lambda_2 n},$
где
$\begin{array}{ll}a_1=-10ie^{+0.5}&\lambda_1=+i\frac{2\pi}{30}\\[0.8ex]a_2=+10ie^{-0.5}&\lambda_2=-i\frac{2\pi}{30}\end{array}$

Правильно так:
$\begin{array}{ll}a_1=-10ie^{+0.5{\color{magenta}i}}\\[0.8ex]a_2=+10ie^{-0.5{\color{magenta}i}}\end{array}$
$\lambda_1$ и $\lambda_2$ правильные.

Ошибку нашёл моментально, как только увидел, что в Ваших файлах сигнал не вещественный, а имеет сравнимую с вещественной мнимую часть. Мнимая часть сигнала должна быть нулевой с точностью до ошибок округления.

[Schrodinger's cat]

Удалось запустить метод по варианту ilghiz. ilghiz спасибо за понятное объяснение алгоритма.
Вопросы:
1) Как выбирать правый сингулярный вектор? Всегда следующий от выбранных? Допустим я выбрал 2 вектора и они описывают допустим 80% дисперсии меня это устраивает, значит для расчета я должен выбрать следующий 3тий вектор даже несмотря на то что его сингулярное число все еще имеет хоть и малый но вес и отличается от нуля. Или нужно выбирать первый около нулевой вектор?
Алгоритм выбора не до конца понятен.
2) Одна из начальных фаз получилась не $0.5$ а $2.642$ понятно, что это $\pi-2.642=0.5$ что не совсем ожидалось.
Это как то можно приводить автоматически к ожидаемым $0.5$?
3) В вычислении частот и фаз какую версию арктангенса нужно использовать? Полагал $\arctg2$, учитывающую знаки, но она дает другие значения, а $\arctg$, не учитывающая знаки, дает верные значения. Это правильно использовать $\arctg$?
В Фурье используется $\arctg2$.
4) Марпл вообще говоря пишет в отношении SVD то что это некий целый процесс, я не очень понял, какие то процедуры форвард, потом бэкворд. А не просто одно SVD разложение. Что это такое и как правильно это сделать?
5) Дальше Мапрл предлагает использовать псевдообратную матрицу Мура-Пенроуза.
Тут я вообще не понял для чего она используется.

А вот классический метод не взлетает, что бы я не делал.
Видимо еще какие то ошибки.
Посмотрите пожалуйста.

Листинги программ
https://disk.yandex.ru/d/TNfQNemTGa9vig

[ilghiz]

Schrodinger's cat в сообщении #1565940 писал(а):

Удалось запустить метод по варианту ilghiz. ilghiz спасибо за понятное объяснение алгоритма.

Рад, что у вас стало получаться!

Schrodinger's cat в сообщении #1565940 писал(а):

1) Как выбирать правый сингулярный вектор? Всегда следующий от выбранных? Допустим я выбрал 2 вектора и они описывают допустим 80% дисперсии меня это устраивает, значит для расчета я должен выбрать следующий 3тий вектор даже несмотря на то что его сингулярное число все еще имеет хоть и малый но вес и отличается от нуля. Или нужно выбирать первый около нулевой вектор?

Это ключевая магия этого алгоритма - если вычислительные ресурсы позволяют, надо взять несколько матриц $(N-R)\times R$, $R=R_1,...,R_2$, для каждой построить сингулярное разложение и искать только его последний правый сингулярный вектор, и для каждого вектора посчитать корни своего полинома, и уже их друг с другом сравнивать.

Как только вы с числом корней переборщили, там будет или дубль, или несколько корней друг с другом замешаются и расползутся, тут надо чувствовать какой у вас шум в задаче и на этой основе принимать правильное решение.

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

Schrodinger's cat в сообщении #1565940 писал(а):

2) Одна из начальных фаз получилась не $0.5$ а $2.642$ понятно, что это $\pi-2.642=0.5$ что не совсем ожидалось.
Это как то можно приводить автоматически к ожидаемым $0.5$?
3) В вычислении частот и фаз какую версию арктангенса нужно использовать? Полагал $\arctg2$, учитывающую знаки, но она дает другие значения, а $\arctg$, не учитывающая знаки, дает верные значения. Это правильно использовать $\arctg$?
В Фурье используется $\arctg2$.

Этот алгоритм может терять знак в фазе или сдвигать на $\pi$, что вы собственно и увидели. Тут ничего не поделаешь. Но фазу постфактум можно довольно точно угадать обычными наименьшими квадратами.

Schrodinger's cat в сообщении #1565940 писал(а):

4) Марпл вообще говоря пишет в отношении SVD то что это некий целый процесс, я не очень понял, какие то процедуры форвард, потом бэкворд. А не просто одно SVD разложение. Что это такое и как правильно это сделать?
5) Дальше Мапрл предлагает использовать псевдообратную матрицу Мура-Пенроуза.
Тут я вообще не понял для чего она используется.

а вы обратите внимание, в каком году эта книга написана - в 1987. И хотя Джек Донгарра свой SVD уже в Lapack за 7 лет до этого засунул, не все еще на лапак начали молиться, и, кстати, даже в середине 90-х в SVD лапака было довольно много ошибок, поэтому каждый ваял свое SVD и не сильно доверял чужому. Ну и были такие наивные методы, которые просто матрицы сингулярного разложения итерационно получали - там куча разного было, многие даже еще в начале 2000-х этим баловались. И ведь только в начале 2000-х была таки доказана численная устойчивость SVD лапака (через дважды Хессенберговы матрицы), и только тогда все начали использовать нормальный SVD, а не какие-то итерационные мур-пенроусовские поделки.

[svv]

Schrodinger's cat в сообщении #1565940 писал(а):

А вот классический метод не взлетает, что бы я не делал.

Спокойно. Взлетит.

Составим СЛАУ для первого этапа. В матричной форме: $Ac=b$ (буква $x$ у меня зарезервирована за сигналом).
Матрица $A$ с элементами $a_{ij}$ имеет $N-p$ строк и только $p$ столбцов.
Столбец неизвестных $c_j$ имеет $p$ строк.
Столбец свободных членов $b_i$ имеет $N-p$ строк.

Система имеет вид
$\begin{bmatrix}x_0&x_1&\cdots&x_{p-1} \\x_1&x_2&\cdots&x_{p\phantom{+0}} \\x_2&x_3&\cdots&x_{p+1} \\ \vdots&\vdots&\ddots&\vdots \\ x_{N-p-1}&x_{N-p}&\cdots&x_{N-2}\end{bmatrix} \begin{bmatrix}c_0\\c_1\\ \vdots\\c_{p-1} \end{bmatrix} = \begin{bmatrix}-x_{p\phantom{+0}}\\-x_{p+1}\\-x_{p+2}\\ \vdots\\-x_{N-1} \end{bmatrix}$

Заполнять по правилам
$\begin{array}{ll|ll}\operatorname{for}& i\in 0\ldots N-p-1&\operatorname{for}& i\in 0\ldots N-p-1\\ \operatorname{for}& j\in 0\ldots p-1&&b_i:=-x_{i+p}\\ &A_{i,j}:=x_{i+j}\end{array}$

Значит, при Вашем выборе $N=100$ и $p=3$ будет $N-p=97$ уравнений и только $3$ неизвестных. То есть система будет сильно переопределённой. Если бы вычисления, с помощью которых получены $A$ и $b$, были абсолютно точными, система была бы совместной. Иначе, с ошибками округления, она будет лишь "почти совместной". Попробуйте сначала решить систему как есть. Возможно, какая-то стандартная функция в Mathcad решает (в некотором смысле) и переопределённые системы. Если нет, я скажу, что делать. Пока меня интересует только решение этой системы, то есть три числа $c_j$.

[Schrodinger's cat]

svv спасибо, начало продвигаться.
СЛАУ я составлял неверно.
Вектор $c = 0.671$, $-0.314$, $-1.285$
Добавляю 1 в начало вектора.
При этих значениях задержки и частоты (которые можно выразить из этих корней) получаются неверными.
Но если методом проб и ошибок поразворачивать вектор ответов и 1 переставлять, то получаются задержки $0.398+ i 3.142$, $- i 0.209$ и $+ i 0.209$
Первое значение непонятно что, но зато второе и третье похоже на заданные.
Та же история с частотами $0$, $-0.209$, $0.209$.
Понятно что это неправильный алгоритм но хотя бы что то стало получаться.

[ilghiz]

Да ктож влоб несовместную систему-то решает???

Вот пусть у вас получилась система, как писал svv $Ac=x$, так и решайте ее классикой, $A = QR$ - то есть классический QR алгоритм (вы же Грамм-Шмита недавно делали? - его для таких маленьких размерностей должно хватить, а будет больше, будете ранг-ревеалинг делать), а деальше $c = R^{-1} Q^T x$. Если совсем лениво, то можно вот так:
$c = (A^T A)^{-1} A^T x$, которое решается опять же через
1. $b=A^T x$,
2. $H=A^T A$,
3. Холецкий для $H$ и решение системы с этой факторизованной матрицей.

[B@R5uk]

ilghiz в сообщении #1565975 писал(а):

Если совсем лениво, то можно вот так:

$$c = (A^T A)^{-1} A^T x$$

Нужно, чтобы столбцы матрицы A были линейно независимы (иначе обратной матрицы не будет существовать). Тогда решение будет удовлетворять условию наименьших квадратов:

$$\lVert x-Ac\rVert^2=\min$$

Почему вы зовёте такой подход ленивым? Разве в большинстве случаев он не единственный?

Алсо, в Матлабе, я смотрю, нет необходимости всё это дело вручную прописывать: mldivide вполне справляется в лоб:

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

c = A \ x;

Надо взять на вооружение и не плодить сущности без необходимости.

ilghiz в сообщении #1565975 писал(а):

Холецкий для $H$

Спасибо за наводку, не знал! Это разложение ещё позволяет посмотреть на величину диагональных элементов и сообразить, когда и как матрица H близка к сингулярной. Плюс в Матлабе есть нативная реализация. Буду теперь пользоваться по мере необходимости. Вот только, интересно, будет ли разница в быстродействии с нативным Матлабовским делением слева?

Алсо, QR-алгоритм применим же только к квадратной обратимой матрица A, разве нет?

[ilghiz]

B@R5uk в сообщении #1565980 писал(а):

ilghiz в сообщении #1565975 писал(а):

Если совсем лениво, то можно вот так:

$$c = (A^T A)^{-1} A^T x$$

Нужно, чтобы столбцы матрицы A были линейно независимы (иначе обратной матрицы не будет существовать). Тогда решение будет удовлетворять условию наименьших квадратов:

$$\lVert x-Ac\rVert^2=\min$$

Почему вы зовёте такой подход ленивым? Разве в большинстве случаев он не единственный?

Не, не единственный, ленивый, и очень не популярный в вычислительной математике. Смотрите, пусть $A$ имеет обусловленность близкую к обратному значению машинной арифметики. Теоретически - матрица $A$ все еще не вырожденная, решение может быть, но вот если построить $(A^T A)^{-1}$, то ее обусловленность будет уще больше обратного машинного эпсилона, и Холецкий, даже с полным пивотированием гарантированно упадет. Да в любом случае, очень не правильно плохо обусловленную задачу приводить к такой, у которой обусловленность в квадрат выше, поэтому я назвал этот метод ленивым.

Правильно делать через $QR$. Конечно можно через SVD, но метод через SVD численно требует больше арифметических операций, поэтому если делать правильно, лучше через QR.

B@R5uk в сообщении #1565980 писал(а):

Алсо, в Матлабе, я смотрю, нет необходимости всё это дело вручную прописывать: mldivide вполне справляется в лоб:

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

c = A \ x;

Надо взять на вооружение и не плодить сущности без необходимости.

А вот ХЗ как матлаб это реализовал, поэтому я бы не советовал бы в случаях с плохой обусловленностью исходной матрицы :)

B@R5uk в сообщении #1565980 писал(а):

Алсо, QR-алгоритм применим же только к квадратной обратимой матрица A, разве нет?

Есть два QR - первый для построения ортогональной $Q$ и верхней треугольной $R$, второй - как один из шагов для итерационного поиска собственных значений.

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

ilghiz в сообщении #1565975 писал(а):

Холецкий для $H$

Спасибо за наводку, не знал! Это разложение ещё позволяет посмотреть на величину диагональных элементов и сообразить, когда и как матрица H близка к сингулярной.
[/quote]

пока вы делаете Холецкого без пивотинга диагонального элемента, вы почти всегда будете видеть в качестве обусловленности погоду на Марсе, то есть оценка по диагонали Холецкого обычно сильно завышена.

Простой пример

$$\left( \begin{tabular}{cc} 0 & 1 \\ 1 & 1 \end{tabular} \right)$$

И даже полное пивотирование тоже часто дает завышенную оценку, смотрим на матрицу с единичной обусловленностью:

$$\left( \begin{tabular}{cc} 0 & 1 \\ 1 & 0 \end{tabular} \right)$$

[B@R5uk]

ilghiz в сообщении #1565986 писал(а):

почти всегда будете видеть в качестве обусловленности погоду на Марсе

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

ilghiz в сообщении #1565986 писал(а):

Правильно делать через $QR$. Конечно можно через SVD

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

-- 02.10.2022, 16:22 --

ilghiz в сообщении #1565986 писал(а):

А вот ХЗ как матлаб это реализовал, поэтому я бы не советовал бы в случаях с плохой обусловленностью исходной матрицы :)

У меня в этом плане выработалось совершенно противоположное практическое правило: не знаешь, как лучше, доверяй разрабам Матлаб. Что, впрочем, не противоречит правилу проверки получившегося результата.

[ilghiz]

B@R5uk в сообщении #1565990 писал(а):

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

пока обусловленность матрицы $A$ маленькая и есть доля везения в Холецком, то да, один и тот же. Хотя я такие задачи обычно советую делать хотя бы через $QR$, по крайней мере есть хоть какая-то гарантия того, что оно не заткнется неожиданно из-за переполнения и сразу видна линейная зависимость столбцов в исходной матрице. Правильнее делать через rank-revealing-QR, так как он в Лапаке по умолчанию, я именно и его и вызываю.

B@R5uk в сообщении #1565990 писал(а):

У меня в этом плане выработалось совершенно противоположное практическое правило: не знаешь, как лучше, доверяй разрабам Матлаб. Что, впрочем, не противоречит правилу проверки получившегося результата.

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

[Schrodinger's cat]

ilghiz в сообщении #1565975 писал(а):

систему-то решает???

Хорошее замечание. Я использую mathCAD для изучения и предварительного написания алгоритмов а потом перекладываю обычно в C#.
Посмотрел справку, функция решающая СЛАУ lsolve использует LU разложение! Именно ее я использовал в прошлый раз.
QR не может вычислиться!
Сделал через SVD, через псевдообратную матрицу. Проверка $A^+A$ показывает почти единичную матрицу.
Соответственно и дальнейший расчет дал несколько иной вектор $c = 0.545$, $-0.055$, $-1.414$
И тут уже как ни комбинируй корни полинома никак не похожи на то что нужно.

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

Если нет, я скажу, что делать.

Давайте попробуем и Ваш вариант. Тем более что в C# у меня не будет волшебных функций mathcad'a которые бы решали все подряд.

ilghiz в сообщении #1565986 писал(а):

А вот ХЗ как матлаб это реализовал

Пока смотрел всякие лекции видел граф этого алгоритма, который рекомендовали смотреть в справке. Граф достаточно разветвленный, куча проверок, куча всяких решателей...
в отличие от matcad'a который просто LU делает

[B@R5uk]

Schrodinger's cat в сообщении #1565999 писал(а):

видел граф этого алгоритма

А можете ссылочку кинуть? Мне как постоянному пользователю весьма интересно, как там у него внутри всё устроено.

[ilghiz]

Schrodinger's cat в сообщении #1565999 писал(а):

QR не может вычислиться!

Выж на днях Грамм-Шмидта обсуждали, так это тот же самый QR, где $Q$ -ортогональная, $R$ - верхняя треугольная. Сделайте, будет вам счастье! LU там заведомо упадет, вернее вероятность, что оно даст адекватное решение реально стремится к нулю. Можно конечно делать LU с пивотингом и останавливаться как в адаптивной кросс-аппроксимации, (ACA - гуглятся на Марио Бебендорфа и Сергея Рязанова, если точнее, то на мозаично-скелетонные аппроксимации Юрия Еремина с Лилией Колотилиной еще в конце 80-х) но, вот внутренний голос мне подсказывает, что их в Маткаде не знают, хотя жаль, хороший метод придумали в конце 90-х.

[Schrodinger's cat]

B@R5uk в сообщении #1566000 писал(а):

А можете ссылочку кинуть?

Видел конкретно вот тут, начиная с 18:50
https://www.youtube.com/watch?v=OiNurcu6OBU

ilghiz в сообщении #1565975 писал(а):

вы же Грамм-Шмита недавно делали?

Да, методами Грамм-Шмидта, модифицированного Грамм-Шмидта, Гивенса и Хаусхолдера.
Попробовал мою функцию модифицированного Грамм-Шмидта, получился вот такой результат $c=0.573$, $-0.121$, $-1.383$
Так же методом проб и ошибок тасую полученные корни и получаю задержки $0.557+ i 3.142$, $+ i 0.209$ и $- i 0.209$
и частоты $0$, $0.209$, $-0.209$

А если покажимать кнопочку принудительного пересчета всего документа, то можно еще получить такие варианты
$c=0.252$, $0.508$, $-1.705$
$c=-0.203$, $1.397$, $-2.159$
$c=1.186$, $-1.321$, $-0.77$
ну я все комбинации не выяснял, видимо как направлен то ортогональный базис мы не задаем... ровно как и в других разложениях например знаки векторов SVD...
При этом значения задержек и частот одни и те же только знаки меняются и положение в векторах.

Гивенс и Хаусхолдер же работают только с квадратными матрицами?
Грамм-Шмидт да, прямоугольные кушает.

ilghiz в сообщении #1566002 писал(а):

Выж на днях Грамм-Шмидта обсуждали

Так это QR mathcad'a не может вычислить.
lsolve на базе LU того же mathcad'a как видите что то вычисляет! Причем результат похож на моего Грамм-Шмидта.

[Schrodinger's cat]

Поэкспериментировав выяснил, что именно в этой задаче, моя реализация метода Грамм-Шмидта выдает почему то разные последний столбец матрицы Q и самое правое нижнее значение матрицы R. Не то чтобы они каждый раз кардинально отличаются, есть некоторый небольшой разброс.
Хотя на тестовых данных алгоритмы отрабатывают стабильно, ничего не меняется. Но на тестовых матрицы мелкие 4х3.
Причем разброс не модифицированного гораздо больше чем модифицированного.
Видимо это и есть численная нестабильность данного метода.
Но, что интересно все остальные значения матриц не меняются.
Чудеса математики