不同的QR分解結果與numpy和CULA

我正在以兩種不同的方式執行QR分解：使用標準numpy方法並使用CULA庫中實現的GEQRF LAPACK函數。這裏是蟒蛇（PyCULA用於訪問CULA）簡單的例子：不同的QR分解結果與numpy和CULA

from PyCULA.cula import culaInitialize,culaShutdown 
from PyCULA.cula import gpu_geqrf, gpu_orgqr 

import numpy as np 
import sys 

def test_numpy(A): 
    Q, R = np.linalg.qr(A) 
    print "Q" 
    print Q 
    print "R" 
    print R 
    print "transpose(Q)*Q" 
    print np.dot(np.transpose(Q), Q) 
    print "Q*R" 
    print np.dot(Q,R) 

def test_cula(A): 
    culaInitialize() 
    QR, TAU = gpu_geqrf(A) 
    R = np.triu(QR) 
    Q = gpu_orgqr(QR, A.shape[0], TAU) 
    culaShutdown() 
    print "Q" 
    print Q 
    print "R" 
    print R 
    print "transpose(Q)*Q" 
    print np.dot(np.transpose(Q), Q) 
    print "Q*R" 
    print np.dot(Q,R) 

def main(): 
    rows = int(sys.argv[1]) 
    cols = int(sys.argv[2]) 
    A = np.array(np.ones((rows,cols)).astype(np.float64)) 
    print "A" 
    print A 
    print "NUMPY" 
    test_numpy(A.copy()) 
    print "CULA" 
    test_cula(A.copy()) 

if __name__ == '__main__': 
    main()

它產生以下的輸出：

A 
[[ 1. 1. 1.] 
[ 1. 1. 1.] 
[ 1. 1. 1.]] 
NUMPY 
Q 
[[-0.57735027 -0.57735027 -0.57735027] 
[-0.57735027 0.78867513 -0.21132487] 
[-0.57735027 -0.21132487 0.78867513]] 
R 
[[-1.73205081 -1.73205081 -1.73205081] 
[ 0.   0.   0.  ] 
[ 0.   0.   0.  ]] 
transpose(Q)*Q 
[[ 1.00000000e+00 2.77555756e-17 0.00000000e+00] 
[ 2.77555756e-17 1.00000000e+00 0.00000000e+00] 
[ 0.00000000e+00 0.00000000e+00 1.00000000e+00]] 
Q*R 
[[ 1. 1. 1.] 
[ 1. 1. 1.] 
[ 1. 1. 1.]] 
CULA 
Q 
[[-0.57735027 -0.57735027 -0.57735027] 
[-0.57735027 0.78867513 -0.21132487] 
[-0.57735027 -0.21132487 0.78867513]] 
R 
[[-1.73205081 0.3660254 0.3660254 ] 
[-0.   0.   0.  ] 
[-0.   0.   0.  ]] 
transpose(Q)*Q 
[[ 1.00000000e+00 2.77555756e-17 0.00000000e+00] 
[ 2.77555756e-17 1.00000000e+00 0.00000000e+00] 
[ 0.00000000e+00 0.00000000e+00 1.00000000e+00]] 
Q*R 
[[ 1.   -0.21132487 -0.21132487] 
[ 1.   -0.21132487 -0.21132487] 
[ 1.   -0.21132487 -0.21132487]]

什麼是錯我的代碼？

來源

2014-04-24 grapescan

QR分解不是唯一的，如果矩陣不可逆。 –

@pv。正如你在我的例子中看到的那樣，CULA產生無效的R矩陣，所以Q * R不等於A.與可逆矩陣相同的問題（例如[[2,2]，[2,3]]）。 – grapescan

我只測試過幾次CULA，但是我發現它會在許多測試中產生不正確的結果（特別是計算矩陣的SVD）。我沒有調查太多，但它看起來像使用32位和64位浮點的問題。 – lmjohns3

我測試了R. CULA您的例子似乎提供相同的結果R.這裏是我的代碼：

