Сингулярное разложение для не полностью заданных данных

[ilghiz]

Добрый день и с наступающим Новым Годом и наступившими и наступающими праздниками!

У меня вопрос в общем виде, как принято решать такие задачи сейчас, а именно, пусть мы ищем набор функций $f_j(x), g_j(y) \in {\mathBB R}, j=1,\dots,r$ пусть даже $r$ - задано, и у нас имеются данные $a_j(x,y) \in {\mathBB R}, U^{(x)}_i, U^{(y)}_i, L^{(x)}_i, L^{(y)}_i, i=1,\dots, N$, что

$$\min_{f_1(x), \dots, f_r(x);~ g_1(y), \dots, g_r(y)} \sum_{i=1}^N \int_{L^{(x)}_i}^{U^{(x)}_i} \int_{L^{(y)}_i}^{U^{(y)}_i} \left| a_i(x,y) - \sum_{j=1}^r f_j(x) g_j(y) \right|^2$$

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

Также обобщение на 3-мерный и любой произвольно больше 2 мерный случай функций $a(x_1,\dots,x_M)$ тоже конечно же интересен.

Если все задано без дырок в двухмерном случае, то вычислительная сложность такой аппроксимации стремиться к $\max(K_1,K_2)r^2$, где $K_1,K_2$ - размерности дискретизации $x$ и $y$.

В общем случае мне пока понятен только один стародавний метод, в котором мы исходные данные вначале дополняем до полных (нулями на первой итерации), решаем обычную малоранговую аппроксимацию, далее начинаем дополнять исходные данные на основе уже вычисленных на текущей итерации $f_i(x), g_i(y)$.

Очевидно, что у такого метода будет существенно больше вычислительная сложность и только монотонная сходимость.

Скажите, пожалуйста, а есть ли какие-то современные методы, про которые бы была бы известна более хорошая и сходимость и меньшая вычислительная сложность?

Спасибо!

[aitap]

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

Также обобщение на 3-мерный и любой произвольно больше 2 мерный случай функций $a(x_1,\dots,x_M)$ тоже конечно же интересен.

В двумерном случае $ A \approx FG^\top = F R R^{-1} G^\top = (F R) (R^{-1} G^\top) = F' {G'}^\top$, т.е. решение задачи можно взять и "повернуть" в другое решение (с такой же невязкой) путём умножения на любую обратимую матрицу $R$ и обратную к ней. В трёхмерном случае выражение выглядит чуть сложнее, а в качестве матрицы $R$ можно подставить только матрицу перестановки и умножения компонентов низкорангового разложения на скаляры; "поворот" в другое решение невозможен (кроме вырожденных случаев).

Такую задачу ещё называют higher-order SVD. В зависимости от того, какого вида данные, можно применять модель PARAFAC (где с "дырками" в тензоре $\mathbf A$ работают при помощи MILES и других подходов) или более общую (но не дающую единственного решения) модель Tucker (где, плюс к матрицам компонентов, появляется тензор $G$, аналогичный $\Sigma$ в сингулярном разложении, содержащий ненулевые элементы не только на главной диагонали).

[ilghiz]

Спасибо большое, aitap за ответ!

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

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

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

поэтому я продолжаю искать методы с надежной (не монотонной, а лучше) сходимостью и существенно меньшей арифметической сложностью.

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

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

Как вариант - вычисляем сингулярное разложение $X=SVD$ через $R=X^TX=D^TV^2D$, а для вычисления элементов R используем только неотсутствующие значения, пропуская пары, где хоть одно значение отсутствуе (деля сумму на фактическое число слагаемых)

[ilghiz]

Спасибо большое, Евгений Машеров, за ответ!

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

Если бы данные генерились по запросу, но их всех создать было бы нельзя (например, много места занимают), то так было бы конечно правильно. Я как раз в той самой статье, что указал выше, аппроксимировал данные, которые должны были занимать тензор 40000х40000х40000 (то есть 512ТБайт), а результат аппроксимировал одномерными объектами размера всего-то в десятки мегабайт.

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

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

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

Ну, книжку Литтл Дж., Рубин Л. Статистический анализ данных с пропусками Вы, наверно, знаете. Но на всякий случай напомню.

[ilghiz]

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

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

В любом случае этих данных - нет. Соответственно варианты:
1. Пренебречь наблюдениями, которые неполны.
1а. Полностью отбрасывать такие наблюдения.
- большой процент потери данных
1б. Использовать те элементы, которые есть в неполном наблюдении (то есть из вектора наблюдённых значений брать имеющиеся для расчёта средних, дисперсий, для корреляций - пары, где оба параметра известны)
- могут возникать несоответствия, скажем, рассчитанная так корреляционная матрица может утратить неотрицательную определённость.
2. Заменить неизвестные наблюдения "правдоподобными"
2а. Условными значениями (нулями, средними, медианами etc.)
- наибольшая недостоверность результата. Можно повторить расчёт, подставив вместо нулей и прочего полученные оценки этих значений. К чему-то сойдётся, но будет ли означать сходимость осмысленность?
2б. Регрессионными оценками, исходя из известных элементов вектора наблюдений.
- но при этом из информации, которая должна была быть в потерянном наблюдении, теряется самое интересное - отклонение от ожидаемого.
2в. Рассматривать наблюдения, как временные ряды, и заменять прогнозом.
- то же, что и в предыдущем, только пилим матрицу данных вдоль, а не поперёк.
3. Рассматривать отсутствующие наблюдения, как параметры модели, подлежащие оценки наряду с прочими.
- размерность пространства параметров возрастает драматически, вместо со сложностью обработки и ошибками оценки параметров. При этом модель должна отражать действительные связи объекта, не просто "подгонять данные".

[ilghiz]

Спасибо большое, Евгений Машеров за советы!

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

Если вконце - то можно забить нулями, так как там и так шум большой. Если в начале - забить сложно, там самые важные данные.

По физике аппаратуры - сдвиг начала всегда скачет вперед-назад на значение порядка нескольких точек оцифровки, но не на целое значение шагов, а на дробное.

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

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

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

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