如何从QR分解输出中获得Q? [英] How to get the Q from the QR factorization output?

问题描述

DGEQRF和SGEQRF 返回QR的Q部分以打包格式分解.拆箱似乎需要O(k^3)步骤(k个低排名产品)，而且似乎不是很简单.另外，对我来说，k顺序乘法的数值稳定性还不清楚.

DGEQRF and SGEQRF from LAPACK return the Q part of the QR factorization in a packed format. Unpacking it seems to require O(k^3) steps (k low-rank products), and doesn't seem to be very straightforward. Plus, the numerical stability of doing k sequential multiplications is unclear to me.

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

Does LAPACK include a subroutine for unpacking Q, and if not, how should I do it?

推荐答案

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

Yes, LAPACK indeed offers a routine to retrieve Q from the elementary reflectors (i.e. the part of data returned by DGEQRF), it is called DORGQR. From the 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.

使用C -wrapper LAPACKE从A完整计算Q和R可能看起来像这样(Fortran改编应该是直接的):

A complete calculation of Q and R from A using the C-wrapper LAPACKE could look like this (a Fortran adaption should be straight forward) :

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被破坏.

It's a piece of my own code, where I removed all checks to improve readability. For productive use you want to check that the input is valid and also mind the return values of the LAPACK calls. Note that the input A is destroyed.

这篇关于如何从QR分解输出中获得Q?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！