如何從QR分解輸出中獲取Q？來自LAPACK的

Question

DGEQRF and SGEQRF以打包格式返回QR分解的Q部分。打開包裝似乎需要O(k^3)步驟（k個低級產品），並且看起來不是非常簡單。另外，做連續乘法的數值穩定性對我來說還不清楚。

LAPACK是否包含解包Q的子程序，如果沒有，我該怎麼做？

Answer 1

是的，LAPACK確實提供了一個例程來從基本反射器中檢索Q（即由DGEQRF返回的數據部分），它被稱爲DORGQR。從describtion：

* DORGQR generates an M-by-N real matrix Q with orthonormal columns, 
* which is defined as the first N columns of a product of K elementary 
* reflectors of order M 
* 
*  Q = H(1) H(2) . . . H(k) 
* as returned by DGEQRF. 

從A使用C -wrapper LAPACKE看起來是這樣的（一個Fortran適應應該是直線前進）的Q和R一個完整的計算：

void qr(double* const _Q, double* const _R, double* const _A, const size_t _m, const size_t _n) { 
    // Maximal rank is used by Lapacke 
    const size_t rank = std::min(_m, _n); 

    // Tmp Array for Lapacke 
    const std::unique_ptr<double[]> tau(new double[rank]); 

    // Calculate QR factorisations 
    LAPACKE_dgeqrf(LAPACK_ROW_MAJOR, (int) _m, (int) _n, _A, (int) _n, tau.get()); 

    // Copy the upper triangular Matrix R (rank x _n) into position 
    for(size_t row =0; row < rank; ++row) { 
     memset(_R+row*_n, 0, row*sizeof(double)); // Set starting zeros 
     memcpy(_R+row*_n+row, _A+row*_n+row, (_n-row)*sizeof(double)); // Copy upper triangular part from Lapack result. 
    } 

    // Create orthogonal matrix Q (in tmpA) 
    LAPACKE_dorgqr(LAPACK_ROW_MAJOR, (int) _m, (int) rank, (int) rank, _A, (int) _n, tau.get()); 

    //Copy Q (_m x rank) into position 
    if(_m == _n) { 
     memcpy(_Q, _A, sizeof(double)*(_m*_n)); 
    } else { 
     for(size_t row =0; row < _m; ++row) { 
      memcpy(_Q+row*rank, _A+row*_n, sizeof(double)*(rank)); 
     } 
    } 
} 

這是一塊我自己的代碼，我刪除了所有檢查以提高可讀性。爲了高效使用，您需要檢查輸入是否有效，並且還要關注LAPACK調用的返回值。請注意，輸入A已被銷燬。

Answer 2

您正在尋找

DORMQR(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO) 

它計算Q * C其中Q = H(1) H(2) . . . H(k)從DGEQRF返回。只需使用C = I即可。

欲瞭解更多信息，請看here。