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

-wrapper LAPACKEを使用した完全な計算はQ、次のようになります (Fortran の適応は簡単なはずです)。RAC

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が破棄されることに注意してください。