В PARAFAC гессиана необратима?

[aitap]

Во многих областях знания выполняется предположение, что воздействие разных факторов на наблюдаемый сигнал разделяется на функции от независимых переменных:

$$ S(a,b,c,d,...) = \sum_r A_r(a) B_r(b) C_r(c) D_r(d) ... $$

Тогда можно измерить сигнал при всех комбинациях всех переменных, получить многомерную ($ \ge 3 $) матрицу наблюдений и оценить профили по каждому из измерений, минимизируя сумму квадратов ошибок восстановления:

$$ \sum_{i,j,k,l} \left( X_{i,j,k,l} - \sum_{r=1}^R A_{i,r} B_{j,r} C_{k,r} D_{l,r} \right)^2 \rightarrow \min_{A, B, C, D} $$

Минимизацию обычно выполняют методом покоординатного спуска (alternating least squares), по очереди фиксируя все, кроме одной, матрицы, вычисляя новые оптимальные значения для оставшейся, затем повторяя процесс, пока не выполнится критерий остановки (относительное изменение параметров меньше заданного порога).

Такой метод называется PARAFAC или CANDECOMP. При этом встаёт вопрос выбора числа компонентов $R$: чем оно больше, тем меньше, естественно, будет ошибка восстановления, но тем менее будут определены сами получающиеся компоненты. Принятым в сообществе способом проверки решения на переобученность является вариант кросс-валидации, когда по одной из осей (например, по той, по которой располагаются отдельные образцы) многомерный массив делят на 2 половины, выполняют PARAFAC для каждой и сравнивают компоненты по оставшимся осям. По идее, хорошо обученная модель должна описывать данные только одним способом (минимум ошибки восстановления в пространстве значений матриц $ A, B, C, D $ только один и компактно расположен), а переобученная - иметь широкую, как сковородка, область, в которой данные описываются с примерно одинаковой ошибкой. Может быть, определить, насколько узка область оптимума, можно аналитически?

В Numerical recipes: the art of scientific computing. W.H. Press, 3rd ed., Cambridge, UK; New York: Cambridge University Press, 2007 §§15.5.1-15.6 рассказывается о том, как оценивать доверительные интервалы параметров, находимых методом наименьших квадратов, если предположение о нормальности распределения остатков выполняется.

(Numerical Recipes о методе наименьших квадратов)

Пусть выполняется поиск параметров $ \mathbf{a} $ для множества точек $ { y_i = y(x_i | \mathbf{a}) \pm \sigma_i, i = 1 ... N } $. Тогда минимизировать нужно функционал:

$$ \chi^2(\mathbf{a}) = \sum_{i=1}^N \left( \frac{y_i - y(x_i | \mathbf{a})}{\sigma_i} \right)^2 $$

Вдали от оптимума решение производится шагами градиентного спуска $ \mathbf{a}_{n+1} = \mathbf{a}_n - c \cdot \nabla \chi^2(\mathbf{a}_n) $; вблизи оптимума функцию можно считать квадратичной формой

$$ \chi^2(\mathbf{a}) \approx c + \mathbf{d} \cdot \mathbf{a} + \frac{1}{2} \mathbf{a} \cdot \mathbf{D} \cdot \mathbf{a} $$

где компоненты вектора $\mathbf{d}$ - первые производные $\chi^2$ по параметрам $ d_k = \frac {\partial \chi^2}{\partial a_k} = -2 \sum_i \frac{y_i - y(x_i | \mathbf{a})}{\sigma_i^2} \frac{\partial y(x_i | \mathbf{a})}{\partial a_k} $, а матрицы $\mathbf{D}$ - вторые производные:

