如何對Matrix包中的「sparseMatrix」類進行QR分解？

Question

我想用我用B<-as(A, "sparseMatrix")創建的矩陣上的Matrix:::qr()函數進行QR分解。我知道我可以用Matrix:::qr.R()得到R矩陣。不過，我也需要Q Matrix。 Matrix包中似乎沒有qr.Q（）函數。我如何獲得Q矩陣？

Answer 1

Q矩陣實際上存儲在V插槽中。看來目前的R Matrix版本包含一個bug ---在進行qr分解之前，它只是神奇地將零行添加到矩陣a中。我希望開發者能夠來解釋它。因此，下面的代碼可以幫助您恢復R和Q：

gx.qr.Q <- function(a){ 
    if(is(a, "qr")){ 
    return(qr.Q(a, complete = T)) 
    }else if(is(a, "sparseQR")){ 
    Q <- diag(nrow = a@V@Dim[1], ncol = a@V@Dim[1]) 
    for(i in rev(1:a@V@Dim[2])){ 
     Q <- Q - (a@V[ ,i] * a@beta[i]) %*% (a@V[ ,i] %*% Q) 
    } 
    return(Q[order(a@p), ][1:a@Dim[1], 1:a@Dim[1]]) 
    }else{ 
    stop(gettextf("gx.qr.Q() fails on class '%s'", class(a)[1])) 
    } 
} 

gx.qr.R <- function(a){ 
    if(is(a, "qr")){ 
    return(qr.R(a, complete = T)[ ,order(a$pivot)]) 
    }else if(is(a, "sparseQR")){ 
    if(length(a@q) == 0){ 
     return(a@R[1:a@Dim[1], ]) 
    }else{ 
     return(a@R[1:a@Dim[1] ,order(a@q)]) 
    } 
    }else{ 
    stop(gettextf("gx.qr.R() fails on class '%s'", class(a)[1])) 
    } 
} 

我已經測試了隨機設置的矩陣大小和稀疏度，它們工作順利。然而，這種風格「只是在不知道原因的情況下才能工作」，並且僅在這裏發佈供討論。因爲我沒有入侵「Matrix」包的實現細節。