Наименьшие квадраты в общем виде для задач, Комп-Томографии

[ilghiz]

Добрый день,

Пусть на некой трехмерной области у нас есть неизвестные искомые функции $g(x) \ge 0, f(x) \ge 0, \in {\mathbb R}$, $x \in {\mathbb R}^3$.

Для их поиска мы имеем множество линейных преобразований в трехмерном пространстве $Z = \{z_i\}_{i=1,\dots,N}$, что $z(x) = y$, где $x, y \in {\mathbb R}^3$, причем эти преобразования содержат только отртогональные матрицы трансформации. То есть $z_i$ описывают только параллельный перенос и вращение, без изменения скалировки.

Для каждой такой матрицы у нас дано соответсвующее значение $h_i$, что

$$\int_{\mathbb R^3} f(x) g(z_i(x)) dx^3 = h_i.$$

Если бы по задаче $g$ было бы задано, то поиск $f(x)$ можно было бы осуществить методом наименьших квадратов

$$\min_{f} \sum_i \left|\int_{\mathbb R^3} f(x) g(z_i(x)) dx^3 - h_i\right|_2^2$$

приводяшей к линейной задаче по дискретным значениям по $f(x)$ - ведь каждое значение $g(z_i(x))$ можно вычислить и дискретизовать. Такая задача очень распространена в компьютерной томографии.

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

Более того, я понимаю, что в общем случае, при произвольных $g(x)$ и $f(x)$ единственного решения может и не быть. В моем случае, я знаю, что $g(x)$ - имеет ненулевые значения на какой-то мне заранее неизвестной поверхности, но, к сожалению, эта поверхность довольно далека от плоскости.

С численной точки зрения - мне надо бы дискретизовать $f(x)$ с точностью около тысячи точек по каждому направлению, то есть наименьшие квадраты в лоб не сойдуться даже если их хорошо регуляризовать даже на хорошем суперкомпьютере. Также мне достаточно искать только $f(x)$, а вид $g(x)$ меня особо не волнует.

Не красивым вариантом решения я вижу попеременную минимизацию

$$\min_{f,g} \sum_i \left|\int_{\mathbb R^3} f(x) g(z_i(x)) dx^3 - h_i\right|_2^2$$

так, что на каждом четном шаге мы фиксируем $f$ и минимизируем $g$, а на каждом нечетном шаге мы фиксируем $g$ и минимизируем $f$, но сходимость будет монотонная, доказать и гарантировать ничего нельзя будет и, скорей всего, решение будет реально кривым.

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

Спасибо!

[B@R5uk]

$g(x)$, я так понимаю, у вас приборная функция? Почему бы не взять какую-нибудь хорошую $f(x)$ и не откалибровать по ней ваш прибор? Это не обязательно делать до измерений, можно и после, главное, чтобы $g(x)$ со временем не плыло. Или вопрос стоит только об обработке чужих данных и возможности что-то изменять нет?

[ilghiz]

Спасибо большое, B@R5uk за комментарий!

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

$g(x)$, я так понимаю, у вас приборная функция?

да, верно.

Плывет ее форма от эксперимента к эксперименту очень сильно. У меня есть точная форма этой функции в идеальных условиях: и промоделированная, и измеренная. Грубо говоря, сама функция $f(x)$ задана в объеме 30х20х20 см, нужна точность хотя бы 1 мм, а форма поверхности этой $g(x)$ плывет на сантиметр, а то и на два. Пока решаю покомпонентной минимизацией, как написал выше, и даже решается, но мне такое решение очень не нравится ибо нет гарантии, что всегда и везде будет хорошо сходиться. А в планах и не за горами бОльшая точность, и, соответсвенно, под миллиард конечных элементов, вместо 12 миллионов как сейчас.

[B@R5uk]

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

[ilghiz]

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

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

ilghiz, что-бы что-то конкретное посоветовать хорошо бы знать, как ваш измерительный прибор устроен.

С этим довольно сложно - постановка - как я написал, аппаратура своя, и из нее выжато все, что можно. По физике функция $g(x)$ может изменяться. Это просто разные эксперименты. Более того, для каждого эксперимента я теоретически эту функцию могу посчитать, но для этого нужно навесить на прибор кучу дополнительны сенсоров, и вогнать все в довольно мощный суперкомпьютер, а это нельзя сделать по причине того, что прибор не должен иметь эти дорогущие сенсоры и суперкомпьютер в кустах.

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

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

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

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

[B@R5uk]

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

Это просто разные эксперименты.