$$ D_{kl} = \frac{\partial^2 \chi^2}{\partial a_k \partial a_l} = 2 \sum_i \frac{1}{\sigma_i^2} \left( \frac{\partial y(x_i | \mathbf{a})}{\partial a_k} \frac{\partial y(x_i | \mathbf{a})}{\partial a_l} - (y_i - y(x_i | \mathbf{a})) \frac{\partial^2 y(x_i | \mathbf{a})}{\partial a_k \partial a_l} \right) $$

Причём Numerical Recipes советует игнорировать слагаемое $ (y_i - y(x_i | \mathbf{a})) \frac{\partial^2 y(x_i | \mathbf{a})}{\partial a_k \partial a_l} $, поскольку вблизи оптимума вторая производная должна быть близка к нулю, а $ y_i - y(x_i | \mathbf{a}) $ - случайная ошибка измерения. Так что если определить матрицу $ \mathbf{J} $ первых производных аппроксимируемой функции $ J_{ik} = \frac{\partial y(x_i)}{\partial a_k} $, то $ \mathbf{D} \approx \mathbf{J}^\intercal \mathbf{J} $.

Имея гессиану $ \mathbf H $, $ H_{k,l} = \frac{\partial^2\chi^2}{\partial a_k \partial a_l} $, посчитанную в найденном оптимуме, мы можем обратить её и получить матрицу ковариации $ \mathbf{C} = \mathbf{H}^{-1} $, после чего использовать диагональные элементы для оценки доверительных интервалов: $ \delta a_i = \sqrt{\Delta \chi^2_\nu} \sqrt{C_{ii}} $ (где $ \Delta \chi^2_\nu $ - обратная функция распределения хи-квадрат с данным числом степеней свободы).

Я попробовал выписать производные:

$$ \frac{\partial X_{i,j,k}}{\partial A_{m,r}} = \begin{cases} {B}_{j,r} {C}_{k,r},& i = m\\ 0,&\text{иначе} \end{cases} $$

и посчитать $ \mathbf{H} $ как $ \mathbf{J}^\intercal \mathbf{J} $, но матрица получилась необратимой. Код на R:

# install.packages('multiway')
library(multiway)

# взято из example(parafac):
# сформируем данные, удовлетворяющие требованиям PARAFAC
set.seed(3)
mydim <- c(50, 20, 5)
nf <- 3
Amat <- matrix(rnorm(mydim[1]*nf), nrow = mydim[1], ncol = nf)
Bmat <- matrix(runif(mydim[2]*nf), nrow = mydim[2], ncol = nf)
Cmat <- matrix(runif(mydim[3]*nf), nrow = mydim[3], ncol = nf)
Xmat <- tcrossprod(Amat, krprod(Cmat, Bmat))
Xmat <- array(Xmat, dim = mydim)
# добавим шум
Emat <- array(rnorm(prod(mydim)), dim = mydim)
Emat <- nscale(Emat, 0, ssnew = sumsq(Xmat)*.1)
X <- Xmat + Emat

# рассчитаем оценки Amat, Bmat, Cmat
pfac <- parafac(X, nfac = nf, nstart = 1)
pfac

# J: матрица производных X[i,j,k] по cbind(A,B,C)[l,r]
J <- array(data = 0, dim = c(mydim, sum(mydim), nf))
for (i in 1:mydim[1]) for (j in 1:mydim[2]) for (k in 1:mydim[3]) for (r in 1:nf) {
   # d X[i,j,k] / d A[i,r] = B[j,r]*C[k,r]
   J[i, j, k, 0+i, r] <- pfac$B[j,r] * pfac$C[k,r]
   # d X[i,j,k] / d B[j,r] = A[i,r]*C[k,r]
   J[i, j, k, mydim[1]+j, r] <- pfac$A[i,r] * pfac$C[k,r]
   # d X[i,j,k] / d C[k,r] = A[i,r]*B[j,r]
   J[i, j, k, mydim[1]+mydim[2]+k, r] <- pfac$A[i,r] * pfac$B[j,r]
}
# развернём первые 3 измерения ("образцы") и последние 2 ("переменные"), получим якобиан
dim(J) <- c(prod(mydim), sum(mydim)*nf)
# посчитаем оценку гессианы и попытаемся её обратить
C <- solve(crossprod(J))

