Мистифицируется qr.Q (): что такое ортонормированная матрица в «компактной» форме? - PullRequest
Новая
       312

Мистифицируется qr.Q (): что такое ортонормированная матрица в «компактной» форме?

18 голосов
/ 13 июня 2010

R имеет функцию qr(), которая выполняет QR-декомпозицию с использованием LINPACK или LAPACK (по моему опыту, последняя на 5% быстрее). Основным возвращаемым объектом является матрица "qr", которая содержится в верхней треугольной матрице R (т.е. R=qr[upper.tri(qr)]). Все идет нормально. Нижняя треугольная часть qr содержит Q "в компактной форме". Можно извлечь Q из разложения qr, используя qr.Q(). Я хотел бы найти обратное значение qr.Q(). Другими словами, у меня есть Q и R, и я хотел бы поместить их в объект "qr". R тривиально, а Q нет. Цель состоит в том, чтобы применить к нему qr.solve(), что намного быстрее, чем solve() в больших системах.

Ответы [ 2 ]

19 голосов
/ 17 сентября 2011

Введение

R использует LINPACK dqrdc по умолчанию или LAPACK DGEQP3 подпрограмма, если указана, для вычисления разложения QR.Обе процедуры вычисляют разложение, используя отражения Домохозяина.Матрица mxn A раскладывается на ортогональную матрицу mxn экономичного размера (Q) и верхнюю треугольную матрицу nxn (R) как A = QR, где Q может быть вычислено как произведение t матриц отражения Домохозяина, причем t является наименьшимm-1 и n: Q = H 1 H 2 ... H t .

Каждая матрица отражений H i может быть представлен вектором длины (m-i + 1).Например, H 1 требует вектора длины m для компактного хранения.Все записи этого вектора, кроме одной, помещаются в первый столбец нижнего треугольника входной матрицы (диагональ используется коэффициентом R).Следовательно, каждому отражению требуется еще один скаляр памяти, и это обеспечивается вспомогательным вектором (называемым $qraux в результате R * qr).

Используемое компактное представление отличается междупроцедуры LINPACK и LAPACK.

Путь LINPACK

Отражение домохозяина вычисляется как H i = I - v i v i T / p i , где I - единичная матрица, p i - соответствующая запись в $qraux, а v i выглядит следующим образом:

  • v i [1..i-1] = 0,
  • v i [i] = p i
  • v i [i + 1: m] = A [i + 1..m, i] (т.е., столбец нижнего треугольника A после вызова qr)

Пример LINPACK

Давайте рассмотрим пример из статьи QR-декомпозиции в Википедиив R.

Разлагаемая матрица имеет вид 

> A <- matrix(c(12, 6, -4, -51, 167, 24, 4, -68, -41), nrow=3)
> A
     [,1] [,2] [,3]
[1,]   12  -51    4
[2,]    6  167  -68
[3,]   -4   24  -41

Мы делаем разложение, и наиболее важныеСоотношения результата показаны ниже: 

> Aqr = qr(A)
> Aqr
$qr
            [,1]         [,2] [,3]
[1,] -14.0000000  -21.0000000   14
[2,]   0.4285714 -175.0000000   70
[3,]  -0.2857143    0.1107692  -35

[snip...]

$qraux
[1]  1.857143  1.993846 35.000000

[snip...]

Это разложение было выполнено (под одеялом) путем вычисления двух отражений Домохозяина и умножения их на A, чтобы получить R. Теперь мы воссоздадим отражения из информации в$qr. 

> p = Aqr$qraux   # for convenience
> v1 <- matrix(c(p[1], Aqr$qr[2:3,1]))
> v1
           [,1]
[1,]  1.8571429
[2,]  0.4285714
[3,] -0.2857143

> v2 <- matrix(c(0, p[2], Aqr$qr[3,2]))
> v2
          [,1]
[1,] 0.0000000
[2,] 1.9938462
[3,] 0.1107692

> I = diag(3)   # identity matrix
> H1 = I - v1 %*% t(v1)/p[1]   # I - v1*v1^T/p[1]
> H2 = I - v2 %*% t(v2)/p[2]   # I - v2*v2^T/p[2]

> Q = H1 %*% H2
> Q
           [,1]       [,2]        [,3]
[1,] -0.8571429  0.3942857  0.33142857
[2,] -0.4285714 -0.9028571 -0.03428571
[3,]  0.2857143 -0.1714286  0.94285714

Теперь давайте проверим правильность вычисленного выше Q: 

> qr.Q(Aqr)
           [,1]       [,2]        [,3]
[1,] -0.8571429  0.3942857  0.33142857
[2,] -0.4285714 -0.9028571 -0.03428571
[3,]  0.2857143 -0.1714286  0.94285714

Выглядит хорошо!Мы также можем убедиться, что QR равно A. 

> R = qr.R(Aqr)   # extract R from Aqr$qr
> Q %*% R
     [,1] [,2] [,3]
[1,]   12  -51    4
[2,]    6  167  -68
[3,]   -4   24  -41

