用 R 代码编写 Householder QR 分解函数
Writing a Householder QR factorization function in R code
我正在编写一段代码来查找 R 中矩阵的 QR 因式分解。
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
。
作为H = I - 2 v v'
,如果我乘以Q
我得到
QH = Q - 2 (Qv) v' // multiplication on the left
HQ = Q - 2 v (Q'v)' // multiplication on the right
现在,这个操作应该在每一步都有效。但是,如果我考虑第一个矩阵 H
而他考虑第二个矩阵 H1
... 这些矩阵将比第一个矩阵小。为了避免我使用了下一行代码:
Q <- Q - beta * Q %*% c(rep(0,k-1),a_r) %*% t(c(rep(0,k-1),a_r))
但是,当我在每一步的第一个 k
个零条目生成新向量 a_r
时,我不确定为什么代码运行良好。
我以为您想要与 qr.default
返回的输出完全相同,它使用紧凑的 QR 存储。但后来我意识到你正在分别存储 Q
和 R
因素。
通常情况下,QR 分解只形成 R
而不是 Q
。下面,我将描述形成两者的 QR 因式分解。对 QR 分解缺乏基本了解的朋友,请先阅读:,其中有用 LaTeX 排列的整齐的数学公式。在下文中,我假设您知道 Householder 反射是什么以及它是如何计算的。
QR 分解程序
首先,Householder 反射向量是 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 次迭代,我们在 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
得到全反射向量,对整个矩阵进行rank-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
来自
a_r <- col; a_r[1] <- a_r[1] + sign(a_r[1]) * norm1
其中 sign
是 R 基函数。
QR 分解的 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
第一个用于薄 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
现在完整的二维码版本:
## 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
所以我们的结果是正确的!
我正在编写一段代码来查找 R 中矩阵的 QR 因式分解。
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
。
作为H = I - 2 v v'
,如果我乘以Q
我得到
QH = Q - 2 (Qv) v' // multiplication on the left
HQ = Q - 2 v (Q'v)' // multiplication on the right
现在,这个操作应该在每一步都有效。但是,如果我考虑第一个矩阵 H
而他考虑第二个矩阵 H1
... 这些矩阵将比第一个矩阵小。为了避免我使用了下一行代码:
Q <- Q - beta * Q %*% c(rep(0,k-1),a_r) %*% t(c(rep(0,k-1),a_r))
但是,当我在每一步的第一个 k
个零条目生成新向量 a_r
时,我不确定为什么代码运行良好。
我以为您想要与 qr.default
返回的输出完全相同,它使用紧凑的 QR 存储。但后来我意识到你正在分别存储 Q
和 R
因素。
通常情况下,QR 分解只形成 R
而不是 Q
。下面,我将描述形成两者的 QR 因式分解。对 QR 分解缺乏基本了解的朋友,请先阅读:
QR 分解程序
首先,Householder 反射向量是 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 次迭代,我们在 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
得到全反射向量,对整个矩阵进行rank-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
来自
a_r <- col; a_r[1] <- a_r[1] + sign(a_r[1]) * norm1
其中
sign
是 R 基函数。
QR 分解的 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
第一个用于薄 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
现在完整的二维码版本:
## 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
所以我们的结果是正确的!