Согласно сообщению об ошибке, число обусловленности у такой матрицы получается порядка $ 3 \cdot 10^{19} $. Обычное замечание в таком случае - это свидетельствует об очень плохо поставленной задаче, не удовлетворяющей модели - не похоже на правду: метод PARAFAC широко используется, о нём публикуются статьи в признанных в сообществе журналах, его продолжают развивать. Где я допускаю ошибку? Почему для данного вида задачи на наименьшие квадраты не получается найти матрицу ковариации?

[vpb]

aitap

Я $R$-язык не знаю, и код Ваш не понял. Но не в коде дело. Напишите пожалуйста, каковы параметры в вашем случае ($R$
и размерности по сомножителям, т.е. до каких значений пробегают $i$, $j$, $k$ ).

[aitap]

vpb

В данном примере: $R = 3$; $i_{max} = 50$; $j_{max} = 20$; $k_{max} = 5$, так что "экспериментальных точек" получается 5000, а переменных - 225.

На практике я хотел бы применить этот метод к задаче с $R \in \{ 2, ..., 6 \}$; $i_{max} = 61$; $j_{max} = 491$; $k_{max} = 58$, так что строить матрицу $\mathbf{J}$, а потом умножать её на себя будет не очень эффективно - такая матрица едва помещается в оперативную память. По-видимому, рано или поздно мне придётся явно выписать выражения для $ \partial^2 \chi^2 \over { \partial a_k \partial a_l } $ и посчитать именно их.

[vpb]

Кстати, позвольте еще сделать замечание, что "гессиан" --- мужского рода. Еще вопрос: у Вас были выписаны частные производные по $A_{m,r}$, а по $B_{j,r}$ ( и по $C$ тоже) тоже присутствуют ? Какого размера сама матрица $J$ в рассматриваемом примере ? (я задаю вопросы, чтоб быть уверенным, что правильно буду в будущем, когда более тщательно вникать, интерпретировать формулы).
Еще два вопроса: 1) Вы собираетесь найти разложение более-менее точно, или только приближенно ? Ибо в этой науке есть задача точного разложения, и задача приближения.
2) Вы знакомы с понятием тензорного произведения ?

[DeBill]

aitap в сообщении #1381593 писал(а):

Я попробовал выписать производные:

$$ \frac{\partial X_{i,j,k}}{\partial A_{m,r}} = \begin{cases} {B}_{j,r} {C}_{k,r},& i = m\\ 0,&\text{иначе} \end{cases} $$

Что бы это значило?
Если Вы считаете производную функционала "хи" - то неправильно сосчитано.
А если дифференцируете ФИКСИРОВАННЫЕ ДАННЫЕ "икс-и-жи-катые" - так производная вааще равна нулю...

[aitap]

vpb в сообщении #1381804 писал(а):

Кстати, позвольте еще сделать замечание, что "гессиан" --- мужского рода.

Спасибо, я буду это знать.

vpb в сообщении #1381804 писал(а):

Еще вопрос: у Вас были выписаны частные производные по $A_{m,r}$, а по $B_{j,r}$ ( и по $C$ тоже) тоже присутствуют ?

Да,

$$ \frac{\partial \hat{X}_{i,j,k}}{\partial B_{l,r}} = \begin{cases} A_{i,r}C_{k,r},&j = l\\ 0,\text{иначе} \end{cases} $$

$$ \frac{\partial \hat{X}_{i,j,k}}{\partial C_{l,r}} = \begin{cases} A_{i,r}B_{j,r},&k = l\\ 0,\text{иначе} \end{cases} $$

vpb в сообщении #1381804 писал(а):

Какого размера сама матрица $J$ в рассматриваемом примере ?