Путь LAPACK

Отражение домохозяина вычисляется как H i = I - p i v i v i T , где I - единичная матрица, p i - соответствующая запись в $qrauxи v i выглядит следующим образом:

  • v i [1..i-1] = 0,
  • v i [i] = 1
  • v i [i + 1: m] = A [i + 1..m, i] (то есть столбецнижний треугольник A после вызова qr)

При использовании подпрограммы LAPACK в R используется другой поворот: используется поворот столбца, поэтому декомпозиция решает другую связанную проблему: AP =QR, где P - матрица перестановок .

Пример LAPACK

В этом разделе делается тот же пример, что и раньше. 

> A <- matrix(c(12, 6, -4, -51, 167, 24, 4, -68, -41), nrow=3)
> Bqr = qr(A, LAPACK=TRUE)
> Bqr
$qr
            [,1]        [,2]       [,3]
[1,] 176.2554964 -71.1694118   1.668033
[2,]  -0.7348557  35.4388886  -2.180855
[3,]  -0.1056080   0.6859203 -13.728129

[snip...]

$qraux
[1] 1.289353 1.360094 0.000000

$pivot
[1] 2 3 1

attr(,"useLAPACK")
[1] TRUE

[snip...]

Обратите внимание на *Поле 1134 *;мы вернемся к этому.Теперь мы генерируем Q из информации Aqr. 

> p = Bqr$qraux   # for convenience
> v1 = matrix(c(1, Bqr$qr[2:3,1]))
> v1
           [,1]
[1,]  1.0000000
[2,] -0.7348557
[3,] -0.1056080


> v2 = matrix(c(0, 1, Bqr$qr[3,2]))
> v2
          [,1]
[1,] 0.0000000
[2,] 1.0000000
[3,] 0.6859203


> H1 = I - p[1]*v1 %*% t(v1)   # I - p[1]*v1*v1^T
> H2 = I - p[2]*v2 %*% t(v2)   # I - p[2]*v2*v2^T
> Q = H1 %*% H2
           [,1]        [,2]       [,3]
[1,] -0.2893527 -0.46821615 -0.8348944
[2,]  0.9474882 -0.01602261 -0.3193891
[3,]  0.1361660 -0.88346868  0.4482655

Еще раз, Q, вычисленный выше, согласуется с предоставленным R Q. 

> qr.Q(Bqr)
           [,1]        [,2]       [,3]
[1,] -0.2893527 -0.46821615 -0.8348944
[2,]  0.9474882 -0.01602261 -0.3193891
[3,]  0.1361660 -0.88346868  0.4482655

Наконец, давайте вычислим QR. 

> R = qr.R(Bqr)
> Q %*% R
     [,1] [,2] [,3]
[1,]  -51    4   12
[2,]  167  -68    6
[3,]   24  -41   -4

Заметили разницу?QR - это A с переставленными столбцами в порядке, указанном в Bqr$pivot выше.

0 голосов
/ 27 марта 2012

Я исследовал эту проблему, как спрашивает ОП, и не думаю, что это возможно. По сути, вопрос OP заключается в том, можно ли восстановить H1 H2 ... Ht, имея явно вычисленный Q. Я не думаю, что это возможно без вычисления QR с нуля, но мне также было бы очень интересно узнать, существует ли такое решение.

У меня проблема, похожая на OP, но в другом контексте мой итеративный алгоритм должен изменять матрицу A, добавляя столбцы и / или строки. В первый раз QR вычисляется с использованием DGEQRF и, следовательно, компактного формата LAPACK. После того, как матрица А мутирует, например, с новыми рядами я могу быстро создать новый набор отражателей или ротаторов, которые уничтожат ненулевые элементы самой низкой диагонали моего существующего R и построят новый R, но теперь у меня есть набор H1_old H2_old ... Hn_old и H1_new H2_new ... Hn_new (и аналогично тау), которые нельзя смешать в одном компактном представлении QR. У меня есть две возможности, и, возможно, у ОП есть те же две возможности:

  1. Всегда поддерживайте явное разделение Q и R, будь то при первом вычислении или после каждого обновления, за счет дополнительных флопов, но с сохранением необходимой ограниченной памяти.
  2. Придерживайтесь компактного формата LAPACK, но каждый раз, когда появляется новое обновление, сохраняйте список всех этих мини-наборов отражателей обновлений. В момент решения системы нужно выполнить большую Q '* c, т.е. H1_u3 * H2_u3 * ... * Hn_u3 * H1_u2 * H2_u2 * ... * Hn_u2 * H1_u1 * H2_u1 ... * Hn_u1 * H1 * H2 * ... * Hn * c, где ui - номер обновления QR, и это потенциально много умножений и памяти для отслеживания, но определенно самый быстрый способ.

Длинный ответ от Дэвида в основном объясняет, что такое компактный формат QR, но не как добраться до этого компактного формата QR, имеющего в качестве входных данных явно вычисленные Q и R.

...