用 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，就是把一个矩阵 $A_{m \times n}$ 分解为三个特殊矩阵 $U_{m \times n}$、$S_{n \times n}$、$V_{n \times n}$ 的乘积。

$$A_{m \times n} = U_{m \times n} \cdot S_{n \times n} \cdot V_{n \times n}^T$$

上面式子中的 $T$ 表示矩阵的转置。分解之后的这三个矩阵还要满足些特殊条件，其中 $U_{m \times n}$ 和 $V_{n \times n}$ 是正交矩阵，也就是满足：

$$U^{T} \cdot U = I\\ V^{T} \cdot V = I$$

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

因为矩阵 $U_{m \times n}$、$S_{n \times n}$、$V_{n \times n}$ 的都具有很好的性质，所以这样的分解可以更好的帮助我们了解原始矩阵 $A$ 的性质。

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

$$A^{-1} = V \cdot S^{-1} \cdot U^{T}$$

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

GSL 中的相关函数

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

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

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

    把这几个函数组合一下就可以合成一个求解线性代数方程组 $A \cdot 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$ 就是我们一般意义上的方程的解。

    下面举一个具体的例子：

    $$\begin{pmatrix} 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 \end{pmatrix} \cdot x = \begin{pmatrix} 1.1\\ 1.6\\ 4.7\\ 9.1\\ 0.1 \end{pmatrix}$$

    下面是测试代码：

        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 \cdot x - b||_2$ 最小的解。下面给个例子：

    $$\begin{pmatrix} 2 & 4 \\ 3 & -5\\ 1 & 2\\ \end{pmatrix} \cdot x = \begin{pmatrix} 11\\ 3\\ 6 \end{pmatrix}$$

    测试代码如下：

        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$ 是不唯一的。这时算出来的其中一个解。 下面给个例子：

    $$\begin{pmatrix} 1 & 2 \\ 2 & 4 \end{pmatrix} \cdot x = \begin{pmatrix} 3\\ 6 \end{pmatrix}$$

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

    $$x = \begin{pmatrix} 1\\ 1 \end{pmatrix} + \begin{pmatrix} -2\\ 1 \end{pmatrix} \cdot 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 = \begin{pmatrix} -3.4\\ 3.2 \end{pmatrix} + \begin{pmatrix} -0.894427\\ 0.447214 \end{pmatrix} \cdot t$$

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

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

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