$ 5000 \times 225 $

vpb в сообщении #1381804 писал(а):

Еще два вопроса: 1) Вы собираетесь найти разложение более-менее точно, или только приближенно ? Ибо в этой науке есть задача точного разложения, и задача приближения.
2) Вы знакомы с понятием тензорного произведения ?

1) Приближённо: экспериментальные данные содержат шум. Мы надеемся, что он нормально и более-менее одинаково распределён.
2) К сожалению, нет. Если это необходимо для продолжения работы, я готов изучить пособия, которые Вы посоветуете.

DeBill в сообщении #1381810 писал(а):

Если Вы считаете производную функционала "хи" - то неправильно сосчитано.
А если дифференцируете ФИКСИРОВАННЫЕ ДАННЫЕ "икс-и-жи-катые" - так производная вааще равна нулю...

Прошу прощения, я неправильно выразился. Это производная одной из $ i_{max} \cdot j_{max} \cdot k_{max} $ модельных функций $ \hat{X}_{i,j,k}(\mathbf{A}, \mathbf{B}, \mathbf{C}) $ по её параметрам. А $ \chi^2 $, соответственно, равен $ \sum_{i,j,k} (X_{i,j,k} - \hat{X}_{i,j,k})^2 $.

[vpb]

aitap в сообщении #1381843 писал(а):

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

Необходимо, необходимо...
Понимаете, какая трагедия: вроде бы это везде написано, и все это знают, но нигде не написано действительно хорошо.

Давайте я Вам посоветую книгу для вычислительных математиков. Там написано весьма дурно, с моей точки зрения, зато у Вас будет мотив дальше грызть. Но учить, конечно, надо по книгам для математиков, а не для вычислительных математиков.
Подумаю, что делать. Книга такая: W.Hackbusch, Tensor spaces and numerical tensor calculus. (есть в либгене). Почитайте там пока из начала, первых трех глав, но не всё подряд, а что поймете. И заглядывайте в тему периодически, я буду писать разные мысли (прежде всего разберемся, почему гессиан вырожденный... хотя это и так, кажется, ясно).

-- 14.03.2019, 16:32 --

Собственно, почему гессиан вырожден, можно объяснить прямо сейчас. Рассмотрим, для примера, такую функцию от двух переменных : $f(x,y)=(xy-1)^2$. Заметим, что она обладает свойством "масштабной инвариантности" : $f(\lambda x,\lambda^{-1}y)=f(x,y)$, для любых $x,y$ и $\lambda\ne0$. Множество точек минимума ее --- это гипербола $xy=1$. Но в любой точке минимума гессиан вырожден. Действительно, если бы он был невырожден, то тогда бы он с необходимостью был положительно определен, и тогда бы данная точка была единственной точкой локального минимума в некоторой своей окрестности. Но это не так, значит он вырожден.

Совершенно аналогичное явление масштабной инвариантности имеет место и в вашем случае.

(Замечание о терминологии. Я выше писал "гессиан" там, где, строго говоря, надо было писать "матрица Гессе"; а гессиан --- это ее определитель. )

[DeBill]

aitap в сообщении #1381843 писал(а):

я неправильно выразился. Это производная одной из

Это я понял. Но - зачем?
В стартовом посте Вы выписали функционал, который надлежит минимизировать (вторая из формул).
Функционал нехорош - неквадратичный он по совокупности переменных. Там же Вы указали и рецепт (применяемый стандартно - от безысходности): фиксировать часть переменных (B,C,D) и минимизировать по остальным (по А), и т.д....
Так таки минимизироваться будет тот же самый функционал! И его производная по А-эм-эртым - вовсе не такая.
И соответствующее уравнение выглядит много приличней. И вроде, его можно записать прямо таки в матричном виде (типовом для МНК). И матрицы будут не таки ужасны, а маленькие. И даже обратимость будет - если, конечно, R меньше всех прочих размеров (и не произойдет случайно несчастный случай...).