#include <Rcpp.h> 
#include <cula.h> 

// [[Rcpp::export]] 
std::vector<float> gpuQR_cula(std::vector<float> x, const uint32_t nRows, const uint32_t nCols) 
{  
    std::vector<float> tau(nCols) ; 

    culaInitialize() ; 
    culaSgeqrf(nRows, nCols, &x.front(), nRows, &tau.front()) ; 
    culaShutdown() ; 

    Rcpp::Rcout << "Tau: " << tau[ 0 ] << ", " << tau[ 1 ] << ", " << tau[ 2 ] << "\n" ; 

    for(uint32_t jj = 0 ; jj < nCols ; ++jj) { 
     for(uint32_t ii = 1 ; ii < nRows ; ++ii) { 
      if(ii > jj) { x[ ii + jj * nRows ] *= tau[ jj ] ; } 
     } 
    } 

    return x ; 
}

你的矩陣：

(A <- matrix(1, 3, 3)) 

    [,1] [,2] [,3] 
[1,] 1 1 1 
[2,] 1 1 1 
[3,] 1 1 1 
n_row <- nrow(A) 
n_col <- ncol(A)

下面是CULA的結果：

matrix(gpuQR_cula(c(A), n_row, n_col), n_row, n_col) 

Tau: 1.57735, 0, 0 
      [,1]  [,2]  [,3] 
[1,] -1.7320509 -1.732051 -1.732051 
[2,] 0.5773503 0.000000 0.000000 
[3,] 0.5773503 0.000000 0.000000

下面是從R中的結果：

(qrA <- qr(A)) 
$qr 
      [,1]  [,2]  [,3] 
[1,] -1.7320508 -1.732051 -1.732051 
[2,] 0.5773503 0.000000 0.000000 
[3,] 0.5773503 0.000000 0.000000 

$qraux 
[1] 1.57735 0.00000 0.00000 

Q <- qr.Q(qrA) 
R <- qr.R(qrA) 
crossprod(Q) 

      [,1]   [,2]   [,3] 
[1,] 1.000000e+00 4.163336e-17 5.551115e-17 
[2,] 4.163336e-17 1.000000e+00 0.000000e+00 
[3,] 5.551115e-17 0.000000e+00 1.000000e+00 

Q %*% R 
    [,1] [,2] [,3] 
[1,] 1 1 1 
[2,] 1 1 1 
[3,] 1 1 1

我希望有幫助！

來源

2014-05-26 03:56:49 chris

這是一個棘手的問題，這裏的問題是Python使用行主要順序，但CULA使用列主要順序作爲R。只需查看CULA文檔以獲取更多詳細信息。

這裏您與scikit-CUDA例如：

import numpy as np 
import pycuda.gpuarray as gpuarray 
import pycuda.autoinit 
from skcuda import linalg 
linalg.init() 


# skcuda 
A = np.ones((3,3)) 
A_gpu = gpuarray.to_gpu(np.array(A, order='F')) 
Q , R = linalg.qr(A_gpu) 
Q, R = Q.get(), R.get() 
print Q.dot(R) #recovers A 
[[ 1. 1. 1.] 
[ 1. 1. 1.] 
[ 1. 1. 1.]] 

print Q.T.dot(Q) # As expected 
[[ 1.00000000e+00 -5.55111512e-17 1.11022302e-16] 
[ -5.55111512e-17 1.00000000e+00 -2.22044605e-16] 
[ 1.11022302e-16 -2.22044605e-16 1.00000000e+00]]

如果您使用，而不是（這是在Python默認）

A_gpu = gpuarray.to_gpu(np.array(A, order='C'))

，你會爲你上面貼得到同樣的錯誤的結果。

這個問題可能會導致一些問題，所以你必須非常小心，並照顧矩陣順序。

乾杯，本

來源

2015-11-25 03:52:10 benli

不同的QR分解結果與numpy和CULA

回答

相關問題