如何从 QR 分解输出中获取 Q

Question

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

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

Answer 1

是的，LAPACK确实提供了一个从基本反射器(即DGEQRF返回的数据部分(中检索Q的例程，它被称为DORGQR。从描述中：

*  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.

使用C包装器LAPACKE从A Q和R的完整计算可能如下所示(Fortran改编应该是直截了当的(：

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)

它计算从DGEQRF返回Q = H(1) H(2) . . . H(k) Q * C。只需使用C = I.

欲了解更多信息，请看这里。