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

文章目录

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

写在前面

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

gsl_linalg.h

完全正交分解

完全正交分解可以看作是QR分解的推广，对于矩阵 A A A ( M × N M\times N M×N)，有如下分解：

A P = Q ( R 11 0 0 0 ) Z T A P=Q\left(\begin{array}{cc}R_{11} & 0 \\ 0 & 0\end{array}\right) Z^{T} AP=Q(R11​0​00​)ZT

其中：

P P P：转置矩阵， N × N N\times N N×N

Q Q Q：正交矩阵， M × M M\times M M×M

R 11 R_{11} R11​：上三角矩阵， r × r , r = r a n k ( A ) r\times r, r=rank(A) r×r,r=rank(A)

Z Z Z：正交矩阵， N × N N\times N N×N

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

对于不满秩 A A A，其对应的最小二乘解 ( min ⁡ x ∥ A x − b ∥ \min_x\|Ax-b\| minx​∥Ax−b∥) 不唯一，如果进一步加上条件 min ⁡ x ∥ x ∥ \min_x\|x\| minx​∥x∥，利用完备正交解，可以得到 x x x

x = P Z ( R 11 − 1 c 1 0 ) x=PZ\left(\begin{array}{c} R_{11}^{-1}c_1\\ 0\end{array}\right) x=PZ(R11−1​c1​0​)

其中 c 1 c_1 c1​ 是 Q T b Q^Tb QTb 的前 r r r 个元素

以吉洪诺夫 (Tikhonov) 正则化的形式表述下列问题

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

该问题的解为：

x = P Z ( y 1 0 ) x=PZ\left(\begin{array}{c} y_1\\ 0\end{array}\right) x=PZ(y1​0​)

其中 y 1 y_1 y1​ 长度为 r r r，可以根据下列方程求出：

( R 11 λ I r ) y 1 = ( c 1 0 ) \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​=(c1​0​)

该方程可以用 QR 分解方法求解

完全正交分解相关函数

完全正交分解

A = Q R Z P T A=QRZP^T A=QRZPT

Q , Z Q,Z Q,Z 一部分存在

A

中，另一部分存在

tau_Q

,

tau_Z

中

R 11 R_{11} R11​ 存在

A

中，

rank

计算方式，参考：GSL 系列 6 — 线性代数 3 — QR 分解

work

为工作空间，长度为 N N N

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 ( ∥ A x − b ∥ 2 + ∥ x ∥ 2 ) \min_x(\|Ax-b\|^2+\|x\|^2) minx​(∥Ax−b∥2+∥x∥2)，当 A A A 满秩时，不考虑正则项

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 ( ∥ A x − 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);
           

计算 A Z AZ AZ

A A A 必须为 N N N 列，

work

长度为 N N N

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