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

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() на больших системах.

Ответы

Ответ 1

Введение

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

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

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

Путь LINPACK

Отражение домохозяйства вычисляется как H i= я - v i v i T/p 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= я - p i v i v i T где я - единичная матрица, 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...]

Обратите внимание на поле $pivot; мы вернемся к этому. Теперь мы генерируем 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 выше.

Ответ 2

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

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

  • Всегда поддерживайте Q и R, явно разделенные, независимо от того, были ли они вычислены в первый раз или после каждого обновления за счет дополнительных провалов, но сохраняя ограниченную память.
  • Придерживайтесь компактного формата LAPACK, но каждый раз, когда появляется новое обновление, сохраняйте список всех этих мини-наборов рефлекторов обновлений. В точке решения системы можно было бы сделать большой Q '* c ie 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. в качестве входных данных.