Тензорный способ нахождения грамоник, прогноз

[Schrodinger's cat]

ilghiz спасибо.
Дочитал и Ваш патент и остальные представленные материалы (Kruskal, Ваши статьи).
Качество понимания прочитанного конечно хромает, я не владею английским свободно и текст сильно отягощен технической спецификой к данной теме по сути не относящейся. Все эти предварительные фильтры, катушки... Уверен математические формулы для моего простого случая могут быть проще. Но в целом понял о чем речь. Было интересно.

Я так понял что представленный там решатель выдает только тензорную декомпозицию. Ну и плюс вспомогательная функция осуществляет наименьшие квадраты в итеративном алгоритме решателя.
Это хорошо, что он есть, но в начале я бы хотел обойтись стандартными средствами. Не углубляясь в математику самой декомпозиции. Алгоритмы декомпозиций рассмотрю отдельно если будет желание. Если конечно Ваша декомпозиция не какая нибудь специальная. Тогда придется рассмотреть.
Вы сказали, что можно использовать любое тензорное разложение (правда ТТ Вы уже забраковали выше )
Из доступных есть Canonical Polyadic, CANDECOMP-PARAFAC, PARAFAC, Tucker.
Читал, что есть еще "скелетная декомпозиция", возможно она представлена в перечисленном но я пока не разобрался. С тензорами работаю впервые. Попробовал парочку, что то разлагают.
Был бы очень признателен если подскажите, что лучше использовать.
Прочитал что декомпозиция Такера более устойчива.

Создал алгоритм преобразующий линейный индекс в многомерный, одномерный массив в тензор. Программа на питоне.
Отладочный вывод программы с основанием 2:

signal = [ 0.  1.  2.  3.  4.  5.  6.  7.  8.  9. 10. 11. 12. 13. 14. 15.]
tensor =
 [[[[ 0.  8.]
   [ 4. 12.]]

  [[ 2. 10.]
   [ 6. 14.]]]

 [[[ 1.  9.]
   [ 5. 13.]]

  [[ 3. 11.]
   [ 7. 15.]]]]

Отладочный вывод программы с основанием 3:

signal = [ 0.  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.]
tensor =
 [[[ 0.  9. 18.]
  [ 3. 12. 21.]
  [ 6. 15. 24.]]

 [[ 1. 10. 19.]
  [ 4. 13. 22.]
  [ 7. 16. 25.]]

 [[ 2. 11. 20.]
  [ 5. 14. 23.]
  [ 8. 17. 26.]]]

Верно ли заполняются тензоры?
Дальше я подставляю туда сигнал и пропускаю через тензорное разложение.

Недостающие точки заполняются нулями. Это не верно?

Подскажите пожалуйста что делать дальше? Я так понимаю метод дает на выходе сумму экспонент. Совершенно непонятна их связь с разложением.

ПС: ilghiz поздравляю Вас с завершением представленной научной работы. Меня всегда позитивно удивляли люди способные выдумать такие вещи.

[ilghiz]

Рад, что Вам, Schrodinger's cat, мои результаты могут быть полезны!

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

Отладочный вывод программы с основанием 2:
...
Отладочный вывод программы с основанием 3:

да, тут все верно.

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

Недостающие точки заполняются нулями. Это не верно?

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

В нашем решателе (NWSPF) там даннные вообще залиты как одномерный массив, но одновременно с этим массивом, там подается массив многомерных индексов где лежит какое данное. Вкрадце об этом формате у меня было даже еще в 10.1002/nla.297 написано, вы ее уже вроде читали, посмотрите там главу 4 - там и блоксемы алгоритмов и структуры данных описаны. Если все-таки это будет не сильно поняттно, спрашивайте что именно вызывает затруднения - с радостью расскажу.

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

Вы сказали, что можно использовать любое тензорное разложение (правда ТТ Вы уже забраковали выше )
Из доступных есть Canonical Polyadic, CANDECOMP-PARAFAC, PARAFAC, Tucker.

тут мы с вами друг друга не поняли.

TT, PARAFAC, Tucker - дают разные разложения, хотя работают с одними и теми же тензорами. Только PARAFAC и несколько его аналогов позволяют получить то разложение, что вам должно подойти. Только sparse PARAFAC способен работать с данными, у которых чего-то не достает (есть еще несколько других методов, но там совсем дичь, и лучше туда пока не лезть).

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

То есть я бы вам посоветовал начать с PARAFAC и строить полный тензор - так будет проще, по первости по крайней мере. Потом, если тут все будет получаться, перейти на sparse PARAFAC и им обрабатывать разреженные данные (тензоры с недостающими элементами). NWSPF - это как раз sparse PARAFAC. Если и это будет с успехом пройдено, можно дальше пробовать - я еще могу чего на эту тему накидать, но не сейчас, иначе точно каша будет у вас в голове.

[Schrodinger's cat]

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

То есть я бы вам посоветовал начать с PARAFAC и строить полный тензор - так будет проще, по первости по крайней мере.

Полностью согласен.

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

Вкрадце об этом формате у меня было даже еще в 10.1002/nla.297 написано, вы ее уже вроде читали, посмотрите там главу 4 - там и блоксемы алгоритмов и структуры данных описаны.

Еще раз перечитал.
Да, там про само разложение. В целом даже понял как этот алгоритм работает.
Решение СЛАУ в нем это и есть ALS как я понял.
Но как составлять матрицы и некоторые детали непонятны. Это пока опустим.
Пока использую готовый PARAFAC.

Разобрался как его запускать. Непонятно какой ранг выставлять, допустим 1. 2 почему то виснет (не может выйти из локального минимума?), 3 работает.
Пока для наглядности за исходные данные беру все тот же увеличивающийся ряд. (Но я уже создал функции создающие гармоники, дальше переключусь на них)
Исходный тензор

signal = [ 0.  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.]
tensor =
 [[[ 0.  9. 18.]
  [ 3. 12. 21.]
  [ 6. 15. 24.]]

 [[ 1. 10. 19.]
  [ 4. 13. 22.]
  [ 7. 16. 25.]]

 [[ 2. 11. 20.]
  [ 5. 14. 23.]
  [ 8. 17. 26.]]]

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

weight = [1.]
factors = 
[array([[-42.65355615], [-45.25686168], [-47.86016722]]), 
array([[0.47106575], [0.57149718], [0.67192862]]), 
array([[-0.16593588], [-0.50613755], [-0.84633922]])]

Реконструкция по разложению

Reconstruction=
 [[[ 3.33408805 10.16963414 17.00518024]
  [ 4.04491715 12.33780482 20.63069248]
  [ 4.75574626 14.50597549 24.25620472]]

 [[ 3.53757987 10.79032483 18.04306978]
  [ 4.29179353 13.09082704 21.88986055]
  [ 5.04600718 15.39132925 25.73665132]]

 [[ 3.74107169 11.41101551 19.08095932]
  [ 4.5386699  13.84384926 23.14902862]
  [ 5.3362681  16.27668301 27.21709791]]]

Остаток (ошибка) разложения

Error =
 [[[ 3.33408805  1.16963414 -0.99481976]
  [ 1.04491715  0.33780482 -0.36930752]
  [-1.24425374 -0.49402451  0.25620472]]

 [[ 2.53757987  0.79032483 -0.95693022]
  [ 0.29179353  0.09082704 -0.11013945]
  [-1.95399282 -0.60867075  0.73665132]]

 [[ 1.74107169  0.41101551 -0.91904068]
  [-0.4613301  -0.15615074  0.14902862]
  [-2.6637319  -0.72331699  1.21709791]]]

Что делать дальше?
Нигде в указанных материалах не увидел как находить коэффициенты экспонент и как точно выглядит модель.
Это упоминается в патенте но там очень мутно написано.

[ilghiz]

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

Решение СЛАУ в нем это и есть ALS как я понял.

Не совсем так. СЛАУ надо решать на каждой итерации ALS, но не хочется это делать с огромной матрицей и, часто плохо обусловленной матрицей, и там как раз я и предлагал эффективные методы решения, чтобы не наступать на очевидные грабли. Итераций ALS, кстати, может быть сказочно много, у метода только монотонная сходимость, что применяли против этого авторы пакета, который вы используете - мне не известно, если у них не было достаточно опыта, могли вообще ничего не применять.

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

Исходный тензор

signal = [ 0.  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.]

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

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

С рангом 3 решение должно быть примерно таким:

$$9\sqrt{15}; (\frac1{\sqrt{3}},\frac1{\sqrt{3}},\frac1{\sqrt{3}}); (\frac1{\sqrt{3}},\frac1{\sqrt{3}},\frac1{\sqrt{3}}); (0,\frac1{\sqrt{5}},\frac2{\sqrt{5}});$$

$$3\sqrt{15}; (\frac1{\sqrt{3}},\frac1{\sqrt{3}},\frac1{\sqrt{3}}); (0,\frac1{\sqrt{5}},\frac2{\sqrt{5}}); (\frac1{\sqrt{3}},\frac1{\sqrt{3}},\frac1{\sqrt{3}});$$

$$\sqrt{15}; (0,\frac1{\sqrt{5}},\frac2{\sqrt{5}}); (\frac1{\sqrt{3}},\frac1{\sqrt{3}},\frac1{\sqrt{3}}); (\frac1{\sqrt{3}},\frac1{\sqrt{3}},\frac1{\sqrt{3}}); $$

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

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

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

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

Что делать дальше?

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

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

Нигде в указанных материалах не увидел как находить коэффициенты экспонент и как точно выглядит модель.

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

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

Это упоминается в патенте но там очень мутно написано.

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

[Schrodinger's cat]

ilghiz Ваше программное обеспечение по адресу ввв Ильгиз точка ком более недоступно. Хотел использовать его для проверки правильности результатов пока учусь.
Пока пробую этот пакет. Насколько я понял он самый популярный.
https://tensorly.org/stable/modules/generated/tensorly.decomposition.parafac.html
Но в нем похоже нет регуляризаций
https://github.com/tensorly/tensorly/issues/169?ysclid=la0p9nhlss826638915
По ссылке что то обсуждается, но я так понял это все не то.

[ilghiz]

Да, верно, давно у меня уже нет домейна ilghiz.com - я его где-то лет так десять назад не продлил и он где-то затерялся. И даже софта с тех времен в том первозданном виде не осталось, реально. Есть конечно новый софт, но он на годовых подписках для узкого круга "своих" заказчиков, и очевидно, этот софт не должен попадать в открытый доступ или чужие руки.

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

Ума не приложу, почему такая простая вещь до сих пор не появилась в таком распространенном пакете, хотя мне конечно это на руку :)

Причина, почему эта регуляризация сильно помогает, кроется в том факте, что в 3+ мерных разложениях имеются решения с бесконечно большими по норме тензорами, например, как у вас в вашем примере выше с рангом равным двум. Для большого класса таких решений эта симметричная матрица становится сингулярной, то есть хотя бы одно ее собственное значение стремится к нулю. Очевидно, в плавающей арифметике это приводит к проблемам с устойчивостью решения, а вот тихоновский регуляризатор просто сдвигает на "чуть-чуть" все собственные значения, избегая численной неустойчивости и практически не коверкая все остальные собственные значения, поэтому с такой регуляризацией сходимость ALS бывает существенно лучшей, чем без него.

[Schrodinger's cat]

Покопавшись в коде я обнаружил, что недокументированная регуляризация все же добавлена но в странном виде. Одним из параметров идет число которое можно задать перед вычислением, далее в коде как Вы и говорите единичный тензор умножается на это число и далее эта конструкция прибавляется к псевдо инверсии на каждой итерации. Очевидно данное значение регуляризатора никак не модифицируется в ходе вычислений. А Ваш Тихоновский должен вычисляться каждый цикл.
Там же есть еще один параметр, ограничитель максимального числа итераций.
Имея эти 2 рычага можно исправить данный недостаток ограничив функцию решателя до 1 итерации и вычислять регуляризатор Тихонова.

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

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

Опа, а вот это я упустил. Ну ничего прочитал сейчас.

Но я до сих пор не могу понять что получается на выходе разложения. Помимо векторов решения есть еще какие то веса, что это за веса? У меня они почему то всегда равны 1.
Давайте рассмотрим пример. Например у нас будут 2 гармоники:

Код генерации и вывод программы

(Оффтоп)

signal1 = np.empty((Length,))
signal2 = np.empty((Length,))
    for i in range(Length):
        signal1[i] = 0.0 + 10.0 * math.sin(i/100 * 2.0*math.pi*1.3  + 0.0)
        signal2[i] = 0.0 + 10.0 * math.sin(i/100 * 2.0*math.pi*10.0 + 0.0)

signal = signal1 + signal2

Дискретные значения:
0.           6.69375864  11.13693681  11.93655747   9.08728862
3.97147891  -1.1708132   -4.09935264  -3.43126219   0.82900324
7.28968627  13.70176063  17.81652415  18.24318971  14.97891223
9.40880769   3.77596387   0.32388689   0.43894501   4.120371
9.98026728  15.77361371  19.25583389  19.04035858  15.12862459
8.91006524   2.6320923   -1.45748631  -1.96805136   1.10380167
6.3742399   11.60217378  14.54679718  13.82512562   9.44197131
2.78991106  -3.88075272  -8.31959356  -9.13366334  -6.31753371
-1.25333234   3.81922644   6.66037254   5.8878115    1.50669486
-5.09041416 -11.65357956 -15.93309169 -16.53706486 -13.46147168
-8.09016994  -2.66492179   0.57215092   0.23611363  -3.67079292
-9.75916762 -15.78246678 -19.49458066 -19.50740706 -15.82086043
-9.82287251  -3.75938426   0.12322659   0.43572092  -2.82398502
-8.27080574 -13.66247554 -16.75709646 -16.17068384 -11.90714794
-5.35826795   1.22634174   5.59682849   6.36069997   3.51286255
-1.56434465  -6.63112058  -9.44773372  -8.6320532   -4.18951808
2.48689887

Выбираю 81 значение, основание системы счисления 3, ранк 5

fac = (weights, factors) : rank-5 CPTensor of shape (3, 3, 3, 3) 
weights = [1. 1. 1. 1. 1.]
factors = [array([
[-1.52870949,  1.33616223,  0.        ,  0.0862158 , -1.54259496], 
[-1.87217963, -0.178804  ,  0.        , -1.35967955, -1.58081122],
[-1.89511359, -2.00860176,  0.        , -2.55644413, -1.43904873]]), 
array([
[ 0.19385824, -0.91555159,  0.        ,  1.76845806,  1.13059629],
[-2.438654  , -0.70153956,  0.        , -0.28094468,  2.1881782 ],
[-1.8571081 ,  2.12635579,  0.        , -2.27710218,  0.93970888]]), 
array([
[ 0.92023131, -1.83597749,  0.        ,  1.68339403, -1.60299346],
[ 1.98287081, -1.48736182,  0.        ,  1.91761919,  0.28275873],
[ 2.15751304, -0.51838487,  0.        ,  1.37130937,  2.07361618]]), 
array([
[ 2.21854656,  1.15246083,  0.        , -2.04440791,  0.07835575],
[-0.16318326, -2.12633716,  0.        ,  1.57984291,  1.7510471 ],
[-2.11775522, -0.04801803,  0.        ,  1.30998508, -1.96904012]])]

Что есть что?

Не понимаю, что из этого есть альфа (из чего Тихоновский регуляризатор считается)?

Там написано, что в разложении каждый вектор представляет точки экспоненты.
Значит, если я возьму первый вектор:

[-1.52870949,  1.33616223,  0.        ,  0.0862158 , -1.54259496]

Это точки одной экспоненты?

Нужно ли до разложения вычитать константную составляющую из данных?

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

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

Как это логарифмировать если есть и нули и отрицательные значения?

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

[ilghiz]

чушь какая-то в выдаче, похоже что-то не так запускается - я этим пакетом не пользовался, только когда-то лет 15 назад, когда этого пакета еще не было, с некоторыми участниками сравнивался. То есть скорей всего вы что-то туда не то подаете. Веса тоже равными не должны быть!!!

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

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

[Schrodinger's cat]

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

Посмотрю и подскажу обязательно

Большое спасибо. Попробовал еще несколько вариантов... но так как у меня опыта 0, не понимаю что правильно а что нет.

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

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

Не проблема. Пока с Прони играюсь.

[Schrodinger's cat]

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

[Schrodinger's cat]

Здравствуйте, с прошедшим новым годом Вас.
Переписал алгоритм из Вашего патента, но у меня 100500 ошибок, без Вас не вижу возможности довести это рабочего состояния.
- откуда взять совместимую с кодом библиотеку комплексных чисел (ни одной не знаю, очень давно последний раз писал на с++)?
- что такое SET, ABS2, DOTU, POSV?
- в других функциях почему то не совпадают по типу enum'ы.
- в некоторых местах понятно какую версию функций (действительную или комплексную) использовать в некоторых нет. Причем там названия длиннее чем указаны в листинге алгоритма.
- в объявлении функций в переменных заданы размерности массивов, синтаксически это ошибка, я так понимаю это просто для понимания что туда пихать?
Как я писал в смежной ветке, удалось подключить блас мкл. Пробую реализовать на нем.

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