0

我用拉帕克。我正在尝试在 C 中对复杂数据进行 QR 分解。为此，我编写了函数（基于 Haatschii 代码How to get the Q from the QR factorization output?）：

// Q - input: matrix that we expand / output: Q matrix
// R - output: R matrix
// rows - input: number of rows of Q
// columns - input: number of columns of Q
// rows >= columns condition is always met
void QR(lapack_complex_double * Q, lapack_complex_double * R, size_t rows, size_t columns){
    size_t i;
    lapack_complex_double* tau = malloc(columns*sizeof(lapack_complex_double));
    LAPACKE_zgeqrf(LAPACK_ROW_MAJOR, (int) rows, (int) columns, Q, (int) columns, tau); // returns the Q, R in a packed format
    // Copy the upper triangular Matrix R (columns x columns).
    for(i = 0; i < columns; ++i)
        memcpy(R+i*columns+i, Q+i*columns+i, (columns-i)*sizeof(lapack_complex_double));
    LAPACKE_zungqr(LAPACK_ROW_MAJOR, (int) rows, (int) columns, (int) columns, Q, (int) columns, tau); // returns the Q
    free(tau);
}

值得注意的是，Alireza 的作者也遇到了函数znugqr的问题，但他切换到函数zunmqr似乎幸福来了（LAPACK QR 分解）。我相信我的问题也与LAPACKE_zungqr矩阵 R 与其他方法相同，因此LAPACKE_zgeqrf可以成功运行。

但是最后，用 Mathematica（QRDecomposition函数）和 Python（numpy.linalg.qr函数）比较相似的结果（QR 分解），我看到矩阵 Q 是不同的，而矩阵 R 是相同的。

输入矩阵，为简单起见 5×5：

1 + 3j  6 + 8j    11 + 13j  16 + 18j  21 + 23j  
2 + 4j  7 + 9j    12 + 14j  17 + 19j  22 + 24j  
3 + 5j  8 + 10j   13 + 15j  18 + 20j  23 + 25j  
4 + 6j  9 + 11j   14 + 16j  19 + 21j  24 + 26j  
5 + 7j  10 + 12j  15 + 17j  20 + 22j  25 + 27j

输出 Q 矩阵（来自我的 C 代码）：前 3 列：

-0.07254 - 0.21764j    -0.61558 - 0.41039j     0.519770 - 0.06712j
-0.14509 - 0.29019j    -0.35909 - 0.25649j    -0.59817 + 0.211099j
-0.21764 - 0.36273j    -0.10259 - 0.10259j     0.035755 - 0.18619j
-0.29019 - 0.43528j     0.153896 + 0.051298j  -0.35605 + 0.007600j
-0.36273 - 0.50783j     0.410391 + 0.205195j   0.398709 + 0.034623j

最后两列：

-0.12316 - 0.06327j    -0.11940 + 0.303152j
 0.221078 + 0.491045j   0.084589 - 0.02148j
-0.06231 - 0.31146j     0.119483 - 0.80553j
-0.04594 - 0.59711j    -0.01512 + 0.462905j
 0.010343 + 0.480803j  -0.06954 + 0.060958j

输出 Q 矩阵（来自 Python 代码）：

-0.11670 - 0.06185j    -0.13105 + 0.301181j
 0.223111 + 0.487988j   0.096454 - 0.02009j
-0.08117 - 0.30874j     0.138515 - 0.80184j
-0.04015 - 0.59906j    -0.04217 + 0.459232j
 0.014923 + 0.481676j  -0.06174 + 0.061519j

（这里我只列出这个矩阵的最后 2 列。前 3 列是相同的）。

我计算了这些矩阵中各列的最大和平均差异。结论是这样的：前三列在水平上10^-11不同，后两列的差异分别是10^-3, 10^-2（肉眼可见的差异）。

随着矩阵大小的增加，观察到差异的增加，并且通常前 2-3 列很好地吻合。

也许有人可以帮助我？

4

1 回答 1

0

输入矩阵 A 是 5x5 但它的秩为 2。前两列是线性独立的，而 A 的后三列是前两列的线性组合。

对于 QR 分解，这意味着 Q 的最后三列不具有唯一性。QR 分解实现（例如 LAPACK、numpy 等）可以返回 (1) 相互正交和 (2) 的任意三列在 A^(perp) 中，这是一个正确的答案。有很多正确答案！The解决方案不是唯一的。

如果您想检查 LAPACK 返回的 Q 和 R（或任何 QR 分解实现），您可以（1）计算 Q'*Q 并检查您是否获得 5x5 单位矩阵，（请使用 BLAS HERK 函数这样做）和（2）计算 Q'*A 并检查您是否具有 5x5 上三角矩阵 R（由 QR 返回）。在您的情况下，您应该看到 R 的最后 3 行全为零，这表明 A 的最后三列是前两列的线性组合。要计算 Q'*A，您可以使用 BLAS GEMM。

我从 LAPACK 中获取了您的输入 A 和输出 Q，并为您检查了 Q'*Q 是身份，Q'*R 是上三角形，最后三行完全为零。所以 LAPACK 的 Q 输出对我来说看起来不错。是的，这个 Q 与另一个实现返回的不同，这是完全可能的。

一般来说，我们通过检查以下内容来检查 QR 分解的质量： (1) || A - 二维码 || / || 一个 || 小，(2) || 我 - Q^TQ || 很小，R 是上三角形。（取任何易于计算的范数。）

检查两个代码是否返回相同的输出不是一个好主意。由于输出 Q 和 R 不唯一性有两个原因：（1）当 A 等级不足时，请参见您给出的示例；(2) 对于任何 j，只要您相应地缩放 R 的第 j 行，对于 Q 的第 j 列使用模数 1 的复数进行任何列缩放都是可能的。这会改变 Q 和 R 因子，但仍然存在有效的 Q 和 R 因子。

此外，假设您强制 R 的所有对角元素为实数非负数（通过适当地重新调整 Q 的列和 R 的列）并且 A 是满秩的，那么您会期望 Q 和 R 的唯一性。但是检查对 (Q,R) 接近另一个与计算 Q 和 R 以进行 QR 分解的条件数有关，因此您可以看到相距很远的对（前向误差），即使它们具有良好的后向误差质量。

于 2022-01-27T18:56:41.740 回答