То есть, во время проведения эксперимента $g(x)$ фиксирована и не плывёт? Если да, то это немного улучшает ситуацию. Можно небольшую часть рабочего пространства выделить под эталонную $\widehat{f}(x)$, специально подобранную так, чтобы приборная функция по ней калибровалась во время эксперимента максимально быстро, максимально точно и максимально удобно в вычислительном плане.

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

Какое-то решение есть, попеременная минимизация, но расстраивает отсутствие доказанной сходимости.

Тут всё гораздо хуже. Даже когда вычислительная часть идеальна, даже когда приборная функция задана точно, задача всё равно является некорректно поставленной. На этот случай есть всякие теории о регуляризации, но я не знаком с этой темой. Как правило в таких ситуациях по физическим процессам в измерениях строят некоторую модель этих явлений, а потом по ней считают. Это будет более грамотно и регуляризация может оказаться не нужна (во всяком случае к этом стоит стремиться). Судя по вашему МНК у вас даже этого нет. Более того, у вас приборная функция тоже не известна, и это ещё больше ухудшает ситуацию. Даже если у вас очень богатая выборка $\{z_i(x),\;h_i\}$ результат может оказаться очень далёк от реальности. Попробуйте оценить погрешность ваших вычислений. Для этого в показания $h_i$ можно внести случайный шум с дисперсией, оцененной по величине невязки результатов МНК и посмотреть как в результате пляшет конечный ответ. Даст вам хотя бы небольшое представление о том, что у вас получается на самом деле. Оно будет смутное, потому что случайный шум случаен и в итоге сам себя гасит, а реальные погрешности как правило коррелируют; но если есть какие-то сингулярности в задаче, шум их подсветит.

Сама же вычислительная задача в том виде в каком вы её поставили очень хорошая. Я бы даже сказал она просто шикарная: минимум не только существует, он единственен (являясь глобальным), и подход к нему квадратичный, даже когда вам надо найти обе функции $f(x)$ и $g(x)$. Нет таких гадостей и кошмаров числовиков, как, например, функция Розенброка. Так что любой алгоритм минимизации типа градиентного спуска или симплекс-метода даст правильный (с вычислительной точки зрения) ответ. Более того, вы фактически минимизируете величину

$$M=\sum_i(F^TZ_iG-h_i)^2$$

где $F$ и $G$ — это искомые столбцы чисел, а $Z_i$ — некоторые матрицы, соответствующие вашим трансформациям, и ответ к задаче можно выписать в явном виде методами линейной алгебры. Для получения ответа надо будет обращать некоторую матрицу, и ваша проблема может оказаться как раз в ней: она будет сингулярной или, во всяком случае, с вычислительной точки зрения сингулярной (так как точность чисел с плавающей точкой ограничена — это другая головная боль числовиков). По этой причине крошечные отклонения в исходных данных приведут к серьёзным искажениям вычисляемых значений, которые окажутся порядка самих значений. Методы регуляризации призваны бороться с такими случаями.

По численной оптимизации познакомьтесь со статьёй, например. (Первая из Гугла, но основные вопросы затронули). И я настоятельно рекомендую оценивать погрешность ваших результатов!

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

Что мне сразу не нравится в этой задаче - что в ней заложена грандиозная неопределённость. Можно G и F домножить на некую ортогональную матрицу, и получим точно такое же по точности совершенно иное решение. Надо вводить какие-то ограничения на G. И то, возможно, ничего не получится.

[B@R5uk]

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

Можно G и F домножить на некую ортогональную матрицу, и получим точно такое же по точности совершенно иное решение.

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

Я понимаю, что вы имеете в виду. Такая проблема бы была, если бы $Z_i$ была только одна, или если этих измерений будет меньше, чем искомых компонент в векторах G и F. Если же измерений много больше, то проблемы нет. Разумеется, весь набор $\{Z_i\}$ должен быть грамотно подобран.

[ilghiz]

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

Позвольте с Вами B@R5uk не согласиться по поводу минимизации функционала четвертой степени методами BFGS и иже с ними.

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

$$\min_{X,Y,Z} \sum_{i,j,k} || a_{i,j,k} - \sum_r x_{ir} y_{jr} z_{kr} ||_2^2$$

Во-первых, начинаются пляски с нормировкой, иначе все квазиньютоны забывают, что они должны сходиться почти квадратично. Я для этого даже решал все так, что аналитически считал, например, $X$, и по вот этому сложному функционалу начинал вычислять BFGS, а аналитически градиенты Бауром-Штрассеном считал. И только при таких ограничениях что-то да начинало сходиться. Свои ссылки или статьи об этом могу в личку послать.

