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);