问题描述
我正在研究一段代码,以查找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
Inp
which 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分解函数的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持!