Я думал, вам нужен точно такой же результат, как и у qr.default
, который использует компактное хранилище QR. Но потом я понял, что вы храните Q
и R
факторы отдельно.
Обычно факторизация QR формирует только R
, но не Q
. Далее я опишу QR-факторизацию, в которой формируются обе. Для тех, кому не хватает базового понимания факторизации QR, сначала прочитайте это: lm(): что такое qraux, возвращаемый декомпозицией QR в LINPACK/LAPACK, где есть аккуратные математические формулы в LaTeX. В дальнейшем я буду предполагать, что кто-то знает, что такое отражение Хаусхолдера и как оно вычисляется.
Процедура QR-факторизации
Прежде всего, вектор отражения Хаусхолдера равен H = I - beta * v v'
(где beta
вычисляется, как в вашем коде), а не H = I - 2 * v v'
.
Затем QR-факторизация A = Q R
продолжается как (Hp ... H2 H1) A = R
, где Q = H1 H2 ... Hp
. Чтобы вычислить Q
, мы инициализируем Q = I
(единичную матрицу), затем итеративно умножаем Hk
справа в цикле. Чтобы вычислить R, мы инициализируем R = A
и итеративно умножаем Hk
слева в цикле.
Теперь, на k-й итерации, у нас есть обновление матрицы ранга 1 для Q
и A
:
Q := Q Hk = Q (I - beta v * v') = Q - (Q v) (beta v)'
A := Hk A = (I - beta v * v') A = A - (beta v) (A' v)'
v = c(rep(0, k-1), a_r)
, где a_r
— сокращенная ненулевая часть полного вектора отражения.
Код у вас делает такое обновление в зверской силе:
Q <- Q - beta * Q %*% c(rep(0,k-1), a_r) %*% t(c(rep(0,k-1),a_r))
Сначала он дополняет a_r
, чтобы получить полный вектор отражения, и выполняет обновление ранга 1 для всей матрицы. Но на самом деле мы можем отбросить эти нули и написать (если неясно, проделать некоторую матричную алгебру):
Q[,k:n] <- Q[,k:n] - tcrossprod(Q[, k:n] %*% a_r, beta * a_r)
A[k:n,k:p] <- A[k:n,k:p] - tcrossprod(beta * a_r, crossprod(A[k:n,k:p], a_r))
так что обновляется только часть Q
и A
.
Несколько других комментариев к вашему коду
- Вы много использовали
t()
и "%*%"
! Но почти все их можно заменить на crossprod()
или tcrossprod()
. Это устраняет явное транспонирование t()
и более эффективно использует память;
Вы инициализируете другую диагональную матрицу Inp
, в которой нет необходимости. Чтобы получить вектор отражения домохозяина a_r
, вы можете заменить
sign <- ifelse(col[1] >= 0, -1, +1)
a_r <- col - sign * Inp[k:n,k] * norm1
by
a_r <- col; a_r[1] <- a_r[1] + sign(a_r[1]) * norm1
где sign
— базовая функция R.
R-код для QR-факторизации
## QR factorization: A = Q %*% R
## if `complete = FALSE` (default), return thin `Q`, `R` factor
## if `complete = TRUE`, return full `Q`, `R` factor
myqr <- function (A, complete = FALSE) {
n <- nrow(A)
p <- ncol(A)
Q <- diag(n)
for(k in 1:p) {
# extract the kth column of the matrix
col <- A[k:n,k]
# calculation of the norm of the column in order to create the vector r
norm1 <- sqrt(drop(crossprod(col)))
# Calculate of the reflection vector a-r
a_r <- col; a_r[1] <- a_r[1] + sign(a_r[1]) * norm1
# beta = 2 / ||a-r||^2
beta <- 2 / drop(crossprod(a_r))
# update matrix Q (trailing matrix only) by Householder reflection
Q[,k:n] <- Q[,k:n] - tcrossprod(Q[, k:n] %*% a_r, beta * a_r)
# update matrix A (trailing matrix only) by Householder reflection
A[k:n, k:p] <- A[k:n, k:p] - tcrossprod(beta * a_r, crossprod(A[k:n,k:p], a_r))
}
if (complete) {
A[lower.tri(A)] <- 0
return(list(Q = Q, R = A))
}
else {
R <- A[1:p, ]; R[lower.tri(R)] <- 0
return(list(Q = Q[,1:p], R = R))
}
}
Теперь давайте проведем тест:
X <- structure(c(0.8147, 0.9058, 0.127, 0.9134, 0.6324, 0.0975, 0.2785,
0.5469, 0.9575, 0.9649, 0.1576, 0.9706, 0.9572, 0.4854, 0.8003
), .Dim = c(5L, 3L))
# [,1] [,2] [,3]
#[1,] 0.8147 0.0975 0.1576
#[2,] 0.9058 0.2785 0.9706
#[3,] 0.1270 0.5469 0.9572
#[4,] 0.9134 0.9575 0.4854
#[5,] 0.6324 0.9649 0.8003
Сначала для тонкой версии QR:
## thin QR factorization
myqr(X)
#$Q
# [,1] [,2] [,3]
#[1,] -0.49266686 -0.4806678 0.17795345
#[2,] -0.54775702 -0.3583492 -0.57774357
#[3,] -0.07679967 0.4754320 -0.63432053
#[4,] -0.55235290 0.3390549 0.48084552
#[5,] -0.38242607 0.5473120 0.03114461
#
#$R
# [,1] [,2] [,3]
#[1,] -1.653653 -1.1404679 -1.2569776
#[2,] 0.000000 0.9660949 0.6341076
#[3,] 0.000000 0.0000000 -0.8815566
Теперь полная QR-версия:
## full QR factorization
myqr(X, complete = TRUE)
#$Q
# [,1] [,2] [,3] [,4] [,5]
#[1,] -0.49266686 -0.4806678 0.17795345 -0.6014653 -0.3644308
#[2,] -0.54775702 -0.3583492 -0.57774357 0.3760348 0.3104164
#[3,] -0.07679967 0.4754320 -0.63432053 -0.1497075 -0.5859107
#[4,] -0.55235290 0.3390549 0.48084552 0.5071050 -0.3026221
#[5,] -0.38242607 0.5473120 0.03114461 -0.4661217 0.5796209
#
#$R
# [,1] [,2] [,3]
#[1,] -1.653653 -1.1404679 -1.2569776
#[2,] 0.000000 0.9660949 0.6341076
#[3,] 0.000000 0.0000000 -0.8815566
#[4,] 0.000000 0.0000000 0.0000000
#[5,] 0.000000 0.0000000 0.0000000
Теперь давайте проверим стандартный результат, возвращаемый qr.default
:
QR <- qr.default(X)
## thin R factor
qr.R(QR)
# [,1] [,2] [,3]
#[1,] -1.653653 -1.1404679 -1.2569776
#[2,] 0.000000 0.9660949 0.6341076
#[3,] 0.000000 0.0000000 -0.8815566
## thin Q factor
qr.Q(QR)
# [,1] [,2] [,3]
#[1,] -0.49266686 -0.4806678 0.17795345
#[2,] -0.54775702 -0.3583492 -0.57774357
#[3,] -0.07679967 0.4754320 -0.63432053
#[4,] -0.55235290 0.3390549 0.48084552
#[5,] -0.38242607 0.5473120 0.03114461
## full Q factor
qr.Q(QR, complete = TRUE)
# [,1] [,2] [,3] [,4] [,5]
#[1,] -0.49266686 -0.4806678 0.17795345 -0.6014653 -0.3644308
#[2,] -0.54775702 -0.3583492 -0.57774357 0.3760348 0.3104164
#[3,] -0.07679967 0.4754320 -0.63432053 -0.1497075 -0.5859107
#[4,] -0.55235290 0.3390549 0.48084552 0.5071050 -0.3026221
#[5,] -0.38242607 0.5473120 0.03114461 -0.4661217 0.5796209
Так что наши результаты верны!
person
Zheyuan Li
schedule
04.10.2016