用R代码编写Householder QR分解函数 [英] Writing a Householder QR factorization function in R code
问题描述
我正在研究一段代码，以查找R中矩阵的QR因式分解.
I am working on a piece of code to find the QR factorization of a matrix in 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))
myqr <- function(A) {
n <- nrow(A)
p <- ncol(A)
Q <- diag(n)
Inp <- diag(nrow = n, ncol = p)
for(k in c(1:ncol(A))) {
# 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
norm1<-sqrt(sum(col^2))
# Define the sign positive if a1 > 0 (-) else a1 < 0(+)
sign <- ifelse(col[1] >= 0, -1, +1)
# Calculate of the vector a_r
a_r <- col - sign * Inp[k:n,k] * norm1
# beta = 2 / ||a-r||^2
beta <- 2 / sum(t(a_r) %*% a_r)
# the next line of code calculates the matrix Q in every step
Q <- Q - beta *Q %*% c(rep(0,k-1),a_r) %*% t(c(rep(0,k-1),a_r))
# calculates the matrix R in each step
A[k:n,k:p] <- A[k:n,k:p] - beta * a_r %*% t(a_r) %*% A[k:n,k:p]
}
list(Q=Q,R=A)
}
但是，这里我没有在每一步中都计算出代表住户反射的矩阵H，也没有在每一步中都计算出矩阵A.
But, Here I have not calculated in every step the matrix H that represents the householder reflection, also I have not calculated the matrix A in every step.
作为H = I - 2 v v'，如果我乘以Q，我将得到
As H = I - 2 v v', if I multiply by Q I obtain
QH = Q - 2 (Qv) v' // multiplication on the left
HQ = Q - 2 v (Q'v)' // multiplication on the right
现在，此操作应该在每个步骤中都可以进行.但是，如果我考虑第一个矩阵H和他的第二个矩阵H1 ...，这些矩阵将小于第一个矩阵.为了避免这种情况，我使用了下一行代码:
Now, this operations should be work in every step. However if I consider the first matrix H and he the second matrix H1.... these matrices will be of smaller that the first one. In order to avoid that I have used the next line of code:
Q <- Q - beta * Q %*% c(rep(0,k-1),a_r) %*% t(c(rep(0,k-1),a_r))
但是，当我生成新矢量a_r时，在每个步骤中前一个k项均为零时，我不确定代码为什么能正常工作.
but, I am not sure why the code is working well, when I generate the new vector a_r with the first k entries of zeros at every step.
推荐答案
我认为您希望获得与qr.default返回的输出完全相同的输出，该输出使用紧凑的QR存储.但是后来我意识到您分别存储了Q和R因素.
I thought you want exactly the same output as returned by qr.default, which uses compact QR storage. But then I realized that you are storing Q and R factors separately.
通常，QR因式分解仅形成R而不形成Q.在下文中，我将描述QR分解的形成方式.对于那些对QR分解基本了解的人，请先阅读以下内容: lm():LINPACK/LAPACK中QR分解返回的qraux是什么，其中LaTeX中安排了整洁的数学公式.在下文中，我将假设一个人知道Householder反射是什么以及如何计算.
Normally, QR factorization only forms R but not Q. In the following, I will describe QR factorization where both are formed. For those who lack basic understanding of QR factorization, please read this first: lm(): What is qraux returned by QR decomposition in LINPACK / LAPACK, where there are neat math formulae arranged in LaTeX. In the following, I will assume that one knows what a Householder reflection is and how it is computed.
首先，Householder复制向量是H = I - beta * v v'(其中beta是根据您的代码计算的)，而不是H = I - 2 * v v'.
First of all, a Householder refection vector is H = I - beta * v v' (where beta is computed as in your code), not 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相乘.
Then, QR factorization A = Q R proceeds as (Hp ... H2 H1) A = R, where Q = H1 H2 ... Hp. To compute Q, we initialize Q = I (identity matrix), then multiply Hk on the right iteratively in the loop. To compute R, we initialize R = A and multiply Hk on the left iteratively in the loop.
现在，在第k次迭代中，我们在Q和A上进行了 rank-1矩阵更新:
Now, at k-th iteration, we have a rank-1 matrix update on Q and 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是全反射矢量的缩减后的非零部分.
v = c(rep(0, k-1), a_r), where a_r is the reduced, non-zero part of a full reflection vector.
您拥有的代码正在以一种残酷的方式进行此类更新:
The code you have is doing such update in a brutal force:
Q <- Q - beta * Q %*% c(rep(0,k-1), a_r) %*% t(c(rep(0,k-1),a_r))
首先填充a_r以获取全反射矢量，然后对整个矩阵执行等级1更新.但实际上，我们可以舍弃这些零并编写(如果不清楚，可以做一些矩阵代数):
It first pads a_r to get the full reflection vector and performs the rank-1 update on the whole matrix. But actually we can drop off those zeros and write (do some matrix algebra if unclear):
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的一部分.
so that only a fraction of Q and A are updated.
- 您已经大量使用了
t()和"%*%"！但是几乎所有它们都可以用crossprod()或tcrossprod()代替.这样就消除了显式转置t()并提高了内存效率； -
您不必初始化另一个对角矩阵
Inp.要获得户主反射向量a_r，您可以替换
- You have used
t()and"%*%"a lot! But almost all of them can be replaced bycrossprod()ortcrossprod(). This eliminates the explicit transposet()and is more memory efficient; You initialize another diagonal matrix
Inpwhich is not necessary. To get householder reflection vectora_r, you can replace
sign <- ifelse(col[1] >= 0, -1, +1)
a_r <- col - sign * Inp[k:n,k] * norm1
作者
a_r <- col; a_r[1] <- a_r[1] + sign(a_r[1]) * norm1
其中sign是R基函数.
## 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
首先是Thin-QR版本:
First for thin-QR version:
## 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版本:
Now the full QR version:
## 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返回的标准结果:
Now let's check standard result returned by 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
所以我们的结果是正确的！
So our results are correct!
这篇关于用R代码编写Householder QR分解函数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
