GSL 系列 6 — 線性代數 2 — LU 分解

文章目錄

• • 0 寫在前面
  • 1 LU 分解相關對象和函數
  • • 1.1 轉置矩陣對象
    • 1.2 LU 分解函數
    • 1.3 LU 求解方程
    • 1.4 LU 求逆
    • 1.5 LU 求行列式

0 寫在前面

關于 LU 分解的背景知識介紹，參見：GSL 系列 6 — 線性代數 1 — 背景知識 1 (LU 分解) 節，本篇隻說明相關函數

1 LU 分解相關對象和函數

1.1 轉置矩陣對象

轉置矩陣對象存儲着一列索引。第 j j j 個數為 k k k ，表示轉置矩陣第 j j j 列是相應機關矩陣的第 k k k 列，定義如下：

// gsl_permutation.h
struct gsl_permutation_struct
{
  size_t size;
  size_t *data;
};
typedef struct gsl_permutation_struct gsl_permutation;
           

1.2 LU 分解函數

将矩陣

A

進行 LU 分解，然後得到轉置矩陣

p

，

LU

矩陣儲存在

A

中，

signum

為 1 或 -1，代表交換的次數，奇數次為 -1，偶數次為 1

因為下三角（梯形）矩陣對角線為 1，不儲存，

A

剛好可以方向

LU

// gsl_linalg.h
int gsl_linalg_LU_decomp(gsl_matrix * A, gsl_permutation * p, int *signum);
int gsl_linalg_complex_LU_decomp(gsl_matrix_complex * A, 
                                  gsl_permutation * p, 
                                  int *signum);
           

1.3 LU 求解方程

solve

版本：給

LU

矩陣，

p

轉置矩陣，向量

b

，得到解向量

x

svx

版本: 給

LU

矩陣，

p

轉置矩陣，向量 b（就是

x

），得到解向量

x

，解向量存在

x

中，即

solve

的置換版本

refine

版本：應用疊代改進的辦法求解

x

，還需要給一個向量工作空間

work

，長度為

x

的長度

// gsl_linalg.h
int gsl_linalg_LU_solve (const gsl_matrix * LU,
                         const gsl_permutation * p,
                         const gsl_vector * b,
                         gsl_vector * x);
int gsl_linalg_complex_LU_solve (const gsl_matrix_complex * LU,
                                 const gsl_permutation * p,
                                 const gsl_vector_complex * b,
                                 gsl_vector_complex * x);
int gsl_linalg_LU_svx (const gsl_matrix * LU,
                       const gsl_permutation * p,
                       gsl_vector * x);
int gsl_linalg_complex_LU_svx (const gsl_matrix_complex * LU,
                               const gsl_permutation * p,
                               gsl_vector_complex * x);   
int gsl_linalg_LU_refine (const gsl_matrix * A,
                          const gsl_matrix * LU,
                          const gsl_permutation * p,
                          const gsl_vector * b,
                          gsl_vector * x,
                          gsl_vector * work);                               
int gsl_linalg_complex_LU_refine (const gsl_matrix_complex * A,
                                  const gsl_matrix_complex * LU,
                                  const gsl_permutation * p,
                                  const gsl_vector_complex * b,
                                  gsl_vector_complex * x,
                                  gsl_vector_complex * work);                                                   
           

1.4 LU 求逆

類似地，分為無置換版本 (

invert

) 和置換版本 (

invx

)，無置換時，存在矩陣

inverse

中，有置換時，存在

LU

中

// gsl_linalg.h
int gsl_linalg_LU_invert (const gsl_matrix * LU,
                          const gsl_permutation * p,
                          gsl_matrix * inverse);
int gsl_linalg_complex_LU_invert (const gsl_matrix_complex * LU,
                                  const gsl_permutation * p,
                                  gsl_matrix_complex * inverse);                          
int gsl_linalg_LU_invx (gsl_matrix * LU, const gsl_permutation * p);
int gsl_linalg_complex_LU_invx (gsl_matrix_complex * LU, const gsl_permutation * p);                                  
           

最好避免通過直接求逆的方法來求解線性方程組，因為線性方程組求解函數更有效率和更準确

1.5 LU 求行列式

求解 det ⁡ ( A ) \det(A) det(A)，需要給

LU

矩陣

// gsl_linalg.h
double gsl_linalg_LU_det (gsl_matrix * LU, int signum);
gsl_complex gsl_linalg_complex_LU_det (gsl_matrix_complex * LU,
                                       int signum);                                
           

求解 ln ⁡ ( ∣ det ⁡ ( A ) ∣ ) \ln(|\det(A)|) ln(∣det(A)∣)，需要給

LU

矩陣

// gsl_linalg.h
double gsl_linalg_LU_lndet (gsl_matrix * LU);
double gsl_linalg_complex_LU_lndet (gsl_matrix_complex * LU);   
           

求解 det ⁡ ( A ) / ∣ det ⁡ ( A ) ∣ \det(A)/|\det(A)| det(A)/∣det(A)∣，需要給

LU

矩陣

// gsl_linalg.h
int gsl_linalg_LU_sgndet (gsl_matrix * lu, int signum);
gsl_complex gsl_linalg_complex_LU_sgndet (gsl_matrix_complex * LU,
                                          int signum);