[vpb]

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

Я думаю, в обсуждении спутались два вопроса: 1) Как найти предполагаемые величины $A=(A_{ri})$, $B$, $C$, и 2) Как найти доверительные интервалы для них. ТС хочет найти доверительные интервалы. Т.е., например, с какой степенью надежности можно считать, что величина $A_{12}$ лежит в интервале $(3,5)$, вроде того. Однако этот вопрос не имеет смысла, именно в силу масштабной инвариантности (или может ее лучше назвать масштабной неопределенностью ) задачи. Что находит отражение в вырожденности гессиана.

Наверное еще, коллега DeBill спутал гессиан по всей совокупности параметров $A,B,C$ с гессианом по группе переменных $A$ (или $B$, или $C$). Первое имеет смысл для оценки доверительных интервалов, а второе --- для минимизации целевой функции.

Чтоб неопределенность снять, надо наложить какие-то дополнительные условия нормировки, скажем $A_{r,1}=B_{r,1}=1$ для всех $r$, или же $||A_r||=||B_r||=||C_r||$ для всех $r$ (где $||x||$ --- евклидова норма вектора).

[aitap]

vpb в сообщении #1381863 писал(а):

Книга такая: W.Hackbusch, Tensor spaces and numerical tensor calculus. (есть в либгене). Почитайте там пока из начала, первых трех глав, но не всё подряд, а что поймете.

Большое спасибо, я изучу.

vpb в сообщении #1381863 писал(а):

Совершенно аналогичное явление масштабной инвариантности имеет место и в вашем случае.

Самое обидное здесь то, что я прекрасно знал, что факторы определены с точностью до перестановки и домножения на константу, осознав это ещё тогда, когда решал задачу их сравнения в ходе кросс-валидации. Это абсолютно точно отвечает на вопрос темы ("почему?"). Теперь мне гораздо понятнее, почему все делают кросс-валидацию и другие вариации на тему bootstrap.

DeBill в сообщении #1381928 писал(а):

Это я понял. Но - зачем?

Последовал рекомендации Numerical Recipes считать $ \frac{\partial^2 \chi^2}{\partial a_k \partial a_l} \approx \sum_i \frac{\partial \hat{y}_i}{\partial a_k}\frac{\partial \hat{y}_i}{\partial a_l} $ в окрестностях оптимума. Согласно книге, стоящий перед и без того малой $ \frac{\partial^2 \hat{y}_i}{\partial a_k \partial a_l} $ множитель $ y_i - \hat{y}_i $ только дестабилизирует расчёт, если в данных есть выбросы (а где их нет?).

DeBill в сообщении #1381928 писал(а):

И даже обратимость будет - если, конечно, R меньше всех прочих размеров (и не произойдет случайно несчастный случай...).

Поправьте меня, пожалуйста, если я ошибаюсь, но по-моему, решение не станет однозначным, поскольку для каждого решения $ \mathbf{A}, \mathbf{B}, \mathbf{C} $ я могу подставить $ \frac{\mathbf{A}}{xy}, \, x\mathbf{B}, \, y\mathbf{C} \; \forall x, y \in \mathbb{R}$ или поменять местами $ A_{i,r} \leftrightarrow A_{i,s}, B_{j,r} \leftrightarrow B_{j,s}, C_{k,r} \leftrightarrow C_{k,s} \, i \ne s; \, i, s \in \{ 1 ... R \}; \, \forall i,j,k $ и получить ещё $R!-1$ решений с точно такой же суммой квадратов ошибок.

vpb в сообщении #1382002 писал(а):

Чтоб неопределенность снять, надо наложить какие-то дополнительные условия нормировки, скажем $A_{r,1}=B_{r,1}=1$ для всех $r$, или же $||A_r||=||B_r||=||C_r||$ для всех $r$ (где $||x||$ --- евклидова норма вектора).

