文章目錄
-
- 0 寫在前面
- 1 QR 分解相關函數
-
- 1.1 QR 分解
- 1.2 QR 解包
- 1.3 QR 求解線性方程組
- 1.4 計算其他
- 1.5 特殊矩陣 QR 分解
- 1.6 帶列轉置的 QR 分解
0 寫在前面
關于 QR 分解的背景知識介紹，參見：GSL 系列 6 — 線性代數 1 — 背景知識 1
(QR 分解) 節，本篇隻說明相關函數
若無特别說明，本篇代碼均來自頭檔案
gsl_linalg.h
1 QR 分解相關函數
1.1 QR 分解
A = Q R A=QR A=QR
Householder 版:
将 A 分解為 Q R, Q 分别存儲在矩陣 A 的下三角部分 ( v v v) 和向量 tau ( τ \tau τ) 中，R 存儲在矩陣 A 的上三角部分（含對角）
推薦此方法用于中小型矩陣，或用于 M < N M<N M<N 時
遞歸版 (帶
_r
)：
将 A 分解為 Q R, Q 分别存儲在矩陣 A 的下三角部分 ( V V V) 和矩陣 T ( T T T) 中，R 存儲在矩陣 A 的上三角部分（含對角）
此方法要求 M ≥ N M\ge N M≥N，并且該方法對于“高瘦”型矩陣 ( M ≫ N M\gg N M≫N) 表現更佳
1.2 QR 解包
根據 QR 和 tau (T) 得到矩陣 Q 和矩陣 R
int gsl_linalg_QR_unpack(const gsl_matrix * QR, const gsl_vector * tau, gsl_matrix * Q, gsl_matrix * R);
int gsl_linalg_QR_unpack_r(const gsl_matrix * QR, const gsl_matrix * T, gsl_matrix * Q, gsl_matrix * R);
1.3 QR 求解線性方程組
- 根據 QR 分解得到的 QR 和 tau (或 T)，求解線性方程組 A x = b Ax=b Ax=b
int gsl_linalg_QR_solve (const gsl_matrix * QR, const gsl_vector * tau, const gsl_vector * b, gsl_vector * x);
int gsl_linalg_QR_solve_r (const gsl_matrix * QR, const gsl_matrix * T, const gsl_vector * b, gsl_vector * x);
// 置換版本
int gsl_linalg_QR_svx (const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * x);
-
求解 A x = b Ax=b Ax=b 的最小二乘解（使用歐幾裡得範數， A ( M > N ) A(M>N) A(M>N)）
Householder 版還輸出殘差向量 (residual)
遞歸版輸出的向量 x 包含兩個部分，x 的前 N N N 個元素為解，後 M − N M-N M−N 個元素的向量範數等于殘差範數，遞歸版還需要一個向量工作空間 (work)，長度為 N N N
int gsl_linalg_QR_lssolve (const gsl_matrix * QR, const gsl_vector * tau, const gsl_vector * b,
gsl_vector * x, gsl_vector * residual);
int gsl_linalg_QR_lssolve_r (const gsl_matrix * QR, const gsl_matrix * T, const gsl_vector * b,
gsl_vector * x, gsl_vector * work);
- 求解 R x = b Rx=b Rx=b
- 使用 QR 矩陣中的 R
int gsl_linalg_QR_Rsolve (const gsl_matrix * QR, const gsl_vector * b, gsl_vector * x);
int gsl_linalg_QR_Rsvx (const gsl_matrix * QR, gsl_vector * x);
- 直接使用 R (上三角矩陣)
int gsl_linalg_R_solve (const gsl_matrix * R, const gsl_vector * b, gsl_vector * x);
int gsl_linalg_R_svx (const gsl_matrix * R, gsl_vector * x);
- 求解 R x = Q T b Rx=Q^Tb Rx=QTb
1.4 計算其他
-
計算 Q T v Q^Tv QTv
遞歸版需要一個向量工作空間，長度 N N N
int gsl_linalg_QR_QTvec(const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * v);
int gsl_linalg_QR_QTvec_r(const gsl_matrix * QR, const gsl_matrix * T, gsl_vector * b, gsl_vector * work);
- 計算 Q v Qv Qv
-
計算 Q T B Q^TB QTB，矩陣 B B B 次元 ( M × K M\times K M×K)
遞歸版需要一個矩陣工作空間，次元為 N × K N\times K N×K
int gsl_linalg_QR_QTmat(const gsl_matrix * QR, const gsl_vector * tau, gsl_matrix * A);
int gsl_linalg_QR_QTmat_r(const gsl_matrix * QR, const gsl_matrix * T, gsl_matrix * B, gsl_matrix * work);
-
QR 分解的秩 1 更新
Q ′ R ′ = Q ( R + w v T ) Q^{'} R^{'}=Q(R+wv^T) Q′R′=Q(R+wvT)
-
計算 R R R 的倒數範數 (使用 1 範數)
1 ∥ R ∥ 1 ∥ R − 1 ∥ 1 \frac{1}{\|R\|_1\|R^{-1}\|_1} ∥R∥1∥R−1∥11
向量工作空間次元為 3 N 3N 3N
1.5 特殊矩陣 QR 分解
背景相關知識參見：GSL 系列 6 — 線性代數 1 — 背景知識 1 (特殊矩陣的 QR 分解) 節
将矩陣 ( S , A ) T (S,A)^T (S,A)T QR 分解，得到矩陣 R, Y, T，矩陣 R 存在矩陣 S 中，矩陣 Y 存在 A 中
1.6 帶列轉置的 QR 分解
Q R P T QRP^T QRPT 分解，
norm
是一個向量工作空間，長度為 N N N
// 分解之後 Q 的 v 部分和 R 存在 A 中
int gsl_linalg_QRPT_decomp (gsl_matrix * A,
gsl_vector * tau,
gsl_permutation * p,
int *signum,
gsl_vector * norm);
// 分解之後Q R 分别儲在 q r 中
int gsl_linalg_QRPT_decomp2 (const gsl_matrix * A,
gsl_matrix * q, gsl_matrix * r,
gsl_vector * tau,
gsl_permutation * p,
int *signum,
gsl_vector * norm);
利用 Q R P T QRP^T QRPT 分解求解 A x = b Ax=b Ax=b
int gsl_linalg_QRPT_solve (const gsl_matrix * QR,
const gsl_vector * tau,
const gsl_permutation * p,
const gsl_vector * b,
gsl_vector * x);
int gsl_linalg_QRPT_svx (const gsl_matrix * QR,
const gsl_vector * tau,
const gsl_permutation * p,
gsl_vector * x);
利用 Q R P T QRP^T QRPT 分解求 A x = b Ax=b Ax=b 的最小二乘解 ( A A A 滿秩， A ( M > N ) A(M>N) A(M>N))
int gsl_linalg_QRPT_lssolve (const gsl_matrix * QR,
const gsl_vector * tau,
const gsl_permutation * p,
const gsl_vector * b,
gsl_vector * x,
gsl_vector * residual);
利用 Q R P T QRP^T QRPT 分解求 A x = b Ax=b Ax=b 的最小二乘解 ( A A A 秩為
rank
， A ( M > N ) A(M>N) A(M>N))
int gsl_linalg_QRPT_lssolve2 (const gsl_matrix * QR,
const gsl_vector * tau,
const gsl_permutation * p,
const gsl_vector * b,
const size_t rank,
gsl_vector * x,
gsl_vector * residual);
求解 R P T x = Q T b RP^Tx=Q^Tb RPTx=QTb
int gsl_linalg_QRPT_QRsolve (const gsl_matrix * Q,
const gsl_matrix * R,
const gsl_permutation * p,
const gsl_vector * b,
gsl_vector * x);
求解 R P T x = b RP^Tx=b RPTx=b
int gsl_linalg_QRPT_Rsolve (const gsl_matrix * QR,
const gsl_permutation * p,
const gsl_vector * b,
gsl_vector * x);
int gsl_linalg_QRPT_Rsvx (const gsl_matrix * QR,
const gsl_permutation * p,
gsl_vector * x);
Q R P T QRP^T QRPT 分解的秩 1 更新
Q ′ R ′ = Q ( R + w v T P ) Q^{'}R^{'}=Q(R+wv^TP) Q′R′=Q(R+wvTP)
int gsl_linalg_QRPT_update (gsl_matrix * Q,
gsl_matrix * R,
const gsl_permutation * p,
gsl_vector * u,
const gsl_vector * v);
R R R 的倒數範數（使用 1 範數），向量工作空間
work
次元為 3 N 3N 3N
1 ∥ R ∥ 1 ∥ R − 1 ∥ 1 \frac{1}{\|R\|_1\|R^{-1}\|_1} ∥R∥1∥R−1∥11
估計 R R R 的秩，統計 R R R 的對角元素絕對值大于
tol
的個數，如果
tol
為負，使用如下預設值：
20 ( M + N ) e p s ( m a x ( ∣ d i a g ( R ) ∣ ) ) 20(M+N)eps(max(|diag(R)|)) 20(M+N)eps(max(∣diag(R)∣))
![C語言第四章自述2第四章 選擇結構程式設計[圖]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=)