用 GSL 求解超定方程组及矩阵的奇异值分解(SVD)

用 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，就是把一个矩阵 Am×n 分解为三个特殊矩阵 Um×n、Sn×n、Vn×n 的乘积。

Am×n=Um×n?Sn×n?VTn×n

上面式子中的 T 表示矩阵的转置。分解之后的这三个矩阵还要满足些特殊条件，其中 Um×n 和 Vn×n 是正交矩阵，也就是满足：

UT?U=IVT?V=I

矩阵 S 是对角矩阵，只有主对角线上的元素非 0。

因为矩阵 Um×n、Sn×n、Vn×n 的都具有很好的性质，所以这样的分解可以更好的帮助我们了解原始矩阵 A 的性质。

举例来说，如果矩阵 A 是个满秩方阵，那么 A 是可逆的。A 的逆可以写为:

A?1=V?S?1?UT

这里 V 和 U 因为是正交矩阵，所以 V?1=VT， U?1=UT。S 是对角矩阵，求逆也很简单，就是把主对角线上每个元素取个倒数而已。

GSL 中的相关函数

gsl 中提供了好几个函数来计算 SVD：

  • gsl_linalg_SV_decomp 这个是最基本的，使用 Golub-Reinsch SVD 算法，一般我们用这个就够了。
  • gsl_linalg_SV_decomp_mod 这个是改进后的 Golub-Reinsch SVD 算法，当 M?N 时比 Golub-Reinsch SVD 算法要快。
  • gsl_linalg_SV_decomp_jacobi 这个算法用到了 Jacobi 正交化，号称计算结果比 Golub-Reinsch SVD 算法要更准确。

除此之外，还有个 gsl_linalg_SV_solve 函数。这个就是利用 SVD 的结果来求解线性代数方程组的。

把这几个函数组合一下就可以合成一个求解线性代数方程组 A?x=b的函数了。

下面是函数代码：

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

当 A 是满秩方阵时，计算出来的 x 就是我们一般意义上的方程的解。

下面举一个具体的例子：

????????1.41.63.84.62.62.11.58.08.22.92.11.19.68.40.17.40.75.40.49.99.65.08.88.07.7?????????x=????????1.11.64.79.10.1????????

下面是测试代码：

    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

可以看出计算结果还是很准确的。

当 A 的行数大于列数时求得的是最小二乘意义下的解，也就是 ||A?x?b||2 最小的解。下面给个例子：

???2314?52????x=???1136???

测试代码如下：

    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

如果 A 不满秩，那么 x 是不唯一的。这时算出来的其中一个解。 下面给个例子：

(1224)?x=(36)

方程很简单，口算就可以出结果，这个方程的解是：

x=(11)+(?21)?t

下面用我们的代码计算一下。

    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

可以验算，(?3.4,3.2)T 确实是方程的一个解。其实用 SVD 我们可以求出方程的全部解的，但是我们需要 S 和 V 的值，所以上面的 linearSolve_SVD 函数就不够用了。

下面我们将 SVD 相关的功能封装成一个类，以方便我们提取 S 和 V 的值。

另外，当我们一个 A 有多组 x 需要求解时，也只需要计算一次 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

我们注意到 S 的第二个元素是 0，这表明 V 的对应列（第二列）是方程解的自由向量。所以我们方程的解可以写为：

x=(?3.43.2)+(?0.8944270.447214)?t

大家可以验证一下，这个解是正确的。

另外，我写的类中还提供了一个 trimVectorS(double abseps) 函数，利用这个函数，可以将 S 所有小于 abseps 的项直接替换为 0。之所以提供了这个函数，是因为由于计算误差等的影响，S 中一些本应该是 0 的项可能计算出的结果不是 0。用这个函数就可以解决这个问题。还有些矩阵，条件数很大，方程呈现病态，用这个函数也能解决些问题。

好了，就先写这么多。希望对大家有用。

时间： 2024-12-20 21:40:50

用 GSL 求解超定方程组及矩阵的奇异值分解(SVD)的相关文章

SVD分解 求解超定方程组

做平差的时候,需要解误差方程组,而 有的书本上说解线性的误差方程组,并不需要初值. 在查阅了测量平差书本之后,书里描述,一般是需要参数的初始值的. 这就产生了疑问. 因为非线性方程的线性化之后,舍掉了二次项之后的值,会造成平差模型的弱化.因此在进行非线性方程的平差过程中,一般是对改正值进行一个迭代计算,使其精化. 而线性化之后的各参数的系数中,包含了其他的未知参数,因此在计算的过程之中,必须使用初值. 原本就是线性方程组的平差模型,也可以直接使用SVD分解来解误差方程组. 1.解最小二乘超定方程

高斯消元法求解异或方程组： cojs.tk 539.//BZOJ 1770 牛棚的灯