Как выразить такие накладываемые на задачу ограничения на языке функции потерь? Чем-то вроде $ \mathfrak{L}(\mathbf{A}, \mathbf{B}, \mathbf{C}) = \chi^2(\mathbf{A}, \mathbf{B}, \mathbf{C}) + \lambda ( | \sum_{r} ||A_r||-||B_r|| | + | \sum_{r} ||B_r||-||C_r|| | ) $ ? А как избавиться от неопределённости, вызванной перестановкой? Например, потребовать, чтобы $ || A_r || > || A_s ||, \; r > s $?

[DeBill]

aitap в сообщении #1382125 писал(а):

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

Нет, конечно, Вы (а также и vpb) здесь правы. Все дело в том, что произошла некая путаница - как и указал vpb, и виноват в ней ровно один я : мой текст относился ТОЛЬКО к задаче отыскания экстремума функционала; причем именно к той метОде, что предлагают "обычно" (т.е., минимизируем по одной переменной при фиксированных других; и вот тут Вы не можете размножать решения, они таки однозначны). Это все, конечно, не поможет при построении дов-х интервалов.

[vpb]

aitap в сообщении #1382125 писал(а):

Как выразить такие накладываемые на задачу ограничения на языке функции потерь? Чем-то вроде

Если речь об ограничениях $A_{r,1}=B_{r,_1}=1$, то функция потерь вообще ни при чем, просто некоторые переменные считаем фиксированными, а по остальным минимизируем. Если про $|A_r|=|B_r|=|C_r|$, то можно поступить так: рассмотреть функцию $\chi^2(A,B,C)+\lambda \sum_{r=1}^R (|A_r|^2+|B_r|^2+|C_r|^2)$, где $\lambda$ --- малое число, и искать ее минимумы. Эти минимумы тогда будут автоматически удовлетворять ограничению $|A_r|=|B_r|=|C_r|$ (следует из того, что минимум величины $x^2+y^2+z^2$ при условии, что $xyz=\operatorname{const}$ фиксировано, достигается тогда, когда $x=y=z$).

-- 16.03.2019, 08:29 --

aitap в сообщении #1382125 писал(а):

А как избавиться от неопределённости, вызванной перестановкой?

Никак. Упорядочить слагаемые, как Вам нравится. Но не между итерациями, а в конце всего вычисления.

-- 16.03.2019, 09:17 --

vpb в сообщении #1381863 писал(а):

Там написано весьма дурно, с моей точки зрения,

Насчет "весьма дурно" --- это я преувеличил. Не так уж плохо. Но как первый учебник не годится явно. Скажем так, средне. Текст там перегружен многоэтажными формулами, а толковых объяснений или хотя бы отсылок к хорошим учебникам (по линейной алгебре или функану) не хватает. В общем, монография, но с претензией (не достаточно основательной), что ее может читать и студент.

[aitap]

vpb в сообщении #1382236 писал(а):

Если про $|A_r|=|B_r|=|C_r|$, то можно поступить так: рассмотреть функцию $\chi^2(A,B,C)+\lambda \sum_{r=1}^R (|A_r|^2+|B_r|^2+|C_r|^2)$, где $\lambda$ --- малое число, и искать ее минимумы. Эти минимумы тогда будут автоматически удовлетворять ограничению $|A_r|=|B_r|=|C_r|$ (следует из того, что минимум величины $x^2+y^2+z^2$ при условии, что $xyz=\operatorname{const}$ фиксировано, достигается тогда, когда $x=y=z$).

Спасибо, я почти пришёл к слагаемому такого вида, но почему-то усомнился в его выводе. Это гораздо удобнее, чем требовать равенства нормы каждой из матриц какой-то определённой величине и нормировать $ \mathbf{X} $ перед выполнением расчёта. Тогда я попробую найти оптимум для модифицированной функции потерь и посмотрю на результаты.