使用LAPACKE_dgeqrf/dorgqr生成行数多于列数的QR分解完整Q矩阵

问题背景 #

我正在使用Intel数学核心库（MKL）的LAPACKE_dgeqrf和LAPACKE_dorgqr函数，对以行优先顺序存储的(m×n)维通用输入矩阵A进行QR分解（其中m为矩阵A的行数，n为列数）。

所用C++代码 #

void qr (double *Q, double *R, double *A, const size_t m, const size_t n) {
    std::size_t k = std::max(std::size_t(1), std::min(m, n)); // 初等反射器的数量

    std::unique_ptr<double[]> tau(new double[k]); // 定义初等反射器的标量

    LAPACKE_dgeqrf(LAPACK_ROW_MAJOR, m, n, A, n, tau.get());

    // 生成R矩阵
    for (std::size_t i(0); i < n; ++i) {
        std::memset(R + i * n    , 0            , i       * sizeof(double));
        std::memcpy(R + i * n + i, A + i * n + i, (n - i) * sizeof(double));
    }

    // 生成Q矩阵
    LAPACKE_dorgqr(LAPACK_ROW_MAJOR, m, k, k, A, n, tau.get());

    if(m == n) {
        std::memcpy(Q, A, sizeof(double) * (m * m));
    } else {
        for(std::size_t i(0); i < m; ++i) {
            std::memcpy(Q + i * k, A + i * n, sizeof(double) * (k));
        }
    }
}

功能现状 #

• LAPACKE_dgeqrf会将(m×n)上三角矩阵R的上三角部分覆盖写入输入矩阵A的上三角区域，可直接提取R矩阵；A的下三角区域存储Householder变换相关信息，tau数组存储反射器的标量因子。
• 使用LAPACKE_dorgqr提取Q矩阵时：
  • 当m ≤ n时功能正常：m=n时A会被完整的(m×m)Q矩阵覆盖；m<n时A的前m列就是(m×m)的Q矩阵，可直接提取。
  • 当m > n时，LAPACKE_dorgqr仅输出Q矩阵的前n列，无法得到完整Q矩阵。

技术问题解答 #

1. 当输入矩阵行数大于列数时，如何生成Q矩阵中缺失的剩余列？ #

完整的Q是(m×m)正交矩阵，前n列是与A列空间张成相同空间的正交向量组（记为Q₁），剩余m-n列是Q₁的正交补（记为Q₂）。可以通过以下高效方式生成：

• 先通过dorgqr得到Q₁；
• 构造一个m×(m-n)的临时矩阵（比如单位矩阵的后m-n列），将其投影到Q₁的正交补空间；
• 对投影后的矩阵做QR分解，得到正交的Q₂；
• 最后将Q₁与Q₂拼接成完整的Q矩阵。

2. 是否是我对API的使用方式有误？ #

参数调用本身没有错误，但你对LAPACKE_dorgqr的功能范围理解有偏差。dorgqr只能生成由dgeqrf输出的Householder反射器对应的正交列——因为dgeqrf仅生成了n个反射器，对应Q的前n列，它不会生成超出这个范围的正交列。

3. 这是否是该API的正常行为？我是否需要通过Gram-Schmidt等方法从现有正交列向量生成剩余的正交列向量？ #

这是LAPACKE_dorgqr的正常行为：它的设计目标是生成与原矩阵列空间相关的正交列，剩余的正交补列不属于这个范畴，因此不会自动生成。
Gram-Schmidt是可选方案，但更推荐用Householder变换构造正交补，后者的数值稳定性更好，误差更小。

4. 由于我已经得到了R矩阵，是否需要通过求解方程组来获取完整的Q矩阵？ #

不需要。QR分解中的Q₂是Q₁的正交补，和R没有直接的方程组关系。求解方程组反而会引入额外数值误差且效率更低，直接构造正交补是更优的选择。

针对m>n场景的代码优化示例 #

// 当m > n时的Q矩阵扩展逻辑
if (m > n) {
    // 提取Q1（前n列）
    for(std::size_t i(0); i < m; ++i) {
        std::memcpy(Q + i * m, A + i * n, sizeof(double) * n);
    }
    // 初始化临时矩阵为单位矩阵的后m-n列
    std::unique_ptr<double[]> temp(new double[m*(m-n)]);
    std::memset(temp.get(), 0, m*(m-n)*sizeof(double));
    for (std::size_t i = 0; i < m-n; ++i) {
        temp.get()[(n+i)*m + i] = 1.0;
    }
    // 投影到Q1的正交补空间：temp = temp - Q1*(Q1^T * temp)
    std::unique_ptr<double[]> qt_temp(new double[n*(m-n)]);
    LAPACKE_dgemm(LAPACK_ROW_MAJOR, 'T', 'N', n, m-n, m, 1.0, Q, m, temp.get(), m-n, 0.0, qt_temp.get(), m-n);
    LAPACKE_dgemm(LAPACK_ROW_MAJOR, 'N', 'N', m, m-n, n, -1.0, Q, m, qt_temp.get(), m-n, 1.0, temp.get(), m-n);
    // 对投影后的矩阵做QR分解得到Q2
    std::unique_ptr<double[]> tau2(new double[m-n]);
    LAPACKE_dgeqrf(LAPACK_ROW_MAJOR, m, m-n, temp.get(), m-n, tau2.get());
    LAPACKE_dorgqr(LAPACK_ROW_MAJOR, m, m-n, m-n, temp.get(), m-n, tau2.get());
    // 将Q2拼接到Q的后m-n列
    for(std::size_t i(0); i < m; ++i) {
        std::memcpy(Q + i * m + n, temp.get() + i*(m-n), sizeof(double)*(m-n));
    }
}

内容的提问来源于stack exchange，提问作者Dimitrios Bouzas