В моем случае, конечно можно выразить аналитически $f(x)$ и минимизировать по $g(x)$, благо я знаю, что $g(x)$ уползает от идеальной не сильно далеко и можно положить очень много значений в $g(x)$ тождественно равных нулю.

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

Еще один момент. Я знаю (исходя из внешних параметров прибора), когда у меня $g(x)$ не меняется. Но я не могу заставить юзера всегда работать с постоянной $g(x)$, и, более того, во время работы я могу по косвенным параметрам определить, что сейчас надо работать с новой $g(x)$.

По поводу устойчивости решения. Да, я уже попробовал - вернее я использую устойчивость решения как некий сигнал юзеру, что он снял еще мало данных и решение еще очень плохое в контектсте этих данных. На самом деле, чтобы повысить устойчивость можно потребовать малые флуктуации в $f(x)$, например, регуляризацией по Тихонову - конечно решение от этого становится менее реалистическим, но сходимость повышается, а величина коэффициента регуляризациии, когда решение не разваливается - как раз показывает юзеру, сколько ему еще надо насобирать данных.

[B@R5uk]

Ваша тензорная задача действительно не очень (хотя на первый взгляд она кажется выпуклой, а это плюс), но ваша исходная задача то совсем другая. (Кстати, зачем вы две двойки после модуля пишите? Одна — квадрат, а нижняя?) И с ней то как раз всё хорошо, ну, или, во всяком случае, значительно лучше, чем это обычно бывает. Градиентные методы должны давать ответ быстро. Нормировка — это вообще всегда хорошая штука, но в вашей задаче она не нужна, если я правильно смысл функций $f(x)$ и $g(x)$ понимаю (первая — это "плотность" объекта, вторая — чувствительность прибора к этой "плотности"). Не нравится считать производные, попробуйте симплекс-метод Нелдера-Мида, отличная штука! Он не использует производные, плюс он от многих проблем градиентных методов свободен от слова "совсем". Как раз будет очень хорошо работать с предложенным вами подходом:

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

В моем случае, конечно можно выразить аналитически $f(x)$ и минимизировать по $g(x)$, благо я знаю, что $g(x)$ уползает от идеальной не сильно далеко и можно положить очень много значений в $g(x)$ тождественно равных нулю.

Количество оптимизируемых переменных уменьшается в разы (на порядок? на два?), плюс удобно контролировать факт того, что матрица для нахождения $f(x)$ оказалась близкой к сингулярной, и надо данных дособирать.

[ilghiz]

Спасибо, B@R5uk за ответ!

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

Я экспериментировал с функциями, которые были дискретизованы на сетке с 12 миллионами неизвестных и у меня было примерно 100 миллионов экспериментальных точек.

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

Стоимость умножения на такую матрицу получалась около минуты, да, она хоть и разреженная, но в ней реально дочерта ненулевых элементов, так как она формируется из $AA^T$, где $A \in {\mathbb R}^{12,000,000 \times 100,000,000}$.

ALS сходились за 20 итераций, итого, где-то сутки считалось.

Тепереь если перейти на квазиньютонов.

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

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

Что дальше? Вот вы сами верите, что квазиньютон с 12 миллионами неизвестными сойдется хотя бы за тысячу итераций? Я - нет, и пробовать не хочу. А тысяча итераций - это тысяча часов счета, то есть больше месяца, против суток в ALS. А сколько симплекс будет это все жевать - вечность???

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

Но, я и сутки ждать не готов, мне надо за час (это когда $g$ - известна), а лучше быстрее, и не градиент, а все решение.

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

[B@R5uk]

С такими размерами данных (их стоило сразу упомянуть) задача выше моей квалификации. Чем мог — помог.

[ilghiz]

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

С такими размерами данных (их стоило сразу упомянуть) задача выше моей квалификации. Чем мог — помог.

Спасибо большое, за советы!!! А про 12 миллионов неизвестных - я сразу написал :)

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

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

[Geen]

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

Более того, я понимаю, что в общем случае, при произвольных $g(x)$ и $f(x)$ единственного решения может и не быть.

Рассмотрим одномерный вариант задачи. Более того, пусть ф-ции $f$ и $g$ заданы в конечном числе точек, например, в двух. Тогда, по условию, нам известны $f_1g_2$, $f_1g_1+f_2g_2$ и $f_2g_1$ - то есть имеется три уравнения на 4 неизвестных (и там ещё и квадратное уравнение в процессе возникает)....