标签: qr-decomposition

我在CSS中有以下转换矩阵

// rotate the element 60deg
element.style.transform = "matrix(0.5,0.866025,-0.866025,0.5,0,0)"

我可以用这个找到旋转...

// where a = [0.710138,0.502055,-0.57735,1,0,0]
var rotation = ((180/Math.PI) * Math.atan2( ((0*a[2])+(1*a[3])),((0*a[0])-(1*a[1]))) - 90
console.log(rotation); // ~60

同样如果......

// skew(30deg,-50deg) 
element.style.transform = "matrix(1,-1.19175,0.57735,1,0,0)"

var skewY = ((180/Math.PI) * Math.atan2( ((0*a[2])+(1*a[3])),((0*a[0])-(1*a[1]))) - 90;
var skewX = (180/Math.PI) * Math.atan2( ((1*a[2])+(0*a[3])),((1*a[0])-(0*a[1])));

console.log([skewX,skewY]); // ~= [30,-50] 

然而,只要我使用倾斜和旋转,一切都变得奇怪,因为旋转公式与倾斜公式相同...所以公式不能正确.

如何确定已应用两个属性的旋转和倾斜,我所知道的是矩阵变换.

同时缩放也搞砸了我的偏斜值,我认为不应该这样.

我试图在R中计算P任意N x J矩阵的投影矩阵S:

P = S (S'S) ^ -1 S'

我一直试图用以下功能执行此操作:

P <- function(S){
  output <- S %*% solve(t(S) %*% S) %*% t(S)
  return(output)
}

但是当我使用它时,我得到的错误看起来像这样:

# Error in solve.default(t(S) %*% S, t(S), tol = 1e-07) : 
#  system is computationally singular: reciprocal condition number = 2.26005e-28

我认为这是数值下溢和/或不稳定的结果在许多地方一样讨论R-帮助和这里,但我没有足够的经验使用SVD或QR分解来解决这个问题,要不然就把这个现有的代码到行动.我也尝试了建议的代码,即将solve写成一个系统:

output <- S %*% solve (t(S) %*% S, t(S), tol=1e-7)

但它仍然无效.任何建议,将不胜感激.

我很确定我的矩阵应该是可逆的并且没有任何共线性,只是因为我尝试用正交虚拟变量矩阵进行测试,但它仍然不起作用.

另外,我想将它应用于相当大的矩阵,所以我正在寻找一个简洁的通用解决方案.

我在处理排名不足的情况时非常依赖该qr()函数，但最近遇到了一些它无法正常工作的示例。考虑下面的矩阵badX：

badX <-
structure(c(-1.641906809157e-10, 0, 0, 0, 0, -0.5, 0, 0, -1.10482935525559e-16, 
            0, -3.06266685765538e-17, 0, -4.83736007092039e-17, 0, -3.14414492582296e-18, 
            -3.06158275836099e-18), dim = c(4L, 4L), dimnames = list(c("(Intercept)", 
            "A2", "A3", "B2"), NULL))

我们不能使用以下方法反转该矩阵solve()：

solve(badX)
## Error in solve.default(badX): system is computationally singular: reciprocal condition number = 5.55308e-18

然而qr()，其相关例程认为该矩阵的秩为 4，并且可以对其求逆：

qr(badX)$rank
## [1] 4

qr.solve(badX)
##             [,1] [,2]          [,3]          [,4]
## [1,] -6090479645    0  2.197085e+10  7.366741e+10
## [2,]           0   -2  0.000000e+00  0.000000e+00
## [3,]           0 …

我在MATLAB中编码QR分解算法,只是为了确保我的机制正确.这是main函数的代码:

    function [Q,R] = QRgivens(A)
        n = length(A(:,1));
        Q = eye(n);
        R = A;

        for j = 1:(n-1)
            for i = n:(-1):(j+1)
                G = eye(n);
                [c,s] = GivensRotation( A(i-1,j),A(i,j) );
                G(i-1,(i-1):i) = [c s];
                G(i,(i-1):i)   = [-s c];
                Q = Q*G';
                R = G*R;
            end
        end
    end

子函数GivensRotation如下:

    function [c,s] = GivensRotation(a,b)
        if b == 0
            c = 1;
            s = 0;
        else
            if abs(b) > abs(a)
                r = -a / b;
                s = 1 / sqrt(1 + r^2);
                c …

避免阵列分配有利于提高性能.但是,我还没有理解什么是最有效的方法,可以执行矩阵的QR分解A.(注意:需要Q和R矩阵)

简单地使用 Q, R = qr(A) 可能不是最好的主意,因为它分配Q和R,两者都可以重新分配.

该函数qrfact允许以打包格式存储因子分解.但是,我之后仍会写: F = qrfact(A); Q = F[:Q]; R = F[:R] 再次为Q和分配新数组R.最后,文档还建议了该qrfact!函数,它通过覆盖输入A来节省空间,而不是创建副本.但是,如果一个人使用F = qrfact!(A) 覆盖写入A是没有用的,因为它不是,Q或者R哪一个(具体来说,我)需要.

所以我的两个问题是:

1. 如果您只关心矩阵Q并且R重新分配它们没有问题,那么执行QR分解的最佳/最有效方法是什么？

2. A当一个人打电话时,矩阵中实际写的是什么qrfact!(A)？

我有一个关于将 matlab 函数转换为 R 的问题，我希望有人能提供帮助。

matlab 和 R 中使用的标准 QR 分解称为 qr()。据我了解，用两种语言执行 qr 分解的标准方法是：

Matlab： [Q,R] = qr(A)满足QR=A

回复：

z <- qr(A)
Q <- qr.Q(z)
R <- qr.R(z)

两者都为我提供了相同的结果，不幸的是，这不是我需要的。我需要的是这样的：

Matlab： [Q,R,e] = qr(A,0)产生经济规模的分解，其中 e 是排列向量，因此 A(:,e) = Q*R。

答：没有任何线索

我尝试将[Q,R,E] = qr(A)与

z <- qr(A);
Q <- qr.Q(z);
R <- qr.R(z);
E <- diag(ncol(A))[z$pivot]

变量 Q 和 E 的结果似乎相同（但 R 的结果不同）。因此，根据定义的输入/输出，将会有不同的结果（这是有道理的）。

所以我的问题是： R中有没有一种方法可以在Matlab中模拟[Q,R,e]=qr(A,0)？

我曾尝试深入研究 matlab 函数，但它导致了一条漫长而曲折的无尽函数定义之路，我希望有一个更好的解决方案。

任何帮助将不胜感激，如果我错过了一些明显的事情，我深表歉意。

我正在尝试扩展lwr()包的功能McSptial，它适合将加权回归作为非参数估计。在函数的核心中lwr()，它使用以下方法反转矩阵solve()，它使用QR 分解而不是 QR 分解来我想更改它，但无法弄清楚如何从 QR 分解中获取帽子矩阵（或其他导数）。

有数据：

set.seed(0); xmat <- matrix(rnorm(500), nrow=50)    ## model matrix
y <- rowSums(rep(2:11,each=50)*xmat)    ## arbitrary values to let `lm.wfit` work
w <- runif(50, 1, 2)    ## weights

该lwr()功能如下：

xmat2 <- w * xmat
xx <- solve(crossprod(xmat, xmat2))
xmat1 <- tcrossprod(xx, xmat2)
vmat <- tcrossprod(xmat1)

我需要的值，例如：

sum((xmat[1,] %*% xmat1)^2)
sqrt(diag(vmat))

目前我使用reg <- lm.wfit(x=xmat, y=y, w=w)但无法设法恢复在我看来是帽子矩阵（xmat1）的东西。

我正在尝试学习QR分解,但无法在不借助传统矩阵计算的情况下弄清楚如何获得beta_hat的方差.我正在练习iris数据集,这是我到目前为止所拥有的:

y<-(iris$Sepal.Length)
x<-(iris$Sepal.Width)
X<-cbind(1,x)
n<-nrow(X)
p<-ncol(X)
qr.X<-qr(X)
b<-(t(qr.Q(qr.X)) %*% y)[1:p]
R<-qr.R(qr.X)
beta<-as.vector(backsolve(R,b))
res<-as.vector(y-X %*% beta)

谢谢你的帮助!

来自LAPACK的DGEQRF和SGEQRF以打包格式返回QR因式分解的Q部分。打开包装似乎需要O(k^3)步骤（k个低档产品），而且似乎不是很简单。另外，k对我来说，进行顺序乘法的数值稳定性还不清楚。

LAPACK是否包括用于解包Q的子例程，如果没有，我应该怎么做？

我正在尝试编写一个使用QR分解来解决多元回归的函数.输入:y向量和X矩阵; 输出:b,e,R ^ 2.到目前为止,我已经得到了这个并且非常困难; 我想我已经让事情变得太复杂了:

QR.regression <- function(y, X) {
X <- as.matrix(X)
y <- as.vector(y)
p <- as.integer(ncol(X))
if (is.na(p)) stop("ncol(X) is invalid")
n <- as.integer(nrow(X))
if (is.na(n)) stop("nrow(X) is invalid")
nr <- length(y)
nc <- NCOL(X)

# Householder 
for (j in seq_len(nc)) {
id <- seq.int(j, nr)
sigma <- sum(X[id, j]^2)
s <- sqrt(sigma)
diag_ej <- X[j, j]
gamma <- 1.0 / (sigma + abs(s * diag_ej))
kappa <- if (diag_ej < 0) s else -s
X[j,j] <- …