高斯消元求解异或方程组: 比较不错的一篇文章:http://blog.sina.com.cn/s/blog_51cea4040100g7hl.html cojs.tk  539. 牛棚的灯 ★★☆   输入文件:lights.in   输出文件:lights.out   简单对比时间限制:1 s   内存限制:128 MB [问题描述] 贝希和她的闺密们在她们的牛棚中玩游戏.但是天不从人愿,突然,牛棚的电源跳闸了,所有的灯都被关闭了.贝希是一个很胆小的女生,在伸手不见拇指的无尽的黑暗中,她感到惊

【poj1830-开关问题】高斯消元求解异或方程组

第一道高斯消元题目~ 题目:有N个相同的开关,每个开关都与某些开关有着联系,每当你打开或者关闭某个开关的时候,其他的与此开关相关联的开关也会相应地发生变化,即这些相联系的开关的状态如果原来为开就变为关,如果为关就变为开.你的目标是经过若干次开关操作后使得最后N个开关达到一个特定的状态.对于任意一个开关,最多只能进行一次开关操作.你的任务是,计算有多少种可以达到指定状态的方法.(不计开关操作的顺序)0<=N<=29 我们用样例来模拟一下: 我的高斯消元求解异或方程组模版: 1 int gauss

小游戏 Lights Out (关灯) 的求解 —— 异或方程组

Author : Evensgn  Blog Link : http://www.cnblogs.com/JoeFan/ Article Link : http://www.cnblogs.com/JoeFan/p/4338003.html   游戏介绍 Lights Out (关灯)是一款据说在20世纪90年代就已经被设计出的小游戏,游戏的玩法十分简单. 首先,给定一个 n 行 m 列的矩形方格阵,每个格子上都有一盏灯. 初始时,有些灯是开着的,有些灯是关着的. 玩家每次进行一次操作,选中一盏

机器学习中的数学(5)-强大的矩阵奇异值分解(SVD)及其应用

机器学习中的数学(5)-强大的矩阵奇异值分解(SVD)及其应用 版权声明: 本文由LeftNotEasy发布于http://leftnoteasy.cnblogs.com, 本文可以被全部的转载或者部分使用,但请注明出处,如果有问题,请联系[email protected] 前言: 上一次写了关于PCA与LDA的文章,PCA的实现一般有两种,一种是用特征值分解去实现的,一种是用奇异值分解去实现的.在上篇文章中便是基于特征值分解的一种解释.特征值和奇异值在大部分人的印象中,往往是停留在纯粹的数学计

POJ 1222-EXTENDED LIGHTS OUT(高斯消元求解异或方程组)

题目地址:POJ 1222 题意:有一个5*6的矩阵,每个位置都表示按钮和灯,1表示亮,0表示灭.每当按下一个位置的按钮,它和它周围灯的状态全部翻转(题目中给出如何影响),问在这样的一个方阵中按下哪些按钮可以把整个方阵都变成灭的,这时1表示按了,0表示没按. #include <stdio.h> #include <math.h> #include <string.h> #include <stdlib.h> #include <iostream>

高斯消元求解异或方程组

POJ1830 开关问题 对于解异或方程组,系数可以采用二进制压缩,如果系数太多可以使用bitset,但是如果少一点就可以使用下述的写法,更加简单快速 使用bitset的写法更正常的没什么区别,只是对应的消除变为异或操作,另外行变换也会更加简单 #include <iostream> #include <cstdio> #include <cstring> #include <algorithm> using namespace std; int a[100

（转）机器学习中的数学(5)-强大的矩阵奇异值分解(SVD)及其应用

版权声明: 本文由LeftNotEasy发布于http://leftnoteasy.cnblogs.com, 本文可以被全部的转载或者部分使用,但请注明出处,如果有问题,请联系[email protected] 前言: 上一次写了关于PCA与LDA的文章,PCA的实现一般有两种,一种是用特征值分解去实现的,一种是用奇异值分解去实现的.在上篇文章中便是基于特征值分解的一种解释.特征值和奇异值在大部分人的印象中,往往是停留在纯粹的数学计算中.而且线性代数或者矩阵论里面,也很少讲任何跟特征值与奇异值有

强大的矩阵奇异值分解(SVD)及其应用

版权声明: 本文由LeftNotEasy发布于http://leftnoteasy.cnblogs.com, 本文可以被全部的转载或者部分使用,但请注明出处,如果有问题,请联系 前言: 上一次写了关于PCA与LDA的文章,PCA的实现一般有两种,一种是用特征值分解去实现的,一种是用奇异值分解去实现的.在上篇文章中便是基于特征值分解的一种解释.特征值和奇异值在大部分人的印象中,往往是停留在纯粹的数学计算中.而且线性代数或者矩阵论里面,也很少讲任何跟特征值与奇异值有关的应用背景.奇异值分解是一个有着