文章目录
-
- 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=)