用 GSL 求解超定方程组及矩阵的奇异值分解(SVD)
2016-11-20 12:58
337 查看
用 GSL 求解超定方程组及矩阵的奇异值分解(SVD)
最近在学习高动态图像(HDR)合成的算法,其中需要求解一个超定方程组,因此花了点时间研究了一下如何用 GSL 来解决这个问题。GSL 里是有最小二乘法拟合(Least-Squares Fitting)的相关算法,这些算法的声明在 gsl_fit.h 中,所以直接用 GSL 提供的 gsl_fit_linear 函数就能解决这个问题。不过我想顺便多学习一些有关 SVD 的知识。所以就没直接使用 gsl_fit_linear 函数。
SVD 分解的一些基本概念
关于 SVD 有两篇不错的科普文:A Singularly Valuable Decomposition: The SVD of a Matrix
We Recommend a Singular Value Decomposition
建议大家找来读读,这两篇文章似乎都已经有人翻译成中文了。
所谓 SVD,就是把一个矩阵 分解为三个特殊矩阵 、、 的乘积。
上面式子中的 表示矩阵的转置。分解之后的这三个矩阵还要满足些特殊条件,其中 和 是正交矩阵,也就是满足:
矩阵 是对角矩阵,只有主对角线上的元素非 。
因为矩阵 、、 的都具有很好的性质,所以这样的分解可以更好的帮助我们了解原始矩阵 的性质。
举例来说,如果矩阵 是个满秩方阵,那么 是可逆的。 的逆可以写为:
这里 和 因为是正交矩阵,所以 , 。 是对角矩阵,求逆也很简单,就是把主对角线上每个元素取个倒数而已。
GSL 中的相关函数
gsl 中提供了好几个函数来计算 SVD:gsl_linalg_SV_decomp 这个是最基本的,使用 Golub-Reinsch SVD 算法,一般我们用这个就够了。
gsl_linalg_SV_decomp_mod 这个是改进后的 Golub-Reinsch SVD 算法,当 时比 Golub-Reinsch SVD 算法要快。
gsl_linalg_SV_decomp_jacobi 这个算法用到了 Jacobi 正交化,号称计算结果比 Golub-Reinsch SVD 算法要更准确。
除此之外,还有个 gsl_linalg_SV_solve 函数。这个就是利用 SVD 的结果来求解线性代数方程组的。
把这几个函数组合一下就可以合成一个求解线性代数方程组 的函数了。
下面是函数代码:
void linearSolve_SVD(const gsl_matrix * A, const gsl_vector * b, gsl_vector * x) { int rows = A->size1; int cols = A->size2; gsl_vector * work = gsl_vector_alloc (cols); gsl_vector * S = gsl_vector_alloc (cols); gsl_matrix * U = gsl_matrix_alloc(rows, cols);; gsl_matrix * V = gsl_matrix_alloc(cols, cols); gsl_matrix_memcpy (U, A); // 为了不破坏 A 中原始的数据,这里全都拷贝到 U 中 gsl_linalg_SV_decomp( U, V, S, work ); gsl_linalg_SV_solve ( U, V, S, b, x ); gsl_vector_free(work); gsl_vector_free(S); gsl_matrix_free(V); gsl_matrix_free(U); }
当 是满秩方阵时,计算出来的 就是我们一般意义上的方程的解。
下面举一个具体的例子:
下面是测试代码:
void test1() { double a_data[] = {1.4, 2.1, 2.1, 7.4, 9.6, 1.6, 1.5, 1.1, 0.7, 5.0, 3.8, 8.0, 9.6, 5.4, 8.8, 4.6, 8.2, 8.4, 0.4, 8.0, 2.6, 2.9, 0.1, 9.9, 7.7}; gsl_matrix_view A = gsl_matrix_view_array (a_data, 5, 5); double b_data[] = {1.1, 1.6, 4.7, 9.1, 0.1}; gsl_vector_view b = gsl_vector_view_array (b_data, 5); gsl_vector * x = gsl_vector_alloc (5); linearSolve_SVD(&A.matrix, &b.vector, x); gsl_vector_fprintf (stdout, x, "%f"); qDebug() << ""; gsl_vector * bb = gsl_vector_alloc (5); gsl_blas_dgemv (CblasNoTrans, 1, &A.matrix, x, 0, bb); gsl_vector_fprintf (stdout, bb, "%f"); }
输出结果如下:
-5.208566
5.736694
-2.537472
-1.029814
0.968151
1.100000
1.600000
4.700000
9.100000
0.100000
可以看出计算结果还是很准确的。
当 的行数大于列数时求得的是最小二乘意义下的解,也就是 最小的解。下面给个例子:
测试代码如下:
void test3() { double a_data[] = {2, 4, 3, -5, 1, 2}; gsl_matrix_view A = gsl_matrix_view_array (a_data, 3, 2); double b_data[] = {11, 3, 6}; gsl_vector_view b = gsl_vector_view_array (b_data, 3); gsl_vector * x = gsl_vector_alloc (2); linearSolve_SVD(&A.matrix, &b.vector, x); gsl_vector_fprintf (stdout, x, "%f"); qDebug() << ""; gsl_vector * bb = gsl_vector_alloc (3); gsl_blas_dgemv (CblasNoTrans, 1, &A.matrix, x, 0, bb); gsl_vector_fprintf (stdout, bb, "%f"); }
计算结果如下:
3.090909
1.254545
11.200000
3.000000
5.600000
如果 不满秩,那么 是不唯一的。这时算出来的其中一个解。 下面给个例子:
方程很简单,口算就可以出结果,这个方程的解是:
下面用我们的代码计算一下。
void test4() { double a_data[] = {1, 2, 2, 4}; gsl_matrix_view A = gsl_matrix_view_array (a_data, 2, 2); double b_data[] = {3, 6}; gsl_vector_view b = gsl_vector_view_array (b_data, 2); gsl_vector * x = gsl_vector_alloc (2); linearSolve_SVD(&A.matrix, &b.vector, x); gsl_vector_fprintf (stdout, x, "%f"); qDebug() << ""; gsl_vector * bb = gsl_vector_alloc (2); gsl_blas_dgemv (CblasNoTrans, 1, &A.matrix, x, 0, bb); gsl_vector_fprintf (stdout, bb, "%f"); }
结果是:
-3.400000
3.200000
3.000000
6.000000
可以验算, 确实是方程的一个解。其实用 SVD 我们可以求出方程的全部解的,但是我们需要 和 的值,所以上面的 linearSolve_SVD 函数就不够用了。
下面我们将 SVD 相关的功能封装成一个类,以方便我们提取 和 的值。
另外,当我们一个 有多组 需要求解时,也只需要计算一次 SVD 分解,用下面的类能减少很多计算量。
头文件如下:
#ifndef GSLSINGULARVALUEDECOMPOSITION_H #define GSLSINGULARVALUEDECOMPOSITION_H #include <gsl/gsl_matrix.h> #include <gsl/gsl_vector.h> #include <gsl/gsl_blas.h> #include <gsl/gsl_linalg.h> #include <gsl/gsl_errno.h> void linearSolve_SVD(const gsl_matrix * A, const gsl_vector * b, gsl_vector * x); class GslSVD { public: GslSVD(); ~GslSVD(); int SV_decomp(const gsl_matrix * A); int SV_decomp_mod(const gsl_matrix * A); int SV_decomp_jacobi (gsl_matrix * A); int SV_solve(const gsl_vector *b, gsl_vector *x); gsl_vector * getVectorS(); gsl_matrix * getMatrixU(); gsl_matrix * getMatrixV(); int trimVectorS(double abseps); private: gsl_vector * S; gsl_matrix * U; gsl_matrix * V; void alloc_suv(int rows, int cols); }; #endif // GSLSINGULARVALUEDECOMPOSITION_H
cpp 文件如下:
#include "gsl_SVD.h" void linearSolve_SVD(const gsl_matrix * A, const gsl_vector * b, gsl_vector * x) { int rows = A->size1; int cols = A->size2; gsl_vector * work = gsl_vector_alloc (cols); gsl_vector * S = gsl_vector_alloc (cols); gsl_matrix * U = gsl_matrix_alloc(rows, cols);; gsl_matrix * V = gsl_matrix_alloc(cols, cols); gsl_matrix_memcpy (U, A); // 为了不破坏 A 中原始的数据,这里全都拷贝到 U 中 gsl_linalg_SV_decomp( U, V, S, work ); gsl_linalg_SV_solve ( U, V, S, b, x ); gsl_vector_free(work); gsl_vector_free(S); gsl_matrix_free(V); gsl_matrix_free(U); } int GslSVD::trimVectorS(double abseps) { int count = 0; for(int i = 0; i < S->size; i++) { if(fabs(gsl_vector_get(S, i)) < abseps) { count ++; gsl_vector_set(S, i, 0); } } return count; } gsl_vector * GslSVD::getVectorS() { if(S == NULL) return NULL; gsl_vector * s = gsl_vector_alloc(S->size); gsl_vector_memcpy(s, S); return s; } gsl_matrix * GslSVD::getMatrixU() { if(U == NULL) return NULL; gsl_matrix * u = gsl_matrix_alloc(U->size1, U->size2); gsl_matrix_memcpy(u, U); return u; } gsl_matrix * GslSVD::getMatrixV() { if(V == NULL) return NULL; gsl_matrix * v = gsl_matrix_alloc(V->size1, V->size2); gsl_matrix_memcpy(v, V); return v; } GslSVD::GslSVD() { S = NULL; U = NULL; V = NULL; } void GslSVD::alloc_suv(int rows, int cols) { if( S != NULL ) { gsl_vector_free(S); gsl_matrix_free(U); gsl_matrix_free(V); } S = gsl_vector_alloc (cols); U = gsl_matrix_alloc(rows, cols); V = gsl_matrix_alloc(cols, cols); } int GslSVD::SV_decomp(const gsl_matrix * A) { int rows = A->size1; int cols = A->size2; gsl_vector * work = gsl_vector_alloc (cols); alloc_suv(rows, cols); gsl_matrix_memcpy (U, A); // 为了不破坏 A 中原始的数据,这里全都拷贝到 U 中 int ret = gsl_linalg_SV_decomp( U, V, S, work ); gsl_vector_free(work); return ret; } int GslSVD::SV_decomp_mod(const gsl_matrix * A) { int rows = A->size1; int cols = A->size2; gsl_vector * work = gsl_vector_alloc (cols); gsl_matrix *X = gsl_matrix_alloc(cols, cols); alloc_suv(rows, cols); gsl_matrix_memcpy (U, A); // 为了不破坏 A 中原始的数据,这里全都拷贝到 U 中 int ret = gsl_linalg_SV_decomp_mod( U, X, V, S, work ); gsl_matrix_free(X); gsl_vector_free(work); return ret; } int GslSVD::SV_decomp_jacobi (gsl_matrix * A) { int rows = A->size1; int cols = A->size2; alloc_suv(rows, cols); gsl_matrix_memcpy (U, A); // 为了不破坏 A 中原始的数据,这里全都拷贝到 U 中 int ret = gsl_linalg_SV_decomp_jacobi( U, V, S ); return ret; } int GslSVD::SV_solve(const gsl_vector *b, gsl_vector *x) { if(U != NULL) { return gsl_linalg_SV_solve (U, V, S, b, x); } return -1; } GslSVD::~GslSVD() { if(S != NULL) { gsl_vector_free(S); gsl_matrix_free(V); gsl_matrix_free(U); } }
下面用这个类来计算一下刚才的问题:
void test5() { double a_data[] = {1, 2, 2, 4}; gsl_matrix_view A = gsl_matrix_view_array (a_data, 2, 2); GslSVD svd; svd.SV_decomp(&A.matrix); puts("S = "); gsl_vector_fprintf (stdout, svd.getVectorS(), "%f"); puts("\nV = "); gsl_matrix_fprintf (stdout, svd.getMatrixV(), "%f"); double b_data[] = {3, 6}; gsl_vector_view b = gsl_vector_view_array (b_data, 2); gsl_vector * x = gsl_vector_alloc (2); svd.SV_solve(b, x); puts("\nx = "); gsl_vector_fprintf (stdout, x, "%f"); }
结果如下:
S =
5.000000
0.000000
V =
-0.447214
-0.894427
-0.894427
0.447214
x =
-3.400000
3.200000
我们注意到 的第二个元素是 ,这表明 的对应列(第二列)是方程解的自由向量。所以我们方程的解可以写为:
大家可以验证一下,这个解是正确的。
另外,我写的类中还提供了一个 trimVectorS(double abseps) 函数,利用这个函数,可以将 所有小于 abseps 的项直接替换为 。之所以提供了这个函数,是因为由于计算误差等的影响, 中一些本应该是 的项可能计算出的结果不是 。用这个函数就可以解决这个问题。还有些矩阵,条件数很大,方程呈现病态,用这个函数也能解决些问题。
好了,就先写这么多。希望对大家有用。
相关文章推荐
- 用 GSL 求解超定方程组及矩阵的奇异值分解(SVD) 2
- 用 GSL 求解超定方程组及矩阵的奇异值分解(SVD)
- 用 GSL 求解超定方程组及矩阵的奇异值分解(SVD) 2
- 人工智能里的数学修炼 | 矩阵的花样分解:特征值分解(EVD)、相似对角化、QR分解、Schur分解、奇异值分解(SVD)的概念纠缠与详解
- 矩阵的奇异值分解(SVD)及其应用
- 矩阵的奇异值分解(SVD)(理论)
- 算法系列-----矩阵(九)-------------矩阵SVD---奇异值分解
- SVD(奇异值分解)及求解最小二乘问题
- 奇异值分解(SVD)原理详解及推导-矩阵相关
- 说文解字----矩阵分析(二)特征值特征向量 奇异值分解(SVD)
- SVD(奇异值分解)算法_计算任意N*M矩阵_C语言代码
- 矩阵分解----奇异值分解(SVD)原理和应用
- 奇异值分解SVD,矩阵范数,函数矩阵
- 矩阵的奇异值分解(SVD)及其应用
- SVD(奇异值分解)及求解最小二乘问题
- 矩阵的奇异值分解(SVD)(理论)
- 奇异值分解 SVD(Singular Value Decomposition)
- 使用矩阵分解(SVD)实现推荐系统
- MIT18.06线性代数课程笔记8:求解Ax=b、矩阵的秩以及矩阵的逆
- 求解线性方程组(SVD,QR,Gauss,LU)