GSL 系列 6 — 线性代数 5 — 完全正交分解

文章目录

写在前面
完全正交分解
完全正交分解相关函数
- 完全正交分解
- 计算带正则的最小二乘解
- 解包
- 计算 AZAZAZ

写在前面

若无特别说明，本篇代码均来自头文件 gsl_linalg.h

完全正交分解

完全正交分解可以看作是QR分解的推广，对于矩阵 AAA (M×NM\times NM×N)，有如下分解：
AP=Q(R11000)ZTA P=Q\left(\begin{array}{cc}R_{11} & 0 \\ 0 & 0\end{array}\right) Z^{T}AP=Q(R11000)ZT
其中：
PPP：转置矩阵，N×NN\times NN×N
QQQ：正交矩阵，M×MM\times MM×M
R11R_{11}R11：上三角矩阵，r×r,r=rank(A)r\times r, r=rank(A)r×r,r=rank(A)
ZZZ：正交矩阵，N×NN\times NN×N

如果 AAA 满秩，R11=R,Z=IR_{11}=R,Z=IR11=R,Z=I，也就是带列转置的 QR 分解

对于不满秩 AAA，其对应的最小二乘解 (min⁡x∥Ax−b∥\min_x\|Ax-b\|minx∥Ax−b∥) 不唯一，如果进一步加上条件 min⁡x∥x∥\min_x\|x\|minx∥x∥，利用完备正交解，可以得到 xxx
x=PZ(R11−1c10)x=PZ\left(\begin{array}{c} R_{11}^{-1}c_1\\ 0\end{array}\right)x=PZ(R11−1c10)
其中 c1c_1c1 是 QTbQ^TbQTb 的前 rrr 个元素

以吉洪诺夫 (Tikhonov) 正则化的形式表述下列问题
min⁡x(∥Ax−b∥2+λ2∥x∥2)\min_x(\|Ax-b\|^2+\lambda^2\|x\|^2)xmin(∥Ax−b∥2+λ2∥x∥2)
该问题的解为：
x=PZ(y10)x=PZ\left(\begin{array}{c} y_1\\ 0\end{array}\right)x=PZ(y10)
其中 y1y_1y1 长度为 rrr，可以根据下列方程求出：
(R11λIr)y1=(c10)\left(\begin{array}{c} R_{11}\\ \lambda I_r\end{array}\right)y_1=\left(\begin{array}{c} c_1\\ 0\end{array}\right)(R11λIr)y1=(c10)
该方程可以用 QR 分解方法求解

完全正交分解相关函数

完全正交分解

A=QRZPTA=QRZP^TA=QRZPT

Q,ZQ,ZQ,Z 一部分存在 A 中，另一部分存在 tau_Q, tau_Z 中
R11R_{11}R11 存在 A 中，
rank 计算方式，参考：GSL 系列 6 — 线性代数 3 — QR 分解
work 为工作空间，长度为 NNN

int gsl_linalg_COD_decomp(gsl_matrix * A, gsl_vector * tau_Q, gsl_vector * tau_Z,gsl_permutation * p, size_t * rank, gsl_vector * work);
int gsl_linalg_COD_decomp_e(gsl_matrix * A, gsl_vector * tau_Q, gsl_vector * tau_Z,gsl_permutation * p, double tol, size_t * rank, gsl_vector * work);

计算带正则的最小二乘解

min⁡x(∥Ax−b∥2+∥x∥2)\min_x(\|Ax-b\|^2+\|x\|^2)minx(∥Ax−b∥2+∥x∥2)，当 AAA 满秩时，不考虑正则项

int gsl_linalg_COD_lssolve (const gsl_matrix * QRZT, const gsl_vector * tau_Q, const gsl_vector * tau_Z,const gsl_permutation * perm, const size_t rank, const gsl_vector * b,gsl_vector * x, gsl_vector * residual);

min⁡x(∥Ax−b∥2+λ∥x∥2)\min_x(\|Ax-b\|^2+\lambda\|x\|^2)minx(∥Ax−b∥2+λ∥x∥2)

工作空间矩阵 S 维度为 rank x rank，工作空间向量 work 维度为 rank

int
gsl_linalg_COD_lssolve2 (const double lambda, const gsl_matrix * QRZT, const gsl_vector * tau_Q, const gsl_vector * tau_Z,const gsl_permutation * perm, const size_t rank, const gsl_vector * b,gsl_vector * x, gsl_vector * residual, gsl_matrix * S, gsl_vector * work);

解包

int gsl_linalg_COD_unpack(const gsl_matrix * QRZT, const gsl_vector * tau_Q,const gsl_vector * tau_Z, const size_t rank, gsl_matrix * Q,gsl_matrix * R, gsl_matrix * Z);

计算 AZAZAZ

AAA 必须为 NNN 列，work 长度为 NNN

int gsl_linalg_COD_matZ(const gsl_matrix * QRZT, const gsl_vector * tau_Z, const size_t rank,gsl_matrix * A, gsl